Overview of SCOPE (Spatial reConstruction via Oligonucleotide Proximity Encoding).

(A) Schematic of structure of oligos tethered to SCOPE beads. Hydrogel beads contain 5’-Acrydite-tethered DNA oligos that are either cleavable (“senders”) or non-cleavable (“receivers”) by USER enzyme, based on the presence or absence of dU in the 5’ stem. Both sender and receiver oligos contain the same barcode structure–4 subsequences each with 96 possible identities–depicted by the 4 colored boxes. Both sender and receiver oligos are also functionalized with a 3’ cap of either poly-dT (“decoders”, designed to capture poly-dA or poly-A tailed molecules) or a sequence encoding bead subtype (“messengers”, designed to chimerize with messenger oligos from proximally located beads). (B) Schematic of 2D SCOPE reaction. Top: Poly-dA hashing oligos are deposited onto a 2D array of SCOPE beads, bearing subsequences as “colors” that encode an image of interest. Decoder oligos form chimeras with proximally located hashing oligos. Bottom: Messenger oligos participate in diffusion-dependent, hybridization-mediated reactions with the messenger oligos of proximally located beads in the 2D array, resulting in sender-receiver messenger chimeras. (C) Top: Spatial mapping of hashing oligos. From sequencing of decoder-hashing oligo chimeras, we obtain a bead-by-molecule count matrix that is informative with respect to which hashing oligos are proximate to each bead barcode. Bottom: Inference of the relative positions of beads. From sequencing of messenger chimeras, we obtain a bead-by-bead count matrix that is informative with respect to which bead barcodes are proximate to which other bead barcodes. (D) Out of 20 designed pairs of complementary messenger sequences, the messenger oligos of each bead subtype are capped with 19 “antisense” sequences and 1 “sense” sequence. For example, the messenger oligos of “A” subtype beads include the “A” hybridization sequence in sense orientation, with the other 19 hybridization sequences in antisense orientation. (E) Illustration of how 3 interaction counts between a proximally located pair of “A” and “T” subtype beads might originate from different kinds of sender-reciever chimeras. (F) Only beads of different subtypes can contribute to the generation of sender-receiver chimeras, which avoids “self-reactions” (i.e. senders chimerizing with receivers from the same bead) but runs the risk of missing interactions between closely located beads of the same subtype.

Bead arrays form local communities that can be computationally reconstructed from pairwise interaction matrices.

