Single molecule molecular inversion probes for targeted, high-accuracy detection of low-frequency variation

Method performance

We first sought to assess performance of the smMIP assay with respect to sensitivity and positive predictive value for clonal variation and uniformity of target enrichment. Because of the heterogeneous nature of the specimens, the number of sequencing reads obtained per sample varied from 1 to 16 million reads (Supplemental Fig. S1). However, 77% of samples (41 of 53) were within a threefold range, and this distribution could likely be improved by automated pooling (Supplemental Note S2). Fr-reads were then aligned to the reference genome and processed using a custom analysis pipeline to yield smc-reads.

We first explored coverage of targeted regions, finding that mean smc-read coverage of the targeted coding bases was 3538× across the HapMap samples and 1051× across the clinical specimens (Fig. 2; Supplemental Note S3). On average, smc-reads were comprised of 2.03 fr-reads in cell line samples, 1.21 fr-reads in fresh clinical samples, and 2.79 read-pairs in FFPE clinical samples (Supplemental Fig. S2). However, some samples exhibited a much higher rate of fr-reads per smc-read, which was likely due to degraded DNA leading to low-complexity capture. Thus, given a certain sequencing depth, the number of fr-reads per smc-read is directly related to capture complexity and can serve as a useful quality control metric.

Open in a separate window

Figure 2.

smMIP capture performance and detection of low-frequency variation. (A) Distributions of minimum coverage in a given percentile of total targeted coding positions, rank-ordered by smc-read coverage, for eight HapMap cell line (red) and 45 clinical cancer (blue and green) samples (box plot center line: median; top and bottom edges: quartiles; whiskers: farthest data point within 150% of interquartile range; dots: outliers). Zeroth-percentile indicates maximum coverage. (B) Distributions of fraction of coding positions above a given smc-read coverage cutoff. (C) Observed versus expected variant frequency in smc-read base-calls from mixtures of HapMap genomic DNA samples at known ratios for positions with at least 100× coverage (R = 0.94). Ideal performance is shown as gray line (y = x).

We then used smc-reads to call clonal genotypes using established tools (McKenna et al. 2010) and, for the HapMap samples, compared our calls to 1000 Genomes (“1KG”) pilot project genotypes (The 1000 Genomes Project Consortium 2010). For NA12892, we detected 24 of 25 1KG variant sites; the remaining position was not adequately covered in our data. After discarding three positions that were systematically misgenotyped by our assay (caused by capture of paralogous sequence and subsequent misalignment), we detected two additional variant positions; these calls were supported by manual inspection of more recent 1KG data. For NA19239, we detected 41 of 44 1KG variant sites; the remaining three positions were not adequately covered in our data. Two additional sites had variant genotypes and were again supported by newer 1-KG data. Therefore, based on this limited comparison, our assay is highly accurate at adequately covered positions. Considering all targeted sites, we estimate the sensitivity of our assay for clonal homozygous or heterozygous variation to be 93%–96%, and the positive predictive value to be near 100%.

Subclonal variant detection

To assess whether the smMIP assay was capable of sensitively detecting and accurately quantifying variants present at subclonal frequencies, we applied it to six synthetic mixtures of genomic DNA from the two HapMap cell lines combined in a twofold serial dilution from 1:8 to 1:256 (resulting in low-abundance genome alternate allele frequencies of ∼11% to ∼0.2%). We adopted a custom variant calling strategy to detect subclonal variation and then compared the expected variant frequency to that observed in smc-reads (Fig. 2C). In general, we observed close agreement between the expected and observed frequency for positions with at least 100× smc-read coverage (R = 0.94), with the deviation from expected frequency largely explained by sampling statistics (Supplemental Fig. S3).

