Mapping human tissues with highly multiplexed RNA in situ hybridization – Nature Communications

Human tissue samples

Human brain

Human Brain tissue was obtained from the University of Washington Biorepository and Integrated Neuropathology (BRaIN) Laboratory under UW School of Medicine and HIPAA compliance. Informed consent was obtained for the use of data and samples. One donor brain with postmortem interval ≤12 h and RIN score ≥7 was selected for DART-FISH assay. Regions were identified and isolated utilizing architectural landmarks, aided by the Allen Brain Human Brain Atlas⁸¹. Multiple parallel 10-μm-thick cryosections were taken from the tissue block and mounted onto vectabond-coated 24 x 60 mm No.1.5 coverslips (Azer Scientific, 1152460). Brain cryosections were stored at −80 °C until use.

Human kidney

Kidney tissue was obtained from the Kidney Translational Research Center (KTRC) biorepository under a protocol approved by the Washington University Institutional Review Board (IRB 201102312). Informed consent was obtained for the use of data and samples. The kidney tissue was dissected from the whole kidney and freshly frozen in Optimal Cutting Temperature embedding media in cryomolds on a liquid nitrogen chilled metal block and stored at −80 ^oC until ready for experimental use⁷⁴. 10-μm-thick sections were cut from the frozen blocks for DART-FISH and flanking sections were used for histopathological assessment by a renal pathologist.

Reagents and enzymes

All reagents were listed as in Supplementary Data 1.

Gene selection

A list of genes was selected based on differential expression analysis of snRNA-seq data from human primary motor cortex^46,48,50 and a few curated marker genes were added manually to target 121 genes in the human M1C. Human kidney gene selection was performed by gpsFISH^82,83 to distinguish subclass level 2 annotation in our kidney reference atlas⁷³. snRNA-seq data from the kidney reference atlas with cell type annotation at subclass level 2 was used as input of gpsFISH. Curated marker genes from prior knowledge were also included as input. The size of the gene panel was set to 300. We ran the optimization for 100 iterations to ensure convergence although the optimization converged around iteration 50.

Probe design and production

DART-FISH probe design

