Ethics statement

The Cold Spring Harbor Laboratory Animal Care and Use Committee approved all animal procedures, which were carried out in accordance with the Institutional Animal Care and Use Committee protocol 19-16-10-07-03-00-4 at Cold Spring Harbor Laboratory.

Data

All the experimental data were obtained either using an improved version of BARseq [13] targeting 79 endogenous mRNAs in mouse visual cortex, or from a recent study [14] targeting 65 endogenous mRNAs in mouse motor cortex. Gene identities were read out using a seven-nucleotide gene identification index (GII), which was designed with a minimal Hamming distance of three nucleotides between each pair of GIIs. Rolonies were prepared as described by [14]. Imaging was performed on an Olympus IX81 inverted scope with a Crest Xlight2 spinning disk confocal, a Photometrics BSI Prime camera, and an 89 North LDI 7-line laser source. All images were acquired using an Olympus UPLFLN 40x 0.75 NA objective. The microscope was controlled by micro-manager [15]. See S1 Appendix (Section A) for the preprocessing steps for this data, and S1 Appendix (Section G) for the process of generating the simulation data.

Notation and observation model

Formally speaking, what is the result of a spatial transcriptomics imaging experiment? For each voxel (m) in the tissue, at each imaging round (r), in each color-channel (c), we record a fluorescence intensity. We will use X_m,r,c to denote this fluorescence intensity. Our task is to use X to uncover the presence and identity of rolonies in the tissue. (Throughout we assume that X is preprocessed, including background removal and image registration, hence that there are no systematic shifts of the image between imaging rounds. See S1 Appendix Section A for more detail.) Below we describe the parameters used to model the physical process that yields these intensities:

1. The rolonies, F. The transcripts in the tissue are amplified in place into a ‘rolony’ structure which is easy for fluorophores to bind to [7]. Each voxel m may contain a different amount of rolony material, and hence a varying level of fluorescence signal. We refer to the amount of material in voxel m for rolonies associated with barcode j as the rolony density. We denote this density by F_m,j. The variable F indicates where rolonies are and how bright we should expect them to be. This density should always be non-negative. Note that F cannot be observed directly—instead, we observe fluorescence signal in different rounds and channels, and must use these signal observations to estimate the rolony densities.
2. The codebook, B. In each imaging round r, the rolonies associated with gene j will bind to specific fluorescently labeled detection probes. We use the binary variable B to indicate which imaging rounds and fluorescent probes each gene is associated with. Specifically, we let B_r,c,j = 1 whenever a rolony with barcode j should bind to a fluorescent probe associated with specific color-channel c in imaging round r, and otherwise we let B_r,c,j = 0. Here we assume B is known. The vector of values of B for a particular gene j is known as the ‘barcode’ for that gene, and the collection B of all the barcodes is known as the ‘codebook.’
3. The probe response functions, K, φ. If a probe centered at a particular voxel is illuminated with a particular wavelength, the probe will emit a certain amount of signal which we can record at the corresponding voxel. We may also observe dimmer responses at neighboring voxels, due to the possible spreading of the single point object in the optical system. We use a non-negative matrix K to denote the point-spread function, i.e., the typical fluorescence signal-levels produced at each voxel in the neighborhood of a probe. We use the matrix φ to represent the responsiveness of each type of fluorescent probe to each wavelength; each element of this matrix lies in the range of [0, 1]. Here we assume that the number of types of fluorescent probes is the same as the number of color-channels measured (though this could be relaxed). We further assume that the voxel-resolution of the rolony density is the same as the voxel-resolution of the original images.
4. Phasing, ρ. A washing process is applied after each round of imaging. However, in practice this washing step may not completely remove all of the reagents from every voxel. This can result in a ‘ghost’ of one round appearing in the next rounds. For each color-channel c, we let ρ_c ∈ [0, 1] indicate the fraction of activity which appears as a ‘ghost’ signal in the next round.
5. Background, a. The images we obtain may also include background fluorescence from the tissue. We assume that the background is constant across rounds. We model this effect using a non-negative per-voxel value a_m for each voxel m.
6. Per-round per-wavelength gain, α, and baseline, b. The brightness observed from all rolonies at a particular color-channel in a particular round may have an associated gain factor. We model this gain factor with a non-negative per-round (r) per-channel (c) multiplier α_r,c and non-negative intercept b_r,c,.

