Research Article

Spatially resolved, highly multiplexed RNA profiling in single cells

See allHide authors and affiliations

Science  24 Apr 2015:
Vol. 348, Issue 6233, aaa6090
DOI: 10.1126/science.aaa6090

Multiplexed RNA imaging in single cells

The basis of cellular function is where and when proteins are expressed and in what quantities. Single-molecule fluorescence in situ hybridization (smFISH) experiments quantify the copy number and location of mRNA molecules; however, the numbers of RNA species that can be simultaneously measured by smFISH has been limited. Using combinatorial labeling with error-robust encoding schemes, Chen et al. simultaneously imaged 100 to 1000 RNA species in a single cell. Such large-scale detection allows regulatory interactions to be analyzed at the transcriptome scale.

Science, this issue p. 10.1126/science.aaa6090

Structured Abstract

INTRODUCTION

The copy number and intracellular localization of RNA are important regulators of gene expression. Measurement of these properties at the transcriptome scale in single cells will give answers to many questions related to gene expression and regulation. Single-molecule RNA imaging approaches, such as single-molecule fluorescence in situ hybridization (smFISH), are powerful tools for counting and mapping RNA; however, the number of RNA species that can be simultaneously imaged in individual cells has been limited. This makes it challenging to perform transcriptomic analysis of single cells in a spatially resolved manner. Here, we report multiplexed error-robust FISH (MERFISH), a single-molecule imaging method that allows thousands of RNA species to be imaged in single cells by using combinatorial FISH labeling with encoding schemes capable of detecting and/or correcting errors.

RATIONALE

We labeled each cellular RNA with a set of encoding probes, which contain targeting sequences that bind the RNA and readout sequences that bind fluorescently labeled readout probes. Each RNA species is encoded with a particular combination of readout sequences. We used successive rounds of hybridization and imaging, each with a different readout probe, to identify the readout sequences bound to each RNA and to decode the RNA. In principle, combinatorial labeling allows the number of detectable RNA species to grow exponentially with the number of imaging rounds, but the detection errors also increase exponentially. To combat such accumulating errors, we exploited error-robust encoding schemes used in digital electronics, such as the extended Hamming code, in the design of our encoding probes but modified these schemes in order to account for the error properties in FISH measurements. We assigned each RNA a binary word in our modified Hamming code and encoded the RNA with a combination of readout sequences according to this binary word.

RESULTS

We first imaged 140 RNA species in human fibroblast cells using MERFISH with 16 rounds of hybridization and a modified Hamming code capable of both error detection and correction. We obtained ~80% detection efficiency and observed excellent correlation of RNA copy numbers determined with MERFISH with both bulk RNA sequencing data and conventional smFISH measurements of individual genes.

Next, we used an alternative MERFISH encoding scheme, which is capable of detecting but not correcting errors, to image 1001 RNA species in individual cells using only 14 rounds of hybridization. The observed RNA copy numbers again correlate well with bulk sequencing data. However, the detection efficiency is only one-third that of the error-correcting encoding scheme.

We performed correlation analysis of the 104 to 106 pairs of measured genes and identified many covarying gene groups that share common regulatory elements. Such grouping allowed us to hypothesize potential functions of ~100 unannotated or partially annotated genes of unknown functions. We further analyzed correlations in the spatial distributions of different RNA species and identified groups of RNAs with different distribution patterns in the cell.

DISCUSSION

This highly multiplexed imaging approach enables analyses based on the variation and correlation of copy numbers and spatial distributions of a large number of RNA species within single cells. Such analyses should facilitate the delineation of regulatory networks and in situ identification of cell types. We envision that this approach will allow spatially resolved transcriptomes to be determined for single cells.

MERFISH for transcriptome imaging.

Numerous RNA species can be identified, counted, and localized in a single cell by using MERFISH, a single-molecule imaging approach that uses combinatorial labeling and sequential imaging with encoding schemes capable of detection and/or correction of errors. This highly multiplexed measurement of individual RNAs can be used to compute the gene expression profile and noise, covariation in expression among different genes, and spatial distribution of RNAs within single cells.

Abstract

Knowledge of the expression profile and spatial landscape of the transcriptome in individual cells is essential for understanding the rich repertoire of cellular behaviors. Here, we report multiplexed error-robust fluorescence in situ hybridization (MERFISH), a single-molecule imaging approach that allows the copy numbers and spatial localizations of thousands of RNA species to be determined in single cells. Using error-robust encoding schemes to combat single-molecule labeling and detection errors, we demonstrated the imaging of 100 to 1000 distinct RNA species in hundreds of individual cells. Correlation analysis of the ~104 to 106 pairs of genes allowed us to constrain gene regulatory networks, predict novel functions for many unannotated genes, and identify distinct spatial distribution patterns of RNAs that correlate with properties of the encoded proteins.

System-wide analyses of the abundance and spatial organization of RNAs in single cells promise to transform our understanding in many areas of cell and developmental biology, such as the mechanisms of gene regulation, the heterogeneous behavior of cells, and the development and maintenance of cell fate (1). Single-molecule fluorescence in situ hybridization (smFISH) has emerged as a powerful tool for studying the copy number and spatial organization of RNAs in single cells either in isolation or in their native tissue context (2, 3). Taking advantage of its ability to map the spatial distributions of specific RNAs with high resolution, smFISH has revealed the importance of subcellular RNA localization in diverse processes such as cell migration, development, and polarization (48). In parallel, the ability of smFISH to precisely measure the copy numbers of specific RNAs without amplification bias has allowed quantitative measurement of the natural fluctuations in gene expression, which has in turn elucidated the regulatory mechanisms that shape such fluctuations and their role in a variety of biological processes (913).

Recent advances in imaging and analysis methods have allowed hundreds of smFISH measurements to be performed in an automated manner, substantially expanding our knowledge of the RNA expression profile and spatial organization in different organisms (14, 15). However, application of the smFISH approach to many systems-level questions remains limited by the number of RNA species that can be simultaneously measured in single cells. State-of-the-art efforts by using combinatorial labeling with either color-based barcodes or sequential hybridization have enabled simultaneous measurements of 10 to 30 different RNA species in individual cells (1619), yet many interesting biological questions would benefit from the measurement of hundreds to thousands of RNAs within a single cell. For example, analysis of how the expression profile of such a large number of RNAs vary from cell to cell and how these variations correlate among different genes could be used to systematically identify coregulated genes and map regulatory networks, knowledge of the subcellular organizations of numerous RNAs and their correlations could help elucidate molecular mechanisms underlying the establishment and maintenance of many local cellular structures, and RNA profiling of individual cells in native tissues could allow in situ identification of cell type.

Here, we report multiplexed error-robust FISH (MERFISH), a highly multiplexed smFISH imaging method that substantially increases the number of RNA species that can be simultaneously imaged in single cells by using combinatorial labeling and sequential imaging with error-robust encoding schemes. We demonstrated this transcriptome imaging approach by simultaneously measuring 140 RNA species with an encoding scheme that can both detect and correct errors and 1001 RNA species with an encoding scheme that can detect but not correct errors. Correlation analyses of the copy number variations and spatial distributions of these genes allowed us to identify groups of genes that are coregulated and groups of genes that share similar spatial distribution patterns inside the cell.

Combinatorial labeling with error-robust encoding schemes

Combinatorial labeling that identifies each RNA species by multiple (N) distinct signals offers a route to rapidly increase the number of RNA species that can be probed simultaneously in individual cells (Fig. 1A). However, this approach to scaling up the throughput of smFISH to the systems scale faces a substantial challenge because not only does the number of addressable RNA species increases exponentially with N, but the detection error rates also grow exponentially with N (Fig. 1, B to D). Imagine a conceptually simple scheme to implement combinatorial labeling, in which each RNA species is encoded with a N-bit binary word, and the sample is probed with N corresponding rounds of hybridization, each round targeting only the subset of RNAs that should read “1” in the corresponding bit (fig. S1). N rounds of hybridization would allow 2N – 1 RNA species to be probed. With just 16 hybridizations, more than 64,000 RNA species—which should cover the entire human transcriptome, including both messenger RNAs (mRNAs) and noncoding RNAs (20)—could be identified (Fig. 1B, black symbols). However, as N increases, the fraction of RNAs properly detected (the calling rate) would rapidly decrease and, more troublingly, the fraction of RNAs that are identified as incorrect species (the misidentification rate) would rapidly increase (Fig. 1, C and D, black symbols). With realistic error rates per hybridization (measured below), the majority of RNA molecules would be misidentified after 16 rounds of hybridizations.

Fig. 1 MERFISH: A highly multiplexed smFISH approach enabled by combinatorial labeling and error-robust encoding.