For short genes (length < 1.5kb), we defined the constitutive exon as the union of all isoforms in GencodeV41. For other genes, the constitutive exons were defined as regions in RefSeq where at least (33% for the brain, 50% for the kidney) of isoforms overlap. We used a modified version of ppDesigner³⁸ (https://github.com/Kiiaan/sppDesigner) to find padlock target sequences along the constitutive exons. ppDesigner was run on two settings: 1) no overlap between probes allowed, 2) overlap of up to 20nt allowed. Individual arms were constrained between 17nt and 22nt long with the total target sequences no longer than 40nt. The resulting target sequences were aligned to GRCh38/hg38 with BWA-MEM⁸⁴ and sequences with MAPQ < 40 or secondary alignment were removed. We further removed probes that have GATC (DpnII recognition site). For the brain, a maximum of 50 probes per gene were selected prioritizing the non-overlapping set. For the kidney, a maximum of 40 probes per gene were selected with no overlap. Finally, the target sequences were concatenated with amplification primer sequences, universal sequence, and gene-specific decoder sequences to produce final padlock probe sequences (Supplementary Fig. 3c) and were ordered as an oligo pool from Twist Bioscience (South San Francisco, CA). Amplification primer pairs pAP1V41U and AP2V4 were used for the kidney probe set, while the brain probe set was amplified with AP1V7U and AP2V7 primer pair (Supplementary Data 1).

To select a set of barcodes, we computationally created all possible barcodes in the compact format: an (n) digit barcode with “1”, “2” and “3” representing each of the three fluorescent channels and “0” indicating off cycles. For example, the barcode for RORB in Fig. 1c is “132000” in the 6-digit format. This amounted to 480 and 840 multi-color barcodes for brain and kidney, respectively. We then used a brute force algorithm to find the largest subset of barcodes, (Q), in which every pair had a Hamming distance > 2. Followed by this, we created a graph, (G), in which every possible barcode is a node, and pairs of nodes are connected with edges if their Hamming distance is 1. We then found a maximal independent set (MIS, networkx v2.6.2) that included the nodes in (Q). This method ensures that every pair of barcodes in the MIS have Hamming distance >1. Because the algorithm for finding MIS is random, we ran it 20,000 times and selected the largest MIS across the runs. For the brain, the MIS consisted of 159 barcodes, 121 of which were randomly assigned to the genes. For the kidney, the MIS had 269 barcodes. We randomly added 31 additional barcodes and counted the number of edges of the induced subgraph of (G) with the selected nodes. We repeated this selection 20,000 times and proceeded with the run with the lowest edge count. 300 genes were randomly assigned to these barcodes.

Large-scale padlock probe production

A step-by-step protocol can be found on protocols.io (dx.doi.org/10.17504/protocols.io.n92ldm3pxl5b/v1) and is illustrated in Supplementary Fig. 3c. Briefly, oligo pools were PCR amplified on a 96-well plate (10pM per reaction) using KAPA SYBR fast and 0.4μM of each amplification primer (pAP1V41U and AP2V4 for kidney, AP1V7U and AP2V7 for brain, Supplementary Data 1, Supplementary Fig. 3c) until plateau. The PCR products were pooled and concentrated with ethanol precipitation and further purified using QIAquick PCR purification kit (Qiagen 28106).

For the brain probe set, the purified amplicons were divided into parallel reactions (about 5ug each) and were digested with Lambda Exonuclease (0.5U/ul) in 1x buffer (NEB M0262L) at 37 °C for 2 h and purified using Zymo ssDNA/RNA clean & concentrator kit following manufacturer’s instructions (Zymo D7011). Next, the single-stranded probes were further digested with 5 units of USER enzyme (NEB M5505L) in 1x DpnII buffer at 37 °C for 3 h. Subsequently, for each reaction we added DpnII guide oligo (Supplementary Data 1) to final concentration of 5uM in 1x DpnII buffer, heated the mix to 94 °C for 2 min, cooled to 37  °C and added 50 units of DpnII in 1x DpnII buffer and incubated for 5 h. Finally, probes were size-selected using a TBE-Urea gel.

For the kidney probe set, DpnII digestion was performed after PCR. In detail, the purified amplicons were divided into parallel reactions (about 5ug each) and were digested with DpnII (1U/ul) in 1x NEBuffer DpnII (NEB R0543L) at 37 °C for 3 h and purified with QIAquick PCR purification kit. The purified products were digested with Lambda Exonuclease (0.5U/ul) in 1x buffer (NEB M0262L) for 2 h and purified with Zymo ssDNA/RNA clean & concentrator kit. Finally, the library was digested with USER (0.0625U/ul, M5505L) in 1x NEBuffer DpnII in parallel reactions (about 2.5ug each) for 6 h at 37 °C followed by 3 h at room temperature and purified with Zymo ssDNA/RNA clean & concentrator kit.

DART-FISH

The overall workflow, including reverse transcription, cDNA crosslinking, padlock probe capture, RCA, rolony crosslinking and image acquisition, is illustrated in Fig. 1. A step-by-step protocol can be found at protocols.io (dx.doi.org/10.17504/protocols.io.e6nvwjxnzlmk/v1).

Reverse transcription and cDNA crosslinking

Tissue sections were fixed in 4% PFA in 1x PBS at 4 °C for 1 h, followed by two 3-minute washes with PBST (1x PBS and 0.1% Tween-20). Then, a series of 50%, 70%, 100%, and 100% ethanol were used to dehydrate the tissue sections at room temperature for 5 min each. Next, tissues were air dried for 5 min and in the meantime silicone isolators (Grace Bio-Labs, 664304) were attached around the tissue sections. Then, the tissue sections were permeabilized with 0.25% Triton X-100 in PBSR (1x PBS, 0.05U/μl Superase In, 0.2U/μl Enzymatics RNase Inhibitor) at room temperature for 10 min, followed by two chilled PBSTR (1x PBS, 0.1% Tween-20, 0.05U/μl Superase In, 0.2U/μl Enzymatics RNase Inhibitor) washes and a water wash. Next, the sections were digested with 0.01% pepsin in 0.1 N HCl (pre-warmed 37 °C for 5 min) at 37 °C for 90 s and washed with chilled PBSTR twice. Afterwards, acrydite-modified dT and N9 primers (Acr_dc7-AF488_dT20 and Acr_dc10-Cy5_N9, Supplementary Data 1) were mixed to a final concentration of 2.5 μM with the reverse-transcription mix (10U/μL SuperScript IV (SSIV) reverse transcriptase, 1x SSIV buffer, 250 μM dNTP, 40 μM aminoallyl-dUTP, 5 mM DTT, 0.05U/ul Superase In and 1U/μL Enzymatics RNase inhibitor). The sections with the mix were incubated at 4 °C for 10 min and then transferred to a humidified 37 °C oven for overnight incubation. After reverse transcription, tissue sections were washed with chilled PBSTR twice and incubated in 0.2 mg/mL Acryloyl-X, SE in 1x PBS at room temperature for 30 min. Then, the tissue sections were washed once with PBSTR, followed by incubation with 4% acrylamide solution (4% acrylamide/bis 37:1, 0.05U/μL Superase-In, and 0.2U/μL RNase inhibitor) at room temperature for 30 min. Subsequently, the acrylamide solution was aspirated and gel polymerization solution (0.16% Ammonium persulfate and 0.2% TEMED in the 4% acrylamide solution) was added. Immediately, the tissues were covered with Gel Slick (Lonza #50640)-treated circular coverslips of 18 mm diameter (Ted Pella, 260369), transferred to an argon-filled chamber at room temperature and incubated for 30 min. After gel formation, the tissue sections were washed with 1x PBS twice and the coverslip was gently removed with a needle. At this point, the cDNA is crosslinked to the polyacrylamide gel.

Padlock probe capture

After cDNA crosslinking in gel, remaining RNA was digested with RNase mix (0.25U/μL RNase H, 2.5% Invitrogen RNase cocktail mix, 1x RNase H buffer) at 37 °C for 1 h followed by two PBST washes, 3 min each. The padlock probe library was mixed with Ampligase buffer. Then, the mixture was heated to 85 °C for 3 min and cooled on ice. Subsequently, the mixture was supplemented with 33.3U/μL Ampligase enzyme such that the final concentration of padlock probe library was 180 nM and 100 nM for the kidney and brain probe set, respectively, in 1x Ampligase buffer. Finally, the samples were incubated with probes at 37 °C for 30 min, and then moved to a 55 °C humidified oven for overnight incubation.

RCA and rolony crosslinking

After padlock probe capture, the tissue sections were washed with 1x PBS three times, 3 min each and hybridized with RCA primer solution (0.5 μM rca_primer, 2x SSC, and 30% formamide) at 37 °C for 1 h. Then, the tissue sections were washed with 2x SSC twice and incubated with Phi29 polymerase solution (0.2 U/μL Phi29 polymerase, 1x Phi29 polymerase buffer, 0.02 mM aminoallyl-dUTP, 1 mg/mL BSA, and 0.25 mM dNTP) at 30  °C in a humidified chamber for 7 h. Afterwards, the tissue sections were washed with 1x PBS twice, 3 min each and the rolonies were crosslinked with 5 mM BS(PEG)9 in 1x PBS at room temperature for 1 h. The crosslinking reaction was terminated with 1M Tris, pH 8.0 solution at room temperature for 30 min. Finally, samples were washed with 1x PBS twice and stored in a 4 °C fridge until image acquisition.

Image acquisition

Human Brain

Human brain tissue sample was stained with 1x TrueBlack in 70% ethanol at room temperature for 2 min to reduce the lipofuscin autofluorescence and washed with 1x PBS three times for 3 min each before imaging. For the anchor round imaging, a mixture of anchor round probes, including DARTFISH_anchor_Cy3, dcProbe10_ATTO647N, and dcProbe7_AF488 probes, were diluted to 500nM in 2x SSC and 30% formamide. Then, the samples were stained with anchor round probes at room temperature for 5 min and washed with 1 mL washing buffer (2x SSC, 10% formamide and 0.1% Tween-20) twice for 2 min each prior to imaging. The samples were immersed in 1 mL imaging buffer (2x SSC and 10% formamide) during imaging. For decoding imaging, each imaging cycle started with incubating samples with stripping buffer (2x SSC, 80% formamide, and 0.1% Tween-20) at room temperature for 5 min, washed with washing buffer twice for 2 min each, stained with a mixture of AlexaFluor488, Cy3, and ATTO647 fluorophore-labeled decoding probes (dcProbe0-AF488, dcProbe0-Cy3, and dcProbe0-ATTO647N as an example for round 1) in 2x SSC and 30% formamide for 10 min, and washed with washing buffer three times for 2 min each. Then, the samples were immersed in 1 mL of imaging buffer while imaging. After the last cycle of decoding imaging, DRAQ5 staining (5 μM, room temperature, 10 min) was performed for nuclei segmentation. Z-stack images were acquired by a resonant-scanning Leica TCS SP8 confocal microscope with 20x oil-immersion objective (NA = 0.75), pinhole size of 1 airy unit, pixel size of 284 nm x 284 nm (zoom=2) with 1024 x 1024 pixels per image, and 2 line averaging with 26 z-stacks (step size 1μm).

Human Kidney

The same fluorescent probes were used as in the human brain imaging in this order: anchor round, decoding rounds 1 to 7, DRAQ5 nuclear staining. All hybridizations were performed with 500nM of each of the fluorescent oligos in 2x SSC and 30% formamide for 15 min. Following hybridization, the unbound probes were washed with 4–5 washes with PBST each 2–3 min. Imaging was performed in PBST on a resonant-scanning Leica SP8 with a 20x oil-immersion objective (NA = 0.75), pinhole size of 2 airy units, pixel size of 366 nm x 366 nm (zoom=1.55) with 1024 x 1024 pixels per image, 3 line averaging, with 24 z-stacks (step size 2.5um). After each imaging round, stripping was performed with 80% formamide in 2x SSC and 0.1% Tween-20, 3 times each 3-5 min, followed by 2 quick washes with PBST to prepare for the next hybridization.

RNAscope

Sample preparation

RNAscope HiPlex 50x probe stocks of human SLC17A7, RELN, CUX2, RORB, CBLN2, FEZF2, GAD2, PVALB, LAMP5, PLP1, AQP4,and APBB1IP with HiPlex12 Reagent Kit v2 (488, 550, 650) Assay (ACD, 324419) were purchased from Advanced Cell Diagnostics (ACD). The 50x probe stocks and RNAscope HiPlex diluent were warmed at 40  °C for 10 min. The pre-warmed 50x probe stocks were pooled and diluted to 1x with pre-warmed RNAscope HiPlex diluent before use. RNAscope experiments were carried out according to the manufacturer’s protocol (document number UM324419) with slight modifications for post-mortem human brain tissue. Briefly, the human brain tissue sections were fixed with 4% PFA in 1x PBS at 4 °C for 1 h and dehydrated with a series of 50%, 70%, 100%, and 100% ethanol at room temperature for 5 min each. Then, silicone isolators of 20 mm in diameter (Grace Bio-Labs, 664304) were applied around the tissue sections and the tissue sections were slightly digested with 5 drops of Protease IV at room temperature for 30 min and washed with 1x PBS for 2 min twice. Subsequently, enough volume of 1x pooled probes was added to cover the tissue sections entirely and the probe hybridization was performed in the 40 °C HybEZ Hybridization System for 2 h. Then, the tissue sections were washed with 1 mL 1x wash buffer at room temperature for 2 min twice. Later, the tissue sections were hybridized with RNAscope HiPlex Amp1, incubated in the 40  °C HybEZ Hybridization System for 30 min, and washed with 1x wash buffer at room temperature for 2 min twice. Afterwards, we followed the same process to hybridize the tissue sections with RNAscope HiPlex Amp2 and RNAscope Hiplex Amp3. Finally, we incubated the tissue section with freshly prepared 5% HiPlex FFPE reagent at room temperature for 30 min and washed the tissue sections with 1 mL 1x wash buffer at room temperature for 2 min twice prior to image acquisition.

Image acquisition

The tissue sections with silicone isolators were mounted on the stage of a Leica SP8 confocal microscope and 4 cycles of imaging were performed to image 12 RNA species. In the first imaging cycle, RNAscope HiPlex Fluoro T1-T3 probes were prewarmed at 40  °C, added to cover the tissue sections entirely, and hybridized with the tissue sections for 5 min thrice. After probe hybridization, the tissue sections were washed with 1 mL 1x wash buffer at room temperature for 2 min twice and immersed in 1 mL 4x SSC buffer. Z-stack images were acquired by Leica TCS SP8 confocal microscope with 63x oil-immersed objective (NA 1.4) and pixel size of 113 nm x 113 nm. Then, the fluorophores were cleaved with freshly prepared 10% cleaving solution (100 μL cleaving solution diluted with 900 μL 4x SSC buffer) at room temperature for 15 min and the tissue sections were washed with 0.5% PBST (1x PBS with 0.5% Tween-20) at room temperature for 2 min twice. The fluorophore cleaving process was repeated once to ensure the fluorophores were removed entirely. This process was repeated 3 more rounds to image RNAscope HiPlex Fluoro T4-T12. An additional “Empty” cycle was performed to image the autofluorescence of the human brain tissue without any probes. After the last imaging cycle, we added 80% formamide in 2x SSC buffer to remove RNAscope probes completely and stained the nuclei with 5 μM DRAQ5 at room temperature for 10 min.

RNAscope data processing

RNAscope data was processed with the DART-FISH pipeline with one modification. The images from the “Empty” cycle were subtracted from all RNAscope images to remove the autofluorescence.

DART-FISH data processing (DF3D)

The DART-FISH datasets were processed by our custom pipeline. The source codes of the pipeline can be found in this Github page (https://github.com/Kiiaan/DF3D). Raw z-stack images with 4 channels (3 fluorescent channels and brightfield) from the microscope were registered to a reference round by affine transformation implemented in SimpleElastix⁸⁵ using the brightfield channel as the anchor. Then, each field of view (FOV) underwent decoding to obtain a list of candidate spots. Spots from all FOVs were pooled and filtered (See Sparse deconvolution (SpD) decoder for more details). To obtain the global position of the rolonies, the FOVs were stitched by applying FIJI’s⁸⁶ Grid/Collection Stitching plugin⁸⁷ (in headless mode) to the registered and maximum-projected brightfield images. Note that the theoretical positions of the FOVs, defined by the microscope, were used as initial positions for stitching.

Cell boundaries were segmented with Cellpose (v2.1.1)^88,89. The “cyto” model in Cellpose was fine tuned on each tissue by manually segmenting a handful of composite images of DRAQ5 (nuclei channel) and N9 cDNA stain (cyto channel) using the package’s graphical user interface.

Sparse deconvolution (SpD) decoding

In DART-FISH, each gene is represented by a barcode that can be read out in (n) rounds of (3)-channel imaging. Each barcode is designed to emit fluorescence (be “on”) in exactly (k) rounds, each time in a single fluorescent channel and stay “off” in other rounds. We concatenate the rounds and channels and represent the barcodes as (3n)-dimensional vectors. In other words, barcode (i) is represented by vector ({x}_{i}) in which (1)’s are placed where “on” signal is expected, and 0’s everywhere else. The codebook matrix (X) (3nx(N)) is then defined as (X=[{{{{{{bf{x}}}}}}}_{1},, {{{{{{bf{x}}}}}}}_{2},…,,{{{{{{bf{x}}}}}}}_{N}]), where (N) is the total number of barcodes. In the same way, for every pixel we concatenate the fluorescent intensity values (scaled between 0 and 1) to create a (3n)-dimensional vector ({{{{{bf{y}}}}}}).

The fluorescence signal at each pixel can be sourced from more than one rolony if the distance between neighboring rolonies is smaller than the optical resolution of the imaging system, or if 3-dimensional stacks are analyzed as maximum-projected 2D images. Nevertheless, because of physical constraints, only a handful of rolonies are expected to be the source of signal to each pixel. In this regard, because of the redundancy in the barcode space, combinations of barcodes in one pixel can be decomposed into their original composing barcodes. We formulated this problem as a regularized linear regression problem where a weighted sum of a few barcodes creates the observed signal intensity, where the vector ({{{{{bf{w}}}}}}={[{w}_{1},,{w}_{2},…,,{w}_{N}]}^{T}) shows the contribution of each barcode (Fig. 1d) with most ({w}_{s}(1le sle N)) elements equal to 0. We initially used lasso to solve this problem ((alpha {^{prime}}=0) in Fig. 1d) to promote the sparsity of ({{{{{bf{w}}}}}}), but later decided to use elastic net with a non-zero value for (alpha {^{prime}}) that is much smaller than (alpha) ((alpha {^{prime}}=alpha /100)) to increase stability. We call the solution to this problem ({widehat{{{{{{bf{w}}}}}}}}_{{{{lasso}}}}). Note that, we constrain the problem to positive weight values (({widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}_{{{{s}}}}}{{{{ge }}}}{{{0}}}) for every (s)). The regression problems are solved for all the foreground pixels (({||}{{{{{bf{y}}}}}}{{{{||}}}}_{{{{2}}}},{{{ > }}},{{{0}}}{{{.}}}{{{25}}})) individually. For every barcode (i), we can construct an image with the estimated weight values as pixels: 0 for background and rejected pixels, and non-zero values from (widehat{{{{{{bf{w}}}}}}}). We call these images weight maps. Figure 1f shows weight maps constructed with ({widehat{{{{{{bf{w}}}}}}}}_{{{{lasso}}}}) which have not been filtered.

With our current barcode space, we can only confidently decompose bi-combinations. Hence, for every instance of the elastic net problem, we applied an elbow filter and accepted the solution only when the top one or two weights were significantly larger than other weights.

In more detail, for every pixel, the weights in ({widehat{{{{{{bf{w}}}}}}}}_{{{{lasso}}}}) are sorted in decreasing order. If the second largest weight is smaller than half of the top weight, then the top weight passes the elbow filter. Otherwise, if the third largest weight is smaller than 30% of the largest weight, the top two weights pass the elbow filter. All the values that do not pass the filter are set to zero. For accepted solutions, we performed an ordinary least square (OLS) regression using the top one or two weights to obtain unbiased weights (({widehat{{{{{{bf{w}}}}}}}}_{{{{OLS}}}})). Supplementary Fig. 2a shows weight maps constructed with ({widehat{{{{{{bf{w}}}}}}}}_{{{{{OLS}}}}}) (OLS maps) after applying a Gaussian smoothing.

Estimating channel-specific coefficients

So far, we have assumed that pixel intensities from different rounds and fluorescent channels all have the same scale and distribution. However, there is usually a variation among rounds and fluorescent channels, with some channel-rounds being brighter than others. To account for this effect, we model the channel-specific variations as a multiplicative factor that connects the weights at each pixel to intensities: ({{{{{bf{y}}}}}}{{{{=}}}}{{{{{bf{c}}}}}},{{odot }},X{{{{{bf{w}}}}}}) where ({{{{{bf{c}}}}}}={[{c}_{1},,{c}_{2},…,,{c}_{3n}]}^{T}) is the channel coefficient vector and ({{{{{boldsymbol{odot }}}}}}) denotes element-wise multiplication. Suppose for a set of pixels ({{{{{{bf{y}}}}}}}^{{{{{(}}}}{{{{1}}}}{{{{)}}}}},{{{{,}}}},{{{{{{bf{y}}}}}}}^{{{{{(}}}}{{{{2}}}}{{{{)}}}}}{{{{,}}}}{{{{.}}}}{{{{.}}}}{{{{.}}}}{{{{,}}}},{{{{{{bf{y}}}}}}}^{{{{{(}}}}{{{{P}}}}{{{{)}}}}}) the true barcode weights ({{{{{{bf{w}}}}}}}^{{{{{{boldsymbol{(}}}}}}{{{{{boldsymbol{1}}}}}}{{{{{boldsymbol{)}}}}}}}{{{{{boldsymbol{,}}}}}},{{{{{{bf{w}}}}}}}^{{{{{{boldsymbol{(}}}}}}{{{{{boldsymbol{2}}}}}}{{{{{boldsymbol{)}}}}}}}{{{{{boldsymbol{,}}}}}}{{{{{boldsymbol{.}}}}}}{{{{{boldsymbol{.}}}}}}{{{{{boldsymbol{.}}}}}}{{{{{boldsymbol{,}}}}}},{{{{{{bf{w}}}}}}}^{{{{{{boldsymbol{(}}}}}}{{{{{boldsymbol{P}}}}}}{{{{{boldsymbol{)}}}}}}},) are given. For pixel (i) and channel (j), we could write: ({y}_{j}^{(i)}={c}_{j}{sum }_{b=1}^{N}{X}_{{jb}}{w}_{b}^{(i)}={c}_{j}{sum }_{b=1}^{N}{({{{{{{bf{x}}}}}}}_{{{{{j}}}}})}_{b}{w}_{b}^{(i)}) where ({({{{{{{bf{x}}}}}}}_{{{{{j}}}}})}_{b}) shows the (b)’s element of the (j)’s barcode. In this case, each ({c}_{j}) can be estimated by solving an OLS problem between ({y}_{j}^{(.)}) and ({sum }_{b=1}^{N}{({{{{{{bf{x}}}}}}}_{{{{{j}}}}})}_{b}{w}_{b}^{(.)}). Conversely, if the channel coefficients are given, we can set up the decoding problem with normalized intensities: (bar{{{{{{bf{y}}}}}}}={{{{{bf{y}}}}}}{{{{{boldsymbol{/}}}}}}{{{{{bf{c}}}}}}=X{{{{{bf{w}}}}}}) with (/) being element-wise division. We estimate the channel coefficients in an iterative manner following the algorithm below:

1. 1.

   Initialize ({{{{{bf{c}}}}}}{{{{=}}}}{{{{{bf{1}}}}}}) (no channel variation)

2. 2.

   Take a random sample of foreground pixels

3. 3.

   Normalize the pixel intensities in the sample with ({{{{{bf{c}}}}}})

4. 4.

   Run SpD on the normalized pixels

5. 5.

   Keep pixels with one dominant unsaturated weight (weight in range 0.1 and 0.5) and obtain unbiased weights through OLS

6. 6.

   Update the values of ({{{{{bf{c}}}}}}) by solving (3n) OLS problems

7. 7.

   Repeat steps 3–6 ({n}_{{iter}}) times

We do this procedure for 2 iterations and apply the obtained values when decoding all fields of view.

Setting the elastic net regularization parameter

Because of physical constraints, the solution to the deconvolution problem must be sparse, i.e., only a few non-zero weights should explain the observed intensities. The sparsity of the solution is directly controlled by the L1 regularization term, (alpha) (Fig. 1d). For a given pixel ({{{{{bf{y}}}}}}), higher values of (alpha) shrink the estimated weights (({{{{{||}}}}widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}}{{{{|}}}}{{{{{|}}}}}_{{{{{1}}}}}{{{{to }}}}{{{{0}}}})). Conversely, lower values of (alpha) allow more weights to be non-zero and ({{{{{||}}}}widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}}{{{{|}}}}{{{{{|}}}}}_{{{{{1}}}}}) to grow larger. In fact, one can show if the L2 regularization term, (alpha {^{prime}}=0), the largest weight to be undetected for a pixel made purely from one barcode is ({w}_{max }=frac{3n}{k}alpha) ⁹⁰. For instance, given (alpha=0.05) and codebook parameters (n=6,,{k}=3), then ({w}_{max }=0.3). This means that a pixel composed of one barcode needs to have an underlying intensity ( > 0.3) to get a non-zero ({widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}}). In other words, setting (alpha) too strictly will result in dimmer pixels to have ({widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}}{{{{=}}}}{{{{{boldsymbol{0}}}}}}), while setting (alpha) too loosely will result in spurious non-zero values in ({widehat{{{{{{bf{w}}}}}}}}_{{{{{lasso}}}}}) for brighter more complex pixels, potentially not passing the elbow filter and thus ({widehat{{{{{{bf{w}}}}}}}}_{{{{{OLS}}}}}{{{{=}}}}{{{{{boldsymbol{0}}}}}}). To accommodate a wide range of rolony intensities, we choose (alpha) adaptively based on the pixel norm ({||}{{{{{bf{y}}}}}}{{{{{||}}}}}_{{{{{2}}}}}). First, we form a training data from a random subset of foreground pixels indexed by (i). For a given pixel norm (u), we find the alpha that maximizes a weighted sum of ({||}{{widehat{{{{{{bf{w}}}}}}}}^{(i)}}_{{{{{OLS}}}}}|{|}_{1}) giving more weights to training pixels with closer norms to (u) (equation ∗):

