Full Text

Turn on search term navigation

ARTICLE

Received 5 Aug 2014 | Accepted 3 Mar 2015 | Published 16 Apr 2015

Cheng-Zhong Zhang1,2,*, Viktor A. Adalsteinsson2,3,4,*, Joshua Francis1,2, Hauke Cornils5,6, Joonil Jung2, Cecile Maire1, Keith L. Ligon1,7,8,9,10, Matthew Meyerson1,2,7,11 & J. Christopher Love2,3,4

Artifacts introduced in whole-genome amplication (WGA) make it difcult to derive accurate genomic information from single-cell genomes and require different analytical strategies from bulk genome analysis. Here, we describe statistical methods to quantitatively assess the amplication bias resulting from whole-genome amplication of single-cell genomic DNA. Analysis of single-cell DNA libraries generated by different technologies revealed universal features of the genome coverage bias predominantly generated at the amplicon level (110 kb). The magnitude of coverage bias can be accurately calibrated from low-pass sequencing (B0.1 ) to predict the depth-of-coverage yield of single-cell DNA

libraries sequenced at arbitrary depths. We further provide a benchmark comparison of single-cell libraries generated by multi-strand displacement amplication (MDA) and multiple annealing and looping-based amplication cycles (MALBAC). Finally, we develop statistical models to calibrate allelic bias in single-cell whole-genome amplication and demonstrate a census-based strategy for efcient and accurate variant detection from low-input biopsy samples.

1 Department of Medical Oncology, Dana-Farber Cancer Institute, Boston, Massachusetts 02215, USA. 2 Cancer Program, Broad Institute of Harvard and MIT, Cambridge, Massachusetts 02142, USA. 3 Department of Chemical Engineering Cambridge, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA. 4 Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA. 5 Department of Pediatric Oncology, Dana-Farber Cancer Institute, Boston, Massachusetts 02215, USA. 6 Department of Cell Biology, Harvard Medical School, Boston, Massachusetts 02115, USA. 7 Department of Pathology, Harvard Medical School, Boston, Massachusetts 02115, USA. 8 Department of Pathology, Brigham and Womens Hospital, Boston, Massachusetts 02115, USA. 9 Department of Pathology, Boston Childrens Hospital, Boston, Massachusetts 02115, USA. 10 Center for Molecular Oncologic Pathology, Dana Farber Cancer Institute, Boston, Massachusetts 02115, USA. 11 Center for Cancer Genome Discovery, Dana Farber Cancer Institute, Boston, Massachusetts 02215, USA. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to M.M. (email: [email protected]) or to J.C.L. (email: [email protected]).

NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications 1

DOI: 10.1038/ncomms7822

Calibrating genomic and allelic coverage bias in single-cell sequencing

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822

Single-cell sequencing has provided unique insights into the genetic diversity of living organisms and among different cells within the same individual13. Recent single-cell

analyses have uncovered different clonal populations within a single tumour4,5, revealed genomic diversity in gametes6,7 and neurons8,9, and resolved historical cellular lineages during development10,11. Single-cell sequencing also has many potential clinical applications, such as characterization of circulating tumour cells12,13 or ne-needle aspirates for clinical diagnostics.

A major drawback of single-cell sequencing, however, is the need to amplify genomic DNA before genomic characterizations1417. Owing to the limited processivity (o100 kb) and strand extension rate (o100 nt per second) of DNA polymerases, the amplication of large genomes requires priming and extension at millions of loci, each amplied 10,000- to 1,000,000-fold. Such a large number of polymerase reactions inevitably generate amplication errors that confound the detection of genetic variants (Supplementary Fig. 1). Furthermore, differential priming efciencies and extension rates result in uneven amplications across the genome18,19 and skewed representations of homologous chromosomes. These variations both compromise variant detection sensitivity and may lead to incorrect genotypes5,12. Although technological innovations may improve the delity of whole-genome amplication (WGA)1517,2023, statistical uctuations in the amplications of millions of different DNA templates will persist.

As genetic variants are detected by the relative abundance of variant-containing DNA templates in the library, non-uniformity in genome coverage directly impacts the sensitivity to detect variants. For example, grossly non-uniform libraries emphasize only overrepresented regions of the genome, and contain little information on other regions. Current methods to assess the uniformity of WGA rely on either direct visual inspection or various statistical measures of the sequencing coverage at the base level18,22 or the allele level5,12. These empirical methods and metrics generally require substantial sequencing (10 or greater)

and only gauge the deviation of amplied DNA from the uniform bulk DNA at a particular sequencing depth. They fail, however, to characterize the intrinsic non-uniformity resulting from WGA that is independent of sequencing depth (Fig. 1a,b). Moreover, the nature of the main sources of bias remains poorly characterized (Fig. 1c).

Here, we report a systematic analysis of the coverage bias in single-cell whole-genome amplication. We show that the structure of individual WGA amplicons imparts a dominant amplication bias on length scales longer than the average size of sequencing fragments. Sequencing at low depths (0.11 ) can

effectively reveal this variation in the amplicon-level coverage and enable accurate predictions of the depth-of-coverage yield when sequencing single-cell libraries to arbitrary depths. We further characterized the amplication bias between homologous chromosomes using analytically solvable models and validated these model predictions of allelic coverage by experimentally

Library complexity and sequencing yield

Bulk DNA libraries Whole-genome amplified single-cell DNA libraries Original genome

DNA library

Depth

Coverage bias at different levels

Representations of DNA sequences

Variation in genome coverage

Original genome

WGA amplicons

Sequencing fragments

Uniform

Over amplified

Under amplified

Amplification bias

Sequencing bias

10 ~ 50 kb for MDA

10~50 kb for MDA

Variation at different scales

~1 kb

Recurrent and random amplification bias

WGA1

WGA2

WGA3

Recurrent amplification bias Random amplification bias

WGA1

WGA2

WGA3

Figure 1 | Non-uniformity in genome coverage and its impact on the sequencing yield. (a) Dependence of the information yield on the sequencing depth. Deeper sequencing of bulk libraries yields information on a larger population of cells; deeper sequencing of whole-genome amplied single-cell libraries reveals information on a larger fraction of the genome (thick lines). (b) Genome coverage bias at different levels. Amplication bias (top): whole-genome amplication generates coverage bias at the amplicon level, which is B1050 kb for multi-strand displacement amplication. Sequencing bias (bottom): non-uniformity in the selection of sequencing fragments can be caused by multiple sources of bias including whole-genome amplication: the variation in sequencing coverage can be observed from 100 bp to multiple megabases. (c) Schematic representations of recurrent and random amplication bias from multiple independent amplications of the same DNA material.

2 NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822 ARTICLE

observed coverage at heterozygous sites. These results provide a framework for quality assurance of single-cell libraries and for estimating the sensitivity to detect local variantssuch as single-nucleotide variants or chromosomal translocationspresent in an individual cell at a given sequencing depth. Finally, we demonstrate that the amplication bias in multi-strand displacement amplication (MDA) is more random than recurrent. Although such random bias cannot be corrected systematically, it suggests an efcient census-based strategy to accurately determine somatic genetic variants in small biopsy samples by sequencing multiple single cells from the same sample at modest depths.