We next sought to quantify the absolute error rate of the smMIP assay and to assess potential sources of error (Supplemental Fig. S4). In the samples derived from HapMap cell lines, smc-read base-calls were ∼13-fold more accurate than fr-read base-calls, with a substitution rate of 8.4 × 10⁻⁶ per base compared to 1.1 × 10⁻⁴ per base for fr-read calls (Table 2). In the hematologic clinical samples, which were not FFPE, the difference was ∼12-fold (9.5 × 10⁻⁶ per base compared to 1.1 × 10⁻⁴ per base), and for FFPE clinical samples, the difference was approximately fivefold (2.7 × 10⁻⁵ per base compared to 1.3 × 10⁻⁴ per base). We then explored substitution rates as a function of the expected nucleotide incorporated into the MIP during the gap-fill versus the observed nucleotide in the fr-read or smc-read base-call (Fig. 3A). We observed substantial variation with respect to different pairs of expected/observed nucleotides, as well as an effect of dinucleotide context (Supplemental Fig. S5). In particular, we observed that two single-nucleotide substitutions, cytosine to adenine and guanine to adenine, exhibited elevated rates compared to all other substitutions across all sample types (Fig. 3B).

Table 2.

Substitution error rates

Open in a separate window

Open in a separate window

Figure 3.

Substitution error rates as a function of expected and observed nucleotide during gap-fill. (A) Schematic illustrating mononucleotide and dinucleotide substitution dependencies being considered. All rates are shown for a given expected gap-fill mono- or dinucleotide, which is the complementary nucleotide(s) to the nucleotide(s) present in the target genomic DNA, considering only ≥Q41 fr-read base-calls and Q60 smc-read base-calls at putative homozygous positions based on GATK calls. (B) Distributions of substitution error rates for eight HapMap cell line and 45 clinical cancer samples, comparing fr-reads and smc-reads, and all substitutions other than C>A or G>A (W>N + N>B, left) to only C>A (middle), or G>A (right). (C) Distributions of substitution error rates comparing fr-reads and smc-reads, and all G>A substitutions occurring in the non-CG dinucleotide context (DG>DA + GN>AN, left) to G>A substitutions occurring only in the CG dinucleotide context (CG>CA, right).

One potential explanation for these patterns is the occurrence of pro-mutagenic chemical processes in individual progenitor molecules, namely, the oxidatively damaged base 8-oxo-guanine (manifested as a cytosine to adenine substitution) (Shibutani et al. 1991) and spontaneous deamination of cytosine and 5-methyl-cytosine (manifested as a guanine to adenine substitution) (Fig. 3B,C). After smc-read formation, elevated substitution rates potentially attributable to 8-oxo-guanine damage are evident across all sample types, but are most markedly elevated in the FFPE clinical samples, suggesting that FFPE treatment may increase 8-oxo-G formation rates. To further investigate the possibility of cytosine deamination, we asked whether the G-to-A substitutions were more frequent in the CpG context, which could be consistent with deamination of 5-methyl-cytosine (Fig. 3C). Indeed, the G-to-A substitution was more frequent in the CpG context across all sample types compared with all other G-to-A substitutions, but the CpG dinucleotide is rare and is not sufficient to fully explain the elevated G-to-A substitution rate (Supplemental Table S3). Finally, we asked whether the corresponding C-to-T mutation rate was also elevated in the CpG context, which would be expected if 5-methyl-C deamination were occurring spontaneously in a replicating population of cells. We observed modest elevation of C-to-T in the CpG context (Supplemental Fig. S5D). Together, these patterns suggest a combined role of in vivo 5-methyl-cytosine deamination as well as ex vivo deamination of 5-methyl-C-to-T and C-to-U. This is consistent with the properties of the polymerase used during the gap-fill step, which, as a proofreading polymerase, is expected to stall at Uracil residues (Lasken et al. 1996), but also bypasses these lesions at some rate (Greagg et al. 1999).