Putting all these pieces together, we obtain an observation model. This model states that the observed brightnesses X_m,r,c should be given by the formulae Here the variable Z is used to incorporate the round-phasing effects; i.e., Z_r,c,j measures the concentration of probes of type c which we would expect at round r, arising from a rolony with barcode j. We will also find it convenient to define This represents the total contribution of fluorescence signal expected to arise in round r and channel c from a rolony of type j. A summary of notation can be found in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Notation.

https://doi.org/10.1371/journal.pcbi.1008256.t001

Overall, the model introduced above could certainly be expanded to model the physical imaging process more accurately, but we found that it was sufficient for our purposes: detecting and demixing rolonies.

Inference

Our task is to uncover the positions and barcodes of rolonies in the tissue. According to the model in the previous section, this information can be obtained from the rolony density variable, F. However, F cannot be directly measured; thus our primary task is to estimate F from the original image data. To do this we must in a sense invert the observation model specified above: the observation model tells us how rolony densities give rise to the fluorescence signal, but we would like to use observations of the fluorescence signal to estimate the rolony densities.

Finding rolonies.

Let us now assume we have used the non-negative regression framework to estimate F (the collection of rolony density images, one for each barcode). These per-barcode density images indicate the positions of rolonies that belong to a particular barcode; see the left side of Fig 1 for a schematic. We can then apply a blob-finding algorithm to these per-barcode images to find the rolonies for each barcode; in practice we simply find local maxima in the per-barcode images.

Finding rolonies, or ‘blobs,’ in the per-barcode images is easier than finding blobs in the original images. See Fig 2A as an example. The per-barcode images include fewer blobs and the blobs are smaller, so there are fewer problems with overlapping blobs. More specifically:

1. The rolony densities are demixed. There are fewer blobs in each rolony density than in the original image stack. In the observed images, the intensity measured for each voxel for each wavelength at each round is a sum of contributions from all nearby rolonies which emit signal at that wavelength in that round. By contrast, the intensity measured at a particular voxel in the per-barcode images is only the sum of contributions from rolonies with that one specific barcode.
2. The rolony densities are deconvolved. The blobs are smaller in the rolony density than in the original image stack. In the observed images, the intensity at a voxel is a contribution from all rolonies which are within the radius of the point-spread function K. Recall that this function smears signal from a single voxel across all nearby voxels. By contrast, the intensity of a per-barcode image at a particular voxel represents the amplification level of rolonies in that one voxel. In this sense, the inference process attempts to invert the point-spread function (i.e., perform deconvolution). On its own, this inversion process would not be numerically stable; however, the sparsity penalty and non-negativity constraint ensure it is numerically well-behaved [16].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. BarDensr accurately estimates sparse, per-gene rolony densities from the proposed observation model.

(A) Rolony densities make it easier to detect rolonies. The left plot shows the max-projection of the original experimental image across all rounds and channels; detecting blob-like structures in this image can be challenging, especially when two rolonies are in close proximity. By contrast, the rolony densities for particular genes are sparser, so it is easier to identify the positions of individual rolonies in the tissue. The middle and right plots show examples of these rolony densities. The orange marks represent rolonies detected by a hand-curated approach. Note that the rolony densities appear to show several rolonies which were missed by the hand-curated approach (see S3 Fig for further details). (B) BarDensr accurately recovers the ground truth in simulated data. The left plot shows the simulated data in all rounds and channels. In the right plot, we applied BarDensr to this simulated data, and found that we were able to largely recover the true rolonies in this simulation (shown on the first column). The final column of plots shows the rolony densities learned by BarDensr, which shows that the algorithm accurately recovers most of the simulated ground truth rolonies, with a few mistakes. The middle column of plots shows a blurred version of the rolony densities and the spots discovered from these rolony densities.