(A) Schematic depiction of the identification of multiple RNA species in N rounds of imaging. Each RNA species is encoded with a N-bit binary word, and during each round of imaging, only the subset of RNAs that should read 1 in the corresponding bit emit signal. (B to D) The number of addressable RNA species (B); the rate at which these RNAs are properly identified—the “calling rate” (C); and the rate at which RNAs are incorrectly identified as a different RNA species—the “misidentification rate” (D); plotted as a function of the number of bits (N) in the binary words encoding RNA. Black indicates a simple binary code that includes all 2N-1 possible binary words. Blue indicates the HD4 code in which the Hamming distance separating words is 4. Purple indicates a modified HD4 (MHD4) code where the number of 1 bits are kept at four. The calling and misidentification rates are calculated with per-bit error rates of 10% for the 1→0 error and 4% for the 0→1 error. (E) Schematic diagram of the implementation of a MHD4 code for RNA identification. Each RNA species is first labeled with ~192 encoding probes that convert the RNA into a specific combination of readout sequences (Encoding hyb). These encoding probes each contain a central RNA-targeting region flanked by two readout sequences, drawn from a pool of N different sequences, each associated with a specific hybridization round. Encoding probes for a specific RNA species contain a particular combination of four of the N readout sequences, which correspond to the four hybridization rounds in which this RNA should read 1. N subsequent rounds of hybridization with the fluorescent readout probes are used to probe the readout sequences (hyb 1, hyb 2, …, hyb N). The bound probes are inactivated by photobleaching between successive rounds of hybridization. For clarity, only one possible pairing of the readout sequences is depicted for the encoding probes; however, all possible pairs of the four readout sequences are used at the same frequency and distributed randomly along each cellular RNA in the actual experiments.

To address this challenge, we designed error-robust encoding schemes in which only a subset of the 2N – 1 words separated by a certain Hamming distance (21) were used to encode RNAs. In a codebook in which the minimum Hamming distance is 4 (HD4 code), at least four bits must be read incorrectly to change one valid code word into another (fig. S2A). As a result, every single-bit error produces a word that is exclusively close to a single RNA-encoding word, allowing such errors to be detected and corrected (fig. S2B). Double-bit errors produce words with an equal Hamming distance of 2 from multiple valid code words and, thus, can be detected but not corrected (fig. S2C). Such a code should substantially increase the calling rate and reduce the misidentification rate (Fig. 1, C and D, blue symbols). To further account for the fact that it is more likely to miss a hybridization event (an 1→0 error) than to misidentify a background spot as an RNA (an 0→1 error) in smFISH measurements, we designed a modified HD4 (MHD4) code, in which the number of 1 bits were kept both constant and relatively low—only four per word in this work—so as to reduce error and avoid biased detection. This MHD4 code should further increase the calling rate and reduce the misidentification rate (Fig. 1, C and D, purple symbols).

In addition to the error considerations, several practical challenges have also made it difficult to probe a large number of RNA species, such as the high cost of the massive number of distinct FISH probes needed and the long time required to complete many rounds of hybridization. An oligopaint approach has been previously developed to generate a large number of oligonucleotide probes to label chromosome DNA and to introduce nontargeting sites for secondary activities (22). Inspired by this approach, we designed a two-step labeling scheme to encode and read out cellular RNAs (Fig. 1E). First, we label cellular RNAs with a set of encoding probes, each probe comprising a RNA targeting sequence and two flanking readout sequences. Four of the N distinct readout sequences were assigned to each RNA species based on the N-bit MHD4 code word of the RNA. Second, we identified these N readout sequences with complementary FISH probes (the readout probes) via N rounds of hybridization and imaging, each round using a different readout probe. To increase the signal-to-background ratio, we labeled every cellular RNA with ~192 encoding probes. Because each encoding probe contained two of the four readout sequences associated with that RNA (Fig. 1E), a maximum of ~96 readout probes can bind to each cellular RNA per hybridization round. To generate the massive number of encoding probes required, we amplified them from array-derived oligonucleotide pools containing tens of thousands of custom sequences using a modified form of the oligopaint protocol comprising in vitro transcription followed by reverse transcription (fig. S3 and supplementary materials, materials and methods, “Probe Synthesis”) (22, 23). This two-step labeling approach dramatically diminished the total hybridization time for an experiment; we found that efficient hybridization to the readout sequences took only 15 min, whereas efficient direct hybridization to cellular RNA required more than 10 hours.

Measuring 140 genes with MERFISH by use of a 16-bit MHD4 code

To test the feasibility of this error-robust, multiplexed imaging approach, we performed a 140-gene measurement on human fibroblast cells (IMR90) using a 16-bit MHD4 code to encode 130 RNA species while leaving 10 code words as misidentification controls (table S1). After each round of hybridization with the fluorescent readout probes, cells were imaged by means of conventional wide-field imaging with an oblique-incidence illumination geometry. Fluorescent spots corresponding to individual RNAs were clearly detected and were then efficiently extinguished via a brief photobleaching step (Fig. 2A). The sample was stable throughout the 16 rounds of iterative labeling and imaging: The change in the number of fluorescent spots from round to round matched the expected change predicted on the basis of the relative abundances of RNA species targeted in each round derived from bulk sequencing, and we did not observe a systematic decreasing trend with increasing number of hybridization rounds (fig. S4A). The average brightness of the spots varied from round to round with a standard deviation of 40%, which is likely due to different binding efficiencies of the readout probes to the different readout sequences on the encoding probes (fig. S4B). We observed only a small, systematic decreasing trend in the spot brightness with increasing hybridization rounds, which was on average 4% per round (fig. S4B).

Fig. 2 Simultaneous measurement of 140 RNA species in single cells by use of MERFISH with a 16-bit MHD4 code.

(A) Images of RNA molecules in an IMR90 cell after each hybridization round (hyb 1 to hyb 16). The images after photobleaching (for example, bleach 1) demonstrate efficient removal of fluorescent signals between hybridizations. (B) The localizations of all detected single molecules in this cell colored according to their measured binary words. (Inset) The composite, false-colored fluorescent image of the 16 hybridization rounds for the boxed subregion with numbered circles indicating potential RNA molecules. A red circle indicates an unidentifiable molecule, the binary word of which does not match any of the 16-bit MHD4 code words even after error correction. (C) Fluorescent images from each round of hybridization for the boxed subregion in (B), with circles indicating potential RNA molecules. (D) Corresponding words for the spots identified in (C). Red crosses represent the corrected bits. (E) The RNA copy number for each gene observed without (green) or with (blue) error correction in this cell. (F) The confidence ratio measured for the 130 RNA species (blue) and the 10 misidentification control words (red) normalized to the maximum value observed from the misidentification controls (dashed line). (G) Scatter plot of the average copy number of each RNA species per cell measured with two shuffled codebooks of the MHD4 code. The Pearson correlation coefficient is 0.94 with a P value of 1 × 10−53. The dashed line corresponds to the y = x line. (H) Scatter plot of the average copy number of each RNA species per cell versus the abundance determined by bulk sequencing in FPKM. The Pearson correlation coefficient between the logarithmic abundances of the two measurements was 0.89 with a P value of 3 × 10−39.

We then constructed binary words from the observed fluorescent spots based on their on-off patterns across the 16 hybridization rounds (Fig. 2, B to D). If the word exactly matched one of the 140 MHD4 code words (exact matches) or differed by only one bit (error-correctable matches), we assigned it to the corresponding RNA species (Fig. 2D). Within the single cell depicted in Fig. 2, A and B, more than 1500 RNA molecules corresponding to 87% of the 130 encoded RNA species were detected after error correction (Fig. 2E). Similar observations were made in ~400 cells from seven independent experiments. On average, ~4 times as many RNA molecules and ~2 times as many RNA species were detected per cell after error correction as compared with the values obtained before error correction (fig. S5).

Two types of errors can occur in the copy number measurement of each RNA species: (i) Some molecules of this RNA species are not detected, leading to a drop in calling rate, and (ii) some molecules from other RNA species are misidentified as this RNA species. To assess the extent of misidentification, we used the 10 misidentification control words—code words that were not associated with any cellular RNA. Although matches to these control words were observed, they occurred far less frequently than did the real RNA-encoding words: 95% of the 130 RNA-encoding words were counted more frequently than the median count for these control words. Moreover, we typically found the ratio of the number of exact matches to the number of matches with one-bit errors for a real RNA-encoding word to be substantially higher than the same ratios observed for the misidentification controls, as expected (fig. S6, A and B). Using this ratio as a measure of the confidence in RNA identification, we found that 91% of the 130 RNA species had a confidence ratio greater than the maximum confidence ratio observed for the misidentification controls (Fig. 2F), demonstrating a high accuracy of RNA identification. Subsequent analyses were conducted only on these 91% of genes.

To estimate the calling rate, we used the error-correction ability of the MHD4 code to determine the 1→0 error rates (10% on average) and 0→1 error rates (4% on average) for each hybridization round (fig. S6, C and D). Using these error rates, we estimated an ~80% calling rate for individual RNA species after error correction—~80% of the fluorescent spots corresponding to a RNA species were decoded correctly (fig. S6E). Although the remaining 20% of spots contributed to a loss in detection efficiency, most of them did not cause species misidentification because they were decoded as double-bit error words and discarded.