where (g(.)) is the Gaussian function. In practice, for the training pixels we solve the sparse decoding problem for every value of (alpha) on a grid from 0.01 to 0.1 with a step size of 0.005, ({{{{{{boldsymbol{alpha }}}}}}}_{{{{{{rm{train}}}}}}}), to obtain estimated weights ({{widehat{{{{{{bf{w}}}}}}}}^{(i)}}_{{{{{OLS}}}}}{{{{{boldsymbol{(}}}}}}alpha {{{{{boldsymbol{)}}}}}}). Then we create a grid of norms ({{{{{{boldsymbol{u}}}}}}}_{{{{{{rm{train}}}}}}}), spanning 0 and 2.8 with 50 steps. For every value of (u) in ({{{{{{boldsymbol{u}}}}}}}_{{{{{{rm{train}}}}}}}), we solve equation ∗ on the ({{{{{{boldsymbol{alpha }}}}}}}_{{{{{{rm{train}}}}}}}) grid. In other words, we create a lookup table connecting values of ({{{{{{boldsymbol{u}}}}}}}_{{{{{{rm{train}}}}}}}) to the best (alpha) in ({{{{{{boldsymbol{alpha }}}}}}}_{{{{{{rm{train}}}}}}}). For new pixels, (alpha) is determined by the closest norm in the lookup table.

