BACKGROUND
evidence that rare genetic variation modulates common, complex phenotypes and diseases continues to mount  <cit> , which has spurred efforts to study larger affected cohorts in order to have adequate statistical power for causative rare variant identification.
however, despite the ever-declining costs of sequencing, recent efforts to sequence entire exomes  <cit>  or use genome-wide array data to inform candidate gene selection  <cit>  result in a huge expansion of the number of candidate genes to be surveyed.
yet, sequencing a large number of candidate loci in thousands of individuals remains largely cost prohibitive in terms of capture costs, sequencing costs and bioinformatic capability for the standard investigator.
to address this limitation, we previously pioneered efforts in cost-effective, population-based rare variant detection using pooled dna sequencing methods and the pooled sequencing analysis tools, snpseeker and splinter  <cit> .
while highly accurate and computationally efficient, these tools were designed for discreet target regions  and their aligners were inadequate for larger reference sequences.
in parallel, commercial capture-by-hybridization products have proven to be powerful tools for the custom enrichment of coding and/or other genomic regions of interest.
these protocols are designed to hybridize complimentary oligonucleotide probes to the target regions of interest from a single dna sample in a single incubation step.
the efficacy of these methods to identify novel, deleterious monogenetic variants in mendelian disorders has been reported in a variety of studies  <cit> .
newer products offer multiplexing of up to eight individuals in a single exome hybridization capture but require individual preparation of dna libraries prior to hybridization, which increases dna handling and reagent cost.
despite modest levels of available multiplexing, performing custom or exome hybridization capture on a large population of individuals followed by sequencing and data analysis remains a costly and time consuming proposition for most investigators.
to address these limitations, a variety of library preparation and multiplexing strategies, with and without individual indexing, have been put forth in the literature  <cit> .
there are a variety of things to consider when selecting a pooled hybridization method.
we have considered these technical aspects when designing and comparing our method: is indexing used or anonymous pooling; how large of a pool can be created; if indexing is used, is possible cross-contamination of indexes reported and quantified; are duplicate reads accounted for in the data analysis; what percentage of raw data is on or near-target; how large of a capture is performed; how uniform is the resultant data for each individual within the pooled or multiplexed sample; and, finally, how accurate are the results.
cost is certainly an important parameter for consideration, but is difficult to quantify given the rapidly changing environment and access to certain platforms at different institutions.
while not all of these studies report their findings in each of these categories, we find that our strategy compares favorably or surpasses reported findings in many instances , yet remains cost-effective and practical for smaller labs requiring manual sample preparation.
obtaining accurate results from pools without individual indexing, even modestly sized pools of 20– <dig> individuals, can be quite challenging  <cit> .
we have previously shown that the splinter algorithm is highly sensitive for rare variant detection in pools of up to  <dig> individuals  <cit> .
previously, splinter had only been used with pcr-amplified target sequences of modest total size.
here, we have expanded the capability of splinter to larger target sequences acquired through solution hybridization capture without compromising its accuracy in a small pool of five exomes and in a larger pool of  <dig> individuals across  <dig>  mb.
our indexing strategy is compatible with low amounts of input dna, less than  <dig> ng, which is often a limiting resource in large populations of dna samples.
this is primarily possible by eliminating all dna sample transfers, including post-sonication, prior to hybridization and integrating a size-selection solid phase reversible immobilization  bead cleanup strategy similar to that described by fisher  <cit> .
we employ a series of unique indexes that are effectively blocked using a single, novel blocking oligonucleotide incorporating deoxyinosine, rather than perfectly complimentary or degenerate sequences.
this reduces cost and library complexity.
we validated our methodology by comparing sequencing results to individual genotyping data and find an excellent concordance between pooled custom or exome hybridization capture results for common and rare variant detection.
overall, we believe this methodology provides a robust, scalable, accurate option for rare variant detection in complex phenotypes.
RESULTS
individual exome analysis
initial reports of sureselect  <cit>  exome sequencing analysis used the program maq  <cit>  to align short reads to the human genome  <cit> .
subsequent steps involved calling substitutions with maq and then using cross_match  <cit>  for indel calling.
however, we have previously demonstrated that maq is inadequate for the analysis of pooled dna samples  <cit> , and wanted to use the same alignment tool for both the individual and pooled sample to minimize differences in variant calling due to variable alignments.
our first experiment was to compare affymetrix  <dig>  genome-wide array , individual exome and pooled exome sequencing from the same five individuals .
because the aligner included in the splinter software package was not designed for large reference sequences, we elected to use novoalign for sequence alignment, which performs gapped alignments enabling simultaneous calling of substitutions and indels.
we then used samtools  <cit>  for variant calling.
we called all variants with ≥5-fold coverage/chromosome  to demonstrate a high degree of accuracy at low coverage thresholds, which facilitates costs savings through reduced sequencing.
there was an average of  <dig>  total variants called/person  and an average of  <dig>  novel variants per person .
additional file 1: table s <dig> lists sequencing metrics for each individual.
to validate our individual and pooled sequencing results, we compared common positions between the exome sequencing and the gwa, as well as individual genotyping by sequenom massarray  <cit>  for  <dig> positions, including six novel variants, not on the gwa and falling within genes of potential biological relevance for the phenotype under investigation within these individuals.
additional file 1: figure s2a shows a highly reproducible average sensitivity across all five individuals of  <dig> % .
similarly, additional file 1: figure s2b shows a highly reproducible average specificity of  <dig> % .
additional file 1: table s <dig> lists results from individual genotyping compared to individual and pooled exome sequencing.
we found that  <dig> % were identical, including all novel variants, between individual exome sequence analysis and genotyping.
these data demonstrate that our bioinformatic pipeline is not only computationally efficient, but also accurate even with as little as 5-fold coverage/chromosome.
pooled exome analysis
after establishing the accuracy of the single exome analysis, we focused on adapting the novoalign alignment pipeline for integration with splinter.
to test this, we prepared a pool of normalized dna from the same five individuals previously analyzed separately .
systematic sequencing platform errors were modeled and subtracted from raw data as previously published  <cit> , additional file 1: figure s1].
additional file 1: figure s2a shows that our pooled sequencing sensitivity slightly outperformed our individual exome sensitivity, demonstrating an accuracy of  <dig> %.
this improvement is likely due to the fact that we only considered called variants from pooled sequencing with ≥20-fold coverage based on prior results  <cit> , supporting the idea that we could improve our individual exome variant calling sensitivity by restricting to positions with higher coverage thresholds.
figure  <dig> shows the correlation for substitution minor allele frequencies  between the individual and pooled sequencing to be r2 =  <dig> , suggesting uniform capture and representation of the minor alleles in the pooled sample.
positions deviating from the idealized correlation are shown to be deficient in sequencing coverage in either the single or pooled sequencing .
additional file 1: figure s2b shows the specificity of the pooled hybridization capture to be  <dig> %.
considering the non-array positions individually genotyped by sequenom , we find a correlation of r2 =  <dig>  between splinter and either samtools or sequenom.
the lower correlation is due to seven no-calls by splinter.
for splinter to call a variant, the complimentary base change must be noted on each strand to a degree that passes the defined p-value cutoff.
in these cases, one strand failed to achieve adequate sequencing that surpasses the p-value threshold, thus leaving these positions uncalled.
pooled custom hybridization capture
for this experiment, we analyzed non-indexed pooled sequencing results of  <dig> candidate genes  from  <dig> individuals enrolled in the long life family study .
each individual had illumina humanomni  <dig> - <dig> beadchip gwa data available for comparison.
using filtering criteria listed in methods under “pooled sequencing analysis of  <dig> individuals at  <dig> genes”, we identified  <dig> substitutions for comparison.
of these,  <dig> were called polymorphic at some frequency by splinter or the gwa, yielding a concordance of r2 =  <dig>  .
there were  <dig> rare substitutions within this dataset  and  <dig>  were called within 1% of the array-based maf by splinter,  <dig> were not called and the remaining  <dig> were called within 2% of the array-based maf.
at least some of the discrepancies were likely due to gwa errors.
when considering only the  <dig> substitutions called homozygous wild type by gwa, we find a specificity of  <dig> %.
however, if we exclude positions from consideration that have an indel within  <dig> bp called by sequencing in at least one individual,  <dig> positions are excluded and our specificity improves to  <dig> %.
this suggests false negative array errors due to nearby indels that adversely affect probe binding.
for indel analysis in this experiment, we aggregated the indel calls from the individual sequencing of all  <dig> individuals and compared against splinter’s calls.
when restricting our analysis to  <dig> bp indels, splinter identified 90%  with a concordance of r2 =  <dig>  .
for rare indel calling , splinter detected  <dig> % .
while these results are good overall, the dropoff in concordance from substitution detection to that of indels further highlights the difficulty of accurate indel detection.
thus, we have expanded the prior capability of the splinter algorithm for use with any size target sequence.
splinter was previously shown to outperform other algorithms, such as crisp  <cit> , for rare variant and indel detection in anonymous pools of up to  <dig> individuals  <cit> , but the aligner was not efficient for survey of large reference sequences.
recent reports of anonymous pooling from smaller pools of 2– <dig> individuals  <cit>  use hybridization capture on smaller capture targets .
both studies removed duplicate reads prior to alignment and report 48-55% of data on target .
however, harakalova failed to validate any called variant predicted at less than 25% maf in pools of  <dig>  <cit> .
day-williams found that solution hybridization capture outperformed pcr or array-based hybridization methods and report a sensitivity and specificity of  <dig> % and  <dig> % for a smaller pool of  <dig> individuals using maq v <dig> .
<dig> for alignment and syzygy for snv calling with no analysis of indels.
bansal and colleagues surveyed  <dig> kb in pools of  <dig> and report a sensitivity of 93% at 30-fold coverage for substitutions and a positive predictive value of 87% for indels compared to a false positive rate of  <dig> %.
in comparison, we find a sensitivity of >96% at 20-fold coverage in a pool of four captures of 22– <dig> each.
indexed custom hybridization capture
while pooling samples without indexing is cost-effective for surveying sequence variation between populations, it’s obviously advantageous to directly link a called variant to a particular individual within the multiplex.
a variety of studies have introduced barcodes or indexing into multiplexes ranging in size from 5– <dig> individuals for target sizes of  <dig> - <dig>  mb  <cit> .
these reports start with anywhere from  <dig> - <dig> μg of dna per individual.
rohland and reich created libraries with as little as  <dig> ng of dna/person but report 50% sample loss from their blunt end adapter ligation strategy and therefore recommend using ≥ <dig> ng of dna/pp  <cit> .
we have also adopted a blunt-end adapter ligation strategy incorporating y-shaped adapters and have not seen significant sample loss, which has enabled comparable outcomes from libraries starting with  <dig> ng/pp up through  <dig> ng/pp.
in addition, we present a novel, simplified adapter blocking strategy using deoxyinosines along with a custom set of  <dig> indexes .
in this experiment, we analyzed sequencing data from the same  <dig> candidate genes  in  <dig> llfs individuals and incorporated individual dna indexing prior to pre-hybridization multiplexing.
to test pre-enrichment multiplexing of various sizes, we created nine separate multiplexes of  <dig> ,  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig>  and using 70– <dig> ng dna/pp.
additional file 1: table s <dig> lists how the multiplexes were then combined for three separate lanes of hiseq  <dig> sequencing with 92– <dig> individuals in each lane.
we first studied the single lane of  <dig> individuals from four multiplexes of 22– <dig> individuals each , and found that >96% of raw data contained an identifiable index and that >97% of the individuals  were within a  <dig> -fold difference in raw read counts .
also, individuals did not cluster predominantly within distinct captures , suggesting that the library preparation and hybridization was highly uniform and did not introduce bias in the sequencing output.
after excluding the two failed samples, aligning unique reads and removing pcr duplicates,  <dig> %  of total reads per person were aligned in a trend that closely followed the number of raw reads per sample .
we observed  <dig> %  pcr duplication .
after removing these duplicates,  <dig> %  of the aligned bases were on or within  <dig> bp  of our targeted intervals and  <dig> %  were directly on targeted regions , demonstrating that the percent duplication and capture efficiency were uniform between all samples and the four independent multiplexes.
in addition, we demonstrated that our candidate regions were enriched uniformly across the four multiplexes, showing an average pairwise r2 =  <dig>  .
thus, our method performs uniformly with respect to both aggregate metrics of an entire capture as well as bait-to-bait capture efficiency.
as discussed by rohland and reich  <cit> , many capture methodology papers do not report the number of duplicate reads observed nor specify if duplicate reads are removed prior to alignment  <cit>  and only report the percentage of their reads that align on or near target, not directly on target, which may inflate the amount of usable data.
rohland reports 24% duplicate reads, and we find a comparable percentage of 14-32% duplicate reads in our various multiplexes.
while our percentage of on/near target reads is comparable to the 18-55% on/near target reads reported  <cit> , these numbers may be inflated by the inclusion of duplicates.
in our initial multiplexed samples, we find an average of  <dig> % of our reads remaining after aligning and removing duplicates, and an average of  <dig> % of the remaining bases aligning to our targeted regions.
thus, approximately  <dig> % of our raw data aligns to our targeted regions in our processed aligned files, yielding an average fold enrichment of  <dig> .
overall, we achieved an average mean target coverage of  when targeting  <dig>  mb with one lane of sequencing and yielding ~ <dig> million reads.
while within range of prior data, we expected to obtain a higher on-target percentage.
the most obvious reason for this seemed to be our modified adapter blocking strategy allowing daisy-chaining, but could also be due to limitations in oligonucleotide hybridization efficiency either through non-specific hybridization to non-target genomic regions or because early versions of agilent sureselect custom captures did not incorporate “bait boosting”, which calibrates the relative amount of bait sequences to normalize the expected binding efficiency to each target.
while we could not alter the baits within our capture set, we did improve our blocking strategy by incorporating  <dig> μl of the illumina pe  <dig>  post-hybridization pcr primer sequence .
this step improved our average on-target percentage to >70% over our next  <dig> multiplexed libraries , which is comparable to the 79% on/near target reads using a shortened adapter strategy published by rohland  <cit> .
analysis of the data from the lane with two 48-individual multiplexes  found very similar results for what we observed from multiplexes of 22– <dig> with respect to the percent of raw reads with an identifiable index, the percentage of aligned data per person, and the uniformity of data between individuals and between captures .
this experiment showed a higher fold enrichment of  <dig>  with  <dig> %  pcr duplication .
after removing these duplicates, our on/near-target percentage improved to  <dig> %  and  <dig> %  were directly on targeted regions .
we attribute the higher pcr duplication rate in the multiplexes of  <dig> individuals over the multiplex of 22– <dig> individuals to the fact that we have reduced the library complexity through two mechanisms.
first, we only pooled  <dig> ng of pre-hybridization capture dna per individual instead of  <dig> ng as done with the multiplexes of 22– <dig>  which would result in a lesser amount of unique dna fragments available for the post-adapter ligation enrichment pcr.
second, with  <dig> individuals pooled for a single bait set, there is approximately half as much dna selected per individual as for captures with 22– <dig> individuals pooled.
thus, with similar overall read counts per person, we would expect to see a greater number of pcr duplicates as a result.
as there were over  <dig>  overlapping positions between the omni  <dig> - <dig> array and our targeted regions, we used these validated positions to accurately assess our sensitivity and specificity at various coverage thresholds.
table  <dig> shows sensitivity, specificity, and coverage metrics for the custom multiplexes at different maf thresholds.
as expected, sensitivity is a function of coverage, starting at >95% with >3-fold coverage for common variants and improving to 100% with >20-fold coverage for rare variants .
specificity was consistently ≥ <dig> % and hardly impacted by coverage regardless of maf, confirming that our method does not systematically introduce false positives .
for rare variants with ≤2% maf in this experiment, we find a sensitivity of  <dig> - <dig> % and  <dig> -100% with coverage ranging from ≥3-20-fold for variants with 2% maf and  <dig> % maf, respectively .
specificity was ≥ <dig> % at all coverage levels for either maf.
these results indicate that allele dropout is not occurring with our method.
the percent of bases in the targeted intervals that reached the specified coverage threshold is listed.
homozygous wild type, heterozygous, and homozygous variant sites surveyed indicates the number of each of those positions as seen by the illumina omni  <dig> - <dig> array that were used to determine sensitivity and specificity.
sensitivity is the percentage of heterozygous and homozygous variant sites correctly called heterozygous and homozygous variant, respectively.
specificity is the percentage of homozygous wild type sites called as homozygous wild type.
for “all allele frequencies” these values were averaged among all non-excluded samples.
for ≤5%, 2% and  <dig> % minor allele frequencies, these values were determined from the cumulative metrics of all non-excluded samples at sites with  <dig> or fewer,  <dig> or fewer, or  <dig> variant allele in the entire pool, respectively.
additional file 1: figure s <dig> compares sequencing coverage from a single individual for all common gwa positions against all bases in our  <dig>  mb capture with  <dig>  m paired-end reads.
once a base position surpasses 3-fold coverage, the representation of gwa positions versus any base position with similar coverage is approximately equivalent .
these results support using high-density gwa positions for validation of sequencing results and extrapolating to the entire dataset.
we then analyzed our non-array variant calls from our multiplexes of  <dig>  for accuracy.
any variant found in the exome variant server and not considered a likely false positive by their support vector machine   <cit>  was considered.
of the  <dig>  non-svm, non-gwa variant calls made,  <dig>  are in dbsnp <dig> .
this degree of overlap would not be observed if variants were being called by random chance.
we have also surveyed all variation in our multiplexes of 22– <dig> individuals  to calculate the ti:tv ratio according to the work of depristo  <cit>   ratio analysis”).
we find additional support that we are calling variants as expected, particularly those positions with ≥15-fold coverage.
often, dna samples are a precious resource, with very little material in which to conduct large scale genomic surveys.
we have been able to use relatively little source dna per person by eliminating all sample transfers prior to pooling, which is enabled by having spri beads present in solution from the beginning of the protocol.
these beads can stay in place and do not appear to alter the quality of the library creation reactions.
our strategy is modified from prior reports by rohland  <cit>  and fisher  <cit> .
many investigators do not have access to full automation of the hybridization capture library creation process.
fisher introduced exclusive bead cleanup for use with fully automated single exome library creation, without multiplexing.
while results were highly uniform when employed with full automation, there was a large dropoff in the percentage and uniformity of captured bases when bead cleanup was performed manually [see ref.
<cit> , additional file  <dig>  rohland performed sonication in a 96-well pcr plate without beads present in solution  and also optimized their bead cleanup for automation.
however, as mentioned above, they report a loss of 50% of source dna after adapter ligation.
we have exclusively performed manual library preparations with multichannel pipettors and observe excellent recovery of source dna, which is a cost effective and accurate alternative that is accessible to any lab.
we investigated whether issues such as bait gc content or strand bias adversely impacted our results.
we observed a large dropoff in sequencing coverage with baits that were high in gc content and targets < <dig> bp , and performed an analysis to centeracterize whether strand bias impacted our variant calling .
from these results, we believe that gc bias in hybridization efficiency can be mitigated through “bait boosting”  and strand bias did not adversely effected our results.
however, one serious potential drawback in integrating large scale indexing would be data contamination due to index switching, or the transfer of an index from one individual to the sequencing read of another individual.
unexpectedly, the multiplexes of 22– <dig> and  <dig> demonstrated superior specificity for rare variants rather than common variants, which seems counterintuitive, but is a consistent finding at all coverage thresholds .
to explore this result further, we considered the validated positions binned by the number of variant alleles at each position called by array genotyping in the entire pool of  <dig> captured in two multiplexes of  <dig> .
the percentage of reads with variant alleles in individuals who were homozygous wild type at these positions increases as a function of the frequency of the variant allele in each bin.
this indicated to us that there was up to 9% of reads with an inappropriate index in our samples.
we use a large excess of adapters to drive ligation of every possible template strand.
we surmised that a large percentage of our misattributed indexes were due to an excess of unligated adapters present in the pcr enrichment, as shown by the bimodal peaks in the green tracing on additional file 1: figure s <dig> and schematically diagrammed in additional file 1: figure s <dig>  we tested this by adding an additional bead purification step to remove to remove unligated adapters from our multiplexes of 30– <dig> individuals .
we further analyzed variants with a maf of ~50%  and see an equivalent specificity for common and rare variants with the inclusion of this step .
in addition, the lesser degree of misattributed indexes in the multiplexes of 30– <dig> improves sensitivity for common and rare variants compared to the prior multiplexes.
our results suggest that index switching is a random and stochastic process resulting in approximately  <dig> % of reads in the entire pool attributed to an incorrect index.
we propose a mechanism for index switching in additional file 1: figure s <dig>  none of the previously mentioned reports on multiplexed hybridization capture discuss this  <cit>  which may further skew their findings, but one possible source of this observation is covered in detail by kircher and colleagues  <cit>  as “jumping pcr”.
kircher reports that over  <dig> % of misattributed indexes are expected when pooling samples prior to pcr for multiplexed capture.
if we subtract this source of error from our result, we are left with  <dig> % of reads with an incorrect index .
the source of this is unclear, but may be caused by further residual unligated indexes in our sample prior to pcr enrichment, a different rate of “jumping pcr” due to innate differences in polymerases used , or perhaps erroneous array calls.
thus, while our percentage of incorrectly indexed reads was not as low as reported by the double indexing strategy of kircher  <cit> , the amount we did have did not appear to impact our accuracy for common or rare variant detection.
CONCLUSIONS
the genomic analysis of rare variants in complex disease will continue to drive pathophysiological research and clinical questions.
sequencing many individuals at any number of candidate loci is required to obtain the statistical power necessary for detecting disease-related, rare functional genetic variants associated with complex phenotypes.
as sequencing capacity continues to increase, the need for additional multiplexing is required to maximize the efficiency of the sequencing output.
pooling will also enable sequencing to move more readily into the diagnostic arena as a more cost effective method for screening multiple genes for rare functional variation.
we report a cost-effective and fully scalable, pre-hybridization multiplexed in-solution capture method, amenable to < <dig> ng of input dna per person due to removal of sample transfers during library preparation enabled by the constant presence of spri beads and resulting in a 50% reduction in the volume of library preparation reagent necessary.
sequencing analysis was performed by bioinformatic pipelines aligning raw data with novoalign and calling variants from anonymous pools using splinter or indexed samples using samtools.
both pipelines showed consistent, coverage-dependent high sensitivity ranging from 95-99% for substitutions with a dropoff in concordance for detecting indels.
for rare variants , we achieved a sensitivity and specificity of  <dig> % and  <dig> % for pooled exome analysis and ≥ <dig> % and ≥ <dig> %, respectively, for custom capture libraries of 22– <dig> individuals.
while within the range of published reports, we found our on-target percentage to generally be in the 40-50% range, but has improved with the addition of another post-enrichment pcr primer to act as an additional blocker.
our blocking strategy utilizes deoxyinosines to span the novel  <dig> bp indexes rather than degenerate nucleotides to reduce the complexity of the dna mixture.
in our system, we found that optimal hybridization and sequencing representation occurred with baits having 25-35% gc content and targets of > <dig> bp.
we also found no negative impact due to strand bias.
however, one serious drawback to larger multiplexing of indexed samples is index switching leading to the attribution of variants to the incorrect individual.
we have quantified the prevalence of this phenomenon and shown it to occur most frequently at wild-type or common variant positions simply due to the abundance of those reads in the sample.
by incorporating an additional purification step during library preparation to remove excess unligated adapters, this effect can be mitigated.
in sum, cost savings can be realized through the use of less pre-enrichment library preparation reagents, use of more standard plasticware for sonication, capturing up to  <dig> individuals with a single hybridization, multiplexing multiple libraries for sequencing, and the ability to call rare variants accurately with as little as 5-fold sequencing as a lower threshold .
while there are a variety of reports describing pooled sequencing or multiplexed solution hybridization capture  <cit> , our method offers a highly uniform, reproducible workflow for smaller labs that may not have access to high-throughput robotics but still wish to quantify rare variation within specific populations.

 evidence that rare genetic variation modulates common, complex phenotypes and diseases continues to mount, which has spurred efforts to study larger affected cohorts in order to have adequate statistical power for causative rare variant identification. however, despite the ever-declining costs of sequencing, recent efforts to sequence entire exomes or use genome-wide array data to inform candidate gene selection result in a huge expansion of the number of candidate genes to be surveyed. yet, sequencing a large number of candidate loci in thousands of individuals remains largely cost prohibitive in terms of capture costs, sequencing costs and bioinformatic capability for the standard investigator. to address this limitation, we previously pioneered efforts in cost-effective, population-based rare variant detection using pooled dna sequencing methods and the pooled sequencing analysis tools, snpseeker and splinter. while highly accurate and computationally efficient, these tools were designed for discreet target regions and their aligners were inadequate for larger reference sequences. in parallel, commercial capture-by-hybridization products have proven to be powerful tools for the custom enrichment of coding and/or other genomic regions of interest. these protocols are designed to hybridize complimentary oligonucleotide probes to the target regions of interest from a single dna sample in a single incubation step. the efficacy of these methods to identify novel, deleterious monogenetic variants in mendelian disorders has been reported in a variety of studies. newer products offer multiplexing of up to eight individuals in a single exome hybridization capture but require individual preparation of dna libraries prior to hybridization, which increases dna handling and reagent cost. despite modest levels of available multiplexing, performing custom or exome hybridization capture on a large population of individuals followed by sequencing and data analysis remains a costly and time consuming proposition for most investigators. to address these limitations, a variety of library preparation and multiplexing strategies, with and without individual indexing, have been put forth in the literature. there are a variety of things to consider when selecting a pooled hybridization method. we have considered these technical aspects when designing and comparing our method: is indexing used or anonymous pooling; how large of a pool can be created; if indexing is used, is possible cross-contamination of indexes reported and quantified; are duplicate reads accounted for in the data analysis; what percentage of raw data is on or near-target; how large of a capture is performed; how uniform is the resultant data for each individual within the pooled or multiplexed sample; and, finally, how accurate are the results. cost is certainly an important parameter for consideration, but is difficult to quantify given the rapidly changing environment and access to certain platforms at different institutions. while not all of these studies report their findings in each of these categories, we find that our strategy compares favorably or surpasses reported findings in many instances, yet remains cost-effective and practical for smaller labs requiring manual sample preparation. obtaining accurate results from pools without individual indexing, even modestly sized pools of 20– individuals, can be quite challenging. we have previously shown that the splinter algorithm is highly sensitive for rare variant detection in pools of up to individuals. previously, splinter had only been used with pcr-amplified target sequences of modest total size. here, we have expanded the capability of splinter to larger target sequences acquired through solution hybridization capture without compromising its accuracy in a small pool of five exomes and in a larger pool of individuals across mb. our indexing strategy is compatible with low amounts of input dna, less than ng, which is often a limiting resource in large populations of dna samples. this is primarily possible by eliminating all dna sample transfers, including post-sonication, prior to hybridization and integrating a size-selection solid phase reversible immobilization bead cleanup strategy similar to that described by fisher. we employ a series of unique indexes that are effectively blocked using a single, novel blocking oligonucleotide incorporating deoxyinosine, rather than perfectly complimentary or degenerate sequences. this reduces cost and library complexity. we validated our methodology by comparing sequencing results to individual genotyping data and find an excellent concordance between pooled custom or exome hybridization capture results for common and rare variant detection. overall, we believe this methodology provides a robust, scalable, accurate option for rare variant detection in complex phenotypes.  individual exome analysis initial reports of sureselect exome sequencing analysis used the program maq to align short reads to the human genome. subsequent steps involved calling substitutions with maq and then using cross_match for indel calling. however, we have previously demonstrated that maq is inadequate for the analysis of pooled dna samples, and wanted to use the same alignment tool for both the individual and pooled sample to minimize differences in variant calling due to variable alignments. our first experiment was to compare affymetrix genome-wide array, individual exome and pooled exome sequencing from the same five individuals. because the aligner included in the splinter software package was not designed for large reference sequences, we elected to use novoalign for sequence alignment, which performs gapped alignments enabling simultaneous calling of substitutions and indels. we then used samtools for variant calling. we called all variants with ≥5-fold coverage/chromosome to demonstrate a high degree of accuracy at low coverage thresholds, which facilitates costs savings through reduced sequencing. there was an average of total variants called/person and an average of novel variants per person. additional file 1: table s lists sequencing metrics for each individual. to validate our individual and pooled sequencing results, we compared common positions between the exome sequencing and the gwa, as well as individual genotyping by sequenom massarray for positions, including six novel variants, not on the gwa and falling within genes of potential biological relevance for the phenotype under investigation within these individuals. additional file 1: figure s2a shows a highly reproducible average sensitivity across all five individuals of %. similarly, additional file 1: figure s2b shows a highly reproducible average specificity of %. additional file 1: table s lists results from individual genotyping compared to individual and pooled exome sequencing. we found that % were identical, including all novel variants, between individual exome sequence analysis and genotyping. these data demonstrate that our bioinformatic pipeline is not only computationally efficient, but also accurate even with as little as 5-fold coverage/chromosome. pooled exome analysis after establishing the accuracy of the single exome analysis, we focused on adapting the novoalign alignment pipeline for integration with splinter. to test this, we prepared a pool of normalized dna from the same five individuals previously analyzed separately. systematic sequencing platform errors were modeled and subtracted from raw data as previously published, additional file 1: figure s1]. additional file 1: figure s2a shows that our pooled sequencing sensitivity slightly outperformed our individual exome sensitivity, demonstrating an accuracy of %. this improvement is likely due to the fact that we only considered called variants from pooled sequencing with ≥20-fold coverage based on prior results, supporting the idea that we could improve our individual exome variant calling sensitivity by restricting to positions with higher coverage thresholds. figure shows the correlation for substitution minor allele frequencies between the individual and pooled sequencing to be r2 = , suggesting uniform capture and representation of the minor alleles in the pooled sample. positions deviating from the idealized correlation are shown to be deficient in sequencing coverage in either the single or pooled sequencing. additional file 1: figure s2b shows the specificity of the pooled hybridization capture to be %. considering the non-array positions individually genotyped by sequenom, we find a correlation of r2 =  between splinter and either samtools or sequenom. the lower correlation is due to seven no-calls by splinter. for splinter to call a variant, the complimentary base change must be noted on each strand to a degree that passes the defined p-value cutoff. in these cases, one strand failed to achieve adequate sequencing that surpasses the p-value threshold, thus leaving these positions uncalled. pooled custom hybridization capture for this experiment, we analyzed non-indexed pooled sequencing results of candidate genes from individuals enrolled in the long life family study. each individual had illumina humanomni - beadchip gwa data available for comparison. using filtering criteria listed in methods under “pooled sequencing analysis of individuals at genes”, we identified substitutions for comparison. of these, were called polymorphic at some frequency by splinter or the gwa, yielding a concordance of r2 = . there were rare substitutions within this dataset and were called within 1% of the array-based maf by splinter, were not called and the remaining were called within 2% of the array-based maf. at least some of the discrepancies were likely due to gwa errors. when considering only the substitutions called homozygous wild type by gwa, we find a specificity of %. however, if we exclude positions from consideration that have an indel within bp called by sequencing in at least one individual, positions are excluded and our specificity improves to %. this suggests false negative array errors due to nearby indels that adversely affect probe binding. for indel analysis in this experiment, we aggregated the indel calls from the individual sequencing of all individuals and compared against splinter’s calls. when restricting our analysis to bp indels, splinter identified 90% with a concordance of r2 = . for rare indel calling, splinter detected %. while these results are good overall, the dropoff in concordance from substitution detection to that of indels further highlights the difficulty of accurate indel detection. thus, we have expanded the prior capability of the splinter algorithm for use with any size target sequence. splinter was previously shown to outperform other algorithms, such as crisp, for rare variant and indel detection in anonymous pools of up to individuals, but the aligner was not efficient for survey of large reference sequences. recent reports of anonymous pooling from smaller pools of 2– individuals use hybridization capture on smaller capture targets. both studies removed duplicate reads prior to alignment and report 48-55% of data on target. however, harakalova failed to validate any called variant predicted at less than 25% maf in pools of. day-williams found that solution hybridization capture outperformed pcr or array-based hybridization methods and report a sensitivity and specificity of % and % for a smaller pool of individuals using maq v. for alignment and syzygy for snv calling with no analysis of indels. bansal and colleagues surveyed kb in pools of and report a sensitivity of 93% at 30-fold coverage for substitutions and a positive predictive value of 87% for indels compared to a false positive rate of %. in comparison, we find a sensitivity of >96% at 20-fold coverage in a pool of four captures of 22– each. indexed custom hybridization capture while pooling samples without indexing is cost-effective for surveying sequence variation between populations, it’s obviously advantageous to directly link a called variant to a particular individual within the multiplex. a variety of studies have introduced barcodes or indexing into multiplexes ranging in size from 5– individuals for target sizes of - mb. these reports start with anywhere from - μg of dna per individual. rohland and reich created libraries with as little as ng of dna/person but report 50% sample loss from their blunt end adapter ligation strategy and therefore recommend using ≥ ng of dna/pp. we have also adopted a blunt-end adapter ligation strategy incorporating y-shaped adapters and have not seen significant sample loss, which has enabled comparable outcomes from libraries starting with ng/pp up through ng/pp. in addition, we present a novel, simplified adapter blocking strategy using deoxyinosines along with a custom set of indexes. in this experiment, we analyzed sequencing data from the same candidate genes in llfs individuals and incorporated individual dna indexing prior to pre-hybridization multiplexing. to test pre-enrichment multiplexing of various sizes, we created nine separate multiplexes of,,,, and and using 70– ng dna/pp. additional file 1: table s lists how the multiplexes were then combined for three separate lanes of hiseq sequencing with 92– individuals in each lane. we first studied the single lane of individuals from four multiplexes of 22– individuals each, and found that >96% of raw data contained an identifiable index and that >97% of the individuals were within a -fold difference in raw read counts. also, individuals did not cluster predominantly within distinct captures, suggesting that the library preparation and hybridization was highly uniform and did not introduce bias in the sequencing output. after excluding the two failed samples, aligning unique reads and removing pcr duplicates, % of total reads per person were aligned in a trend that closely followed the number of raw reads per sample. we observed % pcr duplication. after removing these duplicates, % of the aligned bases were on or within bp of our targeted intervals and % were directly on targeted regions, demonstrating that the percent duplication and capture efficiency were uniform between all samples and the four independent multiplexes. in addition, we demonstrated that our candidate regions were enriched uniformly across the four multiplexes, showing an average pairwise r2 = . thus, our method performs uniformly with respect to both aggregate metrics of an entire capture as well as bait-to-bait capture efficiency. as discussed by rohland and reich, many capture methodology papers do not report the number of duplicate reads observed nor specify if duplicate reads are removed prior to alignment and only report the percentage of their reads that align on or near target, not directly on target, which may inflate the amount of usable data. rohland reports 24% duplicate reads, and we find a comparable percentage of 14-32% duplicate reads in our various multiplexes. while our percentage of on/near target reads is comparable to the 18-55% on/near target reads reported, these numbers may be inflated by the inclusion of duplicates. in our initial multiplexed samples, we find an average of % of our reads remaining after aligning and removing duplicates, and an average of % of the remaining bases aligning to our targeted regions. thus, approximately % of our raw data aligns to our targeted regions in our processed aligned files, yielding an average fold enrichment of. overall, we achieved an average mean target coverage of when targeting mb with one lane of sequencing and yielding ~ million reads. while within range of prior data, we expected to obtain a higher on-target percentage. the most obvious reason for this seemed to be our modified adapter blocking strategy allowing daisy-chaining, but could also be due to limitations in oligonucleotide hybridization efficiency either through non-specific hybridization to non-target genomic regions or because early versions of agilent sureselect custom captures did not incorporate “bait boosting”, which calibrates the relative amount of bait sequences to normalize the expected binding efficiency to each target. while we could not alter the baits within our capture set, we did improve our blocking strategy by incorporating μl of the illumina pe post-hybridization pcr primer sequence. this step improved our average on-target percentage to >70% over our next multiplexed libraries, which is comparable to the 79% on/near target reads using a shortened adapter strategy published by rohland. analysis of the data from the lane with two 48-individual multiplexes found very similar results for what we observed from multiplexes of 22– with respect to the percent of raw reads with an identifiable index, the percentage of aligned data per person, and the uniformity of data between individuals and between captures. this experiment showed a higher fold enrichment of with % pcr duplication. after removing these duplicates, our on/near-target percentage improved to % and % were directly on targeted regions. we attribute the higher pcr duplication rate in the multiplexes of individuals over the multiplex of 22– individuals to the fact that we have reduced the library complexity through two mechanisms. first, we only pooled ng of pre-hybridization capture dna per individual instead of ng as done with the multiplexes of 22– which would result in a lesser amount of unique dna fragments available for the post-adapter ligation enrichment pcr. second, with individuals pooled for a single bait set, there is approximately half as much dna selected per individual as for captures with 22– individuals pooled. thus, with similar overall read counts per person, we would expect to see a greater number of pcr duplicates as a result. as there were over overlapping positions between the omni - array and our targeted regions, we used these validated positions to accurately assess our sensitivity and specificity at various coverage thresholds. table shows sensitivity, specificity, and coverage metrics for the custom multiplexes at different maf thresholds. as expected, sensitivity is a function of coverage, starting at >95% with >3-fold coverage for common variants and improving to 100% with >20-fold coverage for rare variants. specificity was consistently ≥ % and hardly impacted by coverage regardless of maf, confirming that our method does not systematically introduce false positives. for rare variants with ≤2% maf in this experiment, we find a sensitivity of - % and -100% with coverage ranging from ≥3-20-fold for variants with 2% maf and % maf, respectively. specificity was ≥ % at all coverage levels for either maf. these results indicate that allele dropout is not occurring with our method. the percent of bases in the targeted intervals that reached the specified coverage threshold is listed. homozygous wild type, heterozygous, and homozygous variant sites surveyed indicates the number of each of those positions as seen by the illumina omni - array that were used to determine sensitivity and specificity. sensitivity is the percentage of heterozygous and homozygous variant sites correctly called heterozygous and homozygous variant, respectively. specificity is the percentage of homozygous wild type sites called as homozygous wild type. for “all allele frequencies” these values were averaged among all non-excluded samples. for ≤5%, 2% and % minor allele frequencies, these values were determined from the cumulative metrics of all non-excluded samples at sites with or fewer, or fewer, or variant allele in the entire pool, respectively. additional file 1: figure s compares sequencing coverage from a single individual for all common gwa positions against all bases in our mb capture with m paired-end reads. once a base position surpasses 3-fold coverage, the representation of gwa positions versus any base position with similar coverage is approximately equivalent. these results support using high-density gwa positions for validation of sequencing results and extrapolating to the entire dataset. we then analyzed our non-array variant calls from our multiplexes of for accuracy. any variant found in the exome variant server and not considered a likely false positive by their support vector machine was considered. of the non-svm, non-gwa variant calls made, are in dbsnp. this degree of overlap would not be observed if variants were being called by random chance. we have also surveyed all variation in our multiplexes of 22– individuals to calculate the ti:tv ratio according to the work of depristo ratio analysis”). we find additional support that we are calling variants as expected, particularly those positions with ≥15-fold coverage. often, dna samples are a precious resource, with very little material in which to conduct large scale genomic surveys. we have been able to use relatively little source dna per person by eliminating all sample transfers prior to pooling, which is enabled by having spri beads present in solution from the beginning of the protocol. these beads can stay in place and do not appear to alter the quality of the library creation reactions. our strategy is modified from prior reports by rohland and fisher. many investigators do not have access to full automation of the hybridization capture library creation process. fisher introduced exclusive bead cleanup for use with fully automated single exome library creation, without multiplexing. while results were highly uniform when employed with full automation, there was a large dropoff in the percentage and uniformity of captured bases when bead cleanup was performed manually [see ref., additional file rohland performed sonication in a 96-well pcr plate without beads present in solution and also optimized their bead cleanup for automation. however, as mentioned above, they report a loss of 50% of source dna after adapter ligation. we have exclusively performed manual library preparations with multichannel pipettors and observe excellent recovery of source dna, which is a cost effective and accurate alternative that is accessible to any lab. we investigated whether issues such as bait gc content or strand bias adversely impacted our results. we observed a large dropoff in sequencing coverage with baits that were high in gc content and targets < bp, and performed an analysis to centeracterize whether strand bias impacted our variant calling. from these results, we believe that gc bias in hybridization efficiency can be mitigated through “bait boosting” and strand bias did not adversely effected our results. however, one serious potential drawback in integrating large scale indexing would be data contamination due to index switching, or the transfer of an index from one individual to the sequencing read of another individual. unexpectedly, the multiplexes of 22– and demonstrated superior specificity for rare variants rather than common variants, which seems counterintuitive, but is a consistent finding at all coverage thresholds. to explore this result further, we considered the validated positions binned by the number of variant alleles at each position called by array genotyping in the entire pool of captured in two multiplexes of. the percentage of reads with variant alleles in individuals who were homozygous wild type at these positions increases as a function of the frequency of the variant allele in each bin. this indicated to us that there was up to 9% of reads with an inappropriate index in our samples. we use a large excess of adapters to drive ligation of every possible template strand. we surmised that a large percentage of our misattributed indexes were due to an excess of unligated adapters present in the pcr enrichment, as shown by the bimodal peaks in the green tracing on additional file 1: figure s and schematically diagrammed in additional file 1: figure s we tested this by adding an additional bead purification step to remove to remove unligated adapters from our multiplexes of 30– individuals. we further analyzed variants with a maf of ~50% and see an equivalent specificity for common and rare variants with the inclusion of this step. in addition, the lesser degree of misattributed indexes in the multiplexes of 30– improves sensitivity for common and rare variants compared to the prior multiplexes. our results suggest that index switching is a random and stochastic process resulting in approximately % of reads in the entire pool attributed to an incorrect index. we propose a mechanism for index switching in additional file 1: figure s none of the previously mentioned reports on multiplexed hybridization capture discuss this which may further skew their findings, but one possible source of this observation is covered in detail by kircher and colleagues as “jumping pcr”. kircher reports that over % of misattributed indexes are expected when pooling samples prior to pcr for multiplexed capture. if we subtract this source of error from our result, we are left with % of reads with an incorrect index. the source of this is unclear, but may be caused by further residual unligated indexes in our sample prior to pcr enrichment, a different rate of “jumping pcr” due to innate differences in polymerases used, or perhaps erroneous array calls. thus, while our percentage of incorrectly indexed reads was not as low as reported by the double indexing strategy of kircher, the amount we did have did not appear to impact our accuracy for common or rare variant detection.  the genomic analysis of rare variants in complex disease will continue to drive pathophysiological research and clinical questions. sequencing many individuals at any number of candidate loci is required to obtain the statistical power necessary for detecting disease-related, rare functional genetic variants associated with complex phenotypes. as sequencing capacity continues to increase, the need for additional multiplexing is required to maximize the efficiency of the sequencing output. pooling will also enable sequencing to move more readily into the diagnostic arena as a more cost effective method for screening multiple genes for rare functional variation. we report a cost-effective and fully scalable, pre-hybridization multiplexed in-solution capture method, amenable to < ng of input dna per person due to removal of sample transfers during library preparation enabled by the constant presence of spri beads and resulting in a 50% reduction in the volume of library preparation reagent necessary. sequencing analysis was performed by bioinformatic pipelines aligning raw data with novoalign and calling variants from anonymous pools using splinter or indexed samples using samtools. both pipelines showed consistent, coverage-dependent high sensitivity ranging from 95-99% for substitutions with a dropoff in concordance for detecting indels. for rare variants, we achieved a sensitivity and specificity of % and % for pooled exome analysis and ≥ % and ≥ %, respectively, for custom capture libraries of 22– individuals. while within the range of published reports, we found our on-target percentage to generally be in the 40-50% range, but has improved with the addition of another post-enrichment pcr primer to act as an additional blocker. our blocking strategy utilizes deoxyinosines to span the novel bp indexes rather than degenerate nucleotides to reduce the complexity of the dna mixture. in our system, we found that optimal hybridization and sequencing representation occurred with baits having 25-35% gc content and targets of > bp. we also found no negative impact due to strand bias. however, one serious drawback to larger multiplexing of indexed samples is index switching leading to the attribution of variants to the incorrect individual. we have quantified the prevalence of this phenomenon and shown it to occur most frequently at wild-type or common variant positions simply due to the abundance of those reads in the sample. by incorporating an additional purification step during library preparation to remove excess unligated adapters, this effect can be mitigated. in sum, cost savings can be realized through the use of less pre-enrichment library preparation reagents, use of more standard plasticware for sonication, capturing up to individuals with a single hybridization, multiplexing multiple libraries for sequencing, and the ability to call rare variants accurately with as little as 5-fold sequencing as a lower threshold. while there are a variety of reports describing pooled sequencing or multiplexed solution hybridization capture, our method offers a highly uniform, reproducible workflow for smaller labs that may not have access to high-throughput robotics but still wish to quantify rare variation within specific populations. 
these results were obtained aligning raw reads against the entire genome using novoalign followed by variant calling of non-indexed pools using splinter or samtools for indexed samples. we set out to 1) create a scalable, multiplexing method for custom capture with or without individual dna indexing that was amenable to low amounts of input dna and 2) expand the functionality of the splinter algorithm for calling substitutions, insertions and deletions across either candidate genes or the entire exome by integrating the variant calling algorithm with the dynamic programming aligner, novoalign. in addition, dna samples are often precious and hybridization methods typically require large amounts of input dna. however, despite persistent declines in sequencing costs, population-based rare variant detection across large genomic target regions remains cost prohibitive for most investigators. for rare variants, we achieved sensitivity and specificity of ≥ % and ≥ % for custom capture of mb in multiplexed libraries of 22– individuals with only ≥5-fold coverage/chromosome, but these parameters improved to ≥ and 100% with 20-fold coverage/chromosome. rare genetic variation in the human population is a major source of pathophysiological variability and has been implicated in a host of complex phenotypes and diseases. finding disease-related genes harboring disparate functional rare variants requires sequencing of many individuals across many genomic regions and comparing against unaffected cohorts. this highly scalable methodology enables accurate rare variant detection, with or without individual dna sample indexing, while reducing the amount of required source dna and total costs through less hybridization reagent consumption, multi-sample sonication in a standard pcr plate, multiplexed pre-enrichment pooling with a single hybridization and lesser sequencing coverage required to obtain high sensitivity. pooled sample dna sequencing is a cost and time-efficient strategy for surveying populations of individuals for rare variants. custom y-shaped adapters containing novel base pair index sequences with a hamming distance of ≥ were directly ligated onto fragmented source dna eliminating the need for pcr to incorporate indexes, and was followed by a custom blocking strategy using a single oligonucleotide regardless of index sequence. modified solid phase reversible immobilization bead purification strategies enable no sample transfers from sonication in 96-well plates through adapter ligation, resulting in 50% less library preparation reagent consumption.
 rare genetic variation in the human population is a major source of pathophysiological variability and has been implicated in a host of complex phenotypes and diseases. finding disease-related genes harboring disparate functional rare variants requires sequencing of many individuals across many genomic regions and comparing against unaffected cohorts. however, despite persistent declines in sequencing costs, population-based rare variant detection across large genomic target regions remains cost prohibitive for most investigators. in addition, dna samples are often precious and hybridization methods typically require large amounts of input dna. pooled sample dna sequencing is a cost and time-efficient strategy for surveying populations of individuals for rare variants. we set out to 1) create a scalable, multiplexing method for custom capture with or without individual dna indexing that was amenable to low amounts of input dna and 2) expand the functionality of the splinter algorithm for calling substitutions, insertions and deletions across either candidate genes or the entire exome by integrating the variant calling algorithm with the dynamic programming aligner, novoalign.  we report methodology for pooled hybridization capture with pre-enrichment, indexed multiplexing of up to individuals or non-indexed pooled sequencing of up to individuals with as little as ng of dna per person. modified solid phase reversible immobilization bead purification strategies enable no sample transfers from sonication in 96-well plates through adapter ligation, resulting in 50% less library preparation reagent consumption. custom y-shaped adapters containing novel base pair index sequences with a hamming distance of ≥ were directly ligated onto fragmented source dna eliminating the need for pcr to incorporate indexes, and was followed by a custom blocking strategy using a single oligonucleotide regardless of index sequence. these results were obtained aligning raw reads against the entire genome using novoalign followed by variant calling of non-indexed pools using splinter or samtools for indexed samples. with these pipelines, we find sensitivity and specificity of % and % for pooled exome sequencing. sensitivity, and to a lesser degree specificity, proved to be a function of coverage. for rare variants, we achieved sensitivity and specificity of ≥ % and ≥ % for custom capture of mb in multiplexed libraries of 22– individuals with only ≥5-fold coverage/chromosome, but these parameters improved to ≥ and 100% with 20-fold coverage/chromosome.  this highly scalable methodology enables accurate rare variant detection, with or without individual dna sample indexing, while reducing the amount of required source dna and total costs through less hybridization reagent consumption, multi-sample sonication in a standard pcr plate, multiplexed pre-enrichment pooling with a single hybridization and lesser sequencing coverage required to obtain high sensitivity. rare variantsgenomicsexomehybridization capturemultiplexed captureindexed capturesplinter 