To test for potential technical bias in our measurements, we probed the same 130 RNAs species with a different MHD4 codebook by shuffling the code words among different RNA species (table S1) and changing the encoding probe sequences. Measurements with this alternative code gave similar misidentification and calling rates (fig. S7). The copy numbers of individual RNA species per cell measured with these two codebooks showed excellent agreement with a Pearson correlation coefficient of 0.94 (Fig. 2G), indicating that the choice of encoding scheme did not bias the measured counts.

In order to validate the copy numbers derived from our MERFISH experiments, we performed conventional smFISH measurements on 15 of the 130 genes, spanning the full measured abundance range of three orders of magnitude. For each of these genes, both the average copy number and the copy number distribution across many cells agreed quantitatively between our MERFISH and conventional smFISH measurements (fig. S8, A and B). The ratio of the copy numbers determined by these two approaches was 0.82 ± 0.06 (mean ± SEM across the 15 measured RNA species) (fig. S8B), which agrees with the estimated 80% calling rate for our multiplexed imaging approach. The quantitative match between this ratio and our estimated calling rate over the full measured abundance range additionally supports our assessment that the misidentification error was low. Given that the agreement between the MERFISH and conventional smFISH results extended to the genes at the lowest measured abundance (<1 copy per cell) (fig. S8B), we estimate that our measurement sensitivity was better than 1 copy per cell.

As a final validation, we compared the abundance of each RNA species averaged over hundreds of cells to those obtained from a bulk RNA sequencing measurement that we performed on the same cell line. Our imaging results correlated remarkably well with bulk sequencing results, with a Pearson correlation coefficient of 0.89 (Fig. 2H).

High-throughput analysis of cell-to-cell variation in gene expression

The MERFISH approach allows parallelization of measurements of many individual RNA species and covariation analysis between different RNA species. We first illustrated the parallelization aspect by examining the cell-to-cell variation in the expression level of each of the measured genes (Fig. 3A). To quantify the measured variation, we computed the Fano factors, defined as the ratio of the variance to the mean RNA copy number, for all measured RNA species. The Fano factors substantially deviated from 1, the value expected for a simple Poisson process, for many genes and exhibited an increasing trend with the mean RNA abundance (Fig. 3B), which is consistent with a previous observation for other cell types (24). A simple model for promoter regulation—the promoter stochastically switches between on and off states with global constraints on the kinetic rates—has been previously suggested to rationalize such a trend (24, 25). According to this model, this trend of increasing Fano factors with mean RNA abundance can be explained by changes in the transcription rate and/or promoter off-switching rates but not by changes in the promoter on-switching rate.

Fig. 3 Cell-to-cell variations and pairwise correlations for the RNA species determined from the 140-gene measurements.

(A) Comparison of gene expression levels in two individual cells. (B) Fano factors for individual genes. Error bars represent standard error of the mean determined from seven independent data sets. (C) Z-scores of the expression variations of four example pairs of genes showing correlated (top two) or anticorrelated (bottom two) variation for 100 randomly selected cells. Z-score is defined as the difference from the mean normalized by the standard deviation. (D) Matrix of the pairwise correlation coefficients of the cell-to-cell variation in expression for the measured genes, shown together with the hierarchical clustering tree. The seven groups identified by a specific threshold on the cluster tree (dashed line) are indicated by the black boxes in the matrix and colored lines on the tree, with gray lines on the tree indicating ungrouped genes. Different threshold choices on the cluster tree could be made to select either smaller subgroups with tighter correlations or larger super-groups containing more weakly coupled subgroups. Two of the seven groups are enlarged on the right. (E) Enrichment of 30 selected, statistically significantly enriched GO terms in the seven groups. Enrichment refers to the ratio of the fraction of genes within a group that have the specific GO term to the fraction of all measured genes having that term. Top 10 statistically significantly enriched GO terms for each of the seven groups are shown in table S2. Not all of the GO terms presented here are in the top 10 list.

Moreover, we identified several RNA species with substantially larger Fano factors than this average trend. For example, we found that SLC5A3, CENPF, MKI67, TNC, and KIAA1199 displayed Fano factor values substantially higher than those of the other genes expressed at similar abundance levels. The high variability of some of these genes can be explained by their association with the cell cycle. For example, two of these particularly “noisy” genes, MKI67 and CENPF, are both annotated as cell-cycle related genes (26), and based on their bimodal expression (Fig. 3C), we propose that their transcription is strongly regulated by the cell cycle. Other high-variability genes did not show the same bimodal expression patterns and are not known to be associated with the cell cycle. Understanding the origin and implications of noisy gene expression is an active topic of current research (24).

Analysis of expression covariation among different genes

Analysis of covariations in the expression levels of different genes can reveal which genes are coregulated and elucidate gene regulatory pathways. At the population level, such analysis often requires the application of external stimuli to drive gene expression variation; hence, correlated expression changes can be observed among genes that share common regulatory elements influenced by the stimuli (27). At the single-cell level, one can take advantage of the natural stochastic fluctuations in gene expression for such analysis and can thus study multiple regulatory networks without having to stimulate each of them individually. Such covariation analysis can constrain regulatory networks, suggest new regulatory pathways, and predict function for unannotated genes based on associations with covarying genes (11, 28).

We applied this approach to the 140-gene measurements and examined the ~10,000 pairwise correlation coefficients that describe how the expression levels of each pair of genes covaried from cell to cell. Many of the highly variable genes showed tightly correlated or anticorrelated variations (Fig. 3C). To better understand the correlations for all gene pairs, we adopted a hierarchical clustering approach, commonly used in the analysis of both bulk and single-cell expression data (29, 30), to organize these genes on the basis of their correlation coefficients (Fig. 3D). From the cluster tree structure, we identified seven groups of genes with substantially correlated expression patterns (Fig. 3D and table S2). Within each of the seven groups, every gene showed significantly stronger average correlation with other members of the group than with genes outside the group (table S2). To further validate and understand these groups, we identified gene ontology (GO) terms (31) enriched in each of these seven groups. The enriched GO terms within each group shared similar functions and were largely specific to each group (Fig. 3E and table S2), validating the notion that the observed covariation in expression reflects some commonalities in the regulation of these genes.

Here, we describe two of these groups as illustrative examples. The predominant GO terms associated with group 1 were terms associated with the extracellular matrix (ECM) (Fig. 3, D and E, and table S2). Notable members of this group included ECM components—such as FBN1, FBN2, COL5A, COL7A, and TNC—and glycoproteins linking the ECM and cell membranes, such as VCAN and THBS1. The group also included an unannotated gene, KIAA1199, which we would predict to play a role in ECM metabolism on the basis of its association with this cluster. Indeed, this gene has recently been identified as an enzyme involved in the regulation of hyaluronan, which is a major sugar component of the ECM (32).

Group 6 contained many genes that encode vesicle transport proteins and proteins associated with cell motility (Fig. 3, D and E, and table S2). The vesicle transport genes included microtubule motors and related genes DYNC1H, CKAP1, and factors associated with vesicle formation and trafficking, such as DNAJC13 and RAB3B. Again, we found an unannotated gene, KIAA1462, within this cluster. On the basis of its strong correlation with DYNC1H1 and DNAJC13, we predict that this gene may be involved in vesicle transport. The cell motility genes in this group included genes encoding actin-binding proteins such as AFAP1, SPTAN1, SPTBN1, and MYH10, and genes involved in the formation of adhesion complexes, such as FLNA and FLNC. Several guanosine triphosphatase (GTPase)–associated factors involved in the regulation of cell motility, attachment, and contraction also fell into this group, including DOCK7, ROCK2, IQGAP1, PRKCA, and AMOTL1. The observation that some cell motility genes correlated with vesicle transport genes is consistent with the role of vesicle transport in cell migration (33). An additional feature of group 6 is that a subset of these genes—in particular, those related to cell motility—were anticorrelated with members of the ECM group discussed above (Fig. 3D). This anticorrelation may reflect regulatory interactions that mediate the switching of cells between adherent and migratory states.

Mapping spatial distributions of RNAs

As an imaging-based approach, MERFISH also allowed us to investigate the spatial distributions of many RNA species simultaneously. Several patterns emerged from the visual inspection of individual genes, with some RNA transcripts enriched in the perinuclear region, some enriched in the cell periphery, and some scattered throughout the cell (Fig. 4A). To identify genes with similar spatial distributions, we determined the correlation coefficients for the spatial density profiles of all pairs of RNA species and organized these RNAs according to the pairwise correlations again using the hierarchical clustering approach. The correlation coefficient matrix showed groups of genes with correlated spatial organizations, and the two most notable groups with the strongest correlations are indicated in Fig. 4B. Group I RNAs appeared enriched in the perinuclear region, whereas group II RNAs appeared enriched near the cell periphery (Fig. 4C). Quantitative analysis of the distances between each RNA molecule and the cell nucleus or the cell periphery indeed confirmed this visual impression (Fig. 4D).

Fig. 4 Distinct spatial distributions of RNAs observed in the 140-gene measurements.