https://doi.org/10.1371/journal.pcbi.1008256.g002

The rolony density F thus represents a demixed and deconvolved version of the raw data. The original data is mixed, insofar as each raw intensity represents contributions from many barcodes. It is also convolved, insofar as each raw intensity represents contributions from many positions in space via the point-spread function. The non-negative sparse regression allows us to simultaneously demix and deconvolve, yielding per-barcode images which are cleaner and easier to understand.

Although it is easier to find blobs in the rolony densities, there is still one obstacle to be overcome: the threshold. Any blob-finding algorithm must specify an intensity above which a blob is considered real. How can this threshold be chosen? Here we make use of ‘unused barcodes.’ There could be as many as C^R unique barcodes in a codebook for an experiment with R rounds and C channels of measurement (assuming only one channel emits signal in each round, which is the case in the experiments we studied). However, most of these barcodes are not used in the actual experiment. These unused barcodes give us a way to pick a sensible threshold. Along with the real codebook, we additionally include several unused barcodes; we enumerated all possible barcodes such that each round contained exactly one active channel, then selected uniformly at random from the set of barcodes such that each barcode differed from every other barcode in at least three rounds. We then run BarDensr on this augmented codebook. Blobs in the rolony densities associated with the unused barcodes must correspond to noise, since the true data-generating process did not include any signal from such barcodes. We therefore set the threshold to be the smallest value which guarantees that no blobs were detected in the unused barcodes. (In practice, using just two unused barcodes sufficed to estimate a stable and accurate threshold.)

Code availability

The BarDensr Python package is available from https://github.com/jacksonloper/bardensr. The NeuroCAAS implementation of BarDensr can be found at http://www.neurocaas.com/analysis/8. This NeuroCAAS implementation requires no software or hardware installation by the user. As input, the user must upload a stack of images, a codebook, and a configuration file specifying parameters such as the radius of the smallest rolonies of interest (see S1 Appendix, Section C, as well as the NeuroCAAS link above for further details regarding the data format.) We assume that the images have been registered and background-subtracted before input into NeuroCAAS. There are three outputs from BarDensr NeuroCAAS implementation:

1. The first output indicates candidate rolony locations. This output is represented in a comma-separated-value file listing all entries in the rolony density F which have signal greater than zero.
2. The second output is designed to help the user assess the quality of the results; it is an HDF5 file which stores diagnostics about each detected rolony. See the next section for details.
3. The third output is a plot designed to help the user assess the overall quality of the results. For each detected rolony, we estimate the quality of evidence for that rolony by taking a key correlation coefficient (see the next section for details). The plot shows the histogram of these quality measures over all spots.

For further details on the AWS hardware used in the NeuroCAAS implementation, see S1 Appendix (Section D).

The rolony densities estimated by BarDensr provide sparse, single images to detect spots for individual barcodes

As described in the Methods Section, the sparse non-negative regression approach aims to yield per-gene rolony density images which are easy to work with. The cartoon in Fig 1 may help illustrate this idea. Our belief is that the true per-gene rolony densities will be sparse images, so the learned rolony densities should also be sparse images.

To test this belief, we applied BarDensr to the experimental data described in the Methods Section. Fig 2A compares the raw data with the learned rolony densities for Nrgn and Slc17a7 in a small region of the tissue. As hoped, the rolony densities are indeed quite sparse compared to the raw data. This ensures that blob-detection is relatively easy. This figure also shows that many of the bright spots in the rolony density images appear at rolony locations found by a hand-curated method (see S1 Appendix, Section E for details). For visualization purposes, Fig 2A shows the blurred version of the rolony densities (i.e. KF); these make it easier to see the bright spots.

To get a sense for what all the different genes look like, we also examined the rolony densities for all the barcodes (81 in total in this dataset, including two unused barcodes); see S1 and S2 Figs. These sparse images enable us to identify the rolony location easily for each barcode.

BarDensr provides improved demixing and detection accuracy compared to existing approaches