Spot calling

To call spots, Gaussian smoothing is applied to individual OLS maps, followed by peak_local_max filter (scikit-image 0.19.3⁹¹) which returns a binary image with 1’s at the local maxima of the smoothed OLS maps. These peaks are then used as markers for watershed segmentation. From each segmented region, the following features are retained to be used in downstream steps: area, centroid, maximum and average intensity. This formed a list of candidate spots from each FOV.

Spot filtering

To control the specificity of the decoding procedure, we augmented the codebook with a number of barcodes (5-10% of the used barcodes) not used in the probe set (empty barcodes). After spot calling, we record the properties (e.g., area, maximum and average intensity) of spots with an empty barcode. Indeed, we see that empty spots tend to be smaller with lower average/maximum weight (Supplementary Fig. 2c and d). On a small fraction of spots from all fields of view, we train a random forest classifier (scikit-learn v1.1.3) with area, maximum and average weights as features to predict empty/non-empty labels (Supplementary Fig. 2e). We applied the classifier to all spots and obtained emptiness probabilities and set a threshold on these probabilities (0.3–0.35).

Spot assignment to cells

The cell boundaries were computed by applying find_boundaries (scikit-image 0.19.3⁹¹) to the segmentation mask. The distances of all spots were calculated to the closest boundary pixel. The distance was set to 0 if a spot was inside a boundary. A spot was assigned to its closest cell if the distance was less than or equal to 11μm in the kidney, 3μm for non-MBP and 0μm for MBP spots in the brain.

