- Email: [email protected]
Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection
Journal Pre-proof The Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection Samantha N. McNulty, Patrick R. Mann, Joshua A. Robinson, Eric J. Duncavage, John D. Pfeifer PII:
S1525-1578(20)30042-8
DOI:
https://doi.org/10.1016/j.jmoldx.2020.02.003
Reference:
JMDI 895
To appear in:
The Journal of Molecular Diagnostics
Received Date: 2 October 2019 Revised Date:
7 January 2020
Accepted Date: 20 February 2020
Please cite this article as: McNulty SN, Mann PR, Robinson JA, Duncavage EJ, Pfeifer JD, The Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection, The Journal of Molecular Diagnostics (2020), doi: https://doi.org/10.1016/j.jmoldx.2020.02.003. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Copyright © 2020 Published by Elsevier Inc. on behalf of the American Society for Investigative Pathology and the Association for Molecular Pathology.
1
The Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection Samantha N. McNulty,* Patrick R. Mann,* Joshua A. Robinson,† Eric J. Duncavage,* and John D. Pfeifer* From the Department of Pathology and Immunology,* and the Department of Medicine, Division of Oncology,† Washington University School of Medicine, Saint Louis, Missouri Corresponding author: Samantha N. McNulty, Ph.D., Department of Pathology and Immunology, Washington University School of Medicine, Campus Box 8118, 600 S. Euclid Ave, St. Louis, MO 63110. E-mail: [email protected] Running head: Impact of DNA input in clinical NGS. Funding: Supported by the Department of Pathology and Immunology, Washington University School of Medicine. Disclosures: None declared.
2 ABSTRACT PCR amplification, a key step in next-generation sequencing (NGS) library construction, can generate an unlimited amount of product from limited input; however, it cannot create more information than was present in the original template. Thus, NGS libraries can be made from very little DNA, but reducing input may compromise assay sensitivity in ways that are difficult to ascertain unless library complexity (ie, the number of unique DNA molecules represented in the library) and depth of coverage with unique sequence reads (those derived from input DNA molecules) versus duplicate sequence reads (resulting from over-amplification of particular molecules) are discretely measured. We performed a series of experiments to explore the impact of low DNA input on an amplicon-based NGS assay using unique molecular identifiers to track unique versus duplicate reads. At high sequencing depths, unique and total (unique plus duplicate) coverage are not well correlated, so increasing the number of sequenced reads does not necessarily improve sensitivity. Unique coverage depth tends to improve with more input, but improvements are not consistent. Fluctuations in library complexity complicated variant detection using both standardized and clinical specimens, often resulting in technical replicates with vastly different estimates of variant allelic fraction. We conclude that depth of coverage with unique reads must be tracked in clinical NGS to ensure that sensitivity and accuracy are maintained.
3 INTRODUCTION Ideally, a next-generation sequencing (NGS) library would perfectly represent the DNA present in the biological specimen from which it was derived, and the sequencing process would sample that library in a manner that is sufficient to detect any variants present in the original specimen. In this idealized scenario, it is possible to mathematically model the number of reads that would need to be sampled from a given genomic locus (ie, the depth of mapped read coverage) to detect a variant of a specified allelic fraction (VAF) with a specified level of confidence. High depth of coverage is particularly important for the detection of low VAF variants like those routinely encountered in the context of testing for cancer and other somatic diseases (Figure 1)1, 2. In actual clinical practice, however, various sources of error/bias in sample preparation, NGS library construction, amplification, and/or sequencing can enrich or deplete certain DNA sequences in ways that are difficult to predict a priori 3-7. Thus, it is essential to distinguish those sequence reads that trace back to DNA molecules that were present in the original specimen (unique sequence reads) from those reads that arose from over-amplification of particular DNA molecules during the process of library construction (duplicate reads; Figure 2) 4, 8-10
. The inadvertent inclusion of duplicate reads in the bioinformatic process of
variant prediction can skew the sampling of input DNA molecules and cause i) failure to detect sequence variants that were present in the original specimen, ii) over- or under-representation of particular variants, or iii) false positive variant
4 calls resulting from PCR errors that are propagated through library preparation and sequencing. Duplicate sequence reads are easily recognized in the context of whole genome shotgun or hybridization-capture–based NGS assays. Using these methods, input DNA is randomly fragmented during the process of library preparation, so unique sequence reads are unlikely to be identical to one another. Duplicate reads are easily identified based on shared mapping coordinates (for example, using PicardTools, Broad Institute, Cambridge, MA) 11, and disregarded in downstream bioinformatic analyses. Conversely, in the setting of amplicon-based NGS, all reads derived from a given amplicon (those amplified from the same PCR primer set) are expected to be identical, so finding duplicate reads is not straightforward. The best way to assess uniqueness in an amplification-based NGS workflow is to ligate unique molecular identifiers (UMIs) to input DNA molecules prior to library construction (Figure 2); sequence reads with identical UMI can be counted as duplicates in the bioinformatic pipeline since they were derived from the same input molecule 12. There are advantages and disadvantages to both hybridization-capture and amplification-based enrichment strategies, but amplicon-based tests are often favored in situations where DNA input is limited 13, which is precisely when library complexity (ie, the number of unique DNA molecules represented in an NGS library) can become compromised and the number of unique sequence reads mapped to a particular locus can drop to the point that the sensitivity of variant detection suffers. Unfortunately, UMI-based methods add a layer of cost
5 and complexity to both library construction and bioinformatic analyses, so they have not been uniformly adopted, particularly in the clinical laboratory setting. Many product manufacturers (eg, TruSeq Custom Amplicon Low Input Kit, Illumina, Carlsbad, CA; CleanPlex Custom NGS Panels, Paragon Genomics, Hayward, CA; DNA Library Prep Kit for Illumina Sequencing, Applied Biological Materials, Richmond, BC, Canada) and laboratories 10, 14, 15 tout amplicon-based assays with low DNA input, even in the absence of methods that facilitate the tracking of unique versus duplicate reads. Although the important distinction between unique versus total (unique plus duplicate) depth of coverage has been modeled in the past 1, 2,this study experimentally demonstrated the impact of reduced input DNA input on unique read coverage and variant prediction in the setting of an amplicon-based NGS assay. Libraries were generated from as little as 2ng of genomic DNA using a UMI-based preparation protocol, sequenced to high total depth of coverage, and analyzed using a UMI-aware bioinformatic pipeline to track unique versus duplicate sequence reads. The results show that library complexity, and consequently depth of coverage with unique sequence reads, suffers when DNA input is reduced, but not in consistent or predictable ways. Variations in library complexity impaired accurate variant prediction in the context of both standardized and clinical specimens. Our results stress the importance of measuring unique versus duplicate reads in all cases, particularly when DNA input is minimal.
6 MATERIALS AND METHODS
Analytical specimens DNA was isolated from two engineered cell lines, one with a knock-in BRAF p.V600E variant (HD 200-002, Horizon Diagnostics, Cambridge, UK) and another with knock-in KRAS p.G13D and PIK3CA p.H1047R variants (HD 200003, Horizon Diagnostics) as previously described 16. DNA from the two cell lines was mixed to generate variants with VAFs ranging from approximately 1% to 25% (Supplemental Tables S1 and S2).
Clinical Specimens Three de-identified, formalin-fixed, paraffin-embedded tumor specimens were selected for sequencing based on previous findings of clinically significant variants present at VAF <20% when tested in routine practice using a validated, hybridization-capture based NGS assay run in a CAP/CLIA certified laboratory 17. Specimen C1 had an IDH1 p.R132S variant at VAF=9%, C2 had a KRAS p. G12D variant at VAF=12%, and C3 had a PIK3CA p.E542K variant at VAF=3%.
Library preparation and sequencing DNA was isolated from cell lines and clinical specimens using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) and quantified using a Qubit fluorometer with a dsDNA HS Assay Kit (ThermoFisher Scientific, Indianapolis, IN). Sequencing libraries were prepared from 2ng, 5ng, 10ng, or 50ng of pre-mixed
7 cell line DNA (Supplemental Tables S1 and S2) or tumor DNA (Supplemental Table S3) using the HaloPlex HS Target Enrichment System (Agilent Technologies, Santa Clara, CA) as previously described 18-20. Designated quantities of input DNA were digested using a mixture of restriction endonucleases included in the Agilent HaloPlex HS kit, denatured, and hybridized to a biotinylated and UMI-barcoded probe library designed to target the coordinates of the AmpliSeq Cancer Panel v2 for Illumina assay (3.6kb capture space, ThermoFisher Scientific). Biotinylated targets were isolated by hybridization with streptavidin beads and amplified for 23 PCR cycles using primers tailed with Illumina sequencing motifs. Prepared libraries were pooled and sequenced on an Illumina NextSeq or Illumina NovaSeq6000 (Illumina). The number of reads generated from each library is reported in Supplemental Tables S1, S2, and S3.
Bioinformatic analysis Demultiplexed paired-end sequence reads were aligned to the human reference assembly (GRCh37-lite-build37) using BWA-MEM (version 7.9.a; parameters: -M). As previously described 20, variants were predicted from aligned BAM files using Barcrawler, a custom, Java-based GATK locus walker (available from: https://github.com/abelhj/gatk/tree/master/public; parameters: -mmq 20 mbq 20 -minCtBC 1 -dcov 500000 -discardN 1 -minOffset 5 -maxNM 3; accessed 2018-02-07). No further filtering or UMI-based error correction was imposed. Barcrawler outputs the depth of coverage with total and unique sequence reads
8 at each base position in the region of interest (ROI); these values were used to calculate the average coverage of the assay ROI. Aligned bam files were down sampled to desired coverage thresholds using Picard DownsampleSam (http://broadinstitute.github.io/picard, version 2.21.4) and re-analyzed as necessary.
Regulatory Approval This study was approved by the Washington University School of Medicine Human Studies Committee.
RESULTS
Unique Versus Total Read Coverage NGS libraries were constructed from decreasing quantities of DNA using a UMI-based PCR enrichment protocol. 50ng (standard input), 10ng, 5ng, and 2ng input libraries were sequenced to an average total read coverage of 23,685±11,210 (Supplemental Table S1). A UMI-aware bioinformatic pipeline was used to discretely track unique versus duplicate reads and predict variants. During the course of sequencing, the first read encountered with a given UMI sequence is considered unique and all subsequent reads containing that UMI are flagged as duplicates. Thus, at low depth of coverage, a high percentage of reads will be classified as unique. The percent of total reads that are classified as unique is expected to decrease as total depth of coverage increases since there
9 are a finite number of input DNA molecules to survey, but libraries can be sequenced to infinite total depth. At the high depth of coverage achieved in this experiment, there did not appear to be any correlation between depth of unique versus total read coverage (Figure 3A). NGS libraries produced from 50ng input DNA were computationally down sampled to explore the impact of total depth of coverage on the percent of captured reads that are classified as unique. The percent of unique reads detected was highest for libraries down sampled to 1,750x total depth of coverage and decreased as total depth of coverage increased (Figure 3B).
Relationship Between Unique Read Coverage and DNA Input Next, NGS libraries prepared from 50ng, 10ng, 5ng, and 2ng DNA input were compared. The average total depth of coverage of the ROI varied between libraries (range of approximately 10,126 to 44,803x), so alignments were down sampled to approximately 10,000x total depth of coverage (mean: 10,007±11.4x) to facilitate inter-sample comparison. As expected, the depth of coverage with unique reads generally increased with increasing DNA input (Figure 4).
Impact of Low Library Complexity on Variant Detection Using Standard Specimens To explore the impact of the depth of coverage with unique reads on variant prediction, DNA was isolated from two isogenic cell lines containing key heterozygous variants and mixed to produce variants with allelic fractions ranging
10 from 25% to approximately 1% (Supplemental Table S2). NGS libraries were prepared and sequenced in triplicate using either the standard input of 50ng or a reduced input of 2ng, and variants were predicted without UMI-based error correction (Figure 5). Variants from a 1:1 mixture of the two cell lines (approximately 25% VAF) were detected even when the DNA input was reduced. However, the variability between the VAFs calculated from the triplicate libraries was higher in the 2ng input libraries compared to the 50ng input libraries for variants at a range of expected VAFs. For example, the KRAS variant expected to be present at a VAF of approximately 6% in the 1:8 DNA mixture was detected at VAFs ranging from 0.04% to 10.83% in the 2ng libraries. Many of the low VAF variants were consistently detected far below their anticipated fraction in the 2ng input libraries, whereas the VAFs calculated from 50ng libraries were closer to the expected value.
Impact of Low Library Complexity on Variant Detection Using Clinical Specimens The same experiment was repeated, this time using FFPE tumor specimens that had been shown to harbor low VAF, clinically significant variants during routine, clinical testing using a validated hybridization-capture–based NGS assay (C1: IDH1 p.R132S, C2: KRAS p. G12D, and C3: PIK3CA p.E542K). NGS libraries were prepared and sequenced in triplicate using either 50ng or 2ng of input DNA (Supplemental Table S3). These libraries were sequenced to very high total depths of coverage (average: 397,739±234,664x). The average depth
11 of coverage with unique reads was higher for the 50ng input libraries compared to the 2ng input libraries (11,769±841x versus 4,500±1,892x, respectively). Although the variants were reliably detected at a consistent VAF using 50ng of input DNA, the VAFs of the variants identified from 2ng input libraries showed a greater degree of variation (Figure 6). For example, the KRAS p.G12D variant present in specimen C2 was estimated between 21.60% to 25.50% VAF from the 50ng input libraries, but VAF estimates ranged from 0.48% to 47.10% in the 2ng input libraries. This trend was consistent between the three FFPE specimens (Figure 6).
DISCUSSION Library complexity is a direct reflection of the number of independent genomes present in the input material used in library construction. It is certainly possible to prepare a NGS library from extremely low amounts of DNA or RNA template, but the information contained in that library may be limited in ways that are not obvious unless the library preparation, sequencing, and bioinformatic workflow includes steps to facilitate the discrete measurement of unique read coverage. By sequencing a series of NGS libraries prepared from decreasing amounts of input DNA extracted from controlled cell lines and previously characterized clinical specimens, it was shown that the total depth of coverage (unique plus duplicate reads) does not provide an adequate measure of library complexity or depth coverage with unique reads at clinically significant genomic loci.
12 A prepared NGS library can be sequenced to any desired depth of coverage. High total depth of coverage is needed to ensure that a sufficient number of unique molecules have been sampled; however, once the unique molecules have been sampled, further sequencing yields no additional information. This explains the lack of correlation between the depth of coverage with total versus unique reads at the high total depths of coverage achieved in this experiment. The same lack of correlation that can be modelled theoretically 1, 2
holds true for the clinical NGS assay in routine clinical practice 21. Although simply increasing the DNA input may improve library complexity
to some extent, this is not possible in many clinical applications. Amplificationbased enrichment strategies are specifically preferred in instances where input material is limited 22-24, including for small biopsies, evaluation of minimum residual disease, cell-free DNA (cfDNA) testing, etc. For example, the total quantity of cfDNA isolated from a peripheral blood sample is sometimes quite limited (often less than 50ng in 1mL of blood, depending on the nature and extent of the disease 25-27). Numerous vendors market kits for use with nanogram quantities of cfDNA, and many protocols include additional amplification steps to generate a sufficient mass of DNA for library preparation and NGS sequencing. As shown, depth of coverage with unique reads becomes unpredictable when limited amounts of input DNA are used for library preparation, so it must be carefully tracked to ascertain the sensitivity of each individual assay. This is particularly critical in cases where treatment decisions are based on results
13 generated from the analysis of a few nanograms of DNA (for instance, detection of specific mutations or estimate of cancer burden 28, 29). Increasing the DNA input may not adequately boost unique library complexity for all specimens. There is a great deal of variability in the quality of DNA isolated from clinical specimens depending on the tissue source, type, and degree of fixation, age of the specimen, storage conditions, etc 30. Two libraries made from the same amount of DNA may, therefore, have markedly different complexity even if the input nucleic acid mass is the same. Furthermore, biases in the amplification and sequencing of certain genomic loci can lead to regions of poor depth of coverage with unique reads even when the input material is plentiful and library complexity is adequate overall 3, 5-7. Ultimately, the concern for depth of coverage with unique reads relates to the concern for assay sensitivity. Theoretically, the minimum sensitivity of an assay is one divided by the number of unique molecules sampled during sequencing. If this sample size is low, it may not be mathematically possible to achieve the sensitivity that is expected from a given diagnostic test even if every step of library preparation is ideal, every read maps perfectly to the reference, and every variant is correctly called. In reality, library construction does not produce idealized results, and this affects variant detection in ways that are not readily predictable sample-to-sample and assay-to-assay. Even when the expected variants are identified, the estimated VAF is subject to sampling error and that error becomes more significant at lower depth of coverage with unique reads, as we have shown using both model and actual clinical specimens. Thus,
14 library complexity is a metric that should be tracked in every assay for every specimen to ensure that a given test meets or exceeds any validated requirements. Here, we focused on the measurement of library complexity in the context of amplicon-based NGS assays. In this study, we stress the need to track unique versus duplicate reads in the context of amplicon-based NGS not because library complexity is more important in amplicon-based tests compared to shotgun or hybridization-based NGS assays, but rather because it requires specific modifications to the wet-bench protocol that are often not employed (due to complexity, cost, and time constraints). Of course, library complexity is an important metric to measure regardless of the specific library preparation or enrichment method employed. There is no doubt that the measures of the impact of library complexity on the NGS results that were describe (eg, optimal depth of coverage in Figures 3-4, variability in measured VAF by input DNA mass in Figures 5-6, etc.) here reflect the technical details of the NGS assay employed (eg, DNA extraction protocol, library preparation method, bioinformatics pipeline, etc.). It can be expected that different NGS assay designs may cause differences in the magnitude of these impacts. However, since the underlying issues described regarding the effect of library complexity on NGS test results cannot be avoided, laboratories should systematically explore the influence they have on their specific amplicon-based NGS tests. This study was limited to the analysis of single nucleotide variants, which are among the least technically challenging variants to predict. The same
15 issues of overall library complexity and depth of coverage with unique reads can be expected to impact the results of NGS testing for other variant classes, namely small insertions and deletions, copy number variants, and structural variants, as has been anticipated by earlier studies 10, 31. Finally, Figure 1 depicts a mathematical model of the depth of coverage with unique reads that would be required to confidently detect variants at a given VAF under idealized circumstances. Although it would be convenient to use these data to propose specific thresholds as standard requirements for NGS tests, it is not possible to anticipate all of the real-world factors that can affect the sensitivity of variant prediction in different areas of the genome using different NGS assays for different purposes (ie, constitutional versus somatic variant detection). Thus, we recommend that specific thresholds should be established for individual NGS assays by the laboratory during the test validation process.
CONCLUSIONS Library complexity is a vital metric that should be assessed in every instance of NGS-based diagnostic testing. Data produced from libraries with poor depth of coverage with unique reads will suffer from limited sensitivity and high variability due to sampling error. These drawbacks will not be apparent in amplicon-based NGS results unless the assay includes a discrete method (eg, implementation of UMI) for tracking unique versus duplicate reads.
16
17 ACKNOWLEDGEMENTS We thank Xiaopei Zhu and Jessica Hoisington-Lopez for expert technical assistance in sample preparation, NGS library preparation, and sequencing and the Genomics and Pathology Services team at Washington University School of Medicine, who performed initial analyses of the clinical cases described in this study.
18 REFERENCES 1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
Eijkelenboom A, Kamping EJ, Kastner-van Raaij AW, Hendriks-Cornelissen SJ, Neveling K, Kuiper RP, Hoischen A, Nelen MR, Ligtenberg MJ, Tops BB: Reliable Next-Generation Sequencing of Formalin-Fixed, Paraffin-Embedded Tissue Using Single Molecule Tags. J Mol Diagn 2016, 18:851-863. Kerick M, Isau M, Timmermann B, Sultmann H, Herwig R, Krobitsch S, Schaefer G, Verdorfer I, Bartsch G, Klocker H, Lehrach H, Schweiger MR: Targeted high throughput sequencing in clinical cancer settings: formaldehyde fixed-paraffin embedded (FFPE) tumor tissues, input amount and tumor heterogeneity. BMC Med Genomics 2011, 4:68. Aird D, Ross MG, Chen WS, Danielsson M, Fennell T, Russ C, Jaffe DB, Nusbaum C, Gnirke A: Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries. Genome Biol 2011, 12:R18. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell TJ, Kernytsky AM, Sivachenko AY, Cibulskis K, Gabriel SB, Altshuler D, Daly MJ: A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 2011, 43:491-498. Head SR, Komori HK, LaMere SA, Whisenant T, Van Nieuwerburgh F, Salomon DR, Ordoukhanian P: Library construction for next-generation sequencing: overviews and challenges. Biotechniques 2014, 56:61-64, 66, 68, passim. Oyola SO, Otto TD, Gu Y, Maslen G, Manske M, Campino S, Turner DJ, Macinnis B, Kwiatkowski DP, Swerdlow HP, Quail MA: Optimizing Illumina nextgeneration sequencing library preparation for extremely AT-biased genomes. BMC Genomics 2012, 13:1. Tilak MK, Botero-Castro F, Galtier N, Nabholz B: Illumina Library Preparation for Sequencing the GC-Rich Fraction of Heterogeneous Genomic DNA. Genome Biol Evol 2018, 10:616-622. Aigrain L, Gu Y, Quail MA: Quantitation of next generation sequencing library preparation protocol efficiencies using droplet digital PCR assays - a systematic comparison of DNA library preparation kits for Illumina sequencing. BMC Genomics 2016, 17:458. Roy S, Coldren C, Karunamurthy A, Kip NS, Klee EW, Lincoln SE, Leon A, Pullambhatla M, Temple-Smolkin RL, Voelkerding KV, Wang C, Carter AB: Standards and Guidelines for Validating Next-Generation Sequencing Bioinformatics Pipelines: A Joint Recommendation of the Association for Molecular Pathology and the College of American Pathologists. J Mol Diagn 2018, 20:4-27. So AP, Vilborg A, Bouhlal Y, Koehler RT, Grimes SM, Pouliot Y, Mendoza D, Ziegle J, Stein J, Goodsaid F, Lucero MY, De La Vega FM, Ji HP: A robust targeted sequencing approach for low input and variable quality DNA from clinical samples. NPJ Genom Med 2018, 3:2. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009, 25:2078-2079.
19 12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
Kivioja T, Vaharautio A, Karlsson K, Bonke M, Enge M, Linnarsson S, Taipale J: Counting absolute numbers of molecules using unique molecular identifiers. Nat Methods 2011, 9:72-74. Mertes F, Elsharawy A, Sauer S, van Helvoort JM, van der Zaag PJ, Franke A, Nilsson M, Lehrach H, Brookes AJ: Targeted enrichment of genomic DNA regions for next-generation sequencing. Brief Funct Genomics 2011, 10:374386. Bennett CW, Berchem G, Kim YJ, El-Khoury V: Cell-free DNA and nextgeneration sequencing in the service of personalized medicine for lung cancer. Oncotarget 2016, 7:71013-71035. Liu HE, Triboulet M, Zia A, Vuppalapaty M, Kidess-Sigal E, Coller J, Natu VS, Shokoohi V, Che J, Renier C, Chan NH, Hanft VR, Jeffrey SS, Sollier-Christen E: Workflow optimization of whole genome amplification and targeted panel sequencing for CTC mutation detection. NPJ Genom Med 2017, 2:34. Tandon B, Hagemann IS, Maluf HM, Pfeifer JD, Al-Kateb H: Association of LiFraumeni Syndrome With Small Cell Carcinoma of the Ovary, Hypercalcemic Type and Concurrent Pleomorphic Liposarcoma of the Cervix. Int J Gynecol Pathol 2017, 36:593-599. Cottrell CE, Al-Kateb H, Bredemeyer AJ, Duncavage EJ, Spencer DH, Abel HJ, Lockwood CM, Hagemann IS, O’Guin SM, Burcea LC, Sawyer CS, Oschwald DM, Stratman JL, Sher DA, Johnson MR, Brown JT, Cliften PF, George B, McIntosh LD, Shrivastava S, Nguyen TT, Payton JE, Watson MA, Crosby SD, Head RD, Mitra RD, Nagarajan R, Kulkarni S, Seibert K, Virgin HW, Milbrandt J, Pfeifer JD: Validation of a Next-Generation Sequencing Assay for Clinical Molecular Oncology. The Journal of Molecular Diagnostics 2014, 16:89-105. Duncavage EJ, Jacoby MA, Chang GS, Miller CA, Edwin N, Shao J, Elliott K, Robinson J, Abel H, Fulton RS, Fronick CC, O'Laughlin M, Heath SE, Brendel K, Saba R, Wartman LD, Christopher MJ, Pusic I, Welch JS, Uy GL, Link DC, DiPersio JF, Westervelt P, Ley TJ, Trinkaus K, Graubert TA, Walter MJ: Mutation Clearance after Transplantation for Myelodysplastic Syndrome. N Engl J Med 2018, 379:1028-1041. Jacoby MA, Duncavage EJ, Chang GS, Miller CA, Shao J, Elliott K, Robinson J, Fulton RS, Fronick CC, O'Laughlin M, Heath SE, Pusic I, Welch JS, Link DC, DiPersio JF, Westervelt P, Ley TJ, Graubert TA, Walter MJ: Subclones dominate at MDS progression following allogeneic hematopoietic cell transplant. JCI Insight 2018, 3. Uy GL, Duncavage EJ, Chang GS, Jacoby MA, Miller CA, Shao J, Heath S, Elliott K, Reineck T, Fulton RS, Fronick CC, O'Laughlin M, Ganel L, Abboud CN, Cashen AF, DiPersio JF, Wilson RK, Link DC, Welch JS, Ley TJ, Graubert TA, Westervelt P, Walter MJ: Dynamic changes in the clonal structure of MDS and AML in response to epigenetic therapy. Leukemia 2017, 31:872-881. Tabatabaeifar S, Larsen MJ, Thomassen M, Larsen SR, Kruse TA, Sorensen JA: The Optimal Sequencing Depth of Tumor Biopsies for Identifying Clonal Cell Populations. J Mol Diagn 2019, 21:790-795. Paasinen-Sohns A, Koelzer VH, Frank A, Schafroth J, Gisler A, Sachs M, Graber A, Rothschild SI, Wicki A, Cathomas G, Mertz KD: Single-Center Experience
20
23.
24. 25. 26.
27.
28.
29.
30.
31.
with a Targeted Next Generation Sequencing Assay for Assessment of Relevant Somatic Alterations in Solid Tumors. Neoplasia 2017, 19:196-206. Guibert N, Hu Y, Feeney N, Kuang Y, Plagnol V, Jones G, Howarth K, Beeler JF, Paweletz CP, Oxnard GR: Amplicon-based next-generation sequencing of plasma cell-free DNA for detection of driver and resistance mutations in advanced non-small cell lung cancer. Ann Oncol 2018. Chang F, Li MM: Clinical application of amplicon-based next-generation sequencing in cancer. Cancer Genet 2013, 206:413-419. Fleischhacker M, Schmidt B: Circulating nucleic acids (CNAs) and cancer--a survey. Biochim Biophys Acta 2007, 1775:181-232. Johansson G, Andersson D, Filges S, Li J, Muth A, Godfrey TE, Stahlberg A: Considerations and quality controls when analyzing cell-free tumor DNA. Biomol Detect Quantif 2019, 17:100078. Lapin M, Oltedal S, Tjensvoll K, Buhl T, Smaaland R, Garresori H, Javle M, Glenjen NI, Abelseth BK, Gilje B, Nordgard O: Fragment size and level of cellfree DNA provide prognostic information in patients with advanced pancreatic cancer. J Transl Med 2018, 16:300. Newman AM, Bratman SV, To J, Wynne JF, Eclov NC, Modlin LA, Liu CL, Neal JW, Wakelee HA, Merritt RE, Shrager JB, Loo BW, Jr., Alizadeh AA, Diehn M: An ultrasensitive method for quantitating circulating tumor DNA with broad patient coverage. Nat Med 2014, 20:548-554. Zhang Y, Xu Y, Zhong W, Zhao J, Chen M, Zhang L, Li L, Wang M: Total DNA input is a crucial determinant of the sensitivity of plasma cell-free DNA EGFR mutation detection using droplet digital PCR. Oncotarget 2017, 8:5861-5873. Ascierto PA, Bifulco C, Palmieri G, Peters S, Sidiropoulos N: Preanalytic Variables and Tissue Stewardship for Reliable Next-Generation Sequencing (NGS) Clinical Analysis. J Mol Diagn 2019, 21:756-767. Wood HM, Belvedere O, Conway C, Daly C, Chalkley R, Bickerdike M, McKinley C, Egan P, Ross L, Hayward B, Morgan J, Davidson L, MacLennan K, Ong TK, Papagiannopoulos K, Cook I, Adams DJ, Taylor GR, Rabbitts P: Using next-generation sequencing for high resolution multiplex analysis of copy number variation from nanogram quantities of DNA from formalin-fixed paraffin-embedded specimens. Nucleic Acids Res 2010, 38:e151.
21 FIGURE LEGENDS
Figure 1. Modeling theoretical limits of variant detection based on depth of coverage using binomial sampling statistics. The confidence of variant prediction was determined for variants with allelic fractions equal to 2.5%, 5.0%, 10.0%, and 15% for depths of coverage with unique reads ranging from 1 to 500 reads (A). Similarly, the depth coverage was determined with unique reads that would be required to achieve 95% confidence (dashed line) for variants ranging in allelic fraction from 0.1 to 50.0% (B). Power and sample sizes were calculated for a binomial proportion using complimentary log parameterization given a sequencing error rate of 0.01%.
Figure 2. The impact of library complexity on assay sensitivity. Nextgeneration sequencing (NGS) tests are designed to survey the genetic makeup of tumor specimens that are typically comprised of a mixture of normal cells and neoplastic cells. The neoplastic component may be further subdivided into genetically diverse subclones. In this schematic, 10 unique input molecules were isolated and ligated to unique molecular identifiers. When a standard input of five molecules is used in library preparation, the theoretical limit of sensitivity is 1 in 5 and the variant is detected at a fraction similar to its true proportion in the original specimen. When a library is constructed from a minimum input of two molecules, the theoretical limit of sensitivity is reduced (1/2 rather than 1/5), and the variant is missed due to sampling error. In absence of UMIs, the total coverage of the
22 standard input and the minimum input libraries are identical and the vast difference in the sensitivity of variant detection would not be immediately apparent.
Figure 3. Relationship between total and unique read coverage. Total and unique read coverage of the designated analytical space was calculated for 41 NGS libraries generated from cell line DNA inputs of 50ng (green, n=10), 10ng (yellow, n=11), 5ng (red, n=10), or 2ng (blue, n=10). At the high total depths of coverage achieved in this experiment (average: 23,685±11,210), there is no significant correlation between total and unique read coverage (A). Data from NGS libraries prepared from 50ng were randomly down sampled (ds) to ~1,750x and ~10,000x total read coverage, and the percent of reads classified as unique was calculated for down sampled and non–down sampled libraries (B). As total coverage increases, the percent of reads classified as unique decreases (average of 26.7%±2.6%, 20.8%±4.6%, and 10.1%±2.58% for ds-1,750, ds10,000x, and original, respectively). The ends of boxes indicate the upper and lower quartiles. The median is indicated by the middle line. Whiskers extend to 1.5 times the interquartile range.
Figure 4. Average unique coverage at 10,000x total read coverage. Forty-one next-generation sequencing (NGS) libraries prepared from cell line material were sequenced and randomly down-sampled to approximately 10,000x total read coverage of the designated analytical space (average: 10,007±11.4 reads) to
23 facilitate intra-sample comparison. Despite some noise, the average unique depth of coverage tended to decrease with reduced DNA input (average of 776.2±569.2x, 971.0±725.4x, 1,263.6±784.4x, and 281.7±458.6x coverage for 2ng, 5ng, 10ng, and 50ng, respectively). The ends of boxes indicate the upper and lower quartiles, and the median is indicated by the middle line; whiskers extend to 1.5 times the interquartile range.
Figure 5. The impact of reduced DNA input on variant identification using standardized specimens. DNA from an engineered cell line with a heterozygous KRAS p.G13D missense variant was mixed in proportions of 1:1, 1:8, 1:16, and 1:32 to produce variants with approximate variant of a specified allelic fractions (VAFs) of 25%, 5%, 3%, and 1.5%, respectively (marked with black bars on the plot). Triplicate libraries were prepared and sequenced using either 50ng or 2ng of the mixed DNA as input, and variants were predicted. Variance between the triplicate estimates was greater for 2ng libraries (blue) compared to 50ng, standard input libraries (green). Coefficients of variation for VAF estimates were calculated as follows: 74.45% versus 8.50% for 2ng and 50ng of 1:1 mix, respectively; 170.73% versus 40.44% for 2ng and 50ng of 1:8 mix, respectively; 148.38% versus 65.15% for 2ng and 50ng of 1:16 mix, respectively; 20.94% versus 72.05% for 2ng and 50ng of 1:32 mix, respectively. Note that though the variance is actually lower for VAFs estimated from the 2ng libraries of the 1:32 mix compared to the 50ng libraries, the VAFs estimated from the 2ng libraries fell to nearly undetectable levels.
24
Figure 6. The impact of reduced DNA input on variant identification using formalin-fixed, paraffin embedded tumor specimens. Three de-identified, formalin-fixed, paraffin-embedded tumor specimens were selected for testing based on the presence of known, clinically significant variants: C1 with chr2:209113113G>T (IDH1 p.R132S, VAF=9.0%), C2 with chr12:25398284C>T (KRAS p.G12D, VAF=15.0%), and C3 with chr3:178936082G>A (PIK3CA p.E542K, VAF=3.0%). Variant of a specified allelic fractions (VAFs) estimated during previous clinical testing are marked with black bars. Triplicate libraries were prepared and sequenced using either 50ng or 2ng of DNA as input, and variants were predicted. Variance between the triplicate estimates was greater for 2ng libraries (blue) compared to 50ng, standard input libraries (green). Coefficients of variation were calculated as follows: 79.41% versus 10.48% for 2ng and 50ng libraries of C1, respectively; 140.96% versus 8.39% for 2ng and 50ng libraries of C2, respectively; 153.04% and 66.69% for 2ng and 50ng libraries of C3, respectively.
A
1.00 0.95
Confidence
0.75
0.50 Variant Allelic Fraction 2.5% 5.0% 10.0% 15.0%
0.25
0.00 0
100
200
300
400
500
Depth of Read Coverage 1000
Depth of Coverage
B
300
100
30 50%
40%
30%
20%
Variant Allelic Fraction
10%
0%
DNA isolation
Barcode Ligation
PCR amplification
Mixture of tumor and normal tissue
Minimum input
Standard input
Tissue isolation
Read alignment
Variant prediction
A (40%) T (60%)
Sensitivity = 1/5
T (100%)
Sensitivity = 1/2
Average Unique Read Coverage
A
DNA Input 2ng 5ng 10ng 50ng
4000
3000
r2 = 0.0537
2000
1000
0
10000
20000
30000
40000
Average Total Read Coverage
Percent Unique
B
0.30
0.25
0.20
0.15
0.10
ds-1,750x
ds-10,000x
Total Read Coverage
Original
Average Unique Read Coverage
2000
1000
2
5
10
DNA Input (ng)
50
VAF of KRAS p.G13D Variant
DNA Input 2 ng 50 ng
10.0%
1.0%
0.1%
1:1
1:8
1:16
Dilution Ratio
1:32
VAF of Variant
DNA Input 2 ng 50 ng
10.0%
1.0%
0.1%
C1
C2
Specimen
C3
S1525-1578(20)30042-8
DOI:
https://doi.org/10.1016/j.jmoldx.2020.02.003
Reference:
JMDI 895
To appear in:
The Journal of Molecular Diagnostics
Received Date: 2 October 2019 Revised Date:
7 January 2020
Accepted Date: 20 February 2020
Please cite this article as: McNulty SN, Mann PR, Robinson JA, Duncavage EJ, Pfeifer JD, The Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection, The Journal of Molecular Diagnostics (2020), doi: https://doi.org/10.1016/j.jmoldx.2020.02.003. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Copyright © 2020 Published by Elsevier Inc. on behalf of the American Society for Investigative Pathology and the Association for Molecular Pathology.
1
The Impact of Reducing DNA Input on Next-Generation Sequencing Library Complexity and Variant Detection Samantha N. McNulty,* Patrick R. Mann,* Joshua A. Robinson,† Eric J. Duncavage,* and John D. Pfeifer* From the Department of Pathology and Immunology,* and the Department of Medicine, Division of Oncology,† Washington University School of Medicine, Saint Louis, Missouri Corresponding author: Samantha N. McNulty, Ph.D., Department of Pathology and Immunology, Washington University School of Medicine, Campus Box 8118, 600 S. Euclid Ave, St. Louis, MO 63110. E-mail: [email protected] Running head: Impact of DNA input in clinical NGS. Funding: Supported by the Department of Pathology and Immunology, Washington University School of Medicine. Disclosures: None declared.
2 ABSTRACT PCR amplification, a key step in next-generation sequencing (NGS) library construction, can generate an unlimited amount of product from limited input; however, it cannot create more information than was present in the original template. Thus, NGS libraries can be made from very little DNA, but reducing input may compromise assay sensitivity in ways that are difficult to ascertain unless library complexity (ie, the number of unique DNA molecules represented in the library) and depth of coverage with unique sequence reads (those derived from input DNA molecules) versus duplicate sequence reads (resulting from over-amplification of particular molecules) are discretely measured. We performed a series of experiments to explore the impact of low DNA input on an amplicon-based NGS assay using unique molecular identifiers to track unique versus duplicate reads. At high sequencing depths, unique and total (unique plus duplicate) coverage are not well correlated, so increasing the number of sequenced reads does not necessarily improve sensitivity. Unique coverage depth tends to improve with more input, but improvements are not consistent. Fluctuations in library complexity complicated variant detection using both standardized and clinical specimens, often resulting in technical replicates with vastly different estimates of variant allelic fraction. We conclude that depth of coverage with unique reads must be tracked in clinical NGS to ensure that sensitivity and accuracy are maintained.
3 INTRODUCTION Ideally, a next-generation sequencing (NGS) library would perfectly represent the DNA present in the biological specimen from which it was derived, and the sequencing process would sample that library in a manner that is sufficient to detect any variants present in the original specimen. In this idealized scenario, it is possible to mathematically model the number of reads that would need to be sampled from a given genomic locus (ie, the depth of mapped read coverage) to detect a variant of a specified allelic fraction (VAF) with a specified level of confidence. High depth of coverage is particularly important for the detection of low VAF variants like those routinely encountered in the context of testing for cancer and other somatic diseases (Figure 1)1, 2. In actual clinical practice, however, various sources of error/bias in sample preparation, NGS library construction, amplification, and/or sequencing can enrich or deplete certain DNA sequences in ways that are difficult to predict a priori 3-7. Thus, it is essential to distinguish those sequence reads that trace back to DNA molecules that were present in the original specimen (unique sequence reads) from those reads that arose from over-amplification of particular DNA molecules during the process of library construction (duplicate reads; Figure 2) 4, 8-10
. The inadvertent inclusion of duplicate reads in the bioinformatic process of
variant prediction can skew the sampling of input DNA molecules and cause i) failure to detect sequence variants that were present in the original specimen, ii) over- or under-representation of particular variants, or iii) false positive variant
4 calls resulting from PCR errors that are propagated through library preparation and sequencing. Duplicate sequence reads are easily recognized in the context of whole genome shotgun or hybridization-capture–based NGS assays. Using these methods, input DNA is randomly fragmented during the process of library preparation, so unique sequence reads are unlikely to be identical to one another. Duplicate reads are easily identified based on shared mapping coordinates (for example, using PicardTools, Broad Institute, Cambridge, MA) 11, and disregarded in downstream bioinformatic analyses. Conversely, in the setting of amplicon-based NGS, all reads derived from a given amplicon (those amplified from the same PCR primer set) are expected to be identical, so finding duplicate reads is not straightforward. The best way to assess uniqueness in an amplification-based NGS workflow is to ligate unique molecular identifiers (UMIs) to input DNA molecules prior to library construction (Figure 2); sequence reads with identical UMI can be counted as duplicates in the bioinformatic pipeline since they were derived from the same input molecule 12. There are advantages and disadvantages to both hybridization-capture and amplification-based enrichment strategies, but amplicon-based tests are often favored in situations where DNA input is limited 13, which is precisely when library complexity (ie, the number of unique DNA molecules represented in an NGS library) can become compromised and the number of unique sequence reads mapped to a particular locus can drop to the point that the sensitivity of variant detection suffers. Unfortunately, UMI-based methods add a layer of cost
5 and complexity to both library construction and bioinformatic analyses, so they have not been uniformly adopted, particularly in the clinical laboratory setting. Many product manufacturers (eg, TruSeq Custom Amplicon Low Input Kit, Illumina, Carlsbad, CA; CleanPlex Custom NGS Panels, Paragon Genomics, Hayward, CA; DNA Library Prep Kit for Illumina Sequencing, Applied Biological Materials, Richmond, BC, Canada) and laboratories 10, 14, 15 tout amplicon-based assays with low DNA input, even in the absence of methods that facilitate the tracking of unique versus duplicate reads. Although the important distinction between unique versus total (unique plus duplicate) depth of coverage has been modeled in the past 1, 2,this study experimentally demonstrated the impact of reduced input DNA input on unique read coverage and variant prediction in the setting of an amplicon-based NGS assay. Libraries were generated from as little as 2ng of genomic DNA using a UMI-based preparation protocol, sequenced to high total depth of coverage, and analyzed using a UMI-aware bioinformatic pipeline to track unique versus duplicate sequence reads. The results show that library complexity, and consequently depth of coverage with unique sequence reads, suffers when DNA input is reduced, but not in consistent or predictable ways. Variations in library complexity impaired accurate variant prediction in the context of both standardized and clinical specimens. Our results stress the importance of measuring unique versus duplicate reads in all cases, particularly when DNA input is minimal.
6 MATERIALS AND METHODS
Analytical specimens DNA was isolated from two engineered cell lines, one with a knock-in BRAF p.V600E variant (HD 200-002, Horizon Diagnostics, Cambridge, UK) and another with knock-in KRAS p.G13D and PIK3CA p.H1047R variants (HD 200003, Horizon Diagnostics) as previously described 16. DNA from the two cell lines was mixed to generate variants with VAFs ranging from approximately 1% to 25% (Supplemental Tables S1 and S2).
Clinical Specimens Three de-identified, formalin-fixed, paraffin-embedded tumor specimens were selected for sequencing based on previous findings of clinically significant variants present at VAF <20% when tested in routine practice using a validated, hybridization-capture based NGS assay run in a CAP/CLIA certified laboratory 17. Specimen C1 had an IDH1 p.R132S variant at VAF=9%, C2 had a KRAS p. G12D variant at VAF=12%, and C3 had a PIK3CA p.E542K variant at VAF=3%.
Library preparation and sequencing DNA was isolated from cell lines and clinical specimens using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) and quantified using a Qubit fluorometer with a dsDNA HS Assay Kit (ThermoFisher Scientific, Indianapolis, IN). Sequencing libraries were prepared from 2ng, 5ng, 10ng, or 50ng of pre-mixed
7 cell line DNA (Supplemental Tables S1 and S2) or tumor DNA (Supplemental Table S3) using the HaloPlex HS Target Enrichment System (Agilent Technologies, Santa Clara, CA) as previously described 18-20. Designated quantities of input DNA were digested using a mixture of restriction endonucleases included in the Agilent HaloPlex HS kit, denatured, and hybridized to a biotinylated and UMI-barcoded probe library designed to target the coordinates of the AmpliSeq Cancer Panel v2 for Illumina assay (3.6kb capture space, ThermoFisher Scientific). Biotinylated targets were isolated by hybridization with streptavidin beads and amplified for 23 PCR cycles using primers tailed with Illumina sequencing motifs. Prepared libraries were pooled and sequenced on an Illumina NextSeq or Illumina NovaSeq6000 (Illumina). The number of reads generated from each library is reported in Supplemental Tables S1, S2, and S3.
Bioinformatic analysis Demultiplexed paired-end sequence reads were aligned to the human reference assembly (GRCh37-lite-build37) using BWA-MEM (version 7.9.a; parameters: -M). As previously described 20, variants were predicted from aligned BAM files using Barcrawler, a custom, Java-based GATK locus walker (available from: https://github.com/abelhj/gatk/tree/master/public; parameters: -mmq 20 mbq 20 -minCtBC 1 -dcov 500000 -discardN 1 -minOffset 5 -maxNM 3; accessed 2018-02-07). No further filtering or UMI-based error correction was imposed. Barcrawler outputs the depth of coverage with total and unique sequence reads
8 at each base position in the region of interest (ROI); these values were used to calculate the average coverage of the assay ROI. Aligned bam files were down sampled to desired coverage thresholds using Picard DownsampleSam (http://broadinstitute.github.io/picard, version 2.21.4) and re-analyzed as necessary.
Regulatory Approval This study was approved by the Washington University School of Medicine Human Studies Committee.
RESULTS
Unique Versus Total Read Coverage NGS libraries were constructed from decreasing quantities of DNA using a UMI-based PCR enrichment protocol. 50ng (standard input), 10ng, 5ng, and 2ng input libraries were sequenced to an average total read coverage of 23,685±11,210 (Supplemental Table S1). A UMI-aware bioinformatic pipeline was used to discretely track unique versus duplicate reads and predict variants. During the course of sequencing, the first read encountered with a given UMI sequence is considered unique and all subsequent reads containing that UMI are flagged as duplicates. Thus, at low depth of coverage, a high percentage of reads will be classified as unique. The percent of total reads that are classified as unique is expected to decrease as total depth of coverage increases since there
9 are a finite number of input DNA molecules to survey, but libraries can be sequenced to infinite total depth. At the high depth of coverage achieved in this experiment, there did not appear to be any correlation between depth of unique versus total read coverage (Figure 3A). NGS libraries produced from 50ng input DNA were computationally down sampled to explore the impact of total depth of coverage on the percent of captured reads that are classified as unique. The percent of unique reads detected was highest for libraries down sampled to 1,750x total depth of coverage and decreased as total depth of coverage increased (Figure 3B).
Relationship Between Unique Read Coverage and DNA Input Next, NGS libraries prepared from 50ng, 10ng, 5ng, and 2ng DNA input were compared. The average total depth of coverage of the ROI varied between libraries (range of approximately 10,126 to 44,803x), so alignments were down sampled to approximately 10,000x total depth of coverage (mean: 10,007±11.4x) to facilitate inter-sample comparison. As expected, the depth of coverage with unique reads generally increased with increasing DNA input (Figure 4).
Impact of Low Library Complexity on Variant Detection Using Standard Specimens To explore the impact of the depth of coverage with unique reads on variant prediction, DNA was isolated from two isogenic cell lines containing key heterozygous variants and mixed to produce variants with allelic fractions ranging
10 from 25% to approximately 1% (Supplemental Table S2). NGS libraries were prepared and sequenced in triplicate using either the standard input of 50ng or a reduced input of 2ng, and variants were predicted without UMI-based error correction (Figure 5). Variants from a 1:1 mixture of the two cell lines (approximately 25% VAF) were detected even when the DNA input was reduced. However, the variability between the VAFs calculated from the triplicate libraries was higher in the 2ng input libraries compared to the 50ng input libraries for variants at a range of expected VAFs. For example, the KRAS variant expected to be present at a VAF of approximately 6% in the 1:8 DNA mixture was detected at VAFs ranging from 0.04% to 10.83% in the 2ng libraries. Many of the low VAF variants were consistently detected far below their anticipated fraction in the 2ng input libraries, whereas the VAFs calculated from 50ng libraries were closer to the expected value.
Impact of Low Library Complexity on Variant Detection Using Clinical Specimens The same experiment was repeated, this time using FFPE tumor specimens that had been shown to harbor low VAF, clinically significant variants during routine, clinical testing using a validated hybridization-capture–based NGS assay (C1: IDH1 p.R132S, C2: KRAS p. G12D, and C3: PIK3CA p.E542K). NGS libraries were prepared and sequenced in triplicate using either 50ng or 2ng of input DNA (Supplemental Table S3). These libraries were sequenced to very high total depths of coverage (average: 397,739±234,664x). The average depth
11 of coverage with unique reads was higher for the 50ng input libraries compared to the 2ng input libraries (11,769±841x versus 4,500±1,892x, respectively). Although the variants were reliably detected at a consistent VAF using 50ng of input DNA, the VAFs of the variants identified from 2ng input libraries showed a greater degree of variation (Figure 6). For example, the KRAS p.G12D variant present in specimen C2 was estimated between 21.60% to 25.50% VAF from the 50ng input libraries, but VAF estimates ranged from 0.48% to 47.10% in the 2ng input libraries. This trend was consistent between the three FFPE specimens (Figure 6).
DISCUSSION Library complexity is a direct reflection of the number of independent genomes present in the input material used in library construction. It is certainly possible to prepare a NGS library from extremely low amounts of DNA or RNA template, but the information contained in that library may be limited in ways that are not obvious unless the library preparation, sequencing, and bioinformatic workflow includes steps to facilitate the discrete measurement of unique read coverage. By sequencing a series of NGS libraries prepared from decreasing amounts of input DNA extracted from controlled cell lines and previously characterized clinical specimens, it was shown that the total depth of coverage (unique plus duplicate reads) does not provide an adequate measure of library complexity or depth coverage with unique reads at clinically significant genomic loci.
12 A prepared NGS library can be sequenced to any desired depth of coverage. High total depth of coverage is needed to ensure that a sufficient number of unique molecules have been sampled; however, once the unique molecules have been sampled, further sequencing yields no additional information. This explains the lack of correlation between the depth of coverage with total versus unique reads at the high total depths of coverage achieved in this experiment. The same lack of correlation that can be modelled theoretically 1, 2
holds true for the clinical NGS assay in routine clinical practice 21. Although simply increasing the DNA input may improve library complexity
to some extent, this is not possible in many clinical applications. Amplificationbased enrichment strategies are specifically preferred in instances where input material is limited 22-24, including for small biopsies, evaluation of minimum residual disease, cell-free DNA (cfDNA) testing, etc. For example, the total quantity of cfDNA isolated from a peripheral blood sample is sometimes quite limited (often less than 50ng in 1mL of blood, depending on the nature and extent of the disease 25-27). Numerous vendors market kits for use with nanogram quantities of cfDNA, and many protocols include additional amplification steps to generate a sufficient mass of DNA for library preparation and NGS sequencing. As shown, depth of coverage with unique reads becomes unpredictable when limited amounts of input DNA are used for library preparation, so it must be carefully tracked to ascertain the sensitivity of each individual assay. This is particularly critical in cases where treatment decisions are based on results
13 generated from the analysis of a few nanograms of DNA (for instance, detection of specific mutations or estimate of cancer burden 28, 29). Increasing the DNA input may not adequately boost unique library complexity for all specimens. There is a great deal of variability in the quality of DNA isolated from clinical specimens depending on the tissue source, type, and degree of fixation, age of the specimen, storage conditions, etc 30. Two libraries made from the same amount of DNA may, therefore, have markedly different complexity even if the input nucleic acid mass is the same. Furthermore, biases in the amplification and sequencing of certain genomic loci can lead to regions of poor depth of coverage with unique reads even when the input material is plentiful and library complexity is adequate overall 3, 5-7. Ultimately, the concern for depth of coverage with unique reads relates to the concern for assay sensitivity. Theoretically, the minimum sensitivity of an assay is one divided by the number of unique molecules sampled during sequencing. If this sample size is low, it may not be mathematically possible to achieve the sensitivity that is expected from a given diagnostic test even if every step of library preparation is ideal, every read maps perfectly to the reference, and every variant is correctly called. In reality, library construction does not produce idealized results, and this affects variant detection in ways that are not readily predictable sample-to-sample and assay-to-assay. Even when the expected variants are identified, the estimated VAF is subject to sampling error and that error becomes more significant at lower depth of coverage with unique reads, as we have shown using both model and actual clinical specimens. Thus,
14 library complexity is a metric that should be tracked in every assay for every specimen to ensure that a given test meets or exceeds any validated requirements. Here, we focused on the measurement of library complexity in the context of amplicon-based NGS assays. In this study, we stress the need to track unique versus duplicate reads in the context of amplicon-based NGS not because library complexity is more important in amplicon-based tests compared to shotgun or hybridization-based NGS assays, but rather because it requires specific modifications to the wet-bench protocol that are often not employed (due to complexity, cost, and time constraints). Of course, library complexity is an important metric to measure regardless of the specific library preparation or enrichment method employed. There is no doubt that the measures of the impact of library complexity on the NGS results that were describe (eg, optimal depth of coverage in Figures 3-4, variability in measured VAF by input DNA mass in Figures 5-6, etc.) here reflect the technical details of the NGS assay employed (eg, DNA extraction protocol, library preparation method, bioinformatics pipeline, etc.). It can be expected that different NGS assay designs may cause differences in the magnitude of these impacts. However, since the underlying issues described regarding the effect of library complexity on NGS test results cannot be avoided, laboratories should systematically explore the influence they have on their specific amplicon-based NGS tests. This study was limited to the analysis of single nucleotide variants, which are among the least technically challenging variants to predict. The same
15 issues of overall library complexity and depth of coverage with unique reads can be expected to impact the results of NGS testing for other variant classes, namely small insertions and deletions, copy number variants, and structural variants, as has been anticipated by earlier studies 10, 31. Finally, Figure 1 depicts a mathematical model of the depth of coverage with unique reads that would be required to confidently detect variants at a given VAF under idealized circumstances. Although it would be convenient to use these data to propose specific thresholds as standard requirements for NGS tests, it is not possible to anticipate all of the real-world factors that can affect the sensitivity of variant prediction in different areas of the genome using different NGS assays for different purposes (ie, constitutional versus somatic variant detection). Thus, we recommend that specific thresholds should be established for individual NGS assays by the laboratory during the test validation process.
CONCLUSIONS Library complexity is a vital metric that should be assessed in every instance of NGS-based diagnostic testing. Data produced from libraries with poor depth of coverage with unique reads will suffer from limited sensitivity and high variability due to sampling error. These drawbacks will not be apparent in amplicon-based NGS results unless the assay includes a discrete method (eg, implementation of UMI) for tracking unique versus duplicate reads.
16
17 ACKNOWLEDGEMENTS We thank Xiaopei Zhu and Jessica Hoisington-Lopez for expert technical assistance in sample preparation, NGS library preparation, and sequencing and the Genomics and Pathology Services team at Washington University School of Medicine, who performed initial analyses of the clinical cases described in this study.
18 REFERENCES 1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
Eijkelenboom A, Kamping EJ, Kastner-van Raaij AW, Hendriks-Cornelissen SJ, Neveling K, Kuiper RP, Hoischen A, Nelen MR, Ligtenberg MJ, Tops BB: Reliable Next-Generation Sequencing of Formalin-Fixed, Paraffin-Embedded Tissue Using Single Molecule Tags. J Mol Diagn 2016, 18:851-863. Kerick M, Isau M, Timmermann B, Sultmann H, Herwig R, Krobitsch S, Schaefer G, Verdorfer I, Bartsch G, Klocker H, Lehrach H, Schweiger MR: Targeted high throughput sequencing in clinical cancer settings: formaldehyde fixed-paraffin embedded (FFPE) tumor tissues, input amount and tumor heterogeneity. BMC Med Genomics 2011, 4:68. Aird D, Ross MG, Chen WS, Danielsson M, Fennell T, Russ C, Jaffe DB, Nusbaum C, Gnirke A: Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries. Genome Biol 2011, 12:R18. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell TJ, Kernytsky AM, Sivachenko AY, Cibulskis K, Gabriel SB, Altshuler D, Daly MJ: A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 2011, 43:491-498. Head SR, Komori HK, LaMere SA, Whisenant T, Van Nieuwerburgh F, Salomon DR, Ordoukhanian P: Library construction for next-generation sequencing: overviews and challenges. Biotechniques 2014, 56:61-64, 66, 68, passim. Oyola SO, Otto TD, Gu Y, Maslen G, Manske M, Campino S, Turner DJ, Macinnis B, Kwiatkowski DP, Swerdlow HP, Quail MA: Optimizing Illumina nextgeneration sequencing library preparation for extremely AT-biased genomes. BMC Genomics 2012, 13:1. Tilak MK, Botero-Castro F, Galtier N, Nabholz B: Illumina Library Preparation for Sequencing the GC-Rich Fraction of Heterogeneous Genomic DNA. Genome Biol Evol 2018, 10:616-622. Aigrain L, Gu Y, Quail MA: Quantitation of next generation sequencing library preparation protocol efficiencies using droplet digital PCR assays - a systematic comparison of DNA library preparation kits for Illumina sequencing. BMC Genomics 2016, 17:458. Roy S, Coldren C, Karunamurthy A, Kip NS, Klee EW, Lincoln SE, Leon A, Pullambhatla M, Temple-Smolkin RL, Voelkerding KV, Wang C, Carter AB: Standards and Guidelines for Validating Next-Generation Sequencing Bioinformatics Pipelines: A Joint Recommendation of the Association for Molecular Pathology and the College of American Pathologists. J Mol Diagn 2018, 20:4-27. So AP, Vilborg A, Bouhlal Y, Koehler RT, Grimes SM, Pouliot Y, Mendoza D, Ziegle J, Stein J, Goodsaid F, Lucero MY, De La Vega FM, Ji HP: A robust targeted sequencing approach for low input and variable quality DNA from clinical samples. NPJ Genom Med 2018, 3:2. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009, 25:2078-2079.
19 12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
Kivioja T, Vaharautio A, Karlsson K, Bonke M, Enge M, Linnarsson S, Taipale J: Counting absolute numbers of molecules using unique molecular identifiers. Nat Methods 2011, 9:72-74. Mertes F, Elsharawy A, Sauer S, van Helvoort JM, van der Zaag PJ, Franke A, Nilsson M, Lehrach H, Brookes AJ: Targeted enrichment of genomic DNA regions for next-generation sequencing. Brief Funct Genomics 2011, 10:374386. Bennett CW, Berchem G, Kim YJ, El-Khoury V: Cell-free DNA and nextgeneration sequencing in the service of personalized medicine for lung cancer. Oncotarget 2016, 7:71013-71035. Liu HE, Triboulet M, Zia A, Vuppalapaty M, Kidess-Sigal E, Coller J, Natu VS, Shokoohi V, Che J, Renier C, Chan NH, Hanft VR, Jeffrey SS, Sollier-Christen E: Workflow optimization of whole genome amplification and targeted panel sequencing for CTC mutation detection. NPJ Genom Med 2017, 2:34. Tandon B, Hagemann IS, Maluf HM, Pfeifer JD, Al-Kateb H: Association of LiFraumeni Syndrome With Small Cell Carcinoma of the Ovary, Hypercalcemic Type and Concurrent Pleomorphic Liposarcoma of the Cervix. Int J Gynecol Pathol 2017, 36:593-599. Cottrell CE, Al-Kateb H, Bredemeyer AJ, Duncavage EJ, Spencer DH, Abel HJ, Lockwood CM, Hagemann IS, O’Guin SM, Burcea LC, Sawyer CS, Oschwald DM, Stratman JL, Sher DA, Johnson MR, Brown JT, Cliften PF, George B, McIntosh LD, Shrivastava S, Nguyen TT, Payton JE, Watson MA, Crosby SD, Head RD, Mitra RD, Nagarajan R, Kulkarni S, Seibert K, Virgin HW, Milbrandt J, Pfeifer JD: Validation of a Next-Generation Sequencing Assay for Clinical Molecular Oncology. The Journal of Molecular Diagnostics 2014, 16:89-105. Duncavage EJ, Jacoby MA, Chang GS, Miller CA, Edwin N, Shao J, Elliott K, Robinson J, Abel H, Fulton RS, Fronick CC, O'Laughlin M, Heath SE, Brendel K, Saba R, Wartman LD, Christopher MJ, Pusic I, Welch JS, Uy GL, Link DC, DiPersio JF, Westervelt P, Ley TJ, Trinkaus K, Graubert TA, Walter MJ: Mutation Clearance after Transplantation for Myelodysplastic Syndrome. N Engl J Med 2018, 379:1028-1041. Jacoby MA, Duncavage EJ, Chang GS, Miller CA, Shao J, Elliott K, Robinson J, Fulton RS, Fronick CC, O'Laughlin M, Heath SE, Pusic I, Welch JS, Link DC, DiPersio JF, Westervelt P, Ley TJ, Graubert TA, Walter MJ: Subclones dominate at MDS progression following allogeneic hematopoietic cell transplant. JCI Insight 2018, 3. Uy GL, Duncavage EJ, Chang GS, Jacoby MA, Miller CA, Shao J, Heath S, Elliott K, Reineck T, Fulton RS, Fronick CC, O'Laughlin M, Ganel L, Abboud CN, Cashen AF, DiPersio JF, Wilson RK, Link DC, Welch JS, Ley TJ, Graubert TA, Westervelt P, Walter MJ: Dynamic changes in the clonal structure of MDS and AML in response to epigenetic therapy. Leukemia 2017, 31:872-881. Tabatabaeifar S, Larsen MJ, Thomassen M, Larsen SR, Kruse TA, Sorensen JA: The Optimal Sequencing Depth of Tumor Biopsies for Identifying Clonal Cell Populations. J Mol Diagn 2019, 21:790-795. Paasinen-Sohns A, Koelzer VH, Frank A, Schafroth J, Gisler A, Sachs M, Graber A, Rothschild SI, Wicki A, Cathomas G, Mertz KD: Single-Center Experience
20
23.
24. 25. 26.
27.
28.
29.
30.
31.
with a Targeted Next Generation Sequencing Assay for Assessment of Relevant Somatic Alterations in Solid Tumors. Neoplasia 2017, 19:196-206. Guibert N, Hu Y, Feeney N, Kuang Y, Plagnol V, Jones G, Howarth K, Beeler JF, Paweletz CP, Oxnard GR: Amplicon-based next-generation sequencing of plasma cell-free DNA for detection of driver and resistance mutations in advanced non-small cell lung cancer. Ann Oncol 2018. Chang F, Li MM: Clinical application of amplicon-based next-generation sequencing in cancer. Cancer Genet 2013, 206:413-419. Fleischhacker M, Schmidt B: Circulating nucleic acids (CNAs) and cancer--a survey. Biochim Biophys Acta 2007, 1775:181-232. Johansson G, Andersson D, Filges S, Li J, Muth A, Godfrey TE, Stahlberg A: Considerations and quality controls when analyzing cell-free tumor DNA. Biomol Detect Quantif 2019, 17:100078. Lapin M, Oltedal S, Tjensvoll K, Buhl T, Smaaland R, Garresori H, Javle M, Glenjen NI, Abelseth BK, Gilje B, Nordgard O: Fragment size and level of cellfree DNA provide prognostic information in patients with advanced pancreatic cancer. J Transl Med 2018, 16:300. Newman AM, Bratman SV, To J, Wynne JF, Eclov NC, Modlin LA, Liu CL, Neal JW, Wakelee HA, Merritt RE, Shrager JB, Loo BW, Jr., Alizadeh AA, Diehn M: An ultrasensitive method for quantitating circulating tumor DNA with broad patient coverage. Nat Med 2014, 20:548-554. Zhang Y, Xu Y, Zhong W, Zhao J, Chen M, Zhang L, Li L, Wang M: Total DNA input is a crucial determinant of the sensitivity of plasma cell-free DNA EGFR mutation detection using droplet digital PCR. Oncotarget 2017, 8:5861-5873. Ascierto PA, Bifulco C, Palmieri G, Peters S, Sidiropoulos N: Preanalytic Variables and Tissue Stewardship for Reliable Next-Generation Sequencing (NGS) Clinical Analysis. J Mol Diagn 2019, 21:756-767. Wood HM, Belvedere O, Conway C, Daly C, Chalkley R, Bickerdike M, McKinley C, Egan P, Ross L, Hayward B, Morgan J, Davidson L, MacLennan K, Ong TK, Papagiannopoulos K, Cook I, Adams DJ, Taylor GR, Rabbitts P: Using next-generation sequencing for high resolution multiplex analysis of copy number variation from nanogram quantities of DNA from formalin-fixed paraffin-embedded specimens. Nucleic Acids Res 2010, 38:e151.
21 FIGURE LEGENDS
Figure 1. Modeling theoretical limits of variant detection based on depth of coverage using binomial sampling statistics. The confidence of variant prediction was determined for variants with allelic fractions equal to 2.5%, 5.0%, 10.0%, and 15% for depths of coverage with unique reads ranging from 1 to 500 reads (A). Similarly, the depth coverage was determined with unique reads that would be required to achieve 95% confidence (dashed line) for variants ranging in allelic fraction from 0.1 to 50.0% (B). Power and sample sizes were calculated for a binomial proportion using complimentary log parameterization given a sequencing error rate of 0.01%.
Figure 2. The impact of library complexity on assay sensitivity. Nextgeneration sequencing (NGS) tests are designed to survey the genetic makeup of tumor specimens that are typically comprised of a mixture of normal cells and neoplastic cells. The neoplastic component may be further subdivided into genetically diverse subclones. In this schematic, 10 unique input molecules were isolated and ligated to unique molecular identifiers. When a standard input of five molecules is used in library preparation, the theoretical limit of sensitivity is 1 in 5 and the variant is detected at a fraction similar to its true proportion in the original specimen. When a library is constructed from a minimum input of two molecules, the theoretical limit of sensitivity is reduced (1/2 rather than 1/5), and the variant is missed due to sampling error. In absence of UMIs, the total coverage of the
22 standard input and the minimum input libraries are identical and the vast difference in the sensitivity of variant detection would not be immediately apparent.
Figure 3. Relationship between total and unique read coverage. Total and unique read coverage of the designated analytical space was calculated for 41 NGS libraries generated from cell line DNA inputs of 50ng (green, n=10), 10ng (yellow, n=11), 5ng (red, n=10), or 2ng (blue, n=10). At the high total depths of coverage achieved in this experiment (average: 23,685±11,210), there is no significant correlation between total and unique read coverage (A). Data from NGS libraries prepared from 50ng were randomly down sampled (ds) to ~1,750x and ~10,000x total read coverage, and the percent of reads classified as unique was calculated for down sampled and non–down sampled libraries (B). As total coverage increases, the percent of reads classified as unique decreases (average of 26.7%±2.6%, 20.8%±4.6%, and 10.1%±2.58% for ds-1,750, ds10,000x, and original, respectively). The ends of boxes indicate the upper and lower quartiles. The median is indicated by the middle line. Whiskers extend to 1.5 times the interquartile range.
Figure 4. Average unique coverage at 10,000x total read coverage. Forty-one next-generation sequencing (NGS) libraries prepared from cell line material were sequenced and randomly down-sampled to approximately 10,000x total read coverage of the designated analytical space (average: 10,007±11.4 reads) to
23 facilitate intra-sample comparison. Despite some noise, the average unique depth of coverage tended to decrease with reduced DNA input (average of 776.2±569.2x, 971.0±725.4x, 1,263.6±784.4x, and 281.7±458.6x coverage for 2ng, 5ng, 10ng, and 50ng, respectively). The ends of boxes indicate the upper and lower quartiles, and the median is indicated by the middle line; whiskers extend to 1.5 times the interquartile range.
Figure 5. The impact of reduced DNA input on variant identification using standardized specimens. DNA from an engineered cell line with a heterozygous KRAS p.G13D missense variant was mixed in proportions of 1:1, 1:8, 1:16, and 1:32 to produce variants with approximate variant of a specified allelic fractions (VAFs) of 25%, 5%, 3%, and 1.5%, respectively (marked with black bars on the plot). Triplicate libraries were prepared and sequenced using either 50ng or 2ng of the mixed DNA as input, and variants were predicted. Variance between the triplicate estimates was greater for 2ng libraries (blue) compared to 50ng, standard input libraries (green). Coefficients of variation for VAF estimates were calculated as follows: 74.45% versus 8.50% for 2ng and 50ng of 1:1 mix, respectively; 170.73% versus 40.44% for 2ng and 50ng of 1:8 mix, respectively; 148.38% versus 65.15% for 2ng and 50ng of 1:16 mix, respectively; 20.94% versus 72.05% for 2ng and 50ng of 1:32 mix, respectively. Note that though the variance is actually lower for VAFs estimated from the 2ng libraries of the 1:32 mix compared to the 50ng libraries, the VAFs estimated from the 2ng libraries fell to nearly undetectable levels.
24
Figure 6. The impact of reduced DNA input on variant identification using formalin-fixed, paraffin embedded tumor specimens. Three de-identified, formalin-fixed, paraffin-embedded tumor specimens were selected for testing based on the presence of known, clinically significant variants: C1 with chr2:209113113G>T (IDH1 p.R132S, VAF=9.0%), C2 with chr12:25398284C>T (KRAS p.G12D, VAF=15.0%), and C3 with chr3:178936082G>A (PIK3CA p.E542K, VAF=3.0%). Variant of a specified allelic fractions (VAFs) estimated during previous clinical testing are marked with black bars. Triplicate libraries were prepared and sequenced using either 50ng or 2ng of DNA as input, and variants were predicted. Variance between the triplicate estimates was greater for 2ng libraries (blue) compared to 50ng, standard input libraries (green). Coefficients of variation were calculated as follows: 79.41% versus 10.48% for 2ng and 50ng libraries of C1, respectively; 140.96% versus 8.39% for 2ng and 50ng libraries of C2, respectively; 153.04% and 66.69% for 2ng and 50ng libraries of C3, respectively.
A
1.00 0.95
Confidence
0.75
0.50 Variant Allelic Fraction 2.5% 5.0% 10.0% 15.0%
0.25
0.00 0
100
200
300
400
500
Depth of Read Coverage 1000
Depth of Coverage
B
300
100
30 50%
40%
30%
20%
Variant Allelic Fraction
10%
0%
DNA isolation
Barcode Ligation
PCR amplification
Mixture of tumor and normal tissue
Minimum input
Standard input
Tissue isolation
Read alignment
Variant prediction
A (40%) T (60%)
Sensitivity = 1/5
T (100%)
Sensitivity = 1/2
Average Unique Read Coverage
A
DNA Input 2ng 5ng 10ng 50ng
4000
3000
r2 = 0.0537
2000
1000
0
10000
20000
30000
40000
Average Total Read Coverage
Percent Unique
B
0.30
0.25
0.20
0.15
0.10
ds-1,750x
ds-10,000x
Total Read Coverage
Original
Average Unique Read Coverage
2000
1000
2
5
10
DNA Input (ng)
50
VAF of KRAS p.G13D Variant
DNA Input 2 ng 50 ng
10.0%
1.0%
0.1%
1:1
1:8
1:16
Dilution Ratio
1:32
VAF of Variant
DNA Input 2 ng 50 ng
10.0%
1.0%
0.1%
C1
C2
Specimen
C3