Detection of Rare Genomic Variants from Pooled Sequencing Using SPLINTER

As DNA sequencing technology has markedly advanced in recent years2, it has become increasingly evident that the amount of genetic variation between any two individuals is greater than previously thought3. In contrast, array-based genotyping has failed to identify a significant contribution of common sequence variants to the phenotypic variability of common disease4,5. Taken together, these observations have led to the evolution of the Common Disease / Rare Variant hypothesis suggesting that the majority of the "missing heritability" in common and complex phenotypes is instead due to an individual's personal profile of rare or private DNA variants6-8. However, characterizing how rare variation impacts complex phenotypes requires the analysis of many affected individuals at many genomic loci, and is ideally compared to a similar survey in an unaffected cohort. Despite the sequencing power offered by today's platforms, a population-based survey of many genomic loci and the subsequent computational analysis required remains prohibitive for many investigators. To address this need, we have developed a pooled sequencing approach1,9 and a novel software package1 for highly accurate rare variant detection from the resulting data. The ability to pool genomes from entire populations of affected individuals and survey the degree of genetic variation at multiple targeted regions in a single sequencing library provides excellent cost and time savings to traditional single-sample sequencing methodology. With a mean sequencing coverage per allele of 25-fold, our custom algorithm, SPLINTER, uses an internal variant calling control strategy to call insertions, deletions and substitutions up to four base pairs in length with high sensitivity and specificity from pools of up to 1 mutant allele in 500 individuals. Here we describe the method for preparing the pooled sequencing library followed by step-by-step instructions on how to use the SPLINTER package for pooled sequencing analysis (http://www.ibridgenetwork.org/wustl/splinter). We show a comparison between pooled sequencing of 947 individuals, all of whom also underwent genome-wide array, at over 20kb of sequencing per person. Concordance between genotyping of tagged and novel variants called in the pooled sample were excellent. This method can be easily scaled up to any number of genomic loci and any number of individuals. By incorporating the internal positive and negative amplicon controls at ratios that mimic the population under study, the algorithm can be calibrated for optimal performance. This strategy can also be modified for use with hybridization capture or individual-specific barcodes and can be applied to the sequencing of naturally heterogeneous samples, such as tumor DNA.

. In contrast, array-based genotyping has failed to identify a significant contribution of common sequence variants to the phenotypic variability of common disease 4,5 . Taken together, these observations have led to the evolution of the Common Disease / Rare Variant hypothesis suggesting that the majority of the "missing heritability" in common and complex phenotypes is instead due to an individual's personal profile of rare or private DNA variants [6][7][8] . However, characterizing how rare variation impacts complex phenotypes requires the analysis of many affected individuals at many genomic loci, and is ideally compared to a similar survey in an unaffected cohort. Despite the sequencing power offered by today's platforms, a population-based survey of many genomic loci and the subsequent computational analysis required remains prohibitive for many investigators.
To address this need, we have developed a pooled sequencing approach 1,9 and a novel software package 1 for highly accurate rare variant detection from the resulting data. The ability to pool genomes from entire populations of affected individuals and survey the degree of genetic variation at multiple targeted regions in a single sequencing library provides excellent cost and time savings to traditional single-sample sequencing methodology. With a mean sequencing coverage per allele of 25-fold, our custom algorithm, SPLINTER, uses an internal variant calling control strategy to call insertions, deletions and substitutions up to four base pairs in length with high sensitivity and specificity from pools of up to 1 mutant allele in 500 individuals. Here we describe the method for preparing the pooled sequencing library followed by step-by-step instructions on how to use the SPLINTER package for pooled sequencing analysis (http://www.ibridgenetwork.org/wustl/splinter). We show a comparison between pooled sequencing of 947 individuals, all of whom also underwent genome-wide array, at over 20kb of sequencing per person. Concordance between genotyping of tagged and novel variants called in the pooled sample were excellent. This method can be easily scaled up to any number of genomic loci and any number of individuals. By incorporating the internal positive and negative amplicon controls at ratios that mimic the population under study, the algorithm can be calibrated for optimal performance. This strategy can also be modified for use with hybridization capture or individual-specific barcodes and can be applied to the sequencing of naturally heterogeneous samples, such as tumor DNA.

Video Link
The If required, PCR results can typically be improved by: 1) lowering the annealing temperature for small amplicons; 2) raising the annealing temperature for large amplicons; 3. lengthening the extension time for any amplicon. 5. Preparation of SPLINTER controls: Every SPLINTER experiment requires the presence of a negative and positive control to obtain optimal accuracy. A negative control can consist of all homozygous base positions in any individual, bar-coded sample that has been previously sequenced (e.g. a HapMap sample). The positive control would then consist of a mixture of two or more such samples. For this report, the negative control is a 1,934 bp amplified region from the backbone of the M13mp18 ssDNA vector. The PCR product was Sanger sequenced prior to its use in order to confirm that no sequence variation exists from the source material or the PCR amplification. The positive control consists of a panel of pGEM-T Easy vectors with a 72 bp cloned insert engineered with specific insertions, deletions, substitutions ( Table 1).
We mix the vectors together against a wild type background at molar ratios such that the mutations are present at the frequency of a single allele in the pool (i.e. for a 100-allele pool, the frequency of a single allele will be 1%). We then PCR amplify the mixed control template using the M13 PUC primer sites in pGEM-T Easy, generating a final 355bp long PCR product.