ResultsInformation yield from bulk and single-cell sequencing. In bulk DNA libraries, each sequencing fragment represents genomic information from an individual cell; therefore, the information content increases with the sequencing depth until fragments are sequenced to exhaustion. The information content of a DNA library (library complexity) is thus measured by the total number of distinct molecules (sequencing fragments) in the library2426. This measure is essentially determined by the total number of cells (or the total amount of genomic DNA) used to prepare the library (Fig. 1a, left panel). In single-cell DNA sequencing, WGA precedes the construction of a DNA library and introduces non-uniformity across the genome: As sequencing depth increases, more genomic regions are uncovered (Fig. 1a, right panel). Hence the fraction of the single-cells genome uncovered at a given sequencing depth determines the information content of single-cell sequencing. This measure ultimately depends on the uniformity of genome coverage, or the magnitude and spread of whole-genome amplication bias, and is conceptually equivalent to a single-cell DNA library complexity.

Amplicon-level bias dominates coverage variation. Visual inspection of single-cell sequencing coverage suggests that the genome coverage varies at many different length scales (Fig. 1b). To systematically evaluate the amplication bias in single-cell libraries, we sequenced MDA-generated DNA libraries of diploid RPE-1 cells (510 ) and compared the sequencing

coverage to a matched, unamplied bulk DNA library (B12 ).

To eliminate the effects of sequencing depths, we computationally down-sampled the bulk and single-cell DNA libraries and calculated the auto-correlation of base-level coverage in disomic chromosome 1 at various depths to examine coverage correlations at all length scales (Fig. 2a, Supplementary Fig. 2). Both bulk and MDA libraries exhibited a correlation at length scale lcE100 bp, reecting the sequencing read length (101 bp). Looking more closely, we also identied a correlation at lcE250 bp, corresponding to the average size of the paired-end fragments (Supplementary Fig. 2). As expected, the magnitude of such correlations at the fragment scale decays with increasing sequencing depth.

Besides the fragment-level correlations, the bulk DNA sequencing coverage showed minimal correlation between loci separated by more than 1 kb. In contrast, single-cell libraries exhibited a prominent correlation in 1100 kb that is independent of the sequencing depth. Independent sequencing of the same single-cell library to 0.1 on the Illumina MiSeq platform and to

9 on the HiSeq platform revealed the same correlation with a

characteristic length lcE33 kb (Fig. 2a). The sequencing depth-independent correlation reects the intrinsic non-uniformity in the DNA library and suggests a characteristic length scale of amplication bias.

The predominant correlation at lc suggests that adjacent loci within this distance have comparable coverage. This observation implies that the primary source of coverage variation (or amplication bias) is at or above the distance lc. Therefore, statistical variation of coverage at the single-base level should reect coverage variation at the amplicon level. To test this hypothesis, we computed the cumulative distribution of bin-level coverage (bin sizeE17 Kb, half of lc). Normalizing the bin-level coverage by the mean depth-of-coverage, we found the cumulative distribution of bin-level coverage to be nearly identical between independent sequencing at 9 or at 0.1

(Fig. 2b), conrming that the amplicon-level coverage variation is intrinsic to the amplied DNA but independent of the sequencing depth. Furthermore, the cumulative distribution of single-base coverage at 9 sequencing depth aligned with the bin-level

coverage (Fig. 2b, Supplementary Fig. 2), suggesting that the amplicon-level variation was indeed the dominant source of nonuniformity in single-cell libraries.

To further validate this conclusion, we computed the depth-of-coverage curves and the Lorenz curves for the bulk RPE-1 library and a single RPE-1 library by MDA at different bin sizes (Supplementary Fig. 3). For the bulk library, the distribution of single-base-level coverage is indistinguishable from that evaluated at the bin level when the bin size is smaller than the fragment size (B300 bp); above this scale, the bin-level distribution is more uniform than the single-base level distribution, reecting smoothing of coverage non-uniformity.

By contrast, for the MDA-generated library, the distribution of single-base-level coverage remains constant until the bin size exceeds the amplicon size, B10 kb. Characterization of coverage non-uniformity by Lorenz curves22 also conrmed that the same bias was observed for bin sizes less than or comparable to the amplicon size and was independent of the sequencing depth. In particular, at sequencing depths oo1 , the majority of the

genome is uncovered and shows no variation in the single-base-level coverage; amplication bias, however, is manifested in the correlation between covered loci and can be evaluated by low-pass sequencing. For typical MDA-generated libraries, the amplicon size is on the order of 10 kb, hence, at 0.1

sequencing depth, there are 0.1 104/100E10 reads (assuming

100 bp single-end reads) on average for each amplicon. As long as the number of reads per amplicon is much larger than the statistical variation due to random selection in sequencing (for example, assuming poisson distribution, the standard deviation of the observable is given by the square root of the expectation), the percentage of such amplicons can be accurately calculated. At0.1 sequencing, the amplicon-level coverage can accurately

predict the fractional genome coverage down to 0.1 mean

depth, when there is approximately one read for each of these under-represented amplicons; below this depth, low-pass sequencing at 0.1 cannot distinguish between regions that

are severely under-amplied (o0.1 mean depth) and those

that dropped out of amplication.

Magnitude of amplicon-level variation determines coverage. We tested the validity of the correlation analysis by analysing DNA libraries generated from different types of cells and by different amplication technologies. For this purpose, we analysed single-cell sequencing data of additional RPE-1 samples (Supplementary Fig. 2) and data from multiple published studies, including frozen glioblastoma nuclei27 (Supplementary Fig. 4), single diploid lymphoblastoid cells5 (Supplementary Fig. 5), frozen single neuron nuclei8 (Supplementary Fig. 6), single sperms6 (Supplementary Fig. 7) and SW480 tumour cells22 (Supplementary Fig. 8); all samples were amplied by MDA.

NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822

Correlation at thefragment level MDA

Bulk

Bin-level coverage

101

Correlation at the amplicon level

0.10.250.5

4 8

1 0.1 MDA 9 MDA

9x MDA

Locus-level coverage

Normalized auto-correlation

0.8

Fraction of genome

~ 0.17 + 2e /lc

lc 33 kb

0.6

100

0.4

0.2

10 100 102 104 106

102 101 100 101

Separation between loci

Normalized depth

0.6

Glioblastomas (epicentre) Neurons (epicentre)

YH1 (Repli-G 12) Sperms (Repli-G 12)

0.5

Single RPEs (Repli-G 14)

SW480 (MALBAC) SW480 (MDA, GE)

Fraction of genome covered

0.8

RPE single (Repli-G 14)

Glioblastoma (epicentre)

SW480 (MALBAC)

Neuron (epicentre)

0.5

Fraction of genome

0.4

YH1 (Repli-G 12)

0.4

0.6

0.3

0.3 0 1 2 3

0.4

0.2

0.1

0.2

0.861.2 + x

40 50

0 10

y =

20 30

101 100 101

Correlation magnitude

Normalized depth