Cell annotation

We used anndata^92,93,94 (v0.8.0) and scanpy⁹²(v1.9.1) to handle and analyze the data. The data normalization was performed using analytic Pearson residuals⁹⁵ (clipped at 40) with a lower bound placed on gene-level standard deviations⁹⁶. Clustering was done with the Leiden algorithm⁹⁷ implemented in scanpy.

Annotating the Brain data set

Cells with counts less than 5 and more than 300 were removed (2980 out of 26348). The top 100 highly variable genes (scanpy.experimental.pp.highly_variable_gene(., flavor=’pearson_residuals’)) were used for normalization, embedding and annotations. PCA was performed on pearson residuals, and the neighborhood graph was created with this command scanpy.pp.neighbors(., n_neighbors = 20, n_pcs = 15, metric=’cosine’). Single-nucleus RNA-seq reference from Jorstad et al.⁶¹ was subsetted to M1C cells and normalized in the same way as DART-FISH. Pax6 and Scng subclasses were removed since we did not design our probe set to target those. Average normalized counts (centroids) were computed for every other subclass in the “within_area_subclass” slot and all clusters of DART-FISH. To annotate the DART-FISH clusters at the class level (excitatory, inhibitory, non-neuronal), we first correlated each cluster to all single-nucleus subclasses, and assigned that cluster to the class of the most highly correlated subclass. Annotation of each class was done separately.