To benchmark BarDensr against other methods, we generated simulated data with rolony density, gene expression levels, and noise levels matched to the experimental data, as shown in Fig 2B. Then, we examined how well we could recover the ‘true’ rolonies from the simulation data. Qualitative results for several different genes are shown in S4 Fig. Quantitatively, we present a Receiver Operating Characteristic curve (ROC curve) in Fig 3A, which summarizes the percentage of true detected rolonies (also known as ‘1-FNR’, the complement of the False Negative Rate (FNR)). Depending on the False Positive Rate (FPR) we are willing to tolerate, different detection rates can be achieved; the ROC curve summarizes this relationship.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. BarDensr outperforms the existing methods by discovering more rolonies.

(A) Performance of BarDensr on simulated data. What percentage of rolonies are correctly detected? We use the Receiver Operating Characteristic curve (ROC curve) to look at this percentage (the complement of False Negative Rates, or 1-FNR) as a function of the tolerated False Positive Rate (FPR), for BarDensr (red), starfish (orange), Single Round Matching (SRM, green), as well as the correlation-based method (‘corr’, gray); cf. S1 Appendix, Sections E and F for details on these other methods. S4 Fig illustrates these simulation data. In drawing these curves, we consider two qualitatively different kinds of errors: errors because a rolony isn’t detected at all (dotted lines), and errors because a rolony is detected but it is assigned the wrong barcode (solid lines). The left plots show these curves for simulated data. The right plots show these curves for simulated data with ‘dropout’—a form of noise present in some spatial transcriptomic methods (cf. S1 Appendix, Section G for details). For all four kinds of simulations, we found BarDensr is able to find significantly more spots. (B) Performance of BarDensr on the hybrid simulation. Simulated data is always imperfect; to try to measure performance on a more realistic dataset, we used a hybrid method a la [18]. We injected fake rolonies into real data, and quantified how well different methods could recover these fake spots. The plots above show 1-FNR (y-axis) as the function of scale intensity of the fake rolonies (x-axis) and number of fake rolonies injected (S, colored lines), without (top) and with (bottom) dropout, using BarDensr (left) and SRM (right). See S1 Appendix, Section G for details.

https://doi.org/10.1371/journal.pcbi.1008256.g003

We compare BarDensr to several other approaches. Starfish is one package developed for analyzing spatial transcriptomics data. This method has many hyperparameters. To give this method its best chance, we first tried to find the best parameters manually, and additionally used the BayesianOptimization package [17] to find the hyperparameters which allowed it to perform as well as possible on the simulated data. Fig 3A shows that this performance falls short of the detection rates achieved by BarDensr. We also investigated SRM (see S1 Appendix, Section E) and a correlation-based method (‘corr’, see S1 Appendix, Section F) for comparison. These two methods represent ‘blobs-first’ and ‘barcodes-first’ approaches. BarDensr has better recovery prediction than either of these methods.

Our simulated data here do not capture the full biological content of the real observed data. For example, in real data, the tissue often has some regions with dense rolony concentrations (e.g. nuclei) and other regions which are more sparse. In order to quantify performance in more realistic biological contexts, we performed a ‘hybrid’ simulation, a la [18]. We started with the original experimental data and injected varying numbers of spots at random locations in the image with varying peak intensities (cf. S1 Appendix, Section G). To test if the model is able to recover these injected spots with the original image background, we computed 1-FNR (FPR could not be computed here since we do not know the ground truth in the original experimental data). We ran two variants on this simulation: one ordinary simulation and one simulation with ‘dropout,’ in which some rolonies emit a strong bright signal in most of the rounds but simply vanishes in one or more rounds (see S1 Appendix, Section G). The results of the dropout and non-dropout experiments are shown in Fig 3B. As expected, the performance decreases when the intensity of the injected spots is smaller. However, as long as the intensity of injected spots was at least half the maximum intensity of the original image, and the number of injected spots is in a reasonable range, BarDensr was able to find all the spots, even in the simulation with dropout; by contrast, the SRM approach was unable to find all the injected spots in the hybrid experiment, especially in the dropout variant.

