A deep learning method that identifies cellular heterogeneity using nanoscale nuclear features

Cells and culture conditions

The following cells were used in this study (all human, except the monkey cells): A549 cells (lung carcinoma, American Type Culture Collection (ATCC) CRL-185); B lymphocytes (GM12878); BJ fibroblasts (human foreskin skin fibroblasts) (ATCC CRL-2522); IMR90 (fibroblasts isolated from lung tissue, ATCC CCL-186); bone marrow MSCs; HeLa; Müller glia; myocardiocytes; spontaneously arising retinal pigment epithelium (ARPE-19); urine epithelial cells; and Vero (African green monkey) cells. Further, hiPSCs were induced from amniocytes, BJ fibroblasts, dermal fibroblasts, normal lung tissue fibroblasts (hiPS(IMR90)-4: WiCell, #WISCi004) periosteum cells, umbilical cord MSCs and urine epithelial cells. HeLa, BJ fibroblasts, IMR90 fibroblasts and myocardiocytes were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, #C11995500BT) supplemented with 10% fetal bovine serum (FBS) (Gibco, #10099141) and 1% penicillin–streptomycin (pen/strep) (Gibco, #15140122). Urine epithelial cells were cultured in Lonza REGM BulletKit (Lonza, #CC-3190) supplemented with 1% pen/strep. B lymphocytes were cultured in RPMI 1640 with l-glutamine (Corning Cellgro, #10-040-CVR) supplemented with 15% FBS and 1% pen/strep. Müller glia and ARPE-19 cells were cultured in DMEM/F-12, GlutaMAX (Gibco, #10565018) supplemented with 10% FBS and 1% pen/strep. Bone marrow MSCs were cultured in DMEM supplemented with 0.01% human FGF-basic (Peprotech, #100-18B), 10.00% FBS and 1.00% pen/strep. All the hiPSCs were cultured in mTeSR1 (Stem Cell, #85850) supplemented with 1% pen/strep. A549 cells were cultured in F-12K (Kaighn’s) medium (Gibco, Thermo Fisher Scientific, #21127022) supplemented with 10% FBS (Thermo Fisher Scientific, #10270106) and 1% pen/strep. Vero cells were culture in DMEM, high-glucose GlutaMAX Supplement (Thermo Fisher Scientific, #10566016) supplemented with 5% FBS and 1% pen/strep.

HSV-1 sample preparation and cell infection

HSV-1 (KOS strain) was grown in Vero cells. Virus preparations were titrated by a plaque assay, as previously reported³⁸. For imaging purposes, A549 cells were plated in borosilicate-glass-bottom eight-well chambers (μSlide ibidi, #80827) at a confluency of 20,000–30,000 cells cm^–2. A549 cells were cultured with media supplemented with 5 µM ethynil-deoxy-cytdine (EdC) (Sigma-Aldrich, #T511307) for 64 h before HSV-1 infection. After media wash, cells were mock infected or infected with HSV-1 at multiplicity of infection of 3 for 1 h at 37 °C in a humidified incubator with 5% CO₂. The medium was then removed and replaced by fresh DMEM medium (supplemented with 2% FBS). At 1 or 3 hpi, cells were fixed with 4% paraformaldehyde (Alfa Aesar, #43368) for 15 min at room temperature and then rinsed with PBS three times for 5 min each.

Immunolabelling for SR imaging (STORM)

Cells were plated in chambered coverglass (µ-Slide 8 Well Chamber Slide, ibidi, #80826) at a concentration of 40,000 cells cm^–2 in growth medium for 24 h. Cells were fixed with 4% paraformaldehyde (Leagene, #DF10315) for 15 min at room temperature and then rinsed with PBS three times for 5 min each. Cells were permeabilized with 0.4% (v/v) Triton X-100 (Acros Organics, 327371000) in PBS for 15 min and then blocked in blocking buffer (10.00% bovine serum albumin, #B2064-50G; 0.01% (v/v) Triton X-100 in PBS) for 1 h at room temperature. Cells were incubated with primary antibodies (anti-RNA polymerase II CTD repeat YSPTSPS (phospho S5) antibody, Abcam, #ab5131; Histone H3 antibody, Active Motif, #39763) in blocking buffer at 1:50 dilution for STORM of both single- and dual-colour imaging, overnight at 4 °C. On the second day, cells were washed three times for 5 min each with wash buffer (2.00% bovine serum albumin and 0.01% Triton X-100 in PBS) and incubated in secondary antibody. For single-colour imaging of Pol II, the secondary antibody (donkey anti-rabbit IgG H&L (Alexa Fluor 647), Abcam, #ab150075) was added at 1:200 dilution in blocking buffer for 1 h at room temperature. For dual-colour imaging, home-made³⁹ dye pair (Alexa Fluor 405 NHS Ester, Thermo Fisher, #A30000; Cy3 Maleimide Mono-Reactive Dye Pack, Sigma-Aldrich, #GEPA23031; Alexa Fluor 647 NHS Ester (succinimidyl ester), Thermo Scientific, #A20006) was used to label secondary antibodies (AffiniPure donkey anti-rabbit IgG (H + L), Jackson ImmunoResearch, #711005152; Peroxidase AffiniPure goat anti-mouse IgG, Fcγ subclass 1 specific, Jackson ImmunoResearch, #115-005-205), which were added at a 1:50 dilution in blocking buffer and were incubated for 1 h at room temperature. Cells were then washed three times for 5 min each with wash buffer, and once with PBS. Samples were stored in PBS at 4 °C before imaging.

For Extended Data Fig. 4h, anti-RNA polymerase II CTD repeat YSPTSPS (phospho S5; Abcam #ab5408) was used for STORM and anti-NOLC1 (Abcam #ab184550) was used to mark the nucleoli (both diluted 1:50). Secondary antibodies (goat anti-mouse IgG Alexa Fluor 647, Thermo Scientific #A21235; and goat anti-rabbit IgG, Oregon Green 488, Thermo Scientific #O-11038) were used at 1:250 dilution. Immunostaining was performed using the protocol given above. The complete STORM protocol can be found in another work⁴⁰.

STORM imaging

STORM imaging was performed as follows: (1) for proteins, using an N-STORM microscope (Nikon) equipped with a CFI SR HP Apochromat TIRF 100× AC oil objective and an ORCA-Flash4.0 V3 digital CMOS camera; (2) for DNA, using Oxford Nanoimager (ONI); and (3) for Pol II and conventional for NOLC1 (Extended Data Fig. 4h), using an N-STORM 4.0 microscope (Nikon) with an iXon Ultra 897 camera (Andor) and a CFI HP Apochromat TIRF ×100 1.49 oil objective, under the highly inclined and laminated optical sheet illumination mode. NIS element (4.60 and 5.21) software was used for N-STORM image acquisition. Only nuclei that were morphologically coherent with the characteristics expected of an interphase nucleus were imaged.

For single-colour STORM imaging of Pol II, continuous imaging acquisition was performed with simultaneous 405 and 647 nm illumination of the sample and 10 ms exposure time for 60,000 frames. The 647 nm laser was used at constant 3.5 kW cm^–2, and the 405 nm laser power was gradually increased over the imaging period.

For the dual-colour STORM imaging of Pol II and H3, a double activator and single reporter strategy was used by combining AF405_AF647 anti-mouse and Cy3_AF647 anti-rabbit secondary antibodies. Sequential imaging acquisition was performed (one frame of 405 nm activation, followed by four frames of 647 nm reporter; and one frame of 560 nm activation, followed by four frames of 647 nm reporter) with 20 ms exposure time for 120,000 frames. The 647 nm laser was used at a constant 3.5 kW cm^–2 power density, whereas the powers of the 405 and 560 nm lasers were gradually increased over the imaging period.

Single-colour imaging of DNA images was performed on the Nanoimager S Mark II from ONI equipped with an Olympus 1.4-numerical-aperture ×100 oil-immersion super apochromatic objective and a Hamamatsu sCMOS Orca-Flash4.0 V3. Continuous imaging acquisition was performed with simultaneous 647 nm illumination of the sample, 10 ms exposure time for 45,000 frames.

During SR imaging, freshly prepared buffer was added to the sample: 100 mM cysteamine MEA (Sigma-Aldrich, #30070), 1% Glox solution (0.5 mg ml^–1 glucose oxidase, 40 mg ml^–1 catalase; Sigma-Aldrich, #G2133 and #C100), 5% glucose (Sigma-Aldrich, #G8270) in PBS. The buffer was exchanged every hour to avoid enzymatic activity decay and maximize the signal-to-noise ratio of the localizations.

For the ONI images, localization lists were obtained using the integrated simultaneous localization software from ONI (ONI NimOS v. 10.5) and subsequently converted into files compatible with Insight3 (ref. ⁴¹) for post-processing using custom-built software in MATLAB 2016a.

N-STORM and ONI microscopes have a resolution of 160 and 117 nm pixel^–1, respectively. For ONI images, molecule localizations were extracted and exported with ONI NimOS v.10.5 and then converted into insight3-compatible files for post-processing with MATLAB 2016a. For N-STORM images, molecule localizations were extracted using Insight3 (ref. ⁴¹), as described earlier^39,42. Specifically, point spread functions (PSFs) from the emission from single fluorophores were identified at individual frames of the acquired videos, based on restrictive set thresholds of intensity, PSF size and symmetry. They were fit to a two-dimensional Gaussian from which the centroid of the PSF was obtained and the x and y positions for 2D STORM imaging were then obtained from the PSF.

All the images were acquired using HiLo⁴³ to maximize the signal-to-noise ratio in the nucleus. Focus plane was set to z of the maximum nucleus area to ensure the representativity of slices. Cells with uneven intensity profiles or localization artefacts were manually discarded before proceeding with the DL pipeline.

EdC incorporation and DNA labelling

To label DNA, a 48 h incorporation of 5 µM EdC (5-ethynyl-2′-deoxycytidine; Sigma-Aldrich, #T511307) was performed using the following cells: HeLa, BJ fibroblast, myocardiocytes, urine epithelial cells, B lymphocytes, Müller glia, ARPE-19, MSCs and hiPSCs induced from amniocytes, BJ fibroblasts, dermal fibroblasts, periosteum cells, umbilical cord MSCs and urine epithelial cells. Cells were plated in chambered coverglass (µ-Slide 8 Well Chamber Slide, ibidi, #80826) at a concentration of 40,000 cells cm^–2 in growth medium supplemented with 5 µM EdC for 24 h. For somatic cells, the growth medium was replaced on the second day with resting medium supplemented with EdC and incubated for another 24 h. For hiPSCs, 10 µM ROCK Inhibitor (Y-27632, Millipore, SCM075) was added together with 5 mM EdC for the first 24 h and then was replaced with medium supplemented with 5 mM EdC only. Cells were fixed with 4% paraformaldehyde (Leagene, #DF10315) for 15 min at room temperature and then rinsed with PBS three times for 5 min each. Cells were permeabilized with 0.4% Triton X-100 in PBS for 15 min and rinsed with PBS three times for 5 min each. The click chemistry reaction was performed by incubating cells for 45 min at room temperature in click chemistry buffer (150 mM HEPES pH 8.2 (Sigma-Aldrich, #7365-45-9), 50 mM aminoguanidine (Sigma-Aldrich, #396494), 100 mM l-ascorbic acid (Sigma-Aldrich, #A92902), 1 mM CuSO₄ (Sigma-Aldrich, #C1297), 2.0% glucose (Sigma-Aldrich, #G8270), 0.1% Glox solution (see the ‘STORM imaging’ section) and 10 mM AF647 azide (Thermo Fisher Scientific, #A10277)), protected from light. After three washes with PBS, STORM imaging was directly carried out for single-colour DNA imaging experiments (see the ‘STORM imaging’ section).

(ss)RT-qPCR amplification of asincRNAs

RNA extraction

hiPSCs (generated from IMR90, WiCell, no. WISCi004) and IMR90 cells at 70–80% confluency were washed with RNase-free PBS before RNA isolation using a Qiagen RNeasy Mini Kit (Qiagen, #74106). The isolated RNA was treated with DNase I enzyme (Thermo Fisher Scientific, #18068015) and precipitated with 0.8 M trisodium citrate and 1.2 M NaCl.

Strand-specific complementary DNA synthesis

IGS region 28 was chosen based on its identification as one of the regions prominently transcribed by Pol II in the nucleolus³⁵. The complementary DNA synthesis reaction was set up for the transcript of interest using M-MLV reverse transcriptase (Thermo Fisher Scientific, #28025013). For this reaction, we used the following strand-specific primers in concentrations of 10 µM:

The segments were amplified for 10 min at 25 °C, 60 min at 37 °C and 25 min at 70 °C. Samples were diluted 1:10.

qPCR

Real-time qPCR was performed using a ViiA 7 Real-Time PCR System (Thermo Fisher Scientific). The qPCR reactions (10 µl total) each contained LightCycler 480 SYBR Green I Master (Roche Diagnostics, #04887352001), 10 µM of each of the forward (strand-specific) and reverse (ssTag) primers, and 2 ng of diluted complementary DNA, as detailed below:

Amplification of the 7SK region was used for normalization. PCR reactions used the following parameters: 1 cycle of 95 °C for 5 min; 45 cycles of 95 °C for 10 s, 58 °C for 10 s and 72 °C for 19 s; and a final melting curve cycle of 60 °C to 97 °C in 0.05 °C s^–1 steps. The quantification of each IGS was normalized for the 7SK region.

Hardware and software

All the computational operations were conducted using the Center for Genomic Regulation cluster infrastructure, which comprises ten NVIDIA GeForce 2080 Ti with 11 GB of memory per GPU. The DL computations were performed using the PyTorch framework version 1.10, with CUDA drivers version 11.3 and Python version 3.9.7. For ONI images, molecule localizations were extracted and exported with ONI NimOS v. 10.5 and then converted into Insight3-compatible files for post-processing with MATLAB 2016a. For N-STORM images, molecule localizations were extracted using Insight3. The images were segmented by manually selecting the nuclear areas corresponding to each cell nucleus with Fiji⁴⁴. Plots and statistical tests shown in Fig. 4i–k were produced with R version 4.2.2 and GraphPad Prism version 8.01.

Datasets

We collected a total of 349 dual-colour images, including 171 from human somatic cell nuclei and 178 from hiPSC nuclei, featuring H3 and Pol II. Additionally, 384 single-colour images of Pol II were obtained, with 187 from somatic cells and 197 from hiPSCs. Somatic cells comprised different cell types (B lymphocytes, bone marrow MSCs, Müller glia, myocardiocytes, BJ fibroblasts, ARPE-19 and urine epithelial cells) and hiPSCs were generated from some of the somatic cells included in the study (Supplementary Table 1).

We also collected single-colour images of Pol II (33) and dual-colour images of H3 and of Pol II (30) in the nuclei of HeLa cells, single-colour images of DNA in human somatic cells (68), hiPSCs (117) as well as images of DNA in nuclei of mock-infected (81) and HSV-1-infected A549 cells at 1 hpi (98) and 3 hpi (85). To avoid sex chromosome bias during model training, we included images from both male and female cells.

The molecule localization coordinates of the selected region of interest (that is, nucleus) of dual-colour Pol II and H3 were multiplied either by 10 or 30 corresponding to a resolution of 16.0 and 5.3 nm pixel^–1, respectively, and then rendered into images by computing a Gaussian kernel density estimate using Scikit-learn⁴⁵ KDTree function. The images were finally resized with the desired dimensions to be used as an input of the CNN and the aspect ratio was maintained.

To enable the model to differentiate between the background of the image outside the nucleus and regions within the nucleus lacking localizations, we set the background of full-size ×10 rendered images to white.

Both single- and dual-colour images rendered at ×10 magnification (Fig. 1a) were directly used as train/validation and test sets (see the ‘Workflow for identifying the best-performing model’ section), whereas dual-colour ones rendered at ×30 magnification were divided in non-overlapping patches (Extended Data Fig. 2a) of three different sizes (256 × 256, 512 × 512 and 768 × 768 pixels). Only patches with more than 75% of the pixels containing localizations were used as training/validation and test sets.

To avoid changing the aspect ratio, which would have modified the spatial relationship of the imaged molecules, each image was first resized on its longest axis to the desired value (for example, 768 pixels for full-size images) by keeping the aspect ratio and was subsequently padded with the background colour on the shortest axis. For all the other datasets (that is, single-colour Pol II, DNA and dual-colour Pol II and H3 HeLa), we just rendered and used images from ×10 magnification of the molecule localizations.

Workflow for identifying the best-performing model

The dual-colour dataset, consisting of 349 full-size images rendered at ×10 magnification, underwent an 80/20 split into training and validation sets, keeping the proportion of each cell class present in the dataset. Employing a stratified K-fold CV (with five folds)²⁵, we explored 11 different CNN architectures, each trained for 300 epochs using images resized to 512 × 512 pixels. Model hyperparameters remained constant across architectures for fair comparisons (Extended Data Fig. 1a). The best-performing architecture was chosen based on the average accuracy and loss over the five folds. Data augmentation techniques, including flipping, random rotation and colour adjustments, were applied using the Albumentations library⁴⁶. Region dropout and the GridMask technique⁴⁷ were explored to enhance training data diversity (as shown in Extended Data Fig. 1a). The use of the GridMask technique with a ratio of mask holes to unit size set at 0.3 and grid unit sizes ranging from 128 to 384 led to the highest average accuracy and lowest average loss values when assessed through five-fold CV (Extended Data Fig. 1a).

Systematically, we reset the model weights with the best-performing parameters based on the validation set during the training stage. The final architecture with its optimal hyperparameters, named AINU, was used for training, validation and evaluation across various datasets, including single-colour images of Pol II localizations collected in somatic cells and hiPSCs; dual-colour images of Pol II and H3 collected from HeLa cells; and single-colour images of DNA collected from somatic cells, hiPSCs and HSV-1-infected A549 cells.

To explore alternative training approaches, we employed a patch dataset, using patches extracted from the dual-colour STORM images rendered at higher resolutions such as ×30 instead of ×10 (Extended Data Fig. 1b, dataset (ii)). This approach was undertaken to increase the training dataset and to uncover potentially hidden nuclear features that could improve the model performance. To collect patches, the localization coordinates of each blinking fluorophore were multiplied by a factor of 30, effectively rendering images at a ×30 magnification relative to the original camera pixel size. These high-resolution images were then subdivided into smaller, non-overlapping patches, of 256 × 256 or 512 × 512 pixels, effectively increasing the size of the dataset.

The performance of the AINU model trained with patches was evaluated using a separate withheld test set, selected randomly and thus different from those used with full-size images (Extended Data Fig. 1b, grey sector in dataset (ii)). For evaluation, we used various metrics, including the F1 score, validation accuracy and loss, AUC, precision and recall^48,49. Due to the intrinsic differences between the full-size dual-colour images and patch-based images, fine-tuning of AINU was needed when training with patches or simulated images.

When evaluating the model trained with patches, a classification threshold of 50% was applied as the minimum number of patches correctly assigned to the positive class. This means that an image was correctly predicted if more than 50% of its patches were assigned to the true class.

No transformations that altered the aspect ratio of the image were applied to augment the data, as this would have resulted in a distorted spatial relationship of the imaged molecules, which is a key characteristic feature of somatic and hiPSCs.

Because those models were not trained five times with random dataset splits due to high resource consumption, the performance scores were reported with the 95% confidence interval calculated with the normal approximation method.

When training was performed using patches, AINU exhibited optimal validation performance when using non-overlapping patches of 256 × 256 pixels. However, when challenged with patches from the test set of full-size images (Extended Data Fig. 1b, dataset (i)), it showed lower performance than the original full-size image training strategy (see the ‘Results’ section). This suggests that using only partial sections of the nucleus, even if at higher magnification, does not improve the model’s overall performance.

Interpretable AI and cluster analysis

To evaluate the contribution of each input feature (image pixels) on the model’s output and prediction accuracy and to gain an understanding of the model’s decisions, the Captum framework³¹ and its occlusion-based method were used. This approach allowed to visually inspect the areas of the image that exhibited positive or negative attribution scores for each correctly predicted sample (Extended Data Fig. 4a). Positive attribution scores indicate areas of the image that contribute to the model’s prediction, whereas negative attribution scores highlight regions whose removal increases the probability of the predicted class. To apply the occlusion technique, we systematically slid a 64 × 64 pixel window, with a 32 pixel stride, across each image generated from the Pol II localizations. Subsequently, we superimposed the heat map of attribution scores onto the original image, facilitating the identification of the image regions most influential in the model’s classification decision.

These results were further validated using the CAM³³ method, which highlighted specific regions within the nucleoli that are critical for classification. Additionally, the differences in localization quantification, distribution and organization inside the nucleoli of hiPSCs and somatic cells were verified using a distance-based clustering algorithm, as detailed elsewhere⁵ (Fig. 4i,j).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.