For excitatory neurons, all DART-FISH cells that had a class label of “excitatory” and had at least 20 transcripts were kept (5957 cells). We realized that the Leiden clustering was unstable and by mere shuffling of the order of cells, we would obtain very different clusters. We reasoned that by removing some cells that tend to move between clusters, we could get more stable clusters and have more confidence in their annotation. To find cells that don’t stably cluster, we ran clustering 20 times, every time shuffling the order of the cells. For every cell, we calculated the number of times it was co-clustered with every other cell and took the average of the non-zero values as the co-clustering index (CCI). A perfect CCI of 20 means that the cell is clustered with the same partners in every clustering instance, while lower values show deviations from this limit. We removed the cells with a CCI smaller than 6 and repeated this filtering procedure for three more iterations. The final results show a more stable clustering of the remaining 5101 cells. We then constructed a new neighborhood graph using newly computed principal components (n_neighbors=10, n_pcs=15), followed by Leiden clustering. The cluster centroids were calculated and correlated to the reference subclass centroids. We assigned clusters to their maximally correlated reference subclass if we could also see differential expression of their marker genes (scanpy’s rank_genes_groups), otherwise we labeled them as NA. Of note, the DART-FISH population labeled as L6b/CT was highly correlated with reference subclasses L6b and L6 CT (Supplementary Fig. 5b) and showed expression of marker genes from both subclasses.