Errors are mostly mis-identification on the barcodes, not missed detections

We used simulated data to investigate the failures represented by the FPR and FNR described above: are they caused by failure in assigning the rolonies to the correct barcodes (‘barcode mis-identification’), or failure in detecting rolonies? To find out, we computed how the failure rates would change if mis-identified barcodes were not considered ‘errors.’ We denote this the ‘total hit rate’ analysis (cf. Fig 3A, dotted lines); both BarDensr and SRM have very high total hit rates for the simulated data examined here, indicating that both of these methods detect spots well, but sometimes mis-classify the spot identity. See S1 Appendix (Section H) for further details.

BarDensr remains effective on data with low spatial resolution

High-resolution imaging can be expensive and time-consuming. BarDensr can also work on low-resolution images. To show this, we spatially downsampled the experimental images for each frame (each round and each channel). We then fit BarDensr to these lower-resolution images. An example is shown in Fig 4 (additional examples with 5× and 10× lower resolutions can be seen in S5 Fig). These figures suggest that BarDensr correctly detects the overall expression levels of each gene in low-resolution images—even when the downsampling is so extreme that picking out individual rolonies is not feasible.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. BarDensr is the best choice for dense and/or low-resolution experimental data.

(A) The 5× downsampled image is compared to the original ‘fine scale’ image. All these plots show the max-projection across all rounds and channels, with the right two showing the zoomed region indicated by the red rectangles in the left two. Note that it is difficult to visually isolate single spots from the downsampled image. To test the performance of BarDensr on this low-resolution data, we first run the model on the original data, obtain rolony densities, and then finally downsample the rolony densities. Next, we run BarDensr on downsampled data, and examine the estimated rolony densities. (B) The rolony densities for a selected gene (Slc17a7) estimated using the original fine scale (left), as well as these two approaches. For a more complete example, see S5 Fig. (C) The cell-level gene expression quantification, for those genes that have more than four spots in the fine scale in a 1000 × 1000 region. The color of the heatmap indicates the proportion of gene counts (i.e., the total counts of each gene divided by the total counts of all genes detected in the region). The x-axis represents the 24 genes that were chosen, ordered based on the counts in the fine scale. The y-axis represents the cells, ordered based on the hierarchical clustering result from the fine scale, as shown in the dendrogram on the left. A total of 43 cells are segmented from the original image using a seeded watershed algorithm (cf. S1 Appendix, Section I). The two different results yield nearly identical clusterings, indicating that BarDensr recovers gene activity with accuracy sufficient to cluster cells even given low-resolution images. (D) In order to evaluate the performance of BarDensr and compare it with the other state of the art methods, we took a 200 × 200 region of real experimental data from a recent study [14] and created a benchmark dataset (cf. S1 Appendix, Section J). Above we plot the ROC performance on this benchmark for four different methods (BarDensr, correlation-based method, as well as starfish with ‘spot-based’ and ‘pixel-based’ approaches), using the original data and denser / relatively low-resolution data. While all the methods performed quite well on the original images (left), BarDensr has a better performance than the others when the lower-resolution images were used (right).

https://doi.org/10.1371/journal.pcbi.1008256.g004

To more rigorously test if BarDensr can recover the correct gene expression level when applied on the low-resolution data, we quantified the cell-level gene activity on a larger region where 43 cells are detected using a seeded watershed algorithm (see S1 Appendix, Section I for detail). Fig 4C suggests that with 5× downsampled data, the cell-level gene expression, as well as the cell clusters, are preserved with high consistency compared to the results of applying the method to the original fine scale.