(A) Examples of the spatial distributions observed for four different RNA species in a cell. (B) Matrix of the pairwise correlation coefficients describing the degree with which the spatial distributions of each gene pair is correlated, shown together with the hierarchical clustering tree. Two strongly correlating groups are indicated by the black boxes on the matrix and color on the tree. (C) The spatial distributions of all RNAs in the two groups in two example cells. Light blue symbols, group I genes; red symbols, group II genes. (D) Average distances for genes in group I and genes in group II to the cell edge or the nucleus normalized to the average distances for all genes. Error bars represent SEM across seven data sets. (E) Enrichment of GO terms in each of the two groups.

Group I contained genes encoding extracellular proteins such as FBN1, FBN2, and THSB1; secreted proteins such as PAPPA; and integral membrane proteins such as LRP1 and GPR107. These proteins have no obvious commonalities in function. Rather, a GO analysis showed significant enrichment for location terms, such as extracellular region, basement membrane, or perivitelline space (Fig. 4E). To reach these locations, proteins must pass through the secretion pathway, which often requires translation of mRNA at the endoplasmic reticulum (ER) (34, 35). Thus, we propose that the spatial pattern that we observed for these mRNAs reflects their cotranslational enrichment at the ER. The enrichment of these mRNAs in the perinuclear region (Fig. 4, C and D, light blue), where the rough ER resides, supports this conclusion.

Group II contained genes encoding the actin-binding proteins, including filamins FLNA and FLNC, talin TLN1, and spectrins SPTAN1 and SPTBN1; the microtubule-binding protein CKAP5; and the motor proteins MYH10 and DYNC1H1. This group was enriched with GO terms such as cortical actin cytoskeleton, actin filament binding, and cell-cell adherens junction (Fig. 4E). It has been shown previously that β-actin mRNA is enriched near the cell periphery in fibroblasts, as are mRNAs that encode members of the actin-binding Arp2/3 complex (36, 37). The enrichment of group II mRNAs in the peripheral region of the cells (Fig. 4, C and D) suggests that the spatial distribution of the group II genes might be related to the distribution of actin cytoskeleton mRNAs.

Measuring 1001 genes with MERFISH by use of a 14-bit MHD2 code

Last, we sought to further increase the throughput of our MERFISH measurement by simultaneously imaging ~1000 RNA species. This increase could be achieved with our MHD4 code by increasing the number of bits per code word to 32 while maintaining the number of 1 bits per word at four (Fig. 1B). This could be implemented by either increasing the number of hybridization round to 32 or maintaining 16 rounds of hybridization, but using two-color imaging in each round. We pursued an alternative approach that did not require an increase in the number of hybridizations or color channels by relaxing the error correction requirement but keeping the error-detection capability. For example, by reducing the Hamming distance from 4 to 2, we could use all 14-bit words that contain four 1 bits to encode 1001 genes and probe these RNAs with only 14 rounds of hybridization. However, because a single error can produce a word equally close to two different code words, error correction is no longer possible for this modified Hamming-distance-2 (MHD2) code. Hence, we expect the calling rate to be lower and the misidentification rate to be higher with this encoding scheme.

To evaluate the performance of this 14-bit MHD2 code, we set aside 16 of the 1001 possible code words as misidentification controls and used the remaining 985 words to encode cellular RNAs (table S3). Among these 985 RNAs, we included 107 RNA species probed in the 140-gene experiments as an additional control. We performed the 1001-gene experiments in IMR90 cells by using a similar procedure as described above. To allow all encoding probes to be synthesized from a single 100,000-member oligopool, we reduced the number of encoding probes per RNA species to ~94. Fluorescent spots corresponding to individual RNA molecules were again clearly detected in each round of hybridization with the readout probes, and based on their on-off patterns, these spots were decoded into RNA (Fig. 5A and fig. S9, A and B). In the cell shown in Fig. 5A, 430 RNA species were detected, and similar results were obtained in ~200 imaged cells in three independent experiments.

Fig. 5 Simultaneous measurements of 1001 RNA species in single cells by using MERFISH with a 14-bit MHD2 code.

(A) The localizations of all detected single molecules in a cell colored based on their measured binary words. (Inset) The composite, false-colored fluorescent image of the 14 hybridization rounds for the boxed subregion with numbered circles indicating potential RNA molecules. Red circles indicate unidentifiable molecules, the binary words of which do not match any of the 14-bit MHD2 code words. Images of individual hybridization round are shown in fig. S9A. (B) Scatter plot of the average copy number per cell measured in the 1001-gene experiments versus the abundance measured via bulk sequencing. The black symbols are for the 73% of genes detected with confidence ratios higher than the maximum ratio observed for the misidentification controls. The Pearson correlation coefficient is 0.76 with a P value of 3 × 10−133. The red symbols are for the remaining 27% of genes. The Pearson correlation coefficient is 0.65 with a P value of 3 × 10−33. (C) Scatter plot of the average copy number for the 107 genes shared in both the 1001-gene measurement with the MHD2 code and the 140-gene measurement with the MHD4 code. The Pearson correlation coefficient is 0.89 with a P value of 9 × 10−30. The dashed line corresponds to the y = x line.

As expected, the misidentification rate of this scheme was higher than that of the MHD4 code. Of all real RNA words, 77% were detected more frequently than the median count for the misidentification controls, instead of the 95% value observed in the MHD4 measurements. Using the same confidence ratio analysis as described above, we found that 73% (instead of 91% for the MHD4 measurements) of the 985 RNA species were measured with a confidence ratio larger than the maximum value observed for the misidentification controls (fig. S9C). RNA copy numbers measured from these 73% RNA species showed excellent correlation with our bulk RNA sequencing results (Pearson correlation coefficient r = 0.76) (Fig. 5B, black). The remaining 27% of the genes still exhibit good, albeit lower, correlation with the bulk RNA sequencing data (r = 0.65) (Fig. 5B, red), but we took the conservative measure of excluding them from further analysis.

The lack of an error correction capability also decreased the calling rate of each RNA species: When comparing the 107 RNA species common in both the 1001-gene and 140-gene measurements, we found that the copy numbers per cell of these RNA species were lower in the 1001-gene measurements (Fig. 5C and fig. S9D). The total count of these RNAs per cell was ~1/3 of that observed in the 140-gene measurements. Thus, the lack of error correction in the MHD2 code reduced the calling rate to ~30% of that of the MHD4 code, which is consistent with the decrease in calling rate observed for the MHD4 code when error correction was not applied. As expected from the quantitative agreement between 140-gene measurements and conventional smFISH results, comparison of the 1001-gene measurements with conventional smFISH results for 10 RNA species also indicated a calling rate that is ~1/3 of that observed for the MHD4 code (fig. S8C). Despite the expected reduction in calling rate, the good correlations found between the copy numbers observed in the 1001-gene measurements and those observed in the 140-gene measurements, as well as in conventional smFISH and bulk RNA sequencing measurements, indicates that the relative abundance of these RNAs can be quantified with the MHD2 encoding scheme.

Simultaneously imaging ~1000 genes in individual cells substantially expanded our ability to detect coregulated genes. The matrix of pairwise correlation coefficients determined from the cell-to-cell variations in the expression levels of these genes is shown in Fig. 6A. Using the same hierarchical clustering analysis as described above, we identified ~100 groups of genes with correlated expression (table S4). Nearly all of these ~100 groups showed statistically significant enrichment of functionally related GO terms (Fig. 6B and table S4). These included some of the groups identified in the 140-gene measurements, such as the group associated with cell-replication genes and the group associated with cell-motility genes (Fig. 6, A and B, groups 7 and 102), as well as many new groups. The groups identified here included 46 RNA species lacking any previous GO annotations, for which we can now hypothesize function on the basis of their group association (table S4). For example, KIAA1462 is part of the cell motility group, as also shown in the 140-gene experiments, suggesting a potential role of this gene in cell motility (Fig. 6A, group 102). Likewise, KIAA0355 is part of a new group enriched in genes associated with heart development (Fig. 6A, group 79), and C17orf70 is part of a group associated with ribosomal RNA processing (Fig. 6A, group 22). Using these groupings, we can also hypothesize cellular functions for 61 transcription factors and other partially annotated proteins of unknown functions (table S4). For example, the transcription factors Z3CH13 and CHD8 are both members of the cell-motility group, suggesting their potential role in the transcriptional regulation of cell-motility genes. Although these predicted functions based on gene-association analysis require further validation, our covariation data provide a resource for generating hypotheses on gene function and regulation.

Fig. 6 Covariation analysis of the RNA species measured in the 1001-gene measurements.

(A) Matrix of all pairwise correlation coefficients of the cell-to-cell variation in expression for the measured genes shown with the hierarchical clustering tree. The ~100 identified groups of correlated genes are indicated by color on the tree. Zoom-in of four of the groups described in the text are shown on the right. (B) Enrichment of 20 selected, statistically significantly enriched GO terms in the four groups. The statistically most significantly enriched GO terms (maximum 10) for each of the ~100 groups are shown in table S4.