Figure 2 | Statistical analysis of whole-genome amplication bias and coverage uniformity. (a) Autocorrelation in the genome coverage of a two-cell RPE-1 DNA library (RPE#1) amplied by multi-strand displacement amplication (MDA). The same library independently sequenced to 0.1 (open

triangles) and to 8 (solid triangles) and exhibits a correlation above 1 kb that is invariant at intermediate depths (shaded triangles) from downsampling of

the 9 sequencing data. Black-dashed curve represents exponential tting of the autocorrelation in the 1100 kb range as 2 0.17e D/lc with a correlation

length lc 33 kb (95% condence interval: 2742 kb). This correlation is absent in the bulk library sequenced to different depths. Both the bulk and the

MDA-generated libraries show a sequencing-fragment-level correlation (lc 100 bp) that decays with the sequencing depth. (b) The identical normalized

cumulative coverage at bin size 1/2 lc evaluated from the 9 (solid) and from the 0.1 sequencing (dashed) reects the same amplicon-level variation

due to MDA. The agreement between bin-level (dashed and solid lines) and base-level (red dots) depth-of-coverage curves further suggests that the bin-level variation contributes the dominant amplication bias. See Supplementary Figs 2 and 48 for more examples of the correlation (a) and coverage (b) analysis of single-cell sequencing data from different studies. (c) Relationship between genome coverage (% covered at 1 mean sequencing depth) and

amplication bias (measured by the amplitude of the amplicon-level correlation) of single-cell libraries from different studies. Coverage is evaluated at Chr. 1 for both haploid sperms and diploid cells, as well as the SW480 tumour cells (disomic in Chr. 1), and at Chr. 10 (monosomic), Chr. 12 (disomic) and Chr. 13 (disomic) for glioblastoma nuclei. The inverse dependence is tted with an empirical formula, y 0.86/(1.2 Ox) (R2 0.98). (d) Comparison of the

cumulative coverage in the most uniform single-cell library from each study. Data were directly evaluated from high-depth sequencing of all samples except the neuron library for which the curve was interpolated from 0.5 sequencing as in b.

The SW480 cells were also amplied by quasi-linear multiple annealing and looping-based amplication cycles (MALBAC). The amplicon size in MDA-generated libraries ranged from 5 to 50 kb, with the sperm libraries having the lowest lcE5 kb (Supplementary Fig. 7). Interestingly, MDA of hundreds or thousands of neurons exhibited similar amplicon sizes between 10 and 20 kb (Supplementary Fig. 6), consistent with estimates by standard and alkaline gel electrophoresis8. In contrast, MALBAC showed a much shorter correlation length B600 bp (Supplementary Fig. 8), consistent with the reported average amplicon size (5001,500 bp, ref. 22). We also found signicant correlations at the fragment-size level in one single-cell library and the reference bulk library5 that persisted at high sequencing depths (Supplementary Fig. 5); these correlations reected substantial GC bias at the fragment level absent in the other bulk libraries and likely arose during library preparation due to PCR. Despite the vastly different correlation lengths evident in MDA and MALBAC amplications, our analysis accurately predicted the cumulative coverage distribution in all libraries sequenced to above 10 from computationally down-sampled

sequencing data at 1 or less (Supplementary Figs 2 and 48).

To benchmark the performance of different single-cell libraries, we compared the fraction of covered genome (Z1 ) when each

library was sequenced to 1 . This percentage was either

computed directly from down-sampled data (when the original data had higher depths) or inferred from the depth-of-coverage curve when the original data had lower depths. The coverage benchmark was plotted against the magnitude of amplicon-level variation as measured by the plateau correlation strength at the amplicon scale (Methods, Fig. 2c). As expected, smaller amplication bias results in a larger fraction of covered genome. Out of the ve published single-cell DNA sequencing studies analysed here, the single-neuron libraries had the best overall uniformity, followed by the two single YH1 libraries; the MALBAC libraries overall had less amplication bias than MDA, although optimized MDA libraries performed equally well. The frozen glioblastoma libraries (59 total) exhibited a range of variations that can be tted by an empirical relationship

0:86

1:2

p 1

where y is the percentage of covered genome and x is the

4 NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822 ARTICLE

(dimensionless) correlation magnitude. Except for the single-sperm libraries that exhibited substantial bias, all other analysed data closely followed this relationship. This result suggested that the uniformity of genome coverage is solely determined by the amplicon-level variation but not the amplicon size. Therefore, one can directly use this empirical relationship to benchmark the uniformity of single-cell libraries by the correlation magnitude that can be accurately computed from low-pass sequencing B0.1 .

We further selected the best single-cell libraries from each study and compared the fraction of genome covered at different depths as observed in the original high-depth sequencing (Fig. 2d). Owing to the different sequencing depths applied to these libraries, we plotted all cumulative genome coverage against the normalized depth (by the mean depth). The benchmark of amplication uniformity as measured by the depth-of-coverage curve agrees with the computed correlation magnitude (Fig. 2c inset).

Finally, we also analysed the base-level coverage in single-cell libraries amplied by degenerate oligonucleotide primed PCR (DOPPCR)28. The correlation was evident both at the read length level (B50 bp) and on a longer scale of B200 bp (Supplementary Fig. 9) that is consistent with the size of puried DOP-PCR product4. In comparison with MDA- or MALBAC-generated libraries, the smaller overall correlation magnitude (at the amplicon level) explains the better uniformity of DOP-PCR. Interestingly, even for the MDA generated libraries, shorter amplicon size tends to result in better uniformity (Supplementary Fig. 9); the underlying mechanism for this observation requires further characterization.

Genome coverage variation reects allele-level bias. Coverage at the locus level includes contributions from homologous chromosomes (the allele-level coverage). The same non-uniformity in the

Amplification of homologous chromosomes

Allele A

Allele B

Allele coverage predictions for single glioblastoma libraries

0 100 101 102

Allelic distortion in amplification

Uneven coverage at heterozygous sites

Allele A

Allele B

Locus-level (combined)

Mixed template Segregated template

Correlated coverage of allele A and allele B

Independent coverage of allele A and allele B

GBM #1 GBM #2 GBM #3

Obs. (Chr. 13) vs pre. (Chr. 12) Obs. (Chr. 12) vs pre. (Chr. 12)

Obs. (Chr. 5) vs pre. (Chr. 12)

Observation from 9x

Prediction from 1x

MTM STM

Locus-level

Combined

0.8

Ref allele

Alt allele

Percentage of SNPs

0.6

0.4

0.2

0 100 101 102

Depth of coverage

0 100 101 102

Depth of coverage

Figure 3 | Amplication bias of homologous chromosomes. (a) Schematic illustration of the mixed template model and the segregated template model reecting different allele-level contributions to the same locus-level coverage. (Methods, Supplementary Fig. 10). (b) Comparison of the allele coverage predictions (Pre.) from 1 sequencing depth with the observed coverage at heterozygous sites (Obs.) at 9 sequencing depth in three single

glioblastoma libraries. The combined coverage of reference and alternate bases (red dots) at 9 sequencing validates the prediction from 1 sequencing