Finally, we wanted to test if BarDensr can correctly identify individual rolonies on a low-resolution version of the experimental data from a recent study [14]. There was no ground truth available for this real data. We therefore manually inspected every location in a small region (200 × 200 voxels) to identify all of the rolonies and created a benchmark dataset. We compared the performance of BarDensr, the correlation-based method, and starfish with ‘spot-based’ (the method that is recommended for this data based on starfish) and ‘pixel-based’ methods. The left plot in Fig 4D shows that all these four methods perform quite well at the original resolution, with the correlation-based method and starfish performing faster compared to BarDensr (0.59 seconds for the correlation-based method; 1.22 and 8.16 seconds for starfish pixel-based and spot-based method, respectively; BarDensr used 21.77 seconds). When the original images were 5 times blurred and downsampled to emulate denser data, as shown on the right plot of Fig 4D, BarDensr recovered more correct rolonies compared to the other methods. For experiments with dense and/or low-resolution images, these results (and the simulation results above) suggest that BarDensr is the best available analysis tool. More details on how we generated the benchmark can be found in S1 Appendix (Section J).

BarDensr computations can be scaled up to large datasets via sparsifying and coarsening accelerations

In the Methods Section, we described how the barcode sparsity could help us potentially apply the method to a large dataset with more barcodes. To test if we can use a much larger dataset, we considered a simulated example with more unique barcodes (53,000 unique barcodes and 17 sequencing rounds). With so many barcodes, naively running BarDensr is prohibitively expensive (in both compute time and memory) on large datasets. However, we also expect such datasets are extremely sparse in terms of barcodes—any given small region of the image is quite unlikely to include rolonies from all 53,000 barcodes. This is particularly true for some of the latest applications of multiplexed imaging technologies, in which each unique barcode corresponds to dendrites from a unique cell, instead of each barcode corresponding to a gene [13]. In these experiments, a small region of tissue may contain many different transcripts, but it will only contain portions of a small number of different cells. Thus we should be able to take advantage of this sparsity to speed up BarDensr. We simulated a 50 × 80 small region where 40 rolonies were present in total. We then obtained a coarse, downsampled image, and then ran BarDensr and learned the parameters for this low-resolution data. If the learned parameters from the coarse scale indicated a particular barcode did not appear, then we assumed that this barcode should be absent even if we used the data at the original resolution. The result in S6 Fig shows nearly perfect prediction performance. This problem was quite small, so we could also run the method without using any sparsity-based acceleration techniques; we found that the unaccelerated version did not outperform the accelerated version, suggesting that BarDensr can be used for datasets of this kind with a larger number of molecular or cellular barcodes (cf. [13, 19–21]).

When the number of barcodes is small, BarDensr can run without these acceleration techniques—but these accelerations may still be worth applying, to help cut down on computation times and reduce memory usage. As shown in Fig 5A, we found that these techniques reduced runtime by a factor of four in simulated data involving 79 barcodes (see S1 Appendix, Section K for details). Fig 5B further shows that these acceleration techniques do not significantly alter the results, even on real data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. BarDensr can be scaled up with sparsifying and coarse-to-fine approaches.

(A) Coarse-to-fine acceleration. The Area Under the ROC (AUROC) summarizes the performance of a method by calculating the integral of the ROC curve (higher is better). The black curve plots the AUROC performance on the simulated data, against the number of seconds the default iterative algorithm has been allowed to use, up to a maximum of 15 iterations. We can also use a coarse-to-fine strategy, where we first run the algorithm on downsampled data for 20 iterations, and use the results to perform 10 additional iterations on the full high-resolution data; the red curve plots the performance for this strategy. (B) BarDensr can take advantage of gene-sparsity. Here we used two different approaches to analyze a 1000 × 1000 region of the experimental data. The first approach uses BarDensr naively, applying it directly to the image. The second approach, illustrated on the first two plots, accelerates the method using a ‘coarse-to-fine’ method by taking advantage of ‘gene-sparsity.’ Specifically, we split this region into 4 × 4 patches (the borders of these patches are indicated as the white lines on the left plot). After the relatively fast ‘coarse’ step, the barcodes that have very low maximum rolony densities were removed before the following ‘fine’ step. This keeps only a relatively small number of barcodes to consider for each patch (ranging from 38 to 65 out of 81 barcodes, as shown in the middle plot), therefore reducing the computation time and the memory usage for the ‘fine’ step later (cf. S1 Appendix, Section K for more detail). We here show that both methods yield nearly the same result, as shown in the ROC curves on the right plot. In particular, we treated one method as the ‘truth’ and constructed an ROC curve indicating the accuracy of the other method. We can then do the reverse, treating the other method as ‘truth.’ The results suggest strong agreement.