Discussion

We have developed a highly multiplexed detection scheme for transcriptomic-scale RNA imaging in single cells. Using combinatorial labeling, sequential hybridization and imaging, and two different error-robust encoding schemes, we simultaneously imaged either 140 or 1001 genes in hundreds of individual human fibroblast cells. Of the two encoding schemes presented here, the MHD4 code is capable of both error detection and error correction and hence can provide a higher calling rate and a lower misidentification rate than can the MHD2 code, which instead can only detect but cannot correct errors. MHD2, on the other hand, provides a faster scaling of the degree of multiplexing with the number of bits than can MHD4. Other error-robust encoding schemes can also be used for such multiplexed imaging, and experimenters can set the balance between detection accuracy and ease of multiplexing according to the specific requirements of the experiments.

By increasing the number of bits in the code words, it should be possible to further increase the number of detectable RNA species by using MERFISH with either MHD4 or MHD2 codes. Because of their much slower increase in error rates with the number of bits, we expect the error-correcting encoding schemes, such as MHD4, to be more favorable for scaling up the measurements. For example, using the MHD4 code with 32 total bits and four or six 1 bits would increase the number of addressable RNA species to 1240 or 27,776, respectively; the latter is the approximate scale of the human transcriptome. The predicted misidentification and calling rates are still reasonable for the 32-bit MHD4 code (shown in Fig. 1, C and D, purple for the MHD4 code with four 1 bits, and similar rates were calculated for the MHD4 code with six 1 bits). If more accurate measurements are desired, an additional increase in the number of bits would allow the use of encoding schemes with a Hamming distance greater than 4, further enhancing the error detection and correction capability. Although an increase in the number of bits by adding more hybridization rounds would increase the data collection time and potentially lead to sample degradation, these problems could be mitigated by using multiple colors to readout multiple bits in each round of hybridization.

As the degree of multiplexing is increased, it is important to consider the potential increase in the density of RNAs that need to be resolved in each round of imaging. On the basis of our imaging and sequencing results, we estimate that including the whole transcriptome of the IMR90 cells would lead to a total RNA density of ~200 molecules/μm3. Using our current imaging and analysis methods, we could resolve 2 to 3 molecules/μm3 per hybridization round (38), which would reach a total RNA density of ~20 molecules/μm3 after 32 rounds of hybridization. This density should allow all but the top 10% most expressed genes to be imaged simultaneously or a subset of genes with even higher expression levels to be included. By using more advanced image analysis algorithms to better resolve overlapping images of individual molecules, such as compressed sensing (39, 40), it would be possible to extend the resolvable density by approximately fourfold and thus allow nearly the entire transcriptome, except for the top 2% most expressed genes, to be imaged all together. Last, theoretical predictions (17) indicate that the use of superresolution imaging (41, 42) could increase the resolvable density to ~105 molecules/μm3, which should be ample to address the entire transcriptome, even in cell types with RNA densities substantially higher than that of IMR90. However, RNAs in densely packed structures, such as p-bodies and stress granules, may still elude measurement.

We have illustrated the utility of the data derived from highly multiplexed RNA imaging by using covariation and correlation analysis to reveal distinct subcellular distribution patterns of RNAs, to constrain gene regulatory networks, and to predict functions for many previously unannotated or partially annotated genes with unknown functions. We anticipate that many more quantitative analyses could be applied to such data sets that include the spatial localization and copy number information of many RNA species in individual cells. Given its ability to quantify RNAs across a wide range of abundances without amplification bias while preserving native context, we envision that MERFISH will enable many applications of in situ transcriptomic analyses of individual cells in culture or complex tissues.

Materials and Methods

Probe design

Each RNA species in our target set was randomly assigned a binary code word either from all 140 possible code words of the 16-bit MHD4 code or from all 1001 possible code words of the 14-bit MHD2 code, as we describe in the main text. The encoding schemes are provided in tables S1 and S3.

We used array-synthesized oligopools as templates to make the encoding probes (22, 23). The template molecule for each encoding probe contains three components: (i) a central targeting sequence for in situ hybridization to the target RNA, (ii) two flanking readout sequences designed to hybridize each of two distinct readout probes, and (iii) two flanking primer sequences to allow enzymatic amplification of the probes (fig. S3). The readout sequences were taken from the 16 possible readout sequences, each corresponding to one hybridization round. The readout sequences were assigned to the encoding probes so that for any RNA species, each of the four readout sequences were distributed uniformly along the length of the target RNA and appeared at the same frequency. Template molecules for the 140-gene library also included a common 20-nucleotide (nt) priming region between the first polymerase chain reaction (PCR) primer and the first readout sequence. This priming sequence was used for the reverse transcription step described below. All template sequences are provided in table S5.

We embedded multiple experiments in a single array-synthesized oligopool and used PCR to selectively amplify only the oligos required for a specific experiment. Primer sequences for this indexed PCR reaction were generated from a set of orthogonal 25-nt sequences (43). These sequences were trimmed to 20 nt and selected for (i) a narrow melting temperature range (70 to 80°C), (ii) the absence of consecutive repeats of 3 or more identical nucleotides, and (iii) the presence of a GC clamp—one of the two 3′ terminal bases must be G or C. To further improve specificity, these sequences were then screened against the human transcriptome by using Basic Local Alignment Search Tool+ (BLAST+) (44), and primers with 14 or more contiguous bases of homology were eliminated. Last, BLAST+ was again used to identify and exclude primers that had an 11-nt homology region at the 3′ end of any other primer or a 5-nt homology region at the 3′ end of the T7 promoter. The forward primer sequences (primer 1) were determined as described above, whereas the reverse primers each contain a 20-nt sequence as described above plus a 20-nt T7 promoter sequence to facilitate amplification via in vitro transcription (primer 2). The primer sequences used in the 140-gene and 1001-gene experiments are listed in Table 1.

Table 1 Primer sequences used in the 140-gene and 1001-gene experiments.

View this table:

Thirty-nt-long readout sequences were created by concatenating fragments of the same orthogonal primer set generated above by combining one 20-nt primer with a 10-nt fragment of another. These readout sequences were then screened, by using BLAST+, for orthogonality with the index primer sequences and other readout sequences (no more than 11 nt of homology) and for potential off-target binding sites in the human genome (no more than 14 nt of homology). Fluorescently labeled readout probes with sequences complementary to the readout sequences were used to probe these readout sequences, one in each hybridization round. All used readout probes sequences are listed in Table 2.

Table 2 All used readout probes sequences.

View this table:

The readout probes used for the 140-gene libraries were probes 1 through 16. The readout probes used for the 1001-gene experiment were probes 1 through 14. “/3Cy5Sp/” indicates a 3′ Cy5 modification.

To design the central targeting sequences of the encoding probes, we first compiled the abundance of different transcripts in IMR90 cells using Cufflinks v2.1 (45), total RNA data from the Encyclopedia of DNA Elements (ENCODE) project (46), and human genome annotations from Gencode v18 (20). Probes were designed from gene models corresponding to the most abundant isoform by using OligoArray2.1 (47) with the following constraints: The target sequence region is 30-nt long; the melting temperatures of the hybridized region of the probe and cellular RNA target is greater than 70°C; there is no cross hybridization targets with melting temperatures greater than 72°C; there is no predicted internal secondary structures with melting temperatures greater than 76°C; and there is no contiguous repeats of six or more identical nucleotides. Melting temperatures were adjusted to optimize the specificity of these probes and minimize secondary structure while still producing sufficient numbers of probes for our libraries. To decrease computational cost, isoforms were divided into 1-kb regions for probe design. Using BLAST+, all potential probes that mapped to more than one cellular RNA species were rejected. Probes with multiple targets on the same RNA were kept.

For each gene in the 140-gene experiments, we generated 198 putative encoding probe sequences by concatenating the appropriate index primers, readout sequences, and targeting regions as shown in fig S3. To address the possibility that concatenation of these sequences introduced new regions of homology to off-target RNAs, we used BLAST+ to screen these putative sequences against all human ribosomal RNA (rRNA) and transfer RNA (tRNA) sequences as well as highly expressed genes [genes with fragments per kilobase per million reads (FPKM) > 10,000]. Probes with greater than 14 nt of homology to rRNAs or tRNAs or greater than 17 nt of homology to highly expressed genes were removed. After these cuts, we had ~192 (with a standard deviation of 2) probes per gene for both MHD4 codebooks used in the 140-gene experiments. We followed the same protocol for the 1001-gene experiments: Starting with 96 putative targeting sequences per gene, we obtained ~94 (with a standard deviation of 6) encoding probes per gene after these additional homology cuts. We decreased the number of encoding probes per RNA for the 1001-gene experiments so that these probes could be synthesized from a single 100,000-member oligopool as opposed to two separate pools. We designed each encoding probe to contain two of the four readout sequences associated with each code word; hence, only half of the bound encoding probes can bind readout probe during any given hybridization round. We used ~192 or ~94 encoding probes per RNA to obtain high signal-to-background ratios for individual RNA molecules. The number of encoding probes per RNA could be substantially reduced but still allow single RNA molecules to be identified (17, 48, 49). In addition, increasing the number of readout sequences per encoding probe or using optical sectioning methods to reduce the fluorescence background may allow further reduction in the number of the encoding probes per RNA.