For inhibitory neurons and non-neuronal cells, the clustering was more stable to begin with, and we started by constructing the neighborhood matrix (For inhibitory neurons: n_neighbors=20, n_pcs=10. For non-neuronal cells: n_neighbors=25, n_pcs=15) and clustering. Then clusters were assigned to the reference subclass with maximum Pearson’s correlation if the marker genes matched, or otherwise were labeled as NA.

Drawing cortical layer boundaries

Cortical layer boundaries were automatically drawn via Support Vector Machine (SVM) decision boundaries. The Scikit-learn python package (v1.1.3) was used to train a SVM on the following excitatory neuron subtype labels: “L2/3 IT”, “L4 IT”, “L5 IT”, “L6 IT”, “L6b/CT”. First, cells with fewer than 10 total gene counts were filtered out. The x and y coordinates of the cells are standardized via the StandardScaler() function, and the data was fed into a SVM with a radial basis function (RBF) kernel with balanced class weights and one vs. one decision function. The RBF SVM model is then applied to a meshgrid with a fine step size with the same geometric size as the original tissue image. The trained SVM classified the cell type label of each point on the meshgrid to define borders between the cortical layers specified by the excitatory neuron subclasses. We drew contours based on the borders between the various subclasses, and manually superimposed them onto Fig. 3c.