(dashed curve). The allele coverage (reference or alternate) is then predicted from the combined coverage assuming mixed templates (MTM, blue dotted lines) or segregated templates (STM, green dotted lines) and compared with the coverage of reference (blue triangles) or alternate (green triangles) bases at heterozygous sites. The predictions were made from the sequence coverage in disomic Chr. 12 but the agreement with observations in different disomic chromosomes demonstrate that amplication bias is consistent in all chromosomes.

Depth of coverage

NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822

genome coverage, however, may result from different combinations of non-uniformity at the allelic level (Fig. 3a). Although allele coverage determines the sensitivity to detect heterozygous variants, we rarely consider this aspect in bulk sequencing due to the comparable contributions of all alleles and largely uniform coverage of the genome. In single-cell libraries, however, we often observe disproportionately represented alleles and numerous loci may exhibit allelic dropout5,12. Consequently, the detection sensitivity of hemizygous variants is measured by the allele coverage and needs to be derived from the genome coverage.

To predict the allele coverage from the locus-level genome coverage, we considered two limiting scenarios: a segregated template model (STM) assuming completely independent amplication of homologous chromosomes, and a mixed template model (MTM) assuming identical coverage of homologous chromosomes (as expected in bulk sequencing, Fig. 3a). The difference between the two models is most evident in highly amplied regions: STM implies preferential amplication of one allele, whereas MTM suggests that both alleles have been highly amplied. Both models are analytically solvable and can be easily implemented computationally (Methods, Supplementary Fig. 10).

We compared the model predictions for allele-level coverage with the observation at germline heterozygous sites detected from bulk DNA sequencing (Fig. 3b, Supplementary Figs 5 and 11). For glioblastoma libraries (Fig. 3b), both locus- and allele-level coverage was calculated from disomic chromosome 12 at 1 sequencing

depth. Coverage at heterozygous sites was evaluated for different disomic chromosomes (5, 12 and 13) from higher-depth sequencing at 910 . As expected, the total coverage (reference plus

alternate bases) at these sites agreed well with the prediction for locus-level coverage, reecting similar amplication bias for different chromosomes with the same copy number. Meanwhile, coverage of either reference or alternate bases followed the same distribution as predicted by the STM model. These results suggested homologous chromosomes are amplied almost independently during WGA and manifest the same degree of amplication bias. This discovery was further underscored by the agreement between the observed coverage of monosomic chromosome 10 and the STM allele-coverage prediction (Supplementary Fig. 11).

We further veried that coverage of alternate or reference alleles was indeed independent of each other in the glioblastoma samples by looking at the distribution of alternate and reference reads at heterozygous sites in disomic chromosome 5 (Supplementary Fig. 12). Interestingly, the two-cell RPE-1 libraries showed positive correlations between the counts of the reference and of the alternate alleles (Supplementary Fig. 12), consistent with the MTM model (Supplementary Fig. 11). Of the two published single YH1 libraries5, one agreed better with the MTM model and the other agreed with the STM model (Supplementary Fig. 5). Whether this difference resulted from the cells initial condition (frozen versus fresh), the stage of cell cycle, or other factors requires further characterization.