https://doi.org/10.1371/journal.pcbi.1008256.g005

Making use of these acceleration approaches, we applied BarDensr to a large section of mouse motor cortex from a recent study [14]. As shown in Fig 6, BarDensr was able to reveal the striation of the cell type marker genes and identify cell types consistent with known markers. The entire region shown is composed of 35 fields of view (each field of view is of size 2048 × 2048 × 1 voxels), and 65 genes are targeted. The entire process with the usage of acceleration approaches takes approximately 3 hours to compute on a GPU machine (a p3.2xlarge machine on the AWS service).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. BarDensr reveals the laminar distribution of the cell type marker genes, as well as the identified cell types in the motor cortex.

The data was obtained from a recent study [14]. Top: each dot indicates a detected rolony; for each cell-type, we use the same color for all rolonies associated with marker-genes for that cell-type. Bottom: each dot indicates a cell. Dots are colored based on the cell-types; cell-types are estimated using the detected rolonies from the top plot.

https://doi.org/10.1371/journal.pcbi.1008256.g006

BarDensr recovers interpretable parameters

BarDensr uses a data-driven approach to estimate all the relevant features of the physical model: the per-channel phasing factor, the per-round per-channel scale factor, the per-round per-channel offset, the per-voxel background, the per-wavelength response matrix, and the rolony densities (the latter of which have already been described in detail above). In the data analyzed here, we found that the per-channel phasing factor was relatively small, suggesting very little ‘ghosting’ in this data. The wavelength-response matrix was almost diagonal, although we found some slight color-mixing from channel 2 to channel 1, consistent with visual inspection (see the fifth round in S7(B) Fig as an example). This indicates that our model is able to correctly recover the color-mixing effects. We also investigated whether all of the features of our model were necessary for the purposes of finding rolonies. For each feature of the model, we tried removing that aspect of the model and seeing whether the method still performed well. For the data analyzed here, we found that the φ and ρ parameters were not essential (though they did seem to improve the performance, at least qualitatively). By contrast, all of the other parameters were essential; removing any of them yielded nonsensical results.

The BarDensr model correctly predicts the observed signal intensities

This algorithm is based upon a physical model of how this data is generated. Rolonies appear at different positions in the tissue, they emit fluorescence signal in different conditions, the fluorescence signal is smeared by a point-spread function, and finally we observe this signal, together with certain background signal and noise. As long as this model captures all the important features of the physical process, observed intensities should match the predicted intensities at each voxel in each round and in each channel. To think about this more clearly, let’s define these predicted intensities as the ‘reconstruction’: To test our model, we can visually compare the reconstruction to the observed data. If the residual between the two includes significant highly-structured noise, then it is likely that we are missing important aspects of the data. S7–S9 Figs show the results of these comparisons, and the overall reconstruction error is approximately 7.64% (computed as the percentage unexplained) in the 180 × 200 region shown in S7(A)–S9(A) Figs. (Also see the Supplementary video, described in S1 Appendix, Section B.) They appear fairly promising, but certain structured features do appear in the residual. Most strikingly, we have found that a minority of rolonies ‘dropout’ for one or more rounds: a rolony may give a strong bright signal in most of the rounds but simply vanish in one round. Our current physical model does not accommodate this, and this limitation appears in the residual as bright and dark spots. However, as mentioned above and shown in the hybrid simulation data in Fig 3B, our method is robust to these ‘dropout’ effects; it is still able to capture the correct rolony positions when dropout occurs on a small number of rounds.

Diagnostics based on ‘cleaned’ images are useful to check the accuracy of BarDensr