We designed two types of misidentification controls. The first control—blank words—were not represented with encoding probes. The second type of control—no-target words—had encoding probes that were not targeting any cellular RNA. The targeting regions of these probes were composed of random nucleotide sequences subject to the same constraints used to design the RNA targeting sequences described above. Moreover, these random sequences were screened against the human transcriptome to ensure that they contain no substantial homology (>14-nt) to any human RNA. The 140-gene measurements contained five blank words and five no-target words. The 1001-gene measurements contained 11 blank words and five no-target words.

Probe synthesis

The encoding probes were synthesized by using the following four steps, and this synthesis protocol is illustrated in fig. S3.

Step 1: The template oligopool (CustomArray) was amplified via limited-cycle PCR on a Bio-Rad CFX96 by using primer sequences specific to the desired probe set. To facilitate subsequent amplification via in vitro transcription, the reverse primer contained the T7 promoter. All primers were synthesized by Integrated DNA Technologies (IDT). This reaction was column purified (Zymo DNA Clean and Concentrator, D4003).

Step 2: The purified PCR products were then further amplified ~200-fold and converted into RNA via a high yield in vitro transcription according to the manufacturer’s instructions [New England Biolabs (NEB), E2040S]. Each 20 μL reaction contained ~1 μg of template DNA from above, 10 mM of each NTP, 1× reaction buffer, 1× RNase inhibitor (Promega RNasin, N2611) and 2 μL of the T7 polymerase. This reaction was incubated at 37°C for 4 hours to maximize yield. This reaction was not purified before the following steps.

Step 3: The RNA products from the above in vitro transcription reaction were then converted back into DNA via a reverse transcription reaction. Each 50-μL reaction contained the unpurified RNA produce from step 2 supplemented with 1.6 mM of each dNTP, 2 nmol of a reverse transcription primer, 300 units of Maxima H- reverse transcriptase (Thermo Scientific, EP0751), 60 units of RNasin, and a final 1× concentration of the Maxima RT buffer. This reaction was incubated at 50°C for 45 min, and the reverse transcriptase was inactivated at 85°C for 5 min. The templates for the 140-gene libraries contain a common priming region for this reverse transcription step; thus, a single primer was used for this step when creating these probes. Its sequence was CGGGTTTAGCGCCGGAAATG. A common priming region was not included for the 1001-gene library; thus, the reverse transcription was conducted with the forward primer: CGCGGGCTATATGCGAACCG.

Step 4: To remove the template RNA, 20 μL of 0.25 M EDTA and 0.5 N NaOH was added to the above reaction to selectively hydrolyze RNA, and the sample was incubated at 95°C for 10 min. This reaction was then immediately purified by means of column purification using a 100-μg-capacity column (Zymo Research, D4030) and the Zymo Oligo Clean and Concentrator protocol. The final probes were eluted in 100 μL of ribonuclease (RNase)–free deionized water, evaporated in a vacuum concentrator, and then resuspended in 10 μL of encoding hybridization buffer (recipe below). Probes were stored at –20°C. Denaturing polyacrylamide gel electrophoresis and absorption spectroscopy were used to confirm the quality of the probes and revealed that this probe synthesis protocol converts 90 to 100% of the reverse-transcription primer into full-length probe and of the probe that is constructed, 70 to 80% is recovered during the purification step. This protocol is similar to another recently published protocol (23) but provides a substantially larger yield.

Fluorescently labeled readout probes have sequences complementary to the readout sequences described above and a Cy5 dye attached at the 3′ end. These probes were synthesized and purified by means of high-performance liquid chromatography (HPLC) by IDT.

Sample preparation and labeling with encoding probes

Human primary fibroblasts (American Type Culture Collection, IMR90), a commonly used cell line with a previously determined transcriptome (46), were used in this work. These cells are relatively large and flat, facilitating wide-field imaging without the need for optical sectioning. Cells were cultured with Eagle’s Minimum Essential Medium. Cells were plated on 22-mm, #1.5 coverslips (Bioptechs, 0420-0323-2) at 350,000 cells per coverslip and incubated at 37°C with 5% CO2 for 48 to 96 hours within petri dishes. Cells were fixed for 20 min in 4% paraformaldehyde (Electron Microscopy Sciences, 15714) in 1× phosphate buffered saline (PBS; Ambion, AM9625) at room temperature, reduced for 5 min with 0.1% w/v sodium borohydride (Sigma, 480886) in water to reduce background fluorescence, washed three times with ice-cold 1× PBS, permeabilized for 2 min with 0.5% v/v Triton (Sigma, T8787) in 1× PBS at room temperature, and washed three times with ice cold 1× PBS.

Cells were incubated for 5 min in encoding wash buffer comprising 2× saline-sodium citrate buffer (SSC) (Ambion, AM9763), 30% v/v formamide (Ambion, AM9342), and 2 mM vanadyl ribonucleoside complex (NEB, S1402S). Ten microliters of 100 μM (140-gene experiments) or 200 μM (1001-gene experiments) encoding probes in encoding hybridization buffer was added to the cell-containing coverslip and spread uniformly by placing another coverslip on top of the sample. Samples were then incubated in a humid chamber inside a 37°C-hybridization oven for 18 to 36 hours. Encoding hybridization buffer is composed of encoding wash buffer supplemented with 1 mg/mL yeast tRNA (Life technologies, 15401-011) and 10% w/v dextran sulfate (Sigma, D8906-50G).

Cells were then washed with encoding wash buffer, incubated at 47°C for 10 min, and this wash was repeated for a total of three times. A 1:1000 dilution of 0.2-μm-diameter carboxylate-modified orange fluorescent beads (Life Technologies, F-8809) in 2×SSC was sonicated for 3 min and then incubated with the sample for 5 min. The beads were used as fiducial markers to align images obtained from multiple successive rounds of hybridization, as described below. The sample was washed once with 2×SSC, and then post-fixed with 4% v/v paraformaldehyde in 2×SSC at room temperature for 30 min. The sample was then washed three times with 2×SSC and either imaged immediately or stored for no longer than 12 hours at 4°C before imaging. All solutions were prepared as RNase-free.

MERFISH imaging

The sample coverslip was assembled into a Bioptech’s FCS2 flow chamber, and the flow through this chamber was controlled via a home-built fluidics system composed of three computer-controlled eight-way valves (Hamilton, MVP and HVXM 8-5) and a computer-controlled peristaltic pump (Rainin, Dynamax RP-1). The sample was imaged on a home-built microscope constructed around an Olympus IX-71 body and a 1.45 NA, 100× oil immersion objective and configured for oblique incidence excitation. The objective was heated to 37°C with a Bioptechs objective heater. Constant focus was maintained throughout the imaging process with a home-built, autofocusing system. Illumination was provided at 641, 561, and 405 nm by using solid-state lasers (MPB communications, VFL-P500-642; Coherent, 561-200CWCDRH; and Coherent, 1069413/AT) for excitation of our Cy5-labeled readout probes, the fiducial beads, and nuclear counterstains, respectively. These lines were combined with a custom dichroic (Chroma, zy405/488/561/647/752RP-UF1) and the emission was filtered with a custom dichroic (Chroma, ZET405/488/561/647-656/752m). Fluorescence was separated with a QuadView (Photometrics) by using the dichroics T560lpxr, T650lpxr, and 750dcxxr (Chroma) and the emission filters ET525/50m, WT59550m-2f, ET700/75m, and HQ770lp (Chroma) and imaged with an EMCCD camera (Andor, iXon-897). The camera was configured so that a pixel corresponds to 167 nm in the sample plane. The entire system was fully automated, so that imaging and fluid handling were performed for the entire experiment without user intervention.

Sequential hybridization, imaging, and bleaching proceeded as follows. One milliliter of 10 nM of the appropriate fluorescently labeled readout probe in readout hybridization buffer (2×SSC; 10% v/v formamide, 10% w/v dextran sulfate, and 2 mM vanadyl ribonucleoside complex) was flown across the sample, flow was stopped, and the sample was incubated for 15 min. Then 2 mL of readout wash buffer (2×SSC, 20% v/v formamide, and 2 mM vanadyl ribonucleoside complex) was flown across the sample, flow was stopped, and the sample was incubated for 3 min. Two milliliters of imaging buffer comprising 2×SSC, 50 mM TrisHCl pH 8, 10% w/v glucose, 2 mM Trolox (Sigma-Aldrich, 238813), 0.5 mg/mL glucose oxidase (Sigma-Aldrich, G2133), and 40 μg/mL catalase (Sigma-Aldrich, C30) was flown across the sample (50). Flow was then stopped, and then ~75 to 100 regions were exposed to ~25 mW 642-nm and 1 mW of 561-nm light and imaged. Each region was 40 by 40 μm. The laser powers were measured at the microscope backport. Because the imaging buffer is sensitive to oxygen (51), the ~50 mL of imaging buffer used for a single experiment was made fresh at the beginning of the experiment and then stored under a layer of mineral oil throughout the measurement. Buffer stored in this fashion was stable for more than 24 hours.

