Barcoded multiple displacement amplification for high coverage sequencing in spatial genomics

Abstract

Determining mutational landscapes in a spatial context is essential for understanding genetically heterogeneous cell microniches. Current approaches, such as Multiple Displacement Amplification (MDA), offer high genome coverage but limited multiplexing, which hinders large-scale spatial genomic studies. Here, we introduce barcoded MDA (bMDA), a technique that achieves high-coverage genomic analysis of low-input DNA while enhancing the multiplexing capabilities. By incorporating cell barcodes during MDA, bMDA streamlines library preparation in one pot, thereby overcoming a key bottleneck in spatial genomics. We apply bMDA to the integrative spatial analysis of triple-negative breast cancer tissues by examining copy number alterations, single nucleotide variations, structural variations, and kataegis signatures for each spatial microniche. This enables the assessment of subclonal evolutionary relationships within a spatial context. Therefore, bMDA has emerged as a scalable technology with the potential to advance the field of spatial genomics significantly.

Fig. 1. Barcoded MDA (bMDA) enables multiplexed preparation of a single-cell genome sequencing library, realizing cost-effective and higher-throughput single-nucleotide resolution single-cell genome analysis.

a bMDA-seq can expand our understanding of heterogeneous cell populations by allowing a large number of single cells to be analyzed in multiplex and at single-nucleotide resolution. Post bMDA-seq, each single-cell data showed a genome coverage depth that was sufficient to perform integrative spatial genomics. b Conventional multiple displacement amplification (MDA) involves preparing a sequencing library individually for each single-cell because a cell barcode is incorporated at the end of the library preparation step (marked as *). c Schematic of bMDA and its workflow bMDA-seq. Since bMDA uses a barcoded primer instead of the conventional random hexamer. The cell barcode used for sample demultiplexing is incorporated during the bMDA reaction (marked as *). After pooling the barcoded MDA products, a sequencing library can be prepared in a single reaction tube (one-pot), thereby reducing the number of library preparation reactions according to the barcoding capacities. d Schematic illustration of bMDA-seq workflow. After performing MDA using the barcoded primer, bMDA products were pooled in a single reaction tube. Pooled bMDA products were first fragmented to obtain a desired library insert size. Then, only barcoded DNA fragments were enriched using streptavidin–biotin interactions. Finally, ligation-based library preparation was performed in one-pot to obtain the final NGS sequencing library. e, f Detailed procedures of conventional MDA and bMDA show that the multiplexing capability of bMDA remarkably reduces reagent costs and labor required for integrative spatial genomics.

Fig. 2. bMDA-seq was optimized and validated with cell line gDNA.

a Schematic illustration of the equivolume pooling approach, barcode bias, and its correction. b NGS read count before and after barcode bias correction shows that the barcode bias can be successfully corrected by fine-tuning the proportion of barcoded primers. The NGS read count of a specific barcode was normalized to ~1 to obtain the normalized NGS read counts (y-axis). Points and lines represent the mean ± s.e.m. (n = 3 independent experiments). c Barcoding status of the bMDA-seq library was analyzed using NGS, showing that the majority of NGS reads (82.5%) were barcoded as predicted. Other unpredicted NGS reads were classified by their barcode sequences, and are shown as a bar chart. To eliminate the putative confounding effect that may arise after the pooling of differently barcoded bMDA products, a single bMDA product was used for the analysis. The ratios from each bMDA product (n = 8) were averaged to obtain a representative single plot. d The barcode swapping ratio of bMDA-seq was sufficiently lower than the other well-known index swapping ratios, such as index hopping in multiplexed Illumina sequencing. To assess the barcode swapping ratio, 47 out of 48 barcodes were used for amplifying the human genome (HL-60), and the remaining barcode was used for amplifying the mouse genome (NIH3T3). Notably, of the ₄₈C₂ possible barcode swapping cases, the probability of swapping from human to mouse is 47 times higher than the opposite due to its higher number of possible swapping cases. Number 1 vs. number 47 was chosen to simulate the worst-case scenario. Box plot show the median (center line), first and third quartiles (box edges), while the whiskers extend from the box edge to the largest or smallest value no further than 1.5 times the interquartile range (IQR) from the box edge (n = 6 independent experiments). Source data are provided as a Source Data file.

Fig. 3. Technical performance of single-cell bMDA is comparable to that of conventional MDA.