When we removed the contribution of these potential sources of false-positive substitution calls, smc-read base-calls were even more accurate, with substitution rates of 3.5 × 10⁻⁶ per base and 5.1 × 10⁻⁶ per base for the HapMap and clinical samples, respectively, while the substitution rates of the fr-read calls did not decrease substantially (Table 2). These patterns are consistent with gap-fill misincorporations due to DNA damage and actual subclonal heterogeneity constituting the major source of substitutions in smc-read base-calls, and with polymerase errors after the initial gap-fill event constituting a major source of substitutions in fr-reads. Smc-reads are, therefore, at least fivefold and as much as 30-fold more accurate than the most confident base-calls in the fr-reads, with substitution rates as low as 3.5 × 10⁻⁶ per base when ignoring only two of 12 possible substitutions, or 8.4 × 10⁻⁶ per base when considering all possible single nucleotide substitutions.

We also sought to determine the sensitivity and false discovery rate (FDR) of smMIP for low-frequency variation more generally. We used the single-nucleotide substitution error rates calculated as described above to compute P-values for each subclonal variant and adjusted these P-values to account for multiple testing. We then explored sensitivity and FDR in the synthetically mixed HapMap samples. This analysis was limited by the small number of sites that were divergent between the two individuals in the targeted coding regions (n = 18). However, we found that smc-reads were generally more sensitive at a given FDR over fr-reads for variation present at frequencies down to 0.2% (Fig. 4). For example, at an FDR of 20%, sensitivity ranged from 94% for variants at 6%–11% frequency to 44% for variants at 0.2%–0.4% frequency, and variant calls from smc-reads were 3%–22% more sensitive than the fr-reads at the same FDR cutoff across the various mixtures.

Open in a separate window

Figure 4.

Sensitivity and false discovery rates for subclonal variation in synthetic mixtures. Sensitivity versus false discovery rate for low-frequency variants (0.1%–40%) in synthetically mixed HapMap samples for variant calls from fr-reads (red) and smc-reads (blue), for coding positions that were adequately genotyped in both unmixed HapMap samples and for which there was no substantial (binomial adjusted P < 10⁻¹⁰) subclonality in the predominant HapMap sample. Expected subclonal variant frequencies are listed at the top of each panel. Area beneath the curve is shown as an inset in each panel. Candidate subclonal variants occurring in coding sequence and at a frequency of at least 0.1% were prioritized using multiple testing-adjusted binomial P-values that were calculated from substitution error rates.

Detection of somatic variation

Because of the simple experimental workflow, the smMIP assay could also be useful in clinical and high-volume research settings as a replacement for single-gene testing for clinically informative mutations. To assess whether the smMIP method can potentially replace such tests, we performed a blinded comparison of smMIP results to the results of clinical single-gene tests. In particular, a subset of our samples had been previously genotyped for individual actionable substitution and indel mutations in BRAF, EGFR, FLT3, JAK2, KIT, KRAS, NRAS, and PDGFRA. Considering these sites, we detected 25 of 27 (93%) previously identified mutations (Table 3; Supplemental Tables S4, S5). We missed two large (67- and 104-bp) insertions in FLT3, although these could in principle be detected using a more sensitive analysis strategy and/or more densely tiled probes in this region. We further detected two mutations in these sites in two lung cancer samples that had not been previously genotyped at that site (KRAS p.G12C in sample 8 and p.G12V in sample 37); these calls were subsequently confirmed using a melt curve-based assay (data not shown).

Table 3.

Concordance with single mutation tests

Open in a separate window