After imaging, the fluorescence of the readout probes was extinguished via photobleaching. The sample was washed with 2 mL of photobleaching buffer (2×SSC and 2 mM vanadyl ribonucleoside complex), and each imaged region of the sample was exposed to 200 mW of 641-nm light for 3 s. To confirm the efficacy of this photobleaching treatment, imaging buffer was reintroduced, and the sample was imaged as described above.

The above hybridization, imaging, and photobleaching process was repeated either 16 times for the 140-gene measurements by using the MHD4 code or 14 times for the 1001-gene measurements by using the MHD2 code. An entire experiment was typically completed in ~20 hours.

After completion of imaging, 2 mL of a 1:1000 dilution of Hoescht (ENZ-52401) in 2×SSC was flown through the chamber to label the nuclei of the cells. The sample was then washed immediately with 2 mL of 2×SSC followed by 2 mL of imaging buffer. Each region of the sample was then imaged once again with ~1 mW of 405-nm light.

Because we imaged cells using wide-field imaging with oblique-incidence illumination, without optical sectioning and z-scanning, we quantified the fraction of individual RNA species that was outside the axial range of our imaging geometry for six different RNA species using conventional smFISH. For this purpose, we optically sectioned these cells by collecting stacks of images at different focal depths through the entire depth of the cells. We aligned the images in consecutive focal planes and then computed for each cell the fraction of RNAs that were detected in the three-dimensional stack but not in the basal focal plane. We found that only a small fraction, 15 ± 1% (mean ± SEM across six different RNA species) of RNA molecules were outside the imaging range of a fixed focal plane without z-scanning. These measurements also confirmed that our excitation geometry illuminated the full depth of our cells. From an imaging perspective, any optical sectioning technique could be used in MERFISH to allow the imaging of RNAs in thicker cells or tissues.

Construction of measured words

Fluorescent spots were identified and localized in each image by using a multi-Gaussian-fitting algorithm (38) assuming a Gaussian with a uniform width of 167 nm. This algorithm was used to allow partially overlapping spots to be distinguished and individually fit. RNA spots were distinguished from background signal—signal arising from probes bound nonspecifically, by setting the intensity threshold required to fit a spot with this software. Because of variation in the brightness of spots between rounds of hybridization, this threshold was adjusted appropriately for each hybridization round in order to minimize the combined average of the 1→0 and 0→1 error rates across all hybridization rounds (140-gene measurements) or to maximize the ratio of the number of measured words with four 1 bits to those with three or five 1 bits (1001-gene measurements). The location of the fiducial beads was identified in each frame by using a faster single-Gaussian fitting algorithm.

Images of the same sample region in different rounds of hybridization were registered by rotating and translating the image to align the two fiducial beads within the same image that were most similar in location after a coarse initial alignment via image correlation. All images were aligned to a coordinate system established by the images collected in the first round of hybridization. The quality of this alignment was determined from the residual distance between five additional fiducial beads, and alignment error was typically ~20 nm.

Fluorescence spots in different hybridization rounds were connected into a single string, corresponding to a potential RNA molecule, if the distance between spots was smaller than 1 pixel (167 nm). For each string of spots, the on-off sequence of fluorescent signals in all hybridization rounds were used to assign a binary word to the potential RNA molecule, in which 1 was assigned to the hybridization rounds that contained a fluorescent signal above threshold and 0 was assigned to the other hybridization rounds. Measured words were then decoded into RNA species by using the 16-bit MHD4 code or the 14-bit MHD2 code discussed in the main text. In the case of the 16-bit MHD4 code, if the measured binary word matched the code word of a specific RNA perfectly or differed from the code word by one single bit, it was assigned to that RNA. In the case of the 14-bit MHD2 code, only if the measured binary word matched the code word of a specific RNA perfectly was it assigned to that RNA. To determine the copy number per cell, the number of each RNA species was counted in individual cells within each 40- by 40-μm imaging area. This number accounts for the majority but not all RNA molecules within a cell because a fraction of the cell could be outside the imaging area or focal depth. Tiling images of adjacent areas and adjacent focal planes could be used to improve the counting accuracy.

In the 140-gene experiments, some regions of the cell nucleus occasionally contained too much fluorescence signal to properly identify individual RNA spots. In the 1001-gene experiments, the cell nucleus generally contained too much fluorescent signal to allow identification of individual RNA molecules. These bright regions were excluded from all subsequent analysis. This work focuses on mRNAs, which are enriched in the cytoplasm. To estimate the fraction of mRNAs missed by excluding the nucleus region, we used conventional smFISH to quantify the fraction of molecules found inside the nucleus for six different mRNAs species. We found that only 5 ± 2% (mean ± SEM across six RNA species) of these RNA molecules are found in the nucleus. Use of super-resolution imaging and/or optical sectioning could potentially allow individual molecules in these dense nucleus regions to be identified, which will be particularly useful for probing those noncoding RNAs that are enriched in the nucleus.

smFISH measurements of individual genes

Pools of 48 fluorescently labeled (Quasar 670) oligonucleotide probes per RNA were purchased from Biosearch Technologies. Thirty-nucleotide probe sequences were taken directly from a random subset of the targeting regions used for the multiplexed measurements. Cells were fixed and permeabilized as described above. Ten microliters of 250 nM oligonucleotide probes in encoding hybridization buffer (described above) was added to the cell-containing coverslip and spread uniformly by placing another coverslip on top of the sample. Samples were then incubated in a humid chamber inside a 37°C-hybridization oven for 18 hours. Cells were then washed with encoding wash buffer (described above) at 37°C for 10 min, and this wash was repeated for a total of three times. The sample was then washed three times with 2×SSC and imaged in imaging buffer by using the same imaging geometry as described above for MERFISH imaging.

Bulk RNA sequencing

Total RNA was extracted from IMR90 cells cultured as above using the Zymo Quick RNA MiniPrep kit (R1054) according to the manufacturer’s instructions. Polyadenylated [poly(A)] RNA was then selected (NEB, E7490), and a sequencing library was constructed by using the NEBNext Ultra RNA library preparation kit (NEB, E7530), amplified with custom oligonucleotides, and 150–base pair (bp) reads were obtained on a MiSeq. These sequences were aligned to the human genome (Gencode v18) and isoform abundance was computed with cufflinks (45).

Calculation of the predicted scaling and error properties of different encoding schemes

Analytic expressions were derived for the dependence of the number of possible code words, the calling rate, and the misidentification rate on N. The calling rate is defined as the fraction of RNA molecules that are properly identified. The misidentification rate is defined as the fraction of RNA molecules that are misidentified as a wrong RNA species. For encoding schemes with an error-detection capability, the calling rate and misidentification rate does not add up to 1 because a fraction of the molecules not called properly can be detected as errors and discarded and, hence, not misidentified as a wrong species. These calculations assume that the probability of misreading bits is constant for all hybridization rounds but differs for the 1→0 and 0→1 errors. Experimentally measured average 1→0 and 0→1 error rates (10 and 4%, respectively) were used for the estimates shown in Fig. 1, B to D. For simplicity, the word corresponding to all 0s was not removed from calculations.

For the simple binary encoding scheme in which all possible N-bit binary words are assigned to different RNA species, the number of possible code words is 2N. The number of words that could be used to encode RNA is actually 2N − 1 because the code word “00…0” does not contain detectable fluorescence in any hybridization round, but for simplicity the word corresponding to all 0s was not removed from subsequent calculations. The error introduced by this approximation is negligible. For any given word with m 1s and Nm 0s, the probability of measuring that word without error—the fraction of RNAs that is properly called—isEmbedded Image (1)where p1 is 1→0 error rate and p0 is 0→1 error rate per bit. Because different words in this simple binary encoding scheme can have different numbers of 1 bits, the calling rate for different words will differ if p1p0. The average calling rate, reported in Fig. 1C, was determined from the weighted average of the value of Eq. 1 for all words. This weighted average isEmbedded Image (2)where Embedded Image is the binomial coefficient and corresponds to the number of words with m 1 bits in this encoding scheme. Because in this encoding scheme every error produces a binary word that encodes a different RNA, the average misidentification rate for this encoding scheme, reported in Fig. 1D, follows directly from Eq. 2:Embedded Image (3)To calculate the scaling and error properties of the extended Hamming distance 4 (HD4) code, we first created the generator matrix for the desired number of data bits using standard methods (21). The generator matrix determines the specific words that are present in a given encoding scheme and was used to directly determine the number of encoded words as a function of the number of bits. In this encoding scheme, the calling rate corresponds to the fraction of words measured without error as well as the fraction of words measured with a single-bit error. For code words with m 1 bits, this fraction is determined by the following expression:Embedded Image (4)where the first term is the probability of not making any errors, the second term corresponds to the total probability of making one 1→0 error at any of the m 1 bits without making any other 0→1 errors, and the final term corresponds to the total probability of making one 0→1 error at any of the N-m 0 bits without making any 1→0 errors. Because the number of 1 bits can differ between words in this encoding scheme, the average calling rate reported in Fig. 1C was computed from a weighted average over Eq. 4 for different values of m. The weight for each term was determined from the number of words that contain m 1 bits as determined from the generator matrix described above.