Pooled PCR Library Preparation and Sequencing
1. PCR product pooling: Each PCR product should be cleaned of excess primers. We used Qiagen Qiaquick column purification or 96-well filter plates with vacuum manifold for large-scale cleanup. Following purification, each PCR product should be quantified using standard techniques. Combine every PCR product (including the controls) into a pool normalized by molecule number as pooling by concentration will result in overrepresentation of small amplicons over larger products. Concentrations are converted to the absolute number of DNA molecules per volume using the formula: (g / μL) x (1 mol x bp / 660 g) x (1 / # bp in amplicon) x (6 x 10 23 molecules / 1 mol) = molecules / μL. We then determine the volume from each reaction required to pool a normalized number of molecules per amplicon. This number is arbitrary, can be adjusted and really depends upon pipetting volumes large enough to maintain accuracy. We typically pool 1-2 x 10 10 molecules of each amplicon. 2. Ligation of PCR products: This step is necessary to achieve uniform sequencing coverage as sonication of small PCR amplicons will biased their representation toward their ends. To overcome this, we ligate the pooled PCR products into large concatemers (>= 10 Kb) prior to fragmentation. Pfu Ultra HF Polymerase generates blunt ends, leading to efficient ligation (a Taq-based polymerase will add a 3p "A" overhang that will not allow ligation without prior fill-in or blunting). This reaction can be scaled up 2-3 fold if necessary. The ligation reaction contains 10 U T4 polynucleotide kinase, 200 U T4 ligase, 15% w/v polyethylene, 1X T4 ligase buffer, glycol 8000 MW, up to 2 μg of pooled PCR products in a final volume of 50 μL. Reactions are incubated at 22 °C for 16 hours followed by 65 °C for 20 minutes and held at 4 °C thereafter. The success of this step can be checked by loading 50 ng of samples into a 1% agarose gel. Successful ligation will result in a high molecular-weight band present in the lane (see Figure 2, lane 3). 3. DNA fragmentation: At this point you should have large concatemers (>10kb) of PCR products. We have a random sonication strategy using a 24-sample Diagenode Bioruptor sonicator that can fragment these concatemers in 25 minutes (40 sec "on"/20 sec "off" per minute). Sonication is inhibited by the viscosity introduced by the PEG, so this can be overcome by diluting the sample 10:1 in Qiagen PB buffer. Results can be checked on a 2% agarose gel (see Figure 2, lanes 4 & 5). 4. The sample is ready to incorporate directly into the Illumina Genomic Library Sample Preparation protocol beginning with the "End Repair" step. The data reported here are from single-end reads on the Illumina Genome Analyzer IIx, but we have used the HiSeq 2000 and performed single or paired-end reads with comparable results. Given the scale of the library created, we have also used custom barcoded adapters in order to multiplex multiple pooled libraries to accommodate the bandwidth supplied by the HiSeq platform (data not shown). Follow the manufacturer's protocol and recommendations that come with the kit. In order to achieve optimal sensitivity and specificity for variant detection, target coverage of 25-fold or more per allele is recommended (Figure 3). This estimate is independent of pool size and type of variant to be detected. If necessary multiple lanes and runs can be combined to reach adequate coverage. The number of mismatches per read that are allowed compared to the reference sequence is a user-defined parameter. Reads that have an excess number of mismatches will be discarded. We recommend allowing 2 mismatches for 36 bp reads, 4 mismatches for 76 bp reads and 5 mismatches for 101 bp reads. Allowing more mismatches will increase the likelihood of allowing excess sequencing errors into the aligned data. As read lengths continue to become longer, this value can be further increased. 3. Tagging aligned files from the same flowcell: At this point the entire aligned read file should be given a unique identifier ("tag") in order to identify read files belonging to the same sequencing run (i.e. multiple lanes from the same flowcell can be aggregated and given a single tag). The tag is necessary because each machine run generates a unique error profile that can be characterized via the tag. A tag is an alphanumerical string of characters used to distinguish a set of reads (the underscore character "_" should not be used for parsing issues). Different tags should be used for aligned read files generated on different flowcells or machine runs. Tags can be added using the included RAPGAP_alignment_tagger.pl with the following command:

Sequencing Reads Alignment and Analysis
. The EMGENERATOR4 tool will generate 3 files named as the output file name parameter followed by _0, _1 or _2. These files correspond to a 0th, 1st and 2nd order error model respectively. For variant calling with SPLINTER, the 2nd order error model should always be used. 5. For visualizing the error rate profile of a run, the error_model_tabler_v4.pl can be used to generate a PDF error plot on the 0th order error model file (Figure 4): .

/error_model_tabler_v4.pl [Error model 0th order file] [output file name]
The plot file will reveal run-specific error trends and can be used to infer the maximum number of read bases to be used for the analysis, which is explained in the next section. The number of read bases to be used varies and should be evaluated according to each run. We generally recommend using the first 2/3rds of the read as they represent the highest quality data (the first 24 read bases of a 36bp long read, for example). Single read bases can be excluded from the analysis if found to be defective (separated by a comma or N e.g. 5,7,11 or N). The p-value cutoff dictates how stringent the variant calling analysis is going to be. We normally start the analysis by allowing a minimum cutoff of -1.301 (corresponding to a p-value ≤ 0.05 in log10 scale). The pool size option optimizes the algorithms "signal-to-noise" discrimination by eliminating potential variants with minor allele frequencies less than that of a single allele in the actual pool. For example in a pool of 50 individuals, the lowest observed variant can be expected at 0.01 frequency or 1 in 100 alleles. Thus, the pool size option should be set to the closest value that is greater than the actual number of alleles analyzed in the experiment (i.e. if 40 people are surveyed, we expect 80 alleles so the closest option would be a pool size of 100). Variants called at frequencies <0.01 will then be ignored as noise. This file returns all hits that are statistically significant across the sample, with a description of position of the variant, type of variant, p-value per DNA strand, frequency of the variant and total coverage per DNA strand ( where the list file is a list of positive control hits in the form of a tab-delimited file. The first field indicates the amplicon of interest, whereas the second field indicates the position in which the mutation is present. N indicates that the rest of the sequence does not contain any mutation. 3. Determining the optimal p-value thresholds using the positive control data: After normalization, the analysis of the positive control is indispensable for maximizing sensitivity and specificity of a particular sample analysis. This can be achieved by finding the optimal p-value cutoff using the information from the positive control. Most likely, the initial p-value of -1.301 will not be stringent enough, which if so, will result in the calling of false positives from the positive or negative control. Every SPLINTER analysis will show the actual p-value for each called variant (see columns 5 and 6 on Table 2), which could not be predicted a priori. However, the entire analysis can be repeated by using the least stringent p-value displayed on the initial output for the known true positive base positions. This will serve to retain all true positives while excluding most, if not all, false positives and they typically have much less significant p-values compared to true positives. To automate this process, the cutoff_tester.pl can be used. cutoff_tester.pl requires a SPLINTER output file and a list of positive control hits in the form of a tab-delimited file as the one used for normalization: ./cutoff_tester.pl [SPLINTER filtered file] [list file] The resulting output will be a list of cutoffs that progressively reach the optimal one (see Table 3). The last line represents the most optimal cutoff for the run and can therefore be used for data analysis. The optimal result is to achieve sensitivity and specificity of 1. In case this result is not reached, the SPLINTER analysis can be repeated by changing the number of incorporated read bases until the most optimum condition is achieved. 4. Final variant filtering: The final cutoff can be applied to the data using cutoff_cut.pl script, which will filter the SPLINTER output file from hits below the optimal cutoff, ./cutoff_cut.

pl [SPLINTER filtered file] [cutoff] > [SPLINTER final file]
This step will generate the final SPLINTER output file, which will contain SNPs and Indels present in the sample. Please note that the output for insertions is slightly different than for substitutions or deletions ( Table 2).

Representative Results
We pooled a population of 947 individuals and targeted over 20 kb for sequencing. We applied SPLINTER for the detection of rare variants following our standard protocol. Each individual had previously had genotyping performed by genome-wide array genotyping. Concordance between genotyping of tagged and novel variants called in the pooled sample were excellent (Figure 6). Three variants, two of which (rs3822343 and rs3776110) were rare in the population, were called de novo from the sequencing results and were validated by individual pyrosequencing. Minor allele frequencies (MAF) in the pool were similar to the MAF reported in dbSNP build 129. The MAF concordance between pyrosequencing and pooled sequencing was excellent ( Table 3).

Discussion
There is increasing evidence that the incidence and therapeutic response of common, complex phenotypes and diseases such as obesity 8 , hypercholesterolemia 4 , hypertension 7 and others may be moderated by personal profiles of rare variation. Identifying the genes and pathways where these variants aggregate in affected populations will have profound diagnostic and therapeutic implications, but analyzing affected individuals separately can be time and cost prohibitive. Population-based analysis offers a more efficient method for surveying genetic variation at multiple loci.
We present a novel pooled-DNA sequencing protocol paired with the SPLINTER software package designed to identify this type of genetic variation across populations. We demonstrate the accuracy of this method in identifying and quantifying minor alleles within a large pooled population of 947 individuals, including rare variants that were called de novo from the pooled sequencing and validated by individual pyrosequencing. Our strategy mainly differs from other protocols by the incorporation of a positive and a negative control within every experiment. This allows SPLINTER to achieve much higher accuracy and power compared to other approaches 1 . The optimal coverage of 25-fold per allele is fixed independently of the size of the pool, making the analysis of large pools feasible as this requirement only scales linearly with the pool size. Our approach is very flexible and can be applied to any phenotype of interest but also to samples that are naturally heterogeneous, such as mixed cell populations and tumor biopsies. Given the ever-increasing interest in pooled sequencing from large target regions such as the exome or genome, our library prep and SPLINTER analysis is compatible with custom-capture and whole-exome sequencing, but the alignment utility in the SPLINTER package was not designed for large references sequences. Therefore, we have successfully utilized the dynamic programming aligner, Novoalign, for genome-wide alignments followed by variant calling from the pooled sample (Ramos et al., submitted). Thus, our pooled sequencing strategy can scale successfully to larger pools with increasing amounts of target sequence.

Disclosures
No conflicts of interest declared.