Census-based strategy enables efcient variant detection. Our analytical prediction of the allele coverage measures the average probability of capturing a single-variant read in single-cell sequencing. In sequencing analysis, however, more than one observation of the variant is necessary to mitigate sequencing errors. This requirement substantially reduces the percentage of detectable variants at low sequencing depths. In one example (GBM #4, correlation magnitude E4 for disomic chromosomes), the normalized allele coverage implied that only 13.3% of clonal hemizygous variants could be condently detected at a mean sequencing depth of 1 when requiring at least two reads for

each variant (Supplementary Fig. 11). This percentage increased

with sequencing depth to a limit of 79% at 100 . In contrast, the

sensitivity to detect a subclonal mutation with allelic fraction of0.4 in a bulk library at 10 sequencing is B80% and quickly

reaches 495% at a sequencing depth of 20 (ref. 29). The

reduced dependence of detection sensitivity on sequencing depth for single-cell libraries suggested that deep sequencing of an individual library is not an efcient approach to increase power for detecting variants from libraries prepared by WGA.

To overcome this challenge, we devised an approach to sequence a large number of single-cell genomes at only modest depths (B1 ). We simultaneously controlled for errors resulting from

random MDA artifacts or from sequencing by requiring true variants to appear in multiple libraries (census based, Fig. 4a). We expected this population-based approach to be effective only when the amplication bias is random, but not recurrent (Fig. 1c). We thus evaluated the correlation between the coverage of reference and alternate alleles in four independent glioblastoma libraries. The small covariance (B0.01) between the coverage of each given allele in different libraries is consistent with random MDA bias (Table 1). These data contrasted with recurrent locus-specic amplication bias in degenerate-oligonucleotide-primed PCR methods such as GenomePlex30.

We next examined how many single cells sequenced to the same total depth would maximize the total allele coverage by census-based variant detection using a representative library with modest bias (GBM #4, correlation magnitude E4, Fig. 4b). In all the cases, our model predicted maximum allele coverage when each individual cell was sequenced to a modest depth (B1 ).

We repeated this calculation using each of the other libraries as the representative, and found that the optimal depth for detecting clonal and subclonal variants is always t 1 (Fig. 4c).

To test this experimentally, we sequenced each of the following subsets of single glioblastoma libraries to 20 total depth: 59

libraries (B0.33 per library), 22 libraries (B1 per library),

two libraries (B10 each, group A) with minimal bias

(correlation magnitude E0.9 for disomic chromosomes) and two libraries (B10 each, group B) with average bias

(correlation magnitude 2B4). We genotyped germline hetero

zygous single-nucleotide polymorphisms (SNPs) and detected somatic single-nucleotide variants and small insertion/deletions (indels) by the census-based strategy and compared the call sets with results from bulk DNA sequencing. For germline SNPs in disomic chromosome 5, we observed that census-based detection in the two pools of single-cell libraries (59 and 22 each) each uncovered more than 80% of all SNPs detected in bulk, while the two sets of two libraries with minimal and average bias uncovered only B30 and B5% of the heterozygous sites, respectively (Fig. 4d). Even combining all four deeply sequenced libraries together to a total depth of 40 still cannot reach the detection

sensitivity offered by the two larger groups. A similar improvement in sensitivity was observed for the detection of somatic single-nucleotide variants and indels among the single cells sequenced to B0.33 and B1 per library (as opposed to

B10 per library), detecting more somatic variants found in

bulk whole-exome sequencing with fewer private or false positive calls (Fig. 4e, Supplementary Data 15). The false positive calls usually occur at low allele frequencies within each library and likely reect recurrent amplication errors and sequencing errors. Such errors are less frequent when the library is sequenced to a low depth and can be suppressed by requiring more than one read for each variant. Together, these data validate our statistical estimates of the variant detection sensitivity from a population of single cell libraries and demonstrate that a census-based strategy using only modest depths of sequencing for many single cells can substantially improve both sensitivity and specicity for detecting variants compared with deep sequencing of individual libraries.

6 NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822 ARTICLE

Census-based variant calling

Predicted census-based sensitivity

Predicted optimal depth per library

20 total sequencing depth

Library #1

Accept (variant allele detected in 2 libraries) Reject (variant allele detected in 1 library only)

Subclonality

0.60.30.1

Subclonality

0.6

0.3

0.1

0.5

...

0.60.30.1

0 0.5

1 1.5

0.25 0.5 1 2 4 5 10 20 30

Census-based sensitivity = % allele covered in 2 libraries

0.60.30.1

Depth per library (x)

Observed census-based sensitivity (germline/clonal) Observed census-based sensitivity (somatic/subclonal)

Both alleles (correct genotype) At least alt. Allele

22 2

(A)

Also in bulkPrivate / false positive

20 total 40 20 total 40

20 total 40

# Of cells sequenced

SNPs detected (%)

0.8

SSNVs/indels

detected (#)

0.6

0.4

0.2

59 22 2

(A)

2(B)

4 59 22 2(A)

2(B)

# Of cells sequenced

Figure 4 | Variant detection in single-cell genomes. (a) Census-based variant calling requires that acceptable variants be observed in at least two independent single-cell libraries. (b) Estimates of the census-based detection sensitivity for a population of independently amplied single-cell libraries all assumed to have similar amplication bias as GBM #4 (Supplementary Fig. 11). Optimal detection sensitivity is achieved at roughly 0.5 depth-per-library

regardless of the sub-clonal fraction or the total sequencing depth. (c) Optimal depth-per-library for census-based variant detection in a population of independently amplied single-cell libraries assumed to have similar coverage bias. The range of the optimal depths is calculated on the basis of the amplication bias observed in single glioblastoma libraries in Fig. 2b. For libraries with more bias or for the detection of variants with lower clonal fractions, it is optimal to sequence more libraries at modest depths (0.10.5 ). (d) Observed coverage of reference and alternate bases at heterozygous SNP sites in

disomic Chr. 5 as an estimate of the census-based detection sensitivity for clonal variants. A varying number of single glioblastoma nuclei (59, 22 and 2) were sequenced to the same total depth (20 ) and genotyped at germline heterozygous SNP sites. Group (A) included two cells with the best uniformity

and group (B) included two cells with average uniformity. For either heterozygous coverage or the detection of alternate bases, the larger pools offer better sensitivity than the two groups of two cells. (e) Comparison between somatic non-synonymous variants detected in different-sized pools of single cells sequenced to the same total depths (20 ). The truth set (48 variants in total) included 43 variants that were detected in both 30 whole-genome and

120 whole-exome sequencing of bulk tumour DNA, plus ve additional variants detected in bulk whole-genome and single-cell sequencing. At the same

overall sequencing depth, census-based detection from a population of cells (59 and 22) offers higher sensitivity and better specicity over deep sequencing of two libraries. A larger number of private/false positive mutations are observed when individual samples are sequenced to higher depths, and these private calls often arise from sporadic sequencing errors that coincide with amplication errors.

Table 1 | Overlap and correlation between allele coverage in independent single-cell libraries by multi-strand displacement amplication.

(a) Coverage at heterozygous sites in single glioblastoma nuclei libraries.

Depth Total Reference Alternate Allelic % Hets (est.) Hets (obs.)

(i) 9.2 49,457 40,345 40,356 72 28,931 29,336

(ii) 8.1 48,745 39,569 39,521 70 27,787 28,149

(iii) 6.6 35,765 22,163 21,549 39 8,486 7,950

(iv) 9.0 37,507 23,763 23,883 42 10,084 10,144

(b) Overlap between independent single-nuclei libraries (covariance pAB pA pB).

Allele A Allele B Allele A Allele B Allele A Allele B

Cell (i) 40,345 40,356 Cell (i) 39,569 39,521 Cell (i) 40,345 40,356 Cell (ii) 39,569 39,521 Cell (ii) 22,163 21,549 Cell (ii) 23,763 23,883 Overlap 28,912 28,953 Overlap 15,290 15,195 Overlap 17,420 17,521 Covariance 0.010 0.011 Covariance 0.006 0.001 Covariance 0.007 0.007

Total germline heterozygous SNPs in Chr. 5: 56,278 (quality (qual.) Z500, HapMap).

Allele coverage in each library is evaluated by the number of covered HapMap heterozygous SNP sites in disomic chromosome 5 detected in bulk sequencing (combining blood and bulk tumour) by

UniedGenotyper (Qual. Z500). (a) In each single-cell library, coverage of A and B alleles is almost equal and the expected overlap assuming random A or B allele coveragethe estimated coverage of heterozygous sitesis comparable to the observed number of heterozygous sites. (b) The overlap between different single-cell libraries coverage of each allele is also close to the expected overlap based on random allele coverage.

NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822

DiscussionHere we have established a universal method to characterize the amplication bias in single-cell DNA libraries at both locus and allele levels. On the basis of our discovery that intrinsic amplication bias occurs predominantly at the amplicon level, we demonstrated that the cumulative distribution of bin-level coverage (with bin size set to the length scale of dominant amplication bias) directly predicts the depth-of-coverage at any sequencing depth. We further derived a quantitative measure of amplication bias that can directly predict locus-level coverage via an empirical relationship. Our analysis thus provides a statistical description of the relationship between the genomic coverage of single-cell DNA libraries and the intrinsic amplication bias. This metric provides a robust benchmark that enables a quantitative prediction of the complexity of single-cell libraries from low-pass sequencing(0.01B0.1 ).

We demonstrated that amplication of different chromosomes (including different homologous chromosomes) in a single cell is often independent (segregated template model), reecting random priming and amplication. This biophysical feature is fundamentally different from amplication from bulk DNA, where allele-level coverage is strongly correlated31,32 (mixed template model). We proposed analytically solvable models that can quantitatively predict the allele coverage of single-cell libraries at any sequencing depth. These models provide the basic framework for estimating the detection sensitivity of hemizygous genetic variants by single-cell sequencing.

The characteristic length in the coverage autocorrelation also determines the scale at which the source of amplication bias should be characterized. In bulk DNA libraries, a dominant bias at the fragment length level is shown to be associated with the sequence content (GC%), but such bias quickly decays at longer length scales (Supplementary Figs 5 and 6). In MDA-generated libraries, however, we observed substantial variation even in regions with similar GC content (Supplementary Fig. 6). This is in sharp contrast to MDAs from bulk samples18,3133. Such a wide range of variation reects random priming bias17 instead of recurrent polymerase extension bias, and may also depend on the size of DNA templates after cell lysis, which is known to affect displacement efciency21. Our discoveries of the amplicon-level correlation and independent allele amplications are both consistent with the dominant bias being generated in the early stage of amplication of single DNA templates and reect the discrete nature of single-molecule biochemical reaction. As early stage bias can be exponentially amplied during subsequent cycles of amplication, limited amplication should result in better uniformity27,34.

The random nature of single-cell genome amplication further underscores the necessity of single-cell-specic bioinformatic tools and experimental design. Deep sequencing of single-cell libraries to recover measures of variant alleles easily extends the sequencing cost and becomes prohibitive for libraries with extreme bias. Our analyses suggest a more practical approach by (1) preparing individual sequencing libraries from many independent samples and (2) ranking and selecting the best libraries on the basis of the complexity and the allelic coverage predicted based on low-pass whole-genome sequencing of each library (B0.1 ) before extensive sequencing.

For clinical samples with a limited number of cells, such as ne-needle aspirates or circulating tumour cells, the most interesting genetic variants are shared among the cells, including both subclonal and clonal variants. For this purpose, it is most efcient to perform census-based variant detection from multiplexed sequencing of independently amplied single-cell DNA libraries each sequenced to modest depths (B1 ).

The census-based variant detection strategy simultaneously controls random errors due to sequencing (0.11% per sequenced base) or amplication (B1% loci with error reads exceeding 10% allele frequency, Supplementary Fig. 7, refs 27,34) and maximizes the total allele coverage at a given sequencing depth by sampling many independently amplied libraries, thus enabling accurate detection of somatic variants and dissection of clonal heterogeneity.

One technical complication in single-cell sequencing is DNA contamination. Contamination of non-human-genomic DNA before whole-genome amplication will result in a large percentage of sequencing reads that are not mapped to the reference assembly, which can be readily identied and excluded by low-pass sequencing. The census-based strategy also effectively controls human genomic DNA contamination limited to one single-cell library. Contaminations to multiple single-cell libraries are usually present at many more copies than a single-cell genome at the affected loci and should be recognizable as they are substantially amplied after whole-genome amplication.

At the current stage, errors introduced during WGA prohibit an accurate characterization of individual genetic variants within a single cell. (This task can be accomplished through independent amplications of biological replicates after cell division.) It is, however, possible to infer global features of mutagenesis, such as the mutation rates in tumour progenitor cells or circulating tumour cells, by single-cell sequencing after correcting the total number of detected genetic variants by the statistical power for detecting variants in a single-cell library sequenced to a certain depth. Our analyses have laid the foundation for single-cell genetic variant detection by calibrating the amplication bias at both genomic and allelic levels.

Methods

Amplication and sequencing of RPE-1 cells. The hTERT RPE-1 (ATCC) cell line stably expressing GFP-H2B was cultured and treated as follows: Briey, cells were transfected with a pool of siRNAs (Smartpool, Dharmacon) against p53 using RNAiMAX (Invitrogen) according to the manufacturers instructions. Eighteen hours later, cells were treated with Nocodazole (100 ng ml 1; Sigma) for 6 h. The

G2/M arrested cells were collected by mitotic shake-off and replated after three washes with medium. Four hours after replating, G1-released cells were sorted into 384-well tissue culture plates and cultured. Conrmed single cells were allowed to divide once, before being washed twice with PBS and lysed and amplied within the 384-well tissue culture plate as outlined above35.

Amplied DNA from two RPE-1 cells after one round of cell division was subject to standard whole-genome DNA library preparation and assessed by low-pass sequencing B0.1 using the MiSeq platform (Illumina). DNA libraries of

RPE cells (three total) were then sequenced to 49 on the HiSeq2500 platform

(Illumina). Bulk RPE-1 DNA was sequenced to B12 on the HiSeq2500 platform

(Illumina).

Processing of single-cell sequencing data. Sequencing reads from published studies were downloaded from the NCBI Short Read Archive. For the diploid YH genome, we downloaded all sequencing runs of the bulk reference (SRR294761) and two single-cell samples, BGI_YH1 (SRR294759) and BGI_YH2 (SRR294760). For diploid neurons, we downloaded all the data from SRP014781, including sequencing data for the bulk DNA, and for the whole-genome amplied products from single-cell DNA, 100-cell DNA and 50,000-cell DNA. For haploid sperms, we downloaded the deep sequencing data of eight single sperm libraries, Sperm23 (SRS344176), Sperm24 (SRS344190), Sperm 27 (SRS344191), Sperm28 (SRS344192), Sperm101 (SRS344222), Sperm113 (SRS344223), Sperm135 (SRS344224), Sperm136 (SRS344225). For SW480 tumour cells, we obtained data corresponding to the bulk reference (SRS374235), a single-cell MDA library (SRS375060) and ve single-cell MALBAC libraries (SRS373654, SRS374233, SRS375671, SRS375672, SRS375673). The data of the glioblastoma libraries were generated from a previous study and can be accessible from SRP052627.

Reads were aligned to the human genome reference (hg19/GRCh37) using bwa (http://bio-bwa.sourceforge.net/) in the paired-end mode. The RPE and glioblastoma libraries were aligned by bwa aln followed by bwa sampe with default parameters. The remaining data were aligned by bwa mem. PCR duplicates were removed by MarkDuplicates from PICARD (http://picard.sourceforge.net/). Sequencing data of the glioblastoma libraries and the matching blood were re-calibrated and indel-realigned by GATK (http://www.broadinstitute.org/gatk/) before variant detection.

8 NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822 ARTICLE

Down-sampling of deep sequencing data to B1 was done by

DownsampleSam from PICARD. Base-level sequencing coverage was enumerated by the DepthOfCoverage module from GATK with minimum read mapping quality set to 5.

To evaluate the allele coverage in RPE-1 MDA libraries, we detected heterozygous SNPs in Chr. 1 of the RPE-1 cells from the sequencing of bulk RPE-1 DNA (B12 ) and individual MDA libraries by UniedGenotyper from

GATK;

only variants with Qual. Z100 and at least three reference and three alternate reads in the bulk sample were selected to evaluate the allele coverage in MDA libraries.

For other samples, we genotyped HapMap SNPs (v3.3) to estimate the allelic coverage; only variants found to be heterozygous in the matching blood with Qual. Z500 were selected and genotyped in each set of glioblastoma libraries. Somatic single-nucleotide variants and small insertions/deletions were detected by HaplotypeCaller from GATK in each set of glioblastoma libraries and in the bulk library, and by MuTect29 from bulk whole-exome sequencing.

Computation of auto-correlation function of sequence coverage. The dimensionless auto-correlation function of coverage is dened as

G D

C x

C x D

h i C x

f2 l2 2 1 l

2 :

h i2

C x

h i2

: 2

The brackets denote average over all genomic loci x and D measures the spread of correlation. In computing the auto-correlation functions, we only include regions not adjacent to the assembly gaps. (Adjacency is determined by the step D.)

The correlation function is tted to an exponential form to estimate the correlation length lc:

G D

a be D=l : 3 For MDA, the correlation length lc is on the order of 10 kb and the correlation function G(D) is roughly constant above the fragment length (B300 bp) and below the correlation length lc. In this regime, G(D) can be written as

G D

h i

The iteration can be continued to calculate the allele coverage at any depth,

P C M

XM 1k0P m1 k P m2 M k P m1 M 16

or (denoting l0 1, l1 l and so on.)

XM 1 k0

lk1

XM 2 k1

lk1

2 1 l lM lM 1l;

which gives

h i2

: 4

Here

C is the average coverage within each bin [x, x D). It becomes evident that

G(D) measures the standard deviation of bin-level coverage. For convenience, we choose to evaluate G(D) at D 1 kb as a quantitative metric of the magnitude of

amplication bias (correlation magnitude).

Statistical models for predicting allele coverage from genome coverage. The power to detect a genetic variant is given by the probability that this variant locus (usually of one chromosome) is represented in the sequencing data, or the relative abundance of variant-supporting reads. But the direct observable in sequencing data is the total number of reads covering all possible alleles, that is,

C m1 m2 mn; 5 where C is the total observed coverage at a given locus as a sum of contributions from each allele denoted by mi.

In the presence of amplication bias, both C and mis vary across the genome. The distribution of C across different loci can be straightforwardly evaluated from the depth-of-coverage curve; here, we want to infer the statistical distribution of mi when the distribution of C is known. The STM assumes that amplications of homologous chromosomes are independent. As a consequence, the counts of reference and of alternate bases at heterozygous sites are independent, and one highly amplied allele may dominate over the remaining ones. In the MTM, different alleles are assumed to be amplied to the same extent at every individual locus. As a result, the counts of reference and of alternate bases at heterozygous sites follow a symmetric binomial distribution.

In mathematical terms, mis are independent of each other but follow the same distribution in STM. In this scenario, one can numerically compute the distribution of mi from the characteristic functions C(k) and m(k) (that is, the Fourier transforms of the probability distribution for C and m), which satisfy

C k

m k

n: 6 Here, we present an iterative method to calculate the distribution of mi and illustrate this method using a diploid genome (that is, n 2).

At a given sequencing depth, denote the total percentage of loci that are covered Z1 by f,

P C 1

f ; 7 the percentage of loci that are covered in a particular allele is denoted by

P mi 1

l: 8 It is then straightforward to see that

P C 1

1 Yi1 P mi 1 9

f 1 1 l

n: 10

Hence, in a region with n alleles, the probability that a given allele is covered is given by

l 1 1 f

1=n: 11 For diploid genomes, this becomes

l 1 1 f

1=2: 12 We can expand this further to compute the coverage at higher depths. For example,

P C 2

P m1 0

P m2 2

P m1 1

P m2 1

P m1 2

If we denote the percentage of loci where total coverage is at or above two as f2,

and the percentage of loci covered at or above two for each allele as l2, then we have

f2 1 l

l2 l l2

l l2; 14

" #

2 1 l

fM llM

XM 2 k1

lk1

: 18

In the mixed template model, we assume that the local coverage C is a mixture of all alleles randomly sampled at the same frequency. In disomic regions, this implies that m follows a binomial distribution B(C, 0.5) at any total coverage C. Under this model, we have

l P m 1

XMt1P C t 1 0:5t

where the sum runs over all observed local coverage (t 1, 2, y M). The series

converges quickly as both ft and the exponential prefactor decay quickly. Furthermore, one easily veries that when f is small, this result is equal to the segregated template model to the leading order (f/2).

It is also straightforward to calculate the allele coverage at higher depths.

2 P C 1

22 P C 2

2 f

14 f2

2t ft ;

P m k

XMtkP C t 1 2 t

Xk 1s0t !s ! t s !

: 20

Census-based detection sensitivity from a pool of single-cell libraries. As the percentage of genome that is covered at or above 1 at any sequencing depth can

be estimated, we can also predict the census-based detection power for hemizygous variants in a pool of single-cell libraries. Consider a total number of Y libraries having similar amplication bias and the probability of observing a hemizygous variant in any of the Y libraries is given by l, then the probability for observing this variant in a subset of libraries (X out of Y) is given by

P Covered in X libraries

XYmXY !m ! Y m !lm 1 l Y m: 21

We can then compute this for a subclonal variant at clonal fraction y in a total of Z libraries from

P Covered in X libraries

XZYXZ !

Z Y ! Y !yY 1 y Z Y

|{z}

variant is harbored by Z X cells

XYmXY !m! Y m !lm 1 l Y m

|{z}

where random selection of cells containing the subclonal variant follows a binomial distribution B(Z,y).

variant allele is covered inZX libraries

NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7822

References

1. Kalisky, T., Blainey, P. & Quake, S. R. Genomic analysis at the single-cell level. Annu. Rev. Genet. 45, 431445 (2011).

2. Shapiro, E., Biezuner, T. & Linnarsson, S. Single-cell sequencing-based technologies will revolutionize whole-organism science. Nat. Rev. Genet. 14, 618630 (2013).

3. Chi, K. R. Singled out for sequencing. Nat. Methods 11, 1317 (2014).4. Navin, N. et al. Tumour evolution inferred by single-cell sequencing. Nature 472, 9094 (2011).

5. Hou, Y. et al. Single-cell exome sequencing and monoclonal evolution of a JAK2-negative myeloproliferative neoplasm. Cell 148, 873885 (2012).

6. Wang, J., Fan, H. C., Behr, B. & Quake, S. R. Genome-wide single-cell analysis of recombination activity and de novo mutation rates in human sperm. Cell 150, 402412 (2012).

7. Lu, S. et al. Probing meiotic recombination and aneuploidy of single sperm cells by whole-genome sequencing. Science 338, 16271630 (2012).

8. Evrony, G. D. et al. Single-neuron sequencing analysis of L1 retrotransposition and somatic mutation in the human brain. Cell 151, 483496 (2012).

9. McConnell, M. J. et al. Mosaic copy number variation in human neurons. Science 342, 632637 (2013).

10. Shalek, A. K. et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature 498, 236240 (2013).

11. Xue, Z. et al. Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing. Nature 500, 593597 (2013).

12. Lohr, J. G. et al. Whole exome sequencing of circulating tumor cells provides a window into metastatic prostate cancer. Nat. Biotechnol. 32, 479484 (2014).

13. Ni, X. et al. Reproducible copy number variation patterns among single circulating tumor cells of lung cancer patients. Proc. Natl Acad. Sci. USA 110, 2108321088 (2013).

14. Eberwine, J., Sul, J.-Y., Bartfai, T. & Kim, J. The promise of single-cell sequencing. Nat. Methods 11, 2527 (2013).

15. Blainey, P. C. The future is now: single-cell genomics of bacteria and archaea. FEMS Microbiol. Rev. 37, 407427 (2013).

16. Zhang, L. et al. Whole genome amplication from a single cell: implications for genetic analysis. Proc. Natl Acad. Sci. USA 89, 58475851 (1992).

17. Zhang, K. et al. Sequencing genomes from single cells by polymerase cloning. Nat. Biotechnol. 24, 680685 (2006).

18. Pinard, R. et al. Assessment of whole genome amplication-induced bias through high-throughput, massively parallel whole-genome sequencing. BMC Genomics 7, 216 (2006).

19. Geigl, J. B. et al. Identication of small gains and losses in single cells after whole genome amplication on tiling oligo arrays. Nucleic Acids Res. 37, e105 (2009).

20. Dean, F. B. et al. Comprehensive human genome amplication using multiple displacement amplication. Proc. Natl Acad. Sci. USA 99, 52615266 (2002).

21. Lage, J. M. et al. Whole genome analysis of genetic alterations in small DNA samples using hyperbranched strand displacement amplication and array-CGH. Genome Res. 13, 294307 (2003).

22. Zong, C., Lu, S., Chapman, A. R. & Xie, X. S. Genome-wide detection of single-nucleotide and copy-number variations of a single human cell. Science 338, 16221626 (2012).

23. Gole, J. et al. Massively parallel polymerase cloning and genome sequencing of single cells using nanoliter microwells. Nat. Biotechnol. 31, 11261132 (2013).

24. Lander, E. S. & Waterman, M. S. Genomic mapping by ngerprinting random clones: a mathematical analysis. Genomics 2, 231239 (1988).

25. DePristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. 43, 491498 (2011).

26. Daley, T. & Smith, A. D. Predicting the molecular complexity of sequencing libraries. Nat. Methods 10, 325327 (2013).

27. Francis, J. M. et al. EGFR variant heterogeneity in glioblastoma resolved through single-nucleus sequencing. Cancer Discov. 4, 956971 (2014).

28. Wang et al. Clonal evolution in breast cancer revealed by single nucleus genome sequencing. Nature 512, 155160 (2014).

29. Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213219 (2013).

30. Voet, T. et al. Single-cell paired-end genome sequencing reveals structural variation per cell cycle. Nucleic Acids Res. 41, 61196138 (2013).

31. Hosono, S. et al. Unbiased whole-genome amplication directly from clinical samples. Genome Res. 13, 954964 (2003).

32. Paez, J. G. et al. Genome coverage and sequence delity of phi29 polymerase-based multiple strand displacement whole-genome amplication. Nucleic Acids Res. 32, e71 (2004).

33. Pugh, T. J. et al. Impact of whole genome amplication on analysis of copy number variants. Nucleic Acids Res. 36, e80 (2008).

34. De Bourcy et al. A quantitative comparison of single-cell whole genome amplication methods. PLoS ONE 9, e105585 (2014).

35. Ganem, N. J., Godinho, S. A. & Pellman, D. A mechanism linking extra centrosomes to chromosomal instability. Nature 460, 278282 (2009).

Acknowledgements

We thank Dr David Pellman for sharing the sequencing data on RPE cells, M. Leibowitz,X. Cai and G. Evrony for discussions, and the Koch Institute Swanson Biotechnology Center (specically the BioMicro Center) for technical support. C.-Z.Z. was supported by the National Cancer Institute (U24CA143867 to M.M.). V.A.A. was supported, in part, by a graduate fellowship from the National Science Foundation. J.C.L. is a Camille Dreyfus Teacher-Scholar. This work was supported by the Bridge Project, a collaboration between Koch Institute for Integrative Cancer Research at MIT and the Dana-Farber/ Harvard Cancer Center (DF/HCC) to (J.C.L., K.L.L. and M.M.) and the National Brain Tumor Society. This work was also supported in part by Janssen Pharmaceuticals, Inc. and the Koch Institute Support (core) Grant P30-CA14051 from the National Cancer Institute.

Author contributions

C.-Z.Z. and V.A.A. initiated the project and carried out the analysis. C.-Z.Z. performed analysis of amplication bias; V.A.A. performed analysis of census-based detection sensitivity with help from C.-Z.Z. J.F., H.C., C.M. and K.L.L. prepared sequencing libraries for the RPE cell line and glioblastoma samples. C.-Z.Z., V.A.A., J.C.L. and M.M. wrote the manuscript with help from all the authors. M.M. and J.C.L. supervised the study.

Additional information

Accession codes: The sequence data have been deposited in the Short Read Archive from NCBI under the following accession codes: RPE-1 bulk (SRX858057); two-cell RPE libraries (SRX858832, SRR1779331 for RPE #1, SRR1779329 for RPE #2, SRR1779330 for RPE #3); single RPE libraries (SRX858836, SRX858838, SRX858840, SRX858841); glioblastoma bulk whole-genome sequencing (SRX848889); glioblastoma bulk whole-exome sequencing (SRX857666); single-glioblastoma nuclei pool #1 (59 nuclei, SRX858332); single-glioblastoma nuclei pool #2 (22 nuclei, SRR1778915, SRR1779027, SRR1779078, SRR1779079, SRR1779080, SRR1779083, SRR1779085, SRR1779088, SRR1779089, SRR1779091, SRR1779092, SRR1779093, SRR1779095, SRR1779098, SRR1779157, SRR1779161, SRR1779163, SRR1779167, SRR1779172, SRR1779174, SRR1779175, SRR1779177); deeply sequenced single-glioblastoma nuclei (SRX858848, SRR1779345 for GBM #1, SRR1779347 for GBM #2; SRR1779348 for GBM #3; SRR1779350 for GBM #4); whole-genome sequencing of blood reference for the glioblastoma patient (SRX851083); whole-exome sequencing of the blood reference for the glioblastoma patient (SRX857684).

Supplementary Information accompanies this paper at http://www.nature.com/ naturecommunications

Competing nancial interests: M.M. is a founder and equity holder of Foundation Medicine, a for-prot company that provides next-generation sequencing diagnostic services.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/

How to cite this article: Zhang, C.-Z. et al. Calibrating genomic and allelic coverage bias in single-cell sequencing. Nat. Commun. 6:6822 doi: 10.1038/ncomms7822 (2015).

10 NATURE COMMUNICATIONS | 6:6822 | DOI: 10.1038/ncomms7822 | www.nature.com/naturecommunications

Word count: 9244

Show less

Abstract

Translate

Artifacts introduced in whole-genome amplification (WGA) make it difficult to derive accurate genomic information from single-cell genomes and require different analytical strategies from bulk genome analysis. Here, we describe statistical methods to quantitatively assess the amplification bias resulting from whole-genome amplification of single-cell genomic DNA. Analysis of single-cell DNA libraries generated by different technologies revealed universal features of the genome coverage bias predominantly generated at the amplicon level (1-10 kb). The magnitude of coverage bias can be accurately calibrated from low-pass sequencing (∼0.1 × ) to predict the depth-of-coverage yield of single-cell DNA libraries sequenced at arbitrary depths. We further provide a benchmark comparison of single-cell libraries generated by multi-strand displacement amplification (MDA) and multiple annealing and looping-based amplification cycles (MALBAC). Finally, we develop statistical models to calibrate allelic bias in single-cell whole-genome amplification and demonstrate a census-based strategy for efficient and accurate variant detection from low-input biopsy samples.

Details

Title

Calibrating genomic and allelic coverage bias in single-cell sequencing

Author

Zhang, Cheng-zhong; Adalsteinsson, Viktor A; Francis, Joshua; Cornils, Hauke; Jung, Joonil; Maire, Cecile; Ligon, Keith L; Meyerson, Matthew; Love, J Christopher

Pages

6822

Publication year

2015

Publication date

Apr 2015

Publisher

Nature Publishing Group

e-ISSN

20411723

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.1038/ncomms7822

ProQuest document ID

1673379025

Calibrating genomic and allelic coverage bias in single-cell sequencing

Jump to:

Full Text

Abstract

Details

Suggested sources