Because RNA-encoding words are separated by a minimum Hamming distance of 4, at least 4 errors are required to switch one word into another. If error correction is applied, then 3 or 5 errors could also convert one RNA into another. Thus, we estimate the misidentification rate from all possible combinations of 3-bit, 4-bit, and 5-bit errors for code words with m 1 bits. Technically, >5-bit errors could also convert one RNA into another, but the probability of making such errors is negligible because of the small per-bit error rate. We approximate this expression withEmbedded Image (5)The first sum corresponds to all of the ways in which exactly four mistakes can be made. Similarly, the second and third sums correspond to all of the ways in which exactly three or five mistakes can be made. Equation 5 provides an upper bound for the misidentification rate because not all 3-, 4-, or 5-bit errors produce a word that matches or would be corrected to another legitimate word. Again because the number of 1 bits can differ between words, the average misidentification rate reported in Fig. 1D is calculated as a weighted average of Eq. 5 over the number of words that have m 1 bits.

To generate our MHD4 code in which the number of 1 bits for each code word is set to 4, we first generated the HD4 codes as described above, and then removed all code words that did not contain four 1s. The calling rate of this code, reported in Fig. 1C, was directly calculated from Eq. 4, but with m = 4 because all code words in this code have four 1 bits. The misidentification rate of this code, reported in Fig. 1D, was calculated by modifying Eq. 5 with the following considerations: (i) the number of 1 bits, m, was set to 4 and (ii) errors that produce words that do not contain three, four, or five 1 bits were excluded. Thus, the expression in Eq. 5 was simplified toEmbedded Image (6)Again, this expression is an upper bound on the actual misidentification rate because not all words with four 1s are valid code words.

Estimates of the 1→0 and 0→1 error rates for each hybridization round

To compute the probability of misreading a bit at a given hybridization round, we used the error-correcting properties of the MHD4 code. Briefly, the probabilities of 1→0 or 0→1 errors were derived in the following way. Let the probability of making an error at the ith bit, —ith hybridization round—be pi and the actual number of RNA molecules of the given species be A, then the number of exact matches for this RNA will be Embedded Image, and the number of one-bit error corrected matches for this RNA corresponding to errors at the ith bit will be Embedded Image. The pi can be directly derived from the ratio: Embedded Image. This ratio assumes that the 1-bit error-corrected counts were only generated from single-bit errors from the correct word and that multi-error contamination from other RNA words is negligible. Given that our error rate per hybridization round is small and that it takes at least three errors to convert one RNA-encoding word into a word that would be misidentified as another RNA, the above approximation should be a good one.

To compute the average 1→0 or 0→1 error probabilities for each of the 16 hybridization rounds, we use the above approach to calculate the per-bit error rates for each bit of every gene, sort these errors on the basis of whether they correspond to a 1→0 or a 0→1 error, and then take the average of these errors for each bit weighted by the number of counts observed for the corresponding gene.

Estimates of the calling rate for individual RNA species from actual imaging data

With the estimates of the 1→0 or 0→1 error probabilities for each round of hybridization as determined above, it is possible to estimate the calling rate for each RNA according to the specific word used to encoded it. Specifically, the fraction of an RNA species that is called correctly is determined byEmbedded Image (7)where the first term represents the probability of observing an exact match of the code word and the second term represents the probability of observing an error-corrected match (with 1-bit error). The values of the per-bit error rate pi for each RNA species are determined by the specific code word for that RNA and the measured 1→0 or 0→1 error rates for each round of hybridization. If the code word of the RNA contains a 1 in the ith bit, then pi is determined from the 1→0 error rate for the ith hybridization round; if the word contains a 0 in the ith bit, pi is determined from the 0→1 error rate for the ith hybridization round.

Hierarchical clustering analysis of the co-variation in RNA abundance

Hierarchical clustering of the covariation in gene expression for both the 140- and 1001-gene experiments was conducted as follows. First, the distance between every pair of genes was determined as 1 minus the Pearson correlation coefficient of the cell-to-cell variation of the measured copy numbers of these two RNA species, both normalized by the total RNA counted in the cell. Thus, highly correlated genes are “closer” to one another, and highly anticorrelated genes are “further” apart. An agglomerative hierarchical cluster tree was then constructed from these distances using the unweighted pair group method with arithmetic mean (UPGMA). Specifically, starting with individual genes, we constructed hierarchical clusters by identifying the two clusters (or individual genes) that are closest to one another according to the arithmetic mean of the distances between all intercluster gene pairs. The pairs of clusters (or individual genes) with the smallest distance are then grouped together, and the process is repeated. The matrix of pairwise correlations was then sorted according to the order of the genes within these trees.

Groups of genes with substantial covariations were identified by selecting a threshold on the hierarchical cluster tree (indicated by the dashed lines in Figs. 3D and 6A) that produced ~10 groups of genes, each of which contains at least four members for the 140-gene experiments or ~100 groups each of which contains at least three members for the 1001-gene experiments. One can change the threshold in order to identify either more tightly coupled smaller groups or larger groups with relatively loose coupling.

A probability value for the confidence that a gene belongs to a specific group was determined by computing the difference between the average correlation coefficient between that gene and all other members of that group and the average correlation coefficient between that gene and all other measured genes outside that group. The significance (P value) of this difference was determined with the student’s t test and is provided in tables S2 and S4.

Because hierarchical clustering is inherently a one-dimensional analysis—any given genes can only be a member of a single group—this analysis does not allow all correlated gene groups to be identified. Higher-dimension analysis, such as principal component analysis or k-means clustering, could be used to identify more covarying gene clusters (30).

Analysis of RNA spatial distributions

To identify genes that have similar spatial distributions, we subdivided each of the measured cells into 2 by 2 regions and calculated the fraction of each RNA species present in each of these bins. To control for the fact that some regions of the cell naturally contain more RNA than others, we calculate the enrichment for each gene—the ratio of the observed fraction in a given region for a given RNA species to the average fraction observed for all genes in that same region. For each pair of RNA species, we then determined the Pearson correlation coefficient of the region-to-region variation in enrichment of these two RNA species for each cell and averaged the correlation coefficients over ~400 cells imaged in seven independent data sets. We then clustered RNA species on the basis of these average correlation coefficients using the same hierarchical clustering algorithm described above. Because of the large number of cells used for the analysis, we found that the coarse spatial binning (2 by 2 regions per cell) was sufficient to capture the spatial correlation between genes, and finer binning did not produce more significantly correlated groups.

To measure the distances of genes from the nuclei and from the cell edge, we first used brightness thresholds on our cell images to segment the nuclei and identify the cell edge. We then measured the distance from every RNA molecule to the nearest part of the nucleus and nearest part of the cell edge. For each data set, we computed the average distance for each RNA species averaged over all the cells measured. We then averaged these distances for the group I genes, group II genes, or all genes. Only those RNA species with at least 10 counts per cell were used in this analysis to minimize statistical error on the distance values.

GO analysis

Groups of genes were selected from the hierarchical trees as discussed above. A collection of GO terms (31) was determined for all measured RNA species as well as the RNA species associated with each group from the most recent human GO annotations (http://geneontology.org/page/download-annotations) by using both the annotated GO terms and terms immediately upstream or downstream of the found annotations. The enrichment of these annotations was calculated from the ratio of the fraction of genes within each group that have this term to the fraction of all measured genes that have this term and the P value for this enrichment was calculated via the hypergeometric function. Only statistically significantly enriched GO terms with a P value less than 0.05 were considered.

Supplementary Materials

References and Notes

  1. Acknowledgments: We thank H. Babcock for technical advice and help with instrumentation and B. Bintu for aid in error analysis. This work was in part supported by the National Institutes of Health. K.H.C. acknowledges a National Science Scholarship from the Agency for Science, Technology and Research of Singapore. A.N.B. acknowledges support by the Damon Runyon Foundation postdoctoral fellowship. J.R.M. acknowledges support from the Helen Hay Whitney Foundation postdoctoral fellowship. S.W. acknowledges support from Jane Coffins Child Foundation postdoctoral fellowship. X.Z. is a Howard Hughes Medical Institute investigator. RNA-seq data were deposited under the accession number GSE67685. MERFISH data are available at http://zhuang.harvard.edu/merfish. X.Z., K.H.C., A.N.B., J.R.M., and S.W. are inventors on a patent applied for by Harvard University that covers the MERFISH method described here; X.Z., J.R.M., and A.N.B. are inventors on a patent applied for by Harvard University that covers the probe synthesis method described here.
View Abstract

Navigate This Article