The reconstruction is made up of many parts: it has the background component a, the per-round per-channel offset and scale terms (α, b), and rolony contributions arising from F, φ, Z. As shown above, it is straightforward to compare the total reconstruction to the observed data. However, this does not isolate the contributions of individual estimated rolonies.

Therefore we adapted a partial subtraction approach from [22]. We pick one barcode, j*, and focus only on the contributions to the reconstruction from this one barcode. In particular, we assume that all other aspects of the model are exactly correct. We assume that a, α, b, φ and Z are all exactly right. We further assume that F_j is exactly correct for every j ≠ j*. Assuming all these aspects of the model were perfect, we can look at what the data would have looked like if it had only included one type of barcode, namely j*. We call this counterfactual simulation the ‘cleaned image’: (2) This is the data with all aspects of the model subtracted away—except for the contributions from barcode j* (see Fig 7A as an example). The cleaned image for the barcode j* has much in common with the rolony density for j*. However, X^(j*) differs from F_j* in one crucial way. For each voxel m, F_j* gives exactly one value. However, for each voxel m, X^(j*) gives R × C values—one for each round and channel of the experiment. According to our model, however, it should be possible to express all these values in terms of multiplication in the following form: In this formula (known as an ‘outer product’) we see that (which varies across voxels, rounds, and channels) is the product of two objects: the rolony density (which varies across voxels) and the transformed barcode G (which varies across rounds and channels) for j*. This is actually a very strong assumption; most tensors would not exhibit this kind of structure. We can empirically check for this ‘rank-one’ structure by computing the singular value decomposition (SVD) of . If the SVD yields only one strong singular value, then can be well-approximated by this rank-one outer product, and furthermore the SVD yields the correct values for F_j*,m and G_j*,r,c. We can compare the values for these quantities (as returned by the SVD analysis) to the estimated values (as returned by BarDensr). We show some examples in Fig 7B comparing the estimated value of G_j* with SVD results (a similar but more complete set of the spots can be seen in S10 Fig). Note that the match isn’t quite perfect (the temporal singular vector of the corresponding cleaned images varies a bit from our estimate). In future work we hope to investigate whether these differences could be accounted for by a more accurate physical model. For now, we content ourselves that the method is accurate enough to provide a useful diagnostic for the detected rolonies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Using cleaned images and SVD to examine model fit quality and variability.

(A) SVD analysis example using one gene. Spots are identified in F_j* for each barcode j* using local-max-peak-finding. For the gene barcode j* (Deptor) shown here (top left), three spots with the highest accuracy are being analyzed. The right panel shows the zoomed-in R × C plots of the raw image X (top) and ‘cleaned’ image X^(j*) (bottom) at these three spot locations for barcode j*. Note that ‘cleaned’ images are significantly sparser than the raw images, as desired. We then applied SVD to the cleaned image X^(j*) at these three spot locations. The first two columns on the bottom left show the zoomed-in image of the original spot (KF)_j* and the learned weighted barcode matrix (G_j*) corresponding to this gene barcode j*. The top singular vectors are plotted in the last two columns (showing a good match with G_j* and the cropped (KF)_j*). R² is the squared correlation coefficient between X^(j*) and the outer product of these two singular vectors; the high R² values seen here indicate that the model accurately summarizes X^(j*). (B) Results of SVD analysis of cleaned images for the top high-R² spots. This plot summarizes the results of the analysis illustrated in (A). The first column shows (KF)_j* around the brightest spots; the second column shows the top spatial singular vectors for the same region, and the last column shows the top temporal singular vectors for these spots (the top row shows the scaled G_j* learned from the model, and the bottom row shows the corresponding top temporal singular vectors for these spots). Only six barcodes that are most abundant in the selected region are shown here; S10 Fig provides a more complete illustration.

https://doi.org/10.1371/journal.pcbi.1008256.g007

We can also use these cleaned images to help us compare BarDensr with other methods by eye. S11 Fig investigates cleaned images for gene Arpp19, comparing the results of our method to the hand-curated results. In cases where the results of the two approaches disagree, these cleaned images suggest that our results are often reasonable.