To explore other somatic mutations, and because we did not have access to matched normal tissue, we filtered variant sites identified in the clinical cancer specimens against germline variant sites observed in ∼5400 exomes by the Exome Sequencing Project (http://evs.gs.washington.edu/EVS/). We then required at least 30× coverage to remove poorly ascertained sites. Across the 45 clinical samples, this filtering process yielded 134 putative somatic events, of which 74 were found in the Catalogue of Somatic Mutations in Cancer (Supplemental Fig. S6A; Forbes et al. 2011). As expected, several genes were recurrently mutated in specific tumor types (Supplemental Fig. S6B), such as 11 of 16 colon cancer samples harboring at least one APC mutation.

To investigate the possibility of tumor subclones in these specimens, we examined the extent to which putative somatic mutations were observed at similar frequencies within individual samples (Supplemental Fig. S7A). While some clustering of alternate allele frequencies is apparent, we observed substantial variation of alternate allele frequencies within individual samples, which may reflect the presence of multiple, genetically distinct subclones. Copy-number gain is another potential source of alternate allele frequencies substantially different from 0.5 or 1. To explore the possible contribution of copy-number change to allele frequency variation, we also examined alternate allele frequencies for putatively germline variant sites for all clinical samples and the two pure HapMap cell line samples (Supplemental Fig. S7B). Alternate allele frequencies for germline variants appeared substantially more variable in a subset of clinical samples compared with the cell line samples, which is consistent with a contribution of copy-number gain to the observed variation in allele frequencies. We cannot exclude the possibility that other factors, including systematic allelic bias, are also contributing to this phenomenon. However, the general precision observed in the HapMap cell line mixtures decreases the likelihood of widespread systematic allelic bias.

Subclonal somatic variation at clinically informative sites

We next explored the potential for the smMIP assay to detect very low-frequency substitutions at clinically informative sites (BRAF p.V600, EGFR p.L858, JAK2 p.V617, KRAS p.G12/p.G13, NRAS p.Q61, and PDGFRA p.D842) in the clinical samples. These sites were not chosen because of any particular sequence or error rate characteristics, but rather to enable confirmation of subclonal variant calls with established assays and to facilitate interpreting the possible biological significance of putative subclonal variation. Restricting to maximum quality smc-read base-calls but without any filtering for variant frequency or confidence, we detected 17 candidate subclonal variants at these sites across the 45 clinical samples. However, these candidate variants exhibited a highly non-uniform distribution of error probabilities (calculated as described for the HapMap subclonal variants) with seven candidate variants having P < 10⁻⁷ and the remaining 10 having P > 10⁻². We therefore focused further analysis on the seven high-confidence candidates (Table 4).

Table 4.

Low-frequency variation at clinically relevant sites in tumor samples

Open in a separate window

These seven subclonal variants consisted of low-frequency JAK2 p.V617F mutations (n = 4) in two lung cancers, a melanoma, and a colon cancer; a BRAF p.V600E mutation in a melanoma; a KRAS p.G12V mutation in one of the lung tumors that also harbored a low-frequency JAK2 p.V617F mutation; and an NRAS p.Q61R mutation in a colon cancer. To exclude the possibility of artifactual low-frequency variant detection due to sample cross-contamination or index cross-talk, we subjected independent DNA aliquots from the four specimens with low-frequency JAK2 mutations to confirmatory clinical testing using an allele-specific PCR assay; all four mutations were confirmed. Furthermore, JAK2 mutations have been reported at low frequency in a previous study in non-small cell lung cancer samples (Lipson et al. 2012). Additionally, index sequence cross-talk, which would be predicted to give rise to mixed read groups, is not a likely explanation for these low-frequency calls; we required very high-quality smc-read base-calls, and mixed read groups are not a general phenomenon in our data (Supplemental Fig. S8). We note that, prior to our study, the melanoma sample was genotyped clinically for BRAF p.V600, and this mutation was not detected (as expected given that assay's limited sensitivity); additionally, polyclonality in melanomas with respect to BRAF p.V600 mutation status has been previously observed (Lin et al. 2009, 2011; Yancovitz et al. 2012).

We therefore expect, based on experimental and biological evidence, that these seven candidate variants are bona fide. However, based on the FDR analysis performed using the HapMap samples, the P-value cutoff of 10⁻⁷ is expected to yield a nontrivial FDR of ∼40% for variants near 0.1% frequency (Supplemental Fig. S9), so it remains possible that one or more of these variants is artifactual.

Capture and library construction

Genomic DNA for HapMap cell line samples (NA12892 and NA19239) was purchased (Coriell). Clinical specimens consisted of DNA prepared from formalin-fixed paraffin-embedded (FFPE) tissue, peripheral blood, or bone-marrow aspirates from patients with sporadic colorectal cancer (n = 18), melanoma (7), non-small cell lung cancer (11), bladder cancer (1), ovarian cancer (1), gastrointestinal stromal tumor (4), acute myeloid leukemia (2), and myeloproliferative disorders (3). De-identified residual clinical specimens were obtained from the University of Washington molecular diagnostics laboratory in accordance with the declaration of Helsinki and ethics guidelines of the local institutional review board. Specifically, the de-identified samples were determined by the institutional review board to be nonhuman subjects research. Hematoxylin and Eosin-stained slides were used as a guide to manually dissect areas of tumor tissue from unstained slide sections for FFPE tissue samples. Genomic DNA was prepared with the Gentra Puregene DNA Isolation Kit (Qiagen). A 3-h to overnight proteinase K digestion step was included for FFPE samples. Genomic DNA from the HapMap sample NA12892 was serially diluted twofold and added to 500 ng of genomic DNA from HapMap samples NA19239 at six relative ratios ranging from 1:8 to 1:256. Captures of the six cell line mixtures, two pure cell line samples, and 47 clinical samples were then performed using ∼500 ng of each genomic DNA.

Captures were performed as previously described (Shen et al. 2011) with some modifications. A total of 500 ng of genomic DNA, 330 femtomoles of probe mixture (ignoring the contribution of the spiked-in poor performers), and 1 μL of 10× Ampligase DNA ligase buffer (Epicentre) were added to molecular biology-grade water for a total of 10 μL. For the probe hybridization phase, these mixtures were incubated in a thermocycler (Bio-Rad) with a heated lid at 98°C for 3 min, 85°C for 30 min, 60°C for 60 min, and 56°C for 120 min. For the gap-fill and ligation phase, we added 300 picomoles each dNTPs (NEB), 7.5 micromoles betaine (Sigma), 20 nanomoles NAD⁺ (NEB), 1 μL of 10× Ampligase buffer, 5 units of Ampligase DNA ligase (Epicentre), 3.2 units of Phusion DNA polymerase (NEB), and molecular biology grade water to 10 mL for a total reaction volume of 20 μL. The gap-fill and ligation phase was carried out at 56°C for 60 min and 72°C for 20 min.

Following the gap-fill and ligation phase, the reactions were cooled to 37°C, and to each reaction we added 20 units of Exonuclease I (NEB) and 100 units of Exonuclease III (NEB) to degrade uncircularized probe and genomic DNA. The digestion was incubated at 37°C for 45 min, heated to 80°C, and incubated for 20 min to inactivate the exonucleases.

After exonuclease treatment and heat-inactivation, the samples were cooled on ice and, optionally, stored at −20°C. For each capture reaction, two PCR reactions were prepared, each with Phusion HF buffer to 1× (Fermentas), forward primer and indexed reverse PCR primers to 500 nM, SYBR green (Invitrogen) to 0.5×, dNTPs to 200 μm each (NEB), 2 units of Phusion Hot-Start II polymerase, 10 μL of capture reaction, and nuclease-free water to 50 μL. PCR cycling conditions were an initial denaturation step for 2 min at 95°C, followed by 26 cycles of: 15 sec at 98°C, 15 sec at 65°C, and 45 sec at 72°C. A subset of samples was run on a real-time PCR instrument (Bio-Rad MiniOpticon) to estimate the required number of cycles; the remaining samples were run without real-time monitoring (Bio-Rad DNA Engine Tetrad 2).

Library purification and pooling

PCR products were purified individually using Ampure XP beads (Agencourt) at 1.8× according to the manufacturer's instructions. Purified PCR products were then pooled naively (i.e., equal volumes) for initial quality control or based on MiSeq sequence data considering the number of reads mapping on target per sample. To remove a low molecular-weight artifact, the PCR product pool was split across 8 wells of a 10-well pre-cast 6% polyacrylamide TBE gel (Invitrogen), run at 140 V for 50 min, and stained with 5 μL of SYBR Gold (Invitrogen). The capture product band at ∼280 bp was excised, crushed, soaked in 800 μL of Tris-EDTA (pH 8.0), and recovered from the supernatant using 100 μL of Ampure XP beads and 700 μL of home-made Ampure buffer (20% PEG 8000, 2.5M NaCl) according to the manufacturer's instructions. Pool concentration was assessed using Qubit (Invitrogen).

As an alternative strategy to reduce the time, labor, and cost required for library construction, purified PCR products were subjected to an Ampure-based size enrichment and normalization step. Twenty microliters of each purified PCR product was purified using 16 μL of a mixture of one part Ampure XP bead solution and four parts homemade Ampure buffer, and eluted in 20 μL of Buffer EB (Qiagen). Two microliters of each sample was then run on a diagnostic precast 6% polyacrylamide TBE gel as described above to assess relative concentrations of capture product. The gel image was analyzed for band intensity (ImageJ), and the purified PCR products were pooled according to relative band intensity and the pool was quantified via Qubit.

Sequencing and primary analysis

Samples were sequenced using the HiSeq 2000 and MiSeq (Illumina) platforms according to the manufacturer's instructions using custom sequencing primers (Supplemental Table S9). On the HiSeq platform, we collected two 101-nt reads to determine the sequence of the gap-fill and the molecular tag and one 8-nt read to determine the sequence of the sample index. On the MiSeq platform, we collected two 152-nt reads and one 8-nt read. Initial quality control and capture performance of an equivolume, non-size selected pool of all 55 samples was assessed using one run of the MiSeq platform. Purified PCR products were then repooled according to MiSeq data, size-selected using a PAGE gel as described above, and subjected to 2.75 lanes of HiSeq 2000 sequencing (two lanes with no other samples mixed in and one lane with 25% by moles of an unrelated library). For the establishment of a rapid workflow, eight samples were processed as described above and subjected to one run of the MiSeq platform.

Read-pairs were assigned to samples requiring an exact match to the expected 8-nt sample index sequence (Supplemental Table S9) and the first 12 nt of the reverse read (corresponding to the molecular tag sequence) were stripped out and placed in the header. Read-pairs with molecular tags with homopolymers longer than 4 nt were discarded. Overlapping regions of read-pairs were then reconciled to form single “fr-reads” using a custom Smith-Waterman-based strategy. For positions where the read-pairs did not overlap, quality scores from the individual reads were retained. For positions where the read-pairs did overlap, quality scores for the resulting consensus calls were estimated as below for smc-reads. Only successfully overlapped fr-reads were retained for downstream analysis; read-pairs that failed to merge were discarded, although subsequent implementations aimed at greater sensitivity toward large insertions such as those found in FLT3 will retain and analyze these reads, and smMIP does not explicitly require collecting overlapping read-pairs. Fr-reads were aligned to the human reference genome (hg19) using the bwasw alignment mode of the aligner bwa (Li and Durbin 2010) (v0.5.9) with non-default parameters “-r 1”. Based on expected alignment positions according to the probe design, fr-reads were then assigned to individual probes, allowing 1 nt of tolerance in each direction for the beginning of the read, which primarily accommodates insertion and deletion mutations during probe synthesis. Then, for each sample and each probe, fr-reads were grouped by molecular tag sequence to form tag-defined read groups (TDRGs).

Single molecule consensus read calling

Alignments to the reference genome were used to call a consensus sequence (i.e., a single molecule consensus read or “smc-read”) for each TDRG. Positions expected to be derived from probe targeting arms (and therefore synthetic DNA) were excluded from consideration at this step. We adopted a likelihood ratio framework to incorporate both the abundance and associated quality-scores of fr-read base-calls supporting each possible nucleotide at each position. Briefly, we calculated the log-likelihood L_x of consensus nucleotide x as the difference of the log-likelihood of a model in which a given nucleotide x was underlying none of the calls and the log-likelihood of a model in which x was underlying all of the calls. This can be represented as in the equation below, where a given observation is associated with a nucleotide call o and a “phred-like” quality score Q_o.

This value (L_x) is computed for each possible consensus call (A, C, G, T, N, and deletion) at each position in the alignment based on the set of base-calls and associated quality scores in the TDRG at that position in the alignment. The consensus call is then determined as the call with the minimum (i.e., most negative) log-likelihood value, as this indicates the consensus nucleotide where the model assuming that nucleotide did not underlie any of the fr-readcalls was the least likely relative to the model assuming that nucleotide underlay all of the fr-read calls. The final phred-like quality score is calculated as the integer casting of −10*L_x. For example, in the event that a model against a given consensus nucleotide is 10⁻³ times as likely as a model for a given consensus nucleotide, the associated quality score was calculated to be 30. These scores were capped at 60 as we did not observe substitution rates substantially below 1 × 10⁻⁶ in practice. We note that a smc-read base-call that is derived from a single Q60 fr-read base-call (which may have been derived from two Q40 calls at an overlapping position, for example) will be assigned an estimated quality score of 59 because of the integer casting step. In practice, therefore, only smc-read base-calls that were supported by at least two independent fr-read base-calls can attain quality 60.

For reference positions where at least one read in the alignment indicated a deletion and at least one a match, or where at least one read contained an insertion relative to reference and at least one read lacked that insertion, we applied the same framework, assigning deletion calls a phred-like quality score of 40. We note that this framework makes simplifying assumptions: that a given nucleotide, in the event of a sequencing error, is equally likely to give rise to any of the other three substitution nucleotides, and that a single nucleotide truly underlies all calls across all reads in the TDRG. Furthermore, for interpretability, we used three as the denominator to distribute the probability when the observed nucleotide was not the candidate consensus nucleotide, though we computed this value for a total of six possible consensus calls and not four.

This strategy was also used to determine quality scores for consensus calls at overlapping positions in fr-reads, which represents the simpler case of two and only two calls.

J.S. is a member of the scientific advisory board or serves as a consultant for Ariosa Diagnostics, Stratos Genomics, Good Start Genetics, and Adaptive Biotechnologies. A patent application has been filed for aspects of the methods disclosed here (“Subassembly of Short Sequencing Reads”; 12/559124).

We thank C. Lee for assistance with sequencing, and members of the Shendure Lab for helpful discussions. We thank the NHLBI GO Exome Sequencing Project and its ongoing studies that produced and provided exome variant calls for comparison: the Lung GO Sequencing Project (HL-102923), the WHI Sequencing Project (HL-102924), the Broad GO Sequencing Project (HL-102925), the Seattle GO Sequencing Project (HL-102926), and the Heart GO Sequencing Project (HL-103010). We thank C. Smith, K. Koehler, M. Murillo, and H-S. Yi for help genotyping clinical samples. Our work was supported by a gift from the Washington Research Foundation; grant CA160080 from the National Cancer Institute (J.S.); grant AG039173 from the National Institute on Aging (J.B.H.); a fellowship from Achievement Rewards for College Scientists (J.B.H.); and the Department of Laboratory Medicine, University of Washington Medical Center.

Author contributions: J.B.H., C.C.P, S.J.S., and J.S. conceived and designed the study. J.B.H. and B.J.O. designed the smMIPs and developed protocols. C.C.P. and S.J.S. obtained anonymized clinical samples, coordinated or oversaw single mutation genotyping, and aided with interpretation of results. J.B.H. performed all other experiments and all data analysis. J.B.H. and J.S. wrote and revised the manuscript with substantial input and revisions from all other authors. J.S. supervised all aspects of the study.