a Single cells were isolated by the PHLI-seq platform. b Single-cell bMDA showed nearly full genome coverage breadth (0.87×) in contrast to the Tn5-based method (0.22×). bMDA also showed similar coverage breadth compared to the in-house MDA, whose protocol was the same as bMDA except for using N6 as a primer (MDA), or MDA using a commercialized kit (MDA kit). c Targeted deep sequencing confirmed that the library complexity of bMDA is high enough for single-nucleotide resolution genome analysis, covering the single-cell whole genome at a depth of more than 500×. d, e Amplification uniformity of bMDA and conventional MDA was similar in both the 66 pg (p = 0.1 and 95% confidence interval (CI) −0.005, 0.0006, two-sided Wilcoxon rank-sum test) and single-cell groups (p = 0.14, t = −1.6, degree of freedom (df) = 9, 95% CI −0.10, 0.02, and Cohen’s d = 1.00, unpaired two-sided Student’s t test). Box plot show the median (center line), first and third quartiles (box edges), while the whiskers extend from the box edge to the largest or smallest value no further than 1.5 times the interquartile range (IQR) from the box edge. n = 3, 98, 11, 3, 7, 4, 11, and 2 biologically independent samples from the left of the box plot. f Copy number alterations (CNA) plot demonstrates that there is no notable difference between bMDA and other conventional MDA methods in the aspect of resolving CNAs. Each row indicates individual MDA-amplified products obtained by different methods and templates. The values displayed on the heatmap represent the log2-transformed relative changes in copy number compared to the average copy number across the entire genome. g To confirm the single nucleotide variant (SNV) detection performance, the allelic dropout (ADO) rate and false-positive mutation detection rate (FPR) were evaluated. The bMDA, MDA, and MDA kit results showed comparable performance. bMDA barcoded MDA, MDA in-house MDA with random hexamer, MDA kit commercial MDA kit. Source data are provided as a Source Data file.

Fig. 4. Structural variation and kataegis further delineate clonal evolution of T1 triple-negative breast cancer tissue section than when the subclonality is determined by copy number alterations and single nucleotide variations.

a Isolating spatially adjacent cell clusters provides highly confident integrated genomics. b SNV detection sensitivity was higher and less biased for cell clusters compared to single-cell. The SNV sensitivity of cell clusters was calculated by comparing the germline mutations identified in the bMDA-seq data from T1 and T2 tumors to the germline mutations detected in the normal bulk sequencing data (n = 27 biologically independent samples). The value for single-cell was obtained by performing MDA/bMDA (n = 10 biologically independent samples) and Tn5-based methods (n = 6 biologically independent samples) on a cell line. c CNA detection specificity was higher and less biased for cell clusters (n = 10 biologically independent samples) compared to the single cells (n = 26 biologically independent samples). All box plots show the median (center line), first and third quartiles (box edges), while the whiskers extend from the box edge to the largest or smallest value no further than 1.5 times the interquartile range (IQR) from the box edge. d Variant allele frequency distribution of detected heterozygous SNPs was less biased and resembled more to bulk distribution in cell cluster isolation. e Spatial landscape of different microniches of a triple-negative breast cancer tissue section. Subclone information was inferred by analyzing SNVs. The gray color indicates that spatial microniches were not included in the whole exome sequencing analysis (scale bar = 1 mm). f Heatmap illustrating copy number alterations in different microniches. The values displayed on the heatmap represent the log2-transformed relative changes in copy number compared to the average copy number across the entire genome. g Single nucleotide variations of different subclones. h Structural variations of the different subclones within the same tumor. i Genes where kataegis was detected are displayed. Source data are provided as a Source Data file.

Fig. 5. Integrative spatial genomics in T2 breast cancer.

The color and order of each spatial microniche remain consistent throughout the figure. a Spatial landscape of different subclones of a triple-negative breast cancer tissue section. Subclone information was inferred through phylogenetic analysis of SNVs (scale bar = 1 mm). b Copy number alterations of different subclones. The values displayed on the heatmap represent the log2-transformed relative changes in copy number compared to the average copy number across the entire genome. c The phylogenetic tree illustrates the evolutionary relationships among different microniches within the same tumor microenvironment. d Single nucleotide variations of different microniches. e Structural variations of the different microniches within the same tumor. Source data are provided as a Source Data file.

Fig. 6. Integrative spatial genomics of T2 tumor reveals the relationship between kataegis, copy number, and translocation in chromosome 12, thereby providing comprehensive insight into the spatial subclone in cancer.

a Spatial landscape of populations with and without kataegis in chromosome 12. b, c Translocations in two different populations are displayed. d, e Rainfall plots of the two different populations and copy number alterations in chromosome 12. The region where kataegis occurred had copy number amplification. f, g Structural variations in Chromosome 12 and 5. h, i Schematic display of genomic aberrations in chromosomes 12 and 5. Source data are provided as a Source Data file.

Fig. 7. bMDA facilitates integrative spatial genomics, suggesting a more plausible tumor evolutionary history of T2 triple-negative breast cancer.

Source data are provided as a Source Data file.