(A) Segmented mask of a representative image of a SCOPE bead array. (B) Delaunay triangulation for segmented bead centroids. (C) Density plot of bead diameters from the segmented image. (D) Distribution of the number of nearest neighbors observed for each bead. (E) Raw, unique interaction counts (y-axis) of a representative bead with all its interacting neighbors ranked in descending order (x-axis). Left plot (blue) shows counts for a representative bead from a SCOPE reaction performed on beads in suspension. Right plot (red) shows counts for a representative bead from a SCOPE reaction performed on beads in a 2D array. (F) Mean cumulative distribution function (CDF) of proportion of all interaction counts explained by the top neighbors of a given bead, for SCOPE reactions performed on beads in suspension (blue) or 2D array (red). Gray shadings show the range of CDFs for the 1000 randomly sampled beads from which the mean CDFs were derived. (G) Network of interactions of two representative beads (red) from SCOPE reaction performed in 2D array, with and among their top ten neighbors. (H) Proportion of read counts from interactions between each bead of interest and their neighbors are overlaid with 100 simulation instances from the inverse square model. (I) Region of SCOPE array imaged after staining for two specific “fiducial” sub-barcodes at position 2 (red) or position 4 (blue). Scale bar denotes 100 µm. (J) SCOPE reconstruction highlighting the positions of sub-barcodes matching the FISH probes used for imaging. (K) Point matching search performed between segmented imaging fiducials and SCOPE reconstruction fiducials. Left: Purple-colored points (beads in SCOPE reconstruction bearing sub-barcode BC2-#33) were scaled, translated, rotated and point-matched onto red-colored centroids from imaging. Right: The same transformation was applied to green-colored points (beads in SCOPE reconstruction bearing sub-barcode BC4-#91) and these were matched to blue-colored centroids from imaging. (L) Mean point matching error for both sub-barcodes during gradient descent optimization. Rigid transformations were optimized according to the error for beads with sub-barcode BC2-#33 (red), with sub-barcode BC4-#91 positions (green) held out.

Optics-free reconstruction of 2D shapes with SCOPE.

(A) Single channel image of a manually cut circular bead array with a 7 mm diameter subjected to hybridization of a sequencing primer followed by single nucleotide extension with dye-conjugated, reversibly terminated dNTPs. Scale bar denotes 1 mm. (B) Computational reconstruction of the circular bead array shown in panel A with SCOPE. Beads are colored by cluster labels, computed by Leiden clustering on the bead-bead interaction graph. (C) Single channel image of a bead array manually cut to an asymmetric “swoosh” shape and stained with SYBR Gold. Scale bar denotes 1 mm. We estimate that the swoosh shape occupies 44 mm2. (D) Computational reconstruction of the asymmetric swoosh-shaped bead array shown in panel C with SCOPE. Beads are colored by cluster labels, computed by Leiden clustering on the bead-bead interaction graph. (E) Pairwise mean Euclidean distance between each of 10 independent runs of the reconstruction algorithm on the circular bead array after alignment to a common reference frame. Units are in bead diameters. (F) Same as panel E, but for asymmetric swoosh-shaped bead array. (G) Mean Euclidean bead positioning error for reconstruction of the circular bead array with downsampling of messenger UMIs per bead. Error bars correspond to the range, and points to mean, for 5 reconstructions of each downsampling trial. (H) Same as panel G, but for asymmetric swoosh-shaped bead array.

An eye examination for SCOPE.

(A) Composite fluorescent image of a Snellen eye chart for visual acuity printed onto a 2D SCOPE array using a microarray printer. Of the 13 printed oligonucleotides (“paints”), three were fluorescently labeled—either via dye conjugation (Cy3, Cy5) or inclusion of DAPI in solution—allowing direct microscopic visualization. Insets show the “E” (top), “O” (middle), and “D” (bottom) for reference. (B) Maximum-intensity projection of the reconstructed array comprising 1.44 million beads, colored by summed decoder counts across all paints. (C) Maximum intensity projection of a single oligo channel (DAPI channel oligo) from the “E” in Row 1 and the “Z” in Row 3. The region corresponding to each letter was cropped from the full reconstruction, rotated, and displayed at different scales with a common ruler for reference. (D) Locally binned image showing the sum-of-squared oligo counts per bin, with horizontal (top) and vertical (left) line intensity profiles corresponding to the red crosshairs. (E) Dimensional fidelity of the reconstructed “E.” Plotted points show the ratio between each of 10 feature lengths indicated by bars at the right and the length of the base of the E, measured via optics (“imaging”) or SCOPE (“reconstruction”) using ImageJ.

UMAP reconstructions of 3D volumes with SCOPE.

(A) Brightfield image of a bead-embedded polyacrylamide gel matrix from a teddy bear mold. Yellow scale bar denotes 1 mm. (B) Binary masks of brightfield images of the silicone molds used to create bead-embedded 3D gels for star, butterfly and letter E-shaped molds. The depth of each mold is 2 mm. Black scale bars denote 1 mm. (C) Top-views of 3D UMAP reconstructions of the bead-embedded gels from the bear, star, butterfly, and letter E-shaped molds, using n_neighbors=15, min_dist=0.4, and 1,000 training epochs on the interaction matrix. Each barcode was normalized to 10,000 total counts per barcode and then log-transformed. Points are bead barcodes in 3D space overlaid on a concave hull (white) calculated by the Python package “alphashape” and colored by the log10 UMI count of each barcode. Reconstructions have been rotated to show the front view. (D) Same as panel C, but side-views. (E) Z-axis cross-sections of 3D UMAP embeddings from bead-embedded hydrogel reconstructions. Each slice shows reconstructed bead positions at successive depths in the inferred volume. The absence of points in central regions highlights the hollow character of computational reconstructions, consistent with reduced recovery of interior beads due to diffusion or steric limitations during SCOPE reaction.

Barcoded hydrogel bead generation and bead array fabrication.

(A) Image of the flow focusing microfluidics generator while producing polyacrylamide hydrogel beads. (B) 10% polyacrylamide gel of oligos released by USER-mediated cleavage from the initial bead (round 0) vs. after each round of DNA synthesis by ligation of the barcode (rounds 1, 2, 3, 4) or functional cap (Cap). Left: DNA size ladder (bp). In the lane corresponding to the initial bead (round 0), the band labeled “stub” corresponds to what is cleaved off the strand that was initially incorporated into the hydrogel bead, while the the band labeled “splint” is a primer that anneals to the stub to create a 4-bp overhang handle for the ligation of the first of four splint barcodes (36). (C) 96 sub-barcodes were used in each of four rounds of the split-pool procedure for barcode generation (36). Shown in the box plot is the percentage of each sub-barcode incorporated to full barcode at each position. The red line at 1.04% corresponds to uniform incorporation (1/96). (D) Based on amplification and sequencing of UMIs associated with barcodes enzymatically released from individual beads, together with the fact that only ∼50% of the barcodes have a dU at their stem and are expected to be cleaved, we estimate that each bead bears 491,811 (interquartile range: 308,757-833,744) functionalized, barcoded oligos. (E) Schematic of a monolayer of hydrogel beads inside a 6% (v/v) PAA (polyacrylamide) encasing gel. (F) Differential interference contrast microscopy image of beads in an encasing gel. Yellow scale bar: 100 µm.

Complementarity and predicted outcomes of using twenty pairs of messenger sequence designs.

(A) Fractions of oligos predicted to be bound at equilibrium in a pool of 40 designed messenger cap sequences (20 messenger caps x 2 orientations) at 37°C as estimated by Nupack (42). Complementary strands are adjacent in numbering (e.g. 1/2, 3/4 … 39/40). (B) We simulated the random assignment of a given number of bead subtypes to beads arrayed in a hexagonal lattice and computed, for any given bead, the proportion of immediately neighboring beads that are of the same subtype (i.e. homotypic neighbors).

Reconstructing simulated SCOPE bead arrays.

(A) Mean Euclidean error of UMAP (green) or t-SNE (purple) applied to simulated SCOPE bead arrays of different sizes. Mean Euclidean error is measured as the average difference in distance between an inferred position and the ground truth simulated position after a linear point cloud registration and point matching. The line shows the median of 5 trials and the shaded areas show the interquartile range. Units are in bead diameters. (B) A representative simulated rectangular array of 2500 beads reconstructed with UMAP. Lines indicate matching between ground truth and UMAP inferred positions after linear point cloud registration. Bead positioning errors are greater at the boundary of the shape. (C) Mean Euclidean error of a simulated and reconstructed 40,000 bead array as a function of simulated sequencing depth. Each point is the median of 10 trials. The shaded blue area shows the interquartile range for the set of trials.

Doublet detection by community detection.

(A) Doublets are defined as beads that have two spatial positions but share an identical DNA barcode (ground truth – left). In sequencing data these doublets appear as a single bead with collapsed interaction networks. (B) Expected rates of doublet occurrence for arrays of a given size. Estimated by computing the birthday problem for a range of array sizes with a fixed set of 964 barcodes. (C) Doublets detected in simulation using Leiden clustering and a range of cluster resolutions. Recall and precision are shown with error bars indicating standard deviation. (D) Examples of detected doublets during reconstruction of the circular array. Highlighted in red and green are the two communities detected for a single bead (shown in yellow).

Spotted poly-dA fiducials enable orientation of SCOPE reconstructions of symmetric 2D arrays.

Poly-dA–tailed hash oligonucleotides were manually spotted onto the circular bead array at the 12 o’clock, 3 o’clock, and 6 o’clock positions prior to the SCOPE reaction. During decoding, complementary poly-dT–capped decoder molecules captured these fiducials, producing locally concentrated clusters of corresponding barcoded sequences. In the reconstructed array (three versions of which are shown here without rotation or reflection relative to one another), the three fiducial barcodes are locally concentrated at expected relative positions (red dots), confirming correct geometry and providing reference points for orienting symmetric reconstructions. X0 and X1 denote the spatial axes outputted by the reconstruction pipeline.

Ten iterations of computational reconstruction of circular and swoosh arrays for (A) circular bead array or (B) asymmetric “swoosh” array.

The pipeline was re-run independently for each reconstruction with identical inputs and initial hyperparameters each time, only varying the random seed. Each reconstruction has been transformed (including rotation, reflection, scaling) for visualization. X0 and X1 denote the spatial axes outputted by the reconstruction pipeline.

Initial 2D UMAP result for Snellen eye chart experiment prior to 3D flattening.

Maximum-intensity projection of the reconstructed Snellen eye chart generated using the standard 2D UMAP pipeline applied directly to the bead–bead interaction matrix. Each point represents a DNA-barcoded bead (n = 1.44 million), colored by the summed decoder counts across all 13 printed poly(A)-tailed oligonucleotide “paints.” Although major features of the eye chart are discernible, the reconstruction exhibits macroscopic warping and curvature, likely reflecting non-uniform diffusion during the SCOPE reaction across this large array.

Representative results from a parameter grid search performed on the bead–bead interaction matrix derived from the butterfly-shaped hydrogel mold.

The matrix was filtered to include barcodes with at least 1,000 UMIs and reconstructed using a fixed number of training epochs (n = 1,000) and nearest neighbors (n = 15). Columns show different normalization strategies: (left) raw count matrix; (middle) normalized count matrix scaled to 10,000 total counts per barcode and log-transformed; (right) matrix normalized by dividing each entry by its row and column sums, then symmetrized by averaging with its transpose. Rows correspond to varying min_dist values.

Representative results from a parameter grid search performed on the bead–bead interaction matrix derived from the block letter E-shaped hydrogel mold.

The matrix was filtered to include barcodes with at least 1,000 UMIs and reconstructed using a fixed number of training epochs (n = 1,000) and nearest neighbors (n = 15). Columns show different normalization strategies: (left) raw count matrix; (middle) normalized count matrix scaled to 10,000 total counts per barcode and log-transformed; (right) matrix normalized by dividing each entry by its row and column sums, then symmetrized by averaging with its transpose. Rows correspond to varying min_dist values.

Lowering UMI cutoffs for barcode inclusion in 3D reconstructions leads to spurious barcode clusters.

SCOPE sequencing libraries were filtered for barcodes with at least 100 UMIs (left) vs. 1,000 UMIs (right) and reconstructed in 3D using the UMAP algorithm (1,000 training epochs, n_neighbors = 15, min_dist = 0.4) on the normalized count matrix.