Gene concordance analysis

The RNA portion of the SNARE-seq2 (snare) dataset from Bakken et al.⁴⁶ and Plongthongkum et al.⁵⁰ was used in this section. First, the snare data was subsetted to the DART-FISH genes. Then, DART-FISH and snare data were both normalized (scanpy.pp.normalize_total(., target_sum = 1000)) followed by log-normalization (scanpy.pp.log1p(.)). The average normalized gene expression was calculated for all subclasses. For each gene, the concordance was defined as the Pearson’s correlation between the average expressions across the subclasses between the DART-FISH and snare data (top panel of Supplementary Fig. 5c). The same analysis was performed for a MERFISH data set from Fang et al.⁵⁹ (sample H18.06.006.MTG.250.expand.rep1) with the following details: the subclass labels from metadata column “cluster_L2” were renamed to be consistent with DART-FISH annotations. In particular, subclasses L6b and L6 CT were merged, and subclass L5 ET was removed. Note that subclasses Sst Chodl, Chandelier and Lamp5 Lhx6 were not annotated in the MERFISH dataset and were removed from the DART-FISH analysis for consistency. The rest of the analysis was carried out with 242 shared genes between the datasets (bottom panel of Supplementary Fig. 5c).

Annotating the kidney data set

Cells with less than 5 and more than 100 transcripts were filtered (2024 out of 65565). The top 250 highly variable genes were kept for downstream analyses (scanpy.experimental.pp.highly_variable_gene(., flavor=’pearson_residuals’)). PCA was performed on pearson residuals, and the neighborhood graph was constructed using the command scanpy.pp.neighbors(., n_neighbors = 20, n_pcs = 20, metric=’cosine’) followed by Leiden clustering (l1 clustering). From the kidney reference atlas⁷³, degenerative, cycling, transitioning and medullary cell types were removed. The counts were transformed to pearson residuals and the remaining subclass level 1 and level 2 centroids were calculated. We then calculated the Pearson correlations between subclass level 1 centroids and cluster centroids and assigned each l1 cluster to the subclass level 1 with maximum correlation. We then subclustered each of the l1 clusters and assigned those to subclass level 2 identities with maximum correlation, only if the relevant marker genes were expressed. Through this procedure we could not resolve PT-S1 and PT-S2 subtypes separately; thus, we labeled the clusters that were highly correlated with these populations as PT-S1/S2. Similarly, for immune cells, this procedure could confidently resolve MAC-M2 cells and the general myeloid (IMM_Myl) and lymphoid (IMM_Lym) populations. To annotate the immune cells at higher level of granularity, we updated their subclass level 2 labels with the following strategy: Each DART-FISH cell with subclass level 1 label “IMM” was separately correlated with the following immune subtypes in the reference atlas: B, PL, T, MAC-M2, MDC, cDC. The immune subtypes with highest and 2nd highest correlation were kept. If the highest correlation was larger than 0.4 and the ratio of the highest to the 2nd highest correlation was larger than 1.25, the label was updated to that of the highest correlated subtype, otherwise it remained unchanged.

Cell-cell interaction analysis

We used squidpy.gr.co_occurrence function (v1.2.4.dev27+gb644428) with n_splits = 1 and an interval between 7μm and 110μm⁷⁷.

Comparison of decoding methods

Datasets of varying levels of complexity were simulated to compare SpD with StarFish⁴² (pixel-based naive matching), BarDensr⁴⁴ and ISTDECO⁴³ (deconvolution-based methods). The synthetic datasets were constructed using the human brain codebook (3-on-3-off, 121 genes with 10 empty barcodes) with equal abundance of all genes and uniform spatial distribution of spots. The rolonies were modeled as gaussian spots with peak intensity randomly chosen to be between 0.25 and 0.7 and sigma between 2 and 2.5 pixels. To model channel-specific intensity variation, we randomly drew 18 channel-specific coefficients from a uniform distribution between 0.75 and 1.25 to scale their respective images, while clipping the intensity values above 1. We simulated multiple datasets varying the number of spots between ({5 * 10}^{3}) to (4 * {10}^{5}) spots in a field of view of size 1024 x 1024 pixels. Different decoding methods were applied to the synthetic datasets with default settings to the extent possible, with no post-hoc filtering of the spots. The only exception was StarFish for which the distance threshold was set to 0.7 as a fair balance between specificity and sensitivity. Then, the groundtruth spots were matched one-to-one to the decoded spots if the barcodes were identical and the centroids were closer than 6 pixels. Sensitivity is defined as the fraction of groundtruth spots matched with a decoded spot. Specificity is defined as the fraction of matched decoded spots over all decoded spots. Empty rate is the fraction of empty barcodes among all decoded barcodes and is inversely related to specificity.

Reporting summary

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

• SEO Powered Content & PR Distribution. Get Amplified Today.
• PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
• PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
• PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
• PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
• Source: https://www.nature.com/articles/s41467-024-46437-y