BACKGROUND
directed sequencing is important for the discovery of novel polymorphisms in human dna sequences that may play a role in disease susceptibility and progression  <cit> .
in order to efficiently screen large numbers of target regions across a diverse panel of dna samples, a sequencing center needs to apply a high-throughput approach with the goal of minimizing cost.
to achieve this, pcr primers are designed with the aim of eliminating human labor associated with amplification protocol optimization and other manual intervention caused by primer failures.
recent contributions in computational pcr primer design include screening additional thermodynamic attributes  <cit>  and increasing the genomic annotation  <cit>  considered in the design process.
table  <dig> provides a comparison of novel features of our primer design pipeline in relation to recently published primer design software  <cit> .
our primer design pipeline is further enhanced by a higher degree of operational automation and the parameterization of many key attributes in its design.
we have employed this pipeline for over five years at our sequencing center and its use was instrumental in the publication of mutations discovered in a set of  <dig> tyrosine kinase receptor domains consisting of  <dig> exons, across  <dig> tumor samples in glioblastomas  <cit> .
since  <dig>  over  <dig>  mb of sequence have been targeted with a primer design target coverage rate exceeding 95% .
our directed sequencing pipeline has achieved a sequencing success rate of over 90%  and we have detected over  <dig> potential variations in  <dig> genes,  <dig> of which were not previously recorded in ncbi's dbsnp.
∗ feature supported
- feature not supported
‡ tiling pattern is fixed
† searches for the full length of the primer, instead of just the critical region on the 3' end
implementation
computational pcr primer design
software
in the method described here, we constrain the pcr primer design software to generate a set of primers with a high likelihood of success using one of two standardized amplification protocols.
the high-throughput computational pcr primer design pipeline is one of the first computational stages we employ as part of our directed sequencing pipeline .
the computational pcr primer design software is available on the sourceforge.net jcvi primer designer website  <cit> .
the tar file of the latest release is available by navigating through the 'download' tab.
it can be installed and run on any computer with a unix operating system.
inputs
target regions are specified by the investigator through a variety of identifiers.
by using the feature rich ensembl perl application programming interface   <cit> , the following types of identifiers are supported as input: gene, transcript, hugo, exon, and dbsnp.
since the reference genome assembly and annotation may be updated, using identifiers allows us to track features of interest even if their absolute genomic coordinates have changed over time.
in addition, by using these identifiers as anchors, we are able to augment target regions with a specified number of flanking bases upstream and downstream of the region, and, if applicable, around exons.
the specification of flanking bases, which essentially provides padding around features, implicitly allows the ability to target promoter regions, splice sites, and introns, relative to exons.
since there is also an interest in non-coding regions, we support arbitrary genomic regions as input based on chromosome coordinates.
additionally, the investigator has the option of including evolutionary conserved regions  in their target region.
we define ecrs as overlapping conserved regions between the human, mouse, and rat genomes established by a set of alignment criteria .
the ecrs for the individual species are obtained dynamically from the dcode ecr browser website  <cit> .
outputs
the output of the primer design pipeline is a set of primer pair sequences and their theoretical amplicons based on a reference sequence.
in addition, we output a summary of the critique results, such as the primer and amplicon melting temperatures and information about trace phase shifting events.
the primers can be optionally 5' tailed with m <dig>  or any other pre-defined, sequencing primers, and are ready for ordering and synthesis.
a gff file  <cit>  is produced to allow the investigator to review the spatial organization of the amplicons in juxtaposition with sequence annotation.
the gff file also highlights target regions where primer pairs could not be designed.
this allows the investigator to decide whether to proceed with sequencing, despite having areas that lack coverage, or to initiate a more customized primer design, targeting the uncovered regions.
figure  <dig> shows an example of a rendered gff file, in the form of an automatically generated adobe pdf file, incorporating the ensembl annotation and designed amplicons.
in addition to the visual summary, a statistical overview is also generated which reports the number of amplicons that were designed and their average, minimum, maximum and standard deviation of lengths.
for high-throughput primer designs, we have project level tools which summarize the total and percentage of targeted base pairs that were covered.
algorithm
the primer design pipeline is composed of two key software components: the coverage manager  and the primer critiquor .
the coverage manager  is responsible for generating a dynamic amplicon tiling across the supplied target regions.
its parameters include the target, minimum, and maximum amplicon size, the minimum number of overlapping base pairs required between amplicons, and the depth of redundant amplicon coverage required.
overlap between amplicons is typically warranted when using sanger sequencing because approximately  <dig> bp from the 5'end of the read may be unusable due to its low quality.
as a result, an overlap of at least  <dig> bp is required between adjacent amplicons to ensure that the target region can be bi-directionally covered with high-quality sequence.
the cm employs primer <dig>  <cit>  to generate a list of primer pair candidates.
the cm will then invoke the primer critiquor  to determine whether or not to accept a primer pair candidate.
if multiple target regions are proximally located, the cm will subsume the regions with one amplicon.
if the target region is smaller than an amplicon, the cm attempts to center that region within the amplicon.
otherwise, the cm will attempt to dynamically generate a tiling from the 5' to 3' end across the target region.
after the first amplicon has been selected, every proceeding amplicon will depend on the length and position of the prior upstream amplicon and the required minimum amplicon overlap.
in order to assay targeted regions, non-uniform tiling can be justified because of our ability to select primers which will be successful.
this greedy method for amplicon tiling is more efficient than a strictly uniformly spaced tiling methodology  <cit>  because we can maximize coverage while minimizing redundancy.
target regions are left uncovered by amplicons only if all primer pair candidates have been rejected in that region.
the primer critiquor  is responsible for determining whether a primer pair passes or fails our selection criteria.
the parameters for most of the criteria have been determined based on prior directed sequencing experiments, providing us the advantage of avoiding a "fuzzy" logic based approach to evaluating the potential success of primer pairs  <cit> .
the pc was designed to have each criterion specified in a modular fashion.
in the event that a new failure pattern is discovered, a new criterion can easily be incorporated as an additional module.
the pc may be invoked iteratively when the cm needs to test primer pairs during a primer design run, or it may be called upon independently to evaluate an arbitrary set of primer pairs generated by methods that are external to the current pipeline.
the following criteria are evaluated by the pc to determine whether to accept a given primer pair:
alternative amplification detection
alternative amplification occurs during pcr when primers bind to and amplify non-targeted regions.
this is especially problematic when using genomic dna samples.
to check for alternative amplification events, we use blastn  <cit> , with a word size of  <dig>  to search for the forward and reverse primers on a reference genome.
the blast alignments are then filtered and accepted as hits if they meet either of the following two criteria: 1) >80% overall identity and a perfect match of  <dig> bp at the 3' end or 2) a perfect match of  <dig> bp on the 3' end.
the detection of alternative amplification depends on how well characterized the reference genome is or how similar the reference genome is to the target organism.
if the length of a predicted alternative amplicon exceeds a parameterized threshold, its amplification is not considered viable.
we determined this threshold to be  <dig> bp by comparing post-pcr gel results against the computationally predicted alternative products.
we observed multiple instances which indicated that computationally predicted alternative products were false positives when their lengths exceeded this  <dig> bp threshold.
we consider this to be a conservative threshold because we also observed instances where computationally predicted alternative amplicons shorter than  <dig> bp were not observed in the gel.
in these cases, the primer binding for the targeted amplicon was more thermodynamically stable and, therefore, out-competed the alternative amplicons.
a primer pair is rejected if that pair yields at least one viable alternative amplification product.
due to the interest in diseases associated with non-coding regions, such as introns, regulatory regions, as well as intergenic regions, it is important to support the ability to target genomic regions rich in repetitive elements.
we do not mask repeat regions on the reference genome before searching for primer binding sites, since it is possible for a primer pair to be specific even when one or both of the primers binds to a repetitive region.
we employ the following heuristic to expedite alternative amplification detection.
if the primer pair can computationally amplify a product in the repeat library, then the primer pair is immediately failed.
if either primer in the pair can bind to a repetitive element shorter than  <dig> bp, then the primer pair is also eliminated.
primers that bind to longer repeats are passed on to the full genome search.
this heuristic allows primers to flank short repeats while forcing primer design to tackle long repeats where multiple tiled amplicons are necessary.
trace phase shifting event detection
a trace phase shifting event  is a dna sequence motif that may cause sanger dna sequencing of pcr amplified products to exhibit mixed peaks 3' of the motif in the chromatogram.
this results from a population of amplicons with non-homogenous lengths forming during the amplification process.
two known tpses that can be detected result from pcr stutter and insertions/deletions  .
pcr stutter occurs when the polymerase slips during the amplification of mono- or di-nucleotide repeats of a certain length  <cit> .
instead of the amplification resulting in a uniform population of template amplicon copies, a population of amplicons with mixed lengths, varying at the stutter event site, is generated.
we are able to detect sequence potentially prone to slippage by searching for stretches of mono- or di-nucleotide sequence that exceed a specified length threshold.
for mono-nucleotide repeats, a stutter event is expected to occur at a length of  <dig> or greater, n≥  <dig>  for di-nucleotide repeats, stutter events occur when greater than  <dig> pairs are detected ≥  <dig>  where n <dig> ≠ n <dig>
if the donor dna sample genotype contains a heterozygous indel in the amplicon, a tpse will also occur.
this results from the amplification process generating a population of two different sequence lengths, for example, one with the indel and one without.
we incorporate the indel polymorphism information provided by dbsnp to identify loci that might potentially cause a tpse.
we have also developed digital signal processing  software to identify and deconvolve the phase shifted trace files, allowing us to split indels into two synthetic trace files  and to deconvolve the sub-peaks from traces where pcr stutter has occurred .
as a result of our ability to perform these advanced trace manipulation techniques, the primer design process is permitted to accept a single tpse per amplicon.
multiple tpses are split between separate amplicons if their proximity to each other is greater than the specified required overlap between amplicons.
by using this combination of dsp and primer design techniques it is possible to achieve high-quality bi-directional coverage, even in the presence of pcr stutter and indel variants.
figures  <dig> and  <dig> demonstrate the benefits of applying our dsp algorithms to resolving pcr stutter and splitting heterozygous indels.
the figure  <dig> is an example of how a modified multi-path acoustic echo attenuation filter was used to remove the obfuscated signal resulting from a  <dig> bp poly-a stutter.
the filtered output is able to reveal an a/t variation, that opposite strand sequencing supports, which was indiscernible in the original trace.
in figure  <dig> is an example of the results of splitting a trace file with a heterozygous indel into its two allelic forms.
this was accomplished both by discrete mixed base calling and the de-mixing of the signal using a peak model, as well as by using a correlation-based approach in dna signal space.
an insertion of agaaa can be detected in the "long" allele, with respect to the "short" allele, after the dsp analysis has been performed.
stem loop interference detection
stem loop interference at the primer binding site is detected with the help of emboss's palindrome algorithm  <cit> .
thresholds for maximum loop size and minimum stem length were determined based on prior primer pair failures.
empirically derived secondary structure predictions  <cit>  also support the inferred values.
if the loop size, as determined by the distance between paired palindrome elements, is less than  <dig> base pairs and the stem length, as determined by the alignment of the palindrome elements, exceeds  <dig> hydrogen bonds, then stem loop interference is predicted to occur in the primer binding site .
melting temperature constraints
the sequence melting temperature is computed using the same formulas as those used in primer <dig>  for primer length sequences, we use the nearest neighbor thermodynamic calculations  <cit> .
to compute the melting temperature of amplicon length sequences, we follow the formula presented by sambrook et al.
<cit>  which is based on measurements performed by bolton et al.
<cit> .
we have discovered that, by augmenting the application of the formula slightly, we can improve the ability to further discern potential primer failures.
note that the formula: tm =  <dig>  +  <dig> ) +  <dig>  ×  – 600/length, uses a statistical average for the %gc content of the input sequence.
the weakness of this averaging is that the amplicon may have localized islands of high gc content that may prohibit complete melting of the amplicon at the predicted tm.
our simple modification of this formula's application involves computing the standard deviation of the sequence melting temperature by using a sliding window so that we have multiple %gc samples along the length of the amplicon.
this standard deviation is then used in the tm formula to estimate the melting temperature variation one standard deviation above the average.
if this augmented amplicon melting temperature exceeds our maximum melting temperature cutoff, then the primer pair is rejected.
donor specific variation detection
donor specific variation in a primer binding site may interfere with primer annealing and decrease pcr efficiency.
this can lead to differential success rates across dna samples or mask mutations as a symptom of allele-specific amplification  <cit> .
since we are interested in sequencing these polymorphisms instead of inferring donor variation based on pcr success, we design primers to avoid regions of known variability based on information continually updated in dbsnp.
primer dimer detection
primer dimer is calculated using a dynamic programming algorithm.
although primer <dig> performs self and cross primer dimer calculations before generating its list of primer pair candidates, we augment that test by appending the sequencing primer, if necessary, onto the pcr primer before computing primer dimer.
this detects internal primer dimers that could have been introduced by the tailing of the pcr primers with the sequencing primers.
pcr amplification and sequencing protocol
primer design is performed with the expectation that only two protocols will be required to sequence the majority of target regions.
the "standard" protocol best amplifies regions with a percent gc content centered around 50%.
for regions with higher melting temperatures, the "high gc" protocol is used.
depending on the target regions, the ratio of amplicons that are designed for each protocol may vary.
target regions are computationally split into two sets of regions, one for each protocol.
primer design is then executed on each set of regions with parameters calibrated for each protocol.
the "standard" and "high gc" protocols produce amplicons with optimal lengths of  <dig> bp and  <dig> bp, respectively.
we generate amplicons shorter than the length of a sequencing read in order to achieve full bi-directional coverage.
this coverage redundancy increases the accuracy of downstream computational variation detection, therefore, decreasing the amount of manual review labor necessary to the curate indeterminate base calls.
this reduction of labor costs more than offsets the additional costs associated with bi-directional sequencing.
"standard" gc pcr was set up in  <dig> μl reactions containing each forward and reverse pcr primer at  <dig>  μm , 1× geneamp® pcr gold buffer ,  <dig> mm mgcl <dig> ,  <dig>  units amplitaq® gold dna polymerase ,  <dig>  mm mixed dntps , 8% glycerol , molecular biology grade water  and  <dig> ng of genomic dna.
"high gc" regions were amplified in a  <dig> μl reaction containing each forward and reverse pcr primer at  <dig>  μm, 1× geneamp® pcr gold buffer,  <dig> mm mgcl <dig>   <dig>  units of amplitaq® gold dna polymerase,  <dig>  mm mixed dntps, 8% glycerol, molecular biology grade water and  <dig> ng of genomic dna.
amplification was done in 384-well format on dual 384-well geneamp® pcr system  <dig> thermal cyclers  with the following cycling conditions:  <dig> cycle of 96°c ;  <dig> cycles of 94°c , 60°c , and 72°c  and  <dig> cycle of 72°c .
after amplification, unincorporated dntps were dephosphorylated and excess primers were removed from "standard" products using a sap/exo i reaction containing  <dig>  units of shrimp alkaline phosphatase ,  <dig>  unit of exonuclease i  and molecular biology grade water.
"high gc" products were treated with  <dig>  units of shrimp alkaline phosphatase,  <dig>  units of exonuclease i and molecular biology grade water.
the sap/exo clean up was carried out in a total volume of  <dig> μl, in 384-well plate format.
the cycling program for "standard" products was:  <dig> cycle of 37°c  and  <dig> cycle of 47°c  and 72°c ; for "high gc" products it was:  <dig> cycle of 37°c  and  <dig> cycle of 50°c  and 72°c .
all pcr products were sequenced in  <dig>  μl reactions containing  <dig>  μm m <dig> sequencing primer ,  <dig> × sequencing buffer with sucrose, molecular biology grade water and  <dig>  μl of bigdye® terminator v <dig>  ready reaction mix .
sequencing reactions were carried out in 384-well format and thermal cycled using dual 384-well geneamp® pcr system 9700s using the following sequencing program:  <dig> cycle of 96°c  and  <dig> cycles of 96°c , 50°c , and 60°c .
reaction products were precipitated with a sodium acetate/ethanol mixture followed by a final ethanol rinse and dried at room temperature.
biomek fx robots  and pixsys  <dig> nanoliter liquid handling systems  were employed in setting up pcr, pcr clean up, sequencing reactions and precipitation.
resuspended sequencing reactions were sequenced on 3730xl dna analyzers .
immediately before sequencing, the precipitated sequencing reactions were resuspended in  <dig> μl of  <dig>  mm edta  using a μfill™  and shaken using a microplate shaker .
we used either dna samples supplied by the investigators or obtained from the coriell institute to sequence the target regions for all projects.
coriell samples were acquired from the coriell cell repositories  <cit>  and were employed as a readily accessible source of anonymous dna samples approved for biomedical research.
we have also used coriell dna samples for both positive controls and for the preliminary testing of primers across an ethnically diverse dna panel.
since primers are selected based on a reference genome, preliminary testing is useful to receive feedback on the possibility of unexpected primer failures due to variations yet undocumented in dbsnp.
RESULTS
the primer design pipeline was used to generate primer sets targeting  <dig> , <dig> bp of human genome sequence across  <dig> genes in  <dig> separate jcvi projects.
jcvi projects # <dig> and # <dig> consisted of  <dig> genes where, for each gene, the regions targeted for primer design consisted of exons, intronic ecrs and  <dig> kbp upstream and downstream of the first and last exon, respectively.
jcvi project # <dig> targeted one gene in full and included ~ <dig>  kbp upstream of the first exon and  <dig> bp downstream of the last exon.
for jcvi project # <dig>  which consisted of two genes, primers were designed to target only the promoter region of one gene and only the coding region of the other gene.
primers pairs for the last project, jcvi # <dig>  were designed to target only the exons and splice sites of  <dig> genes.
primer design is a massively parallel process.
since noncontiguous target regions are independent, primer pairs for each individual target region can be designed on separate processors.
in  <dig> and  <dig>  our production pipeline generated  <dig>  independent jobs, designing primer pairs across  <dig> , <dig> bp of target region.
we utilized  <dig> compute nodes, each of which were comprised of four  <dig>  ghz intel® xeon™ processors with  <dig> kb l <dig> caches and  <dig> gb of shared physical memory.
overall, 42% of the jobs completed in less than  <dig> hour, with a median run time of  <dig>  hours .
because of the variability of repeat density and/or other genomic characteristics that tend to increase run time, we observed a positive skew.
any job exceeding  <dig> hours was typically terminated.
for the analysis presented here, only  <dig> jobs were halted after exceeding this limit.
jobs that require these longer run times are typically a result of target regions densely flanked with repetitive sequence.
since the tiling algorithm will not abandon a portion of the target region until a successful candidate primer pair has been found or until all candidates have been assessed, a region that is impossible to select successful primers for may need to undergo the slow primer specificity search for thousands of unsuccessful primer pairs.
discussion
we have developed an automated pcr primer design pipeline that incorporates novel computational methods for improving the prediction of primer pair sequencing success.
these computational methods have been parameterized and finely tuned based on the results of our high-throughput directed sequencing experiments, using genomic dna with sanger sequencing.
the fine tuning of primer pair selection has lead to our high sequencing success rate and our ability to hone our pipeline for the optimal coverage of targeted regions.
the success rate of any primer design pipeline should not only include the amplification success rate of the designed primers, but should also take into account the target regions that could have been covered, but were missed due to overly stringent design criteria.
since we have laboratory data to support the parameters of primer design, we were able to more precisely predict potential amplification errors, thus improving our success rate without unnecessarily compromising target region coverage.
it is commonly observed that, as the density of repetitive elements and low complexity sequence increase in intronic and intergenic regions, the ability to design successful primers decreases, therefore reducing coverage.
as a result, it tends to be easier to design primers for exonic regions.
our ability to visualize the designed amplicons aligned with sequence annotation helps to attribute uncovered regions to genomic landmarks, such as repeats .
this can minimize the concern of leaving regions uncovered.
our primer design pipeline has achieved an overall target region coverage rate of 95% .
the exonic subset of this target region has a coverage rate of 97% .
compared to existing publications detailing exon coverage results, our corresponding coverage rate is on par with the best of other primer design software  <cit> .
the success rate of designed primers will also vary depending on the chosen criteria for success.
for the results presented, we consider a primer pair successful if the expected sequencing result is generated from both sequencing directions by at least two donors in a given a dna panel.
this criterion for success is more stringent than those used by other similar tools during primer validation  <cit> .
in table  <dig>  we present the results of the recent large scale studies performed on  <dig> genes.
previously published primer pair success rates focused primarily on the directed sequencing of exons.
in order to make our results more comparable, we also calculated our success rate for amplicons associated with exons of the genes in our large scale study .
we consider an amplicon to be exonic if more than half of the amplicon overlaps with an exon or at least one exon is fully contained within the amplicon.
even with this strict criterion for success, we still achieve an overall exonic sequencing success rate of greater than 95%, which is comparable to the gel-only quantified success rate in  <cit> , or an improvement over other  <cit>  primer design studies.
the modular and parameterized implementation of the primer design pipeline has not only made the computational parallelization of the code possible, but has also allowed us to apply its componentized functionality to other projects.
in addition to targeted sanger sequencing on human subjects, we have designed successful primer pairs in the genome sequences of mouse, dog and staphylococcus aureus.
we have also successfully designed primer pairs, with this pipeline, that were subsequently employed in  <dig> pyrosequencing.
future work
the step that requires the largest amount of compute time is the alternative amplification detection, which is done with blastn.
our ongoing research suggests that our primer binding criteria may be more stringent than necessary.
currently, our criteria allows for mismatches within the first  <dig> bases of the primers' 3' end, which requires a more computationally expensive algorithm than perfect string matching.
however, crystal structures of taq dna polymerase i in a binary complex with primer and template dna  <cit>  may indicate that the polymerase only needs to interact with the first  <dig> bases on the 3' end of the primer.
stronger base pairing on the 5' end of the primer will increase the thermodynamic stability of the primer-template complex, leading to more efficient priming and a target amplicon which can out-compete alternative amplicons.
however, it is not required to initiate the dna synthesis that can also be responsible for non-specific amplification.
we have, however, recorded instances of alternative products that resulted despite mismatches within the first  <dig> bases of the primer's 3' end.
analysis of these alternative products revealed that they did not occur because of nonspecific binding of these primers, but because the sequence of the target amplicon could be found elsewhere on the genome with a percent identity greater than 96%.
we hypothesize that a "megaprimer"  <cit>  could be responsible for amplifying the non-target sequence.
megapriming occurs when the dna copy strand is not completely synthesized before melting from the target template.
the partial copy acts as a primer in the subsequent annealing step of the pcr cycle and if the megaprimer anneals to a homologous position on the genome, then this will result in the amplification of the alternative product.
additional experiments will be performed to determine whether the computational detection of megapriming coupled with an exact match primer binding search will improve the compute time of the primer specificity check while maintaining our high degree of primer specificity.
CONCLUSIONS
we have developed a high-throughput pcr primer design pipeline providing high-quality directed sequencing results while minimizing costs associated with labor and reprocessing.
the consistent performance of our laboratory protocols has made it possible to isolate and analyze failures, and to systematically increase our exonic sequencing success rate to 95%.
we can readily incorporate new primer critique tests based on iterative feedback from the laboratory, and integrate with downstream processes as a result of the modular software design.
this work bridges the gap between computational theory and laboratory practice, and discloses the novel set of computational methods that made it possible.
availability and requirements
project name: jcvi primer designer
project home page:
operating system: tested and in production on linux.
programming language: perl
license: gnu gpl
any restrictions to use by non-academics: none
authors' contributions
sm, sl, and kl conceived of the study and participated in its design.
kl wrote the software.
ab, kl, and ts analyzed data and evaluated results.
kb, tcm, db, and sf designed laboratory protocols, ran experiments, and generated data.
kl and ab wrote the manuscript.
all authors read and approved the submitted manuscript.

 directed sequencing is important for the discovery of novel polymorphisms in human dna sequences that may play a role in disease susceptibility and progression. in order to efficiently screen large numbers of target regions across a diverse panel of dna samples, a sequencing center needs to apply a high-throughput approach with the goal of minimizing cost. to achieve this, pcr primers are designed with the aim of eliminating human labor associated with amplification protocol optimization and other manual intervention caused by primer failures. recent contributions in computational pcr primer design include screening additional thermodynamic attributes and increasing the genomic annotation considered in the design process. table provides a comparison of novel features of our primer design pipeline in relation to recently published primer design software. our primer design pipeline is further enhanced by a higher degree of operational automation and the parameterization of many key attributes in its design. we have employed this pipeline for over five years at our sequencing center and its use was instrumental in the publication of mutations discovered in a set of tyrosine kinase receptor domains consisting of exons, across tumor samples in glioblastomas. since over mb of sequence have been targeted with a primer design target coverage rate exceeding 95%. our directed sequencing pipeline has achieved a sequencing success rate of over 90% and we have detected over potential variations in genes, of which were not previously recorded in ncbi's dbsnp. ∗ feature supported - feature not supported ‡ tiling pattern is fixed † searches for the full length of the primer, instead of just the critical region on the 3' end implementation computational pcr primer design software in the method described here, we constrain the pcr primer design software to generate a set of primers with a high likelihood of success using one of two standardized amplification protocols. the high-throughput computational pcr primer design pipeline is one of the first computational stages we employ as part of our directed sequencing pipeline. the computational pcr primer design software is available on the sourceforge.net jcvi primer designer website. the tar file of the latest release is available by navigating through the'download' tab. it can be installed and run on any computer with a unix operating system. inputs target regions are specified by the investigator through a variety of identifiers. by using the feature rich ensembl perl application programming interface, the following types of identifiers are supported as input: gene, transcript, hugo, exon, and dbsnp. since the reference genome assembly and annotation may be updated, using identifiers allows us to track features of interest even if their absolute genomic coordinates have changed over time. in addition, by using these identifiers as anchors, we are able to augment target regions with a specified number of flanking bases upstream and downstream of the region, and, if applicable, around exons. the specification of flanking bases, which essentially provides padding around features, implicitly allows the ability to target promoter regions, splice sites, and introns, relative to exons. since there is also an interest in non-coding regions, we support arbitrary genomic regions as input based on chromosome coordinates. additionally, the investigator has the option of including evolutionary conserved regions in their target region. we define ecrs as overlapping conserved regions between the human, mouse, and rat genomes established by a set of alignment criteria. the ecrs for the individual species are obtained dynamically from the dcode ecr browser website. outputs the output of the primer design pipeline is a set of primer pair sequences and their theoretical amplicons based on a reference sequence. in addition, we output a summary of the critique results, such as the primer and amplicon melting temperatures and information about trace phase shifting events. the primers can be optionally 5' tailed with m or any other pre-defined, sequencing primers, and are ready for ordering and synthesis. a gff file is produced to allow the investigator to review the spatial organization of the amplicons in juxtaposition with sequence annotation. the gff file also highlights target regions where primer pairs could not be designed. this allows the investigator to decide whether to proceed with sequencing, despite having areas that lack coverage, or to initiate a more customized primer design, targeting the uncovered regions. figure shows an example of a rendered gff file, in the form of an automatically generated adobe pdf file, incorporating the ensembl annotation and designed amplicons. in addition to the visual summary, a statistical overview is also generated which reports the number of amplicons that were designed and their average, minimum, maximum and standard deviation of lengths. for high-throughput primer designs, we have project level tools which summarize the total and percentage of targeted base pairs that were covered. algorithm the primer design pipeline is composed of two key software components: the coverage manager and the primer critiquor. the coverage manager is responsible for generating a dynamic amplicon tiling across the supplied target regions. its parameters include the target, minimum, and maximum amplicon size, the minimum number of overlapping base pairs required between amplicons, and the depth of redundant amplicon coverage required. overlap between amplicons is typically warranted when using sanger sequencing because approximately bp from the 5'end of the read may be unusable due to its low quality. as a result, an overlap of at least bp is required between adjacent amplicons to ensure that the target region can be bi-directionally covered with high-quality sequence. the cm employs primer to generate a list of primer pair candidates. the cm will then invoke the primer critiquor to determine whether or not to accept a primer pair candidate. if multiple target regions are proximally located, the cm will subsume the regions with one amplicon. if the target region is smaller than an amplicon, the cm attempts to center that region within the amplicon. otherwise, the cm will attempt to dynamically generate a tiling from the 5' to 3' end across the target region. after the first amplicon has been selected, every proceeding amplicon will depend on the length and position of the prior upstream amplicon and the required minimum amplicon overlap. in order to assay targeted regions, non-uniform tiling can be justified because of our ability to select primers which will be successful. this greedy method for amplicon tiling is more efficient than a strictly uniformly spaced tiling methodology because we can maximize coverage while minimizing redundancy. target regions are left uncovered by amplicons only if all primer pair candidates have been rejected in that region. the primer critiquor is responsible for determining whether a primer pair passes or fails our selection criteria. the parameters for most of the criteria have been determined based on prior directed sequencing experiments, providing us the advantage of avoiding a "fuzzy" logic based approach to evaluating the potential success of primer pairs. the pc was designed to have each criterion specified in a modular fashion. in the event that a new failure pattern is discovered, a new criterion can easily be incorporated as an additional module. the pc may be invoked iteratively when the cm needs to test primer pairs during a primer design run, or it may be called upon independently to evaluate an arbitrary set of primer pairs generated by methods that are external to the current pipeline. the following criteria are evaluated by the pc to determine whether to accept a given primer pair: alternative amplification detection alternative amplification occurs during pcr when primers bind to and amplify non-targeted regions. this is especially problematic when using genomic dna samples. to check for alternative amplification events, we use blastn, with a word size of to search for the forward and reverse primers on a reference genome. the blast alignments are then filtered and accepted as hits if they meet either of the following two criteria: 1) >80% overall identity and a perfect match of bp at the 3' end or 2) a perfect match of bp on the 3' end. the detection of alternative amplification depends on how well characterized the reference genome is or how similar the reference genome is to the target organism. if the length of a predicted alternative amplicon exceeds a parameterized threshold, its amplification is not considered viable. we determined this threshold to be bp by comparing post-pcr gel results against the computationally predicted alternative products. we observed multiple instances which indicated that computationally predicted alternative products were false positives when their lengths exceeded this bp threshold. we consider this to be a conservative threshold because we also observed instances where computationally predicted alternative amplicons shorter than bp were not observed in the gel. in these cases, the primer binding for the targeted amplicon was more thermodynamically stable and, therefore, out-competed the alternative amplicons. a primer pair is rejected if that pair yields at least one viable alternative amplification product. due to the interest in diseases associated with non-coding regions, such as introns, regulatory regions, as well as intergenic regions, it is important to support the ability to target genomic regions rich in repetitive elements. we do not mask repeat regions on the reference genome before searching for primer binding sites, since it is possible for a primer pair to be specific even when one or both of the primers binds to a repetitive region. we employ the following heuristic to expedite alternative amplification detection. if the primer pair can computationally amplify a product in the repeat library, then the primer pair is immediately failed. if either primer in the pair can bind to a repetitive element shorter than bp, then the primer pair is also eliminated. primers that bind to longer repeats are passed on to the full genome search. this heuristic allows primers to flank short repeats while forcing primer design to tackle long repeats where multiple tiled amplicons are necessary. trace phase shifting event detection a trace phase shifting event is a dna sequence motif that may cause sanger dna sequencing of pcr amplified products to exhibit mixed peaks 3' of the motif in the chromatogram. this results from a population of amplicons with non-homogenous lengths forming during the amplification process. two known tpses that can be detected result from pcr stutter and insertions/deletions. pcr stutter occurs when the polymerase slips during the amplification of mono- or di-nucleotide repeats of a certain length. instead of the amplification resulting in a uniform population of template amplicon copies, a population of amplicons with mixed lengths, varying at the stutter event site, is generated. we are able to detect sequence potentially prone to slippage by searching for stretches of mono- or di-nucleotide sequence that exceed a specified length threshold. for mono-nucleotide repeats, a stutter event is expected to occur at a length of or greater, n≥ for di-nucleotide repeats, stutter events occur when greater than pairs are detected ≥ where n ≠ n if the donor dna sample genotype contains a heterozygous indel in the amplicon, a tpse will also occur. this results from the amplification process generating a population of two different sequence lengths, for example, one with the indel and one without. we incorporate the indel polymorphism information provided by dbsnp to identify loci that might potentially cause a tpse. we have also developed digital signal processing software to identify and deconvolve the phase shifted trace files, allowing us to split indels into two synthetic trace files and to deconvolve the sub-peaks from traces where pcr stutter has occurred. as a result of our ability to perform these advanced trace manipulation techniques, the primer design process is permitted to accept a single tpse per amplicon. multiple tpses are split between separate amplicons if their proximity to each other is greater than the specified required overlap between amplicons. by using this combination of dsp and primer design techniques it is possible to achieve high-quality bi-directional coverage, even in the presence of pcr stutter and indel variants. figures and demonstrate the benefits of applying our dsp algorithms to resolving pcr stutter and splitting heterozygous indels. the figure is an example of how a modified multi-path acoustic echo attenuation filter was used to remove the obfuscated signal resulting from a bp poly-a stutter. the filtered output is able to reveal an a/t variation, that opposite strand sequencing supports, which was indiscernible in the original trace. in figure is an example of the results of splitting a trace file with a heterozygous indel into its two allelic forms. this was accomplished both by discrete mixed base calling and the de-mixing of the signal using a peak model, as well as by using a correlation-based approach in dna signal space. an insertion of agaaa can be detected in the "long" allele, with respect to the "short" allele, after the dsp analysis has been performed. stem loop interference detection stem loop interference at the primer binding site is detected with the help of emboss's palindrome algorithm. thresholds for maximum loop size and minimum stem length were determined based on prior primer pair failures. empirically derived secondary structure predictions also support the inferred values. if the loop size, as determined by the distance between paired palindrome elements, is less than base pairs and the stem length, as determined by the alignment of the palindrome elements, exceeds hydrogen bonds, then stem loop interference is predicted to occur in the primer binding site. melting temperature constraints the sequence melting temperature is computed using the same formulas as those used in primer for primer length sequences, we use the nearest neighbor thermodynamic calculations. to compute the melting temperature of amplicon length sequences, we follow the formula presented by sambrook et al. which is based on measurements performed by bolton et al.. we have discovered that, by augmenting the application of the formula slightly, we can improve the ability to further discern potential primer failures. note that the formula: tm = + ) + × – 600/length, uses a statistical average for the %gc content of the input sequence. the weakness of this averaging is that the amplicon may have localized islands of high gc content that may prohibit complete melting of the amplicon at the predicted tm. our simple modification of this formula's application involves computing the standard deviation of the sequence melting temperature by using a sliding window so that we have multiple %gc samples along the length of the amplicon. this standard deviation is then used in the tm formula to estimate the melting temperature variation one standard deviation above the average. if this augmented amplicon melting temperature exceeds our maximum melting temperature cutoff, then the primer pair is rejected. donor specific variation detection donor specific variation in a primer binding site may interfere with primer annealing and decrease pcr efficiency. this can lead to differential success rates across dna samples or mask mutations as a symptom of allele-specific amplification. since we are interested in sequencing these polymorphisms instead of inferring donor variation based on pcr success, we design primers to avoid regions of known variability based on information continually updated in dbsnp. primer dimer detection primer dimer is calculated using a dynamic programming algorithm. although primer performs self and cross primer dimer calculations before generating its list of primer pair candidates, we augment that test by appending the sequencing primer, if necessary, onto the pcr primer before computing primer dimer. this detects internal primer dimers that could have been introduced by the tailing of the pcr primers with the sequencing primers. pcr amplification and sequencing protocol primer design is performed with the expectation that only two protocols will be required to sequence the majority of target regions. the "standard" protocol best amplifies regions with a percent gc content centered around 50%. for regions with higher melting temperatures, the "high gc" protocol is used. depending on the target regions, the ratio of amplicons that are designed for each protocol may vary. target regions are computationally split into two sets of regions, one for each protocol. primer design is then executed on each set of regions with parameters calibrated for each protocol. the "standard" and "high gc" protocols produce amplicons with optimal lengths of bp and bp, respectively. we generate amplicons shorter than the length of a sequencing read in order to achieve full bi-directional coverage. this coverage redundancy increases the accuracy of downstream computational variation detection, therefore, decreasing the amount of manual review labor necessary to the curate indeterminate base calls. this reduction of labor costs more than offsets the additional costs associated with bi-directional sequencing. "standard" gc pcr was set up in μl reactions containing each forward and reverse pcr primer at μm, 1× geneamp® pcr gold buffer, mm mgcl, units amplitaq® gold dna polymerase, mm mixed dntps, 8% glycerol, molecular biology grade water and ng of genomic dna. "high gc" regions were amplified in a μl reaction containing each forward and reverse pcr primer at μm, 1× geneamp® pcr gold buffer, mm mgcl units of amplitaq® gold dna polymerase, mm mixed dntps, 8% glycerol, molecular biology grade water and ng of genomic dna. amplification was done in 384-well format on dual 384-well geneamp® pcr system thermal cyclers with the following cycling conditions: cycle of 96°c; cycles of 94°c, 60°c, and 72°c and cycle of 72°c. after amplification, unincorporated dntps were dephosphorylated and excess primers were removed from "standard" products using a sap/exo i reaction containing units of shrimp alkaline phosphatase, unit of exonuclease i and molecular biology grade water. "high gc" products were treated with units of shrimp alkaline phosphatase, units of exonuclease i and molecular biology grade water. the sap/exo clean up was carried out in a total volume of μl, in 384-well plate format. the cycling program for "standard" products was: cycle of 37°c and cycle of 47°c and 72°c; for "high gc" products it was: cycle of 37°c and cycle of 50°c and 72°c. all pcr products were sequenced in μl reactions containing μm m sequencing primer, × sequencing buffer with sucrose, molecular biology grade water and μl of bigdye® terminator v ready reaction mix. sequencing reactions were carried out in 384-well format and thermal cycled using dual 384-well geneamp® pcr system 9700s using the following sequencing program: cycle of 96°c and cycles of 96°c, 50°c, and 60°c. reaction products were precipitated with a sodium acetate/ethanol mixture followed by a final ethanol rinse and dried at room temperature. biomek fx robots and pixsys nanoliter liquid handling systems were employed in setting up pcr, pcr clean up, sequencing reactions and precipitation. resuspended sequencing reactions were sequenced on 3730xl dna analyzers. immediately before sequencing, the precipitated sequencing reactions were resuspended in μl of mm edta using a μfill™ and shaken using a microplate shaker. we used either dna samples supplied by the investigators or obtained from the coriell institute to sequence the target regions for all projects. coriell samples were acquired from the coriell cell repositories and were employed as a readily accessible source of anonymous dna samples approved for biomedical research. we have also used coriell dna samples for both positive controls and for the preliminary testing of primers across an ethnically diverse dna panel. since primers are selected based on a reference genome, preliminary testing is useful to receive feedback on the possibility of unexpected primer failures due to variations yet undocumented in dbsnp.  the primer design pipeline was used to generate primer sets targeting, bp of human genome sequence across genes in separate jcvi projects. jcvi projects # and # consisted of genes where, for each gene, the regions targeted for primer design consisted of exons, intronic ecrs and kbp upstream and downstream of the first and last exon, respectively. jcvi project # targeted one gene in full and included ~ kbp upstream of the first exon and bp downstream of the last exon. for jcvi project # which consisted of two genes, primers were designed to target only the promoter region of one gene and only the coding region of the other gene. primers pairs for the last project, jcvi # were designed to target only the exons and splice sites of genes. primer design is a massively parallel process. since noncontiguous target regions are independent, primer pairs for each individual target region can be designed on separate processors. in and our production pipeline generated independent jobs, designing primer pairs across, bp of target region. we utilized compute nodes, each of which were comprised of four ghz intel® xeon™ processors with kb l caches and gb of shared physical memory. overall, 42% of the jobs completed in less than hour, with a median run time of hours. because of the variability of repeat density and/or other genomic characteristics that tend to increase run time, we observed a positive skew. any job exceeding hours was typically terminated. for the analysis presented here, only jobs were halted after exceeding this limit. jobs that require these longer run times are typically a result of target regions densely flanked with repetitive sequence. since the tiling algorithm will not abandon a portion of the target region until a successful candidate primer pair has been found or until all candidates have been assessed, a region that is impossible to select successful primers for may need to undergo the slow primer specificity search for thousands of unsuccessful primer pairs. discussion we have developed an automated pcr primer design pipeline that incorporates novel computational methods for improving the prediction of primer pair sequencing success. these computational methods have been parameterized and finely tuned based on the results of our high-throughput directed sequencing experiments, using genomic dna with sanger sequencing. the fine tuning of primer pair selection has lead to our high sequencing success rate and our ability to hone our pipeline for the optimal coverage of targeted regions. the success rate of any primer design pipeline should not only include the amplification success rate of the designed primers, but should also take into account the target regions that could have been covered, but were missed due to overly stringent design criteria. since we have laboratory data to support the parameters of primer design, we were able to more precisely predict potential amplification errors, thus improving our success rate without unnecessarily compromising target region coverage. it is commonly observed that, as the density of repetitive elements and low complexity sequence increase in intronic and intergenic regions, the ability to design successful primers decreases, therefore reducing coverage. as a result, it tends to be easier to design primers for exonic regions. our ability to visualize the designed amplicons aligned with sequence annotation helps to attribute uncovered regions to genomic landmarks, such as repeats. this can minimize the concern of leaving regions uncovered. our primer design pipeline has achieved an overall target region coverage rate of 95%. the exonic subset of this target region has a coverage rate of 97%. compared to existing publications detailing exon coverage results, our corresponding coverage rate is on par with the best of other primer design software. the success rate of designed primers will also vary depending on the chosen criteria for success. for the results presented, we consider a primer pair successful if the expected sequencing result is generated from both sequencing directions by at least two donors in a given a dna panel. this criterion for success is more stringent than those used by other similar tools during primer validation. in table we present the results of the recent large scale studies performed on genes. previously published primer pair success rates focused primarily on the directed sequencing of exons. in order to make our results more comparable, we also calculated our success rate for amplicons associated with exons of the genes in our large scale study. we consider an amplicon to be exonic if more than half of the amplicon overlaps with an exon or at least one exon is fully contained within the amplicon. even with this strict criterion for success, we still achieve an overall exonic sequencing success rate of greater than 95%, which is comparable to the gel-only quantified success rate in, or an improvement over other primer design studies. the modular and parameterized implementation of the primer design pipeline has not only made the computational parallelization of the code possible, but has also allowed us to apply its componentized functionality to other projects. in addition to targeted sanger sequencing on human subjects, we have designed successful primer pairs in the genome sequences of mouse, dog and staphylococcus aureus. we have also successfully designed primer pairs, with this pipeline, that were subsequently employed in pyrosequencing. future work the step that requires the largest amount of compute time is the alternative amplification detection, which is done with blastn. our ongoing research suggests that our primer binding criteria may be more stringent than necessary. currently, our criteria allows for mismatches within the first bases of the primers' 3' end, which requires a more computationally expensive algorithm than perfect string matching. however, crystal structures of taq dna polymerase i in a binary complex with primer and template dna may indicate that the polymerase only needs to interact with the first bases on the 3' end of the primer. stronger base pairing on the 5' end of the primer will increase the thermodynamic stability of the primer-template complex, leading to more efficient priming and a target amplicon which can out-compete alternative amplicons. however, it is not required to initiate the dna synthesis that can also be responsible for non-specific amplification. we have, however, recorded instances of alternative products that resulted despite mismatches within the first bases of the primer's 3' end. analysis of these alternative products revealed that they did not occur because of nonspecific binding of these primers, but because the sequence of the target amplicon could be found elsewhere on the genome with a percent identity greater than 96%. we hypothesize that a "megaprimer" could be responsible for amplifying the non-target sequence. megapriming occurs when the dna copy strand is not completely synthesized before melting from the target template. the partial copy acts as a primer in the subsequent annealing step of the pcr cycle and if the megaprimer anneals to a homologous position on the genome, then this will result in the amplification of the alternative product. additional experiments will be performed to determine whether the computational detection of megapriming coupled with an exact match primer binding search will improve the compute time of the primer specificity check while maintaining our high degree of primer specificity.  we have developed a high-throughput pcr primer design pipeline providing high-quality directed sequencing results while minimizing costs associated with labor and reprocessing. the consistent performance of our laboratory protocols has made it possible to isolate and analyze failures, and to systematically increase our exonic sequencing success rate to 95%. we can readily incorporate new primer critique tests based on iterative feedback from the laboratory, and integrate with downstream processes as a result of the modular software design. this work bridges the gap between computational theory and laboratory practice, and discloses the novel set of computational methods that made it possible. availability and requirements project name: jcvi primer designer project home page: operating system: tested and in production on linux. programming language: perl license: gnu gpl any restrictions to use by non-academics: none authors' contributions sm, sl, and kl conceived of the study and participated in its design. kl wrote the software. ab, kl, and ts analyzed data and evaluated results. kb, tcm, db, and sf designed laboratory protocols, ran experiments, and generated data. kl and ab wrote the manuscript. all authors read and approved the submitted manuscript. 
we reveal the laboratory protocols and their associated, empirically determined computational parameters, as well as describe the novel computational methods which may benefit others in future primer design research. as the first step in pcr directed sequencing, effective pcr primer design is crucial for obtaining high-quality sequence data for target regions. since current computational primer design tools are not fully tuned with stable underlying laboratory protocols, researchers may still be forced to iteratively optimize protocols for failed amplifications after the primers have been ordered. this inefficient approach to primer design is further intensified in a high-throughput laboratory, where hundreds of genes may be targeted in one experiment. we have developed a fully integrated computational pcr primer design pipeline that plays a key role in our high-throughput directed sequencing pipeline. the high-throughput pcr primer design pipeline has been very successful in providing the basis for high-quality directed sequencing results and for minimizing costs associated with labor and reprocessing. the modular architecture of the primer design software has made it possible to readily integrate additional primer critique tests based on iterative feedback from the laboratory. as a result, the primer design software, coupled with the laboratory protocols, serves as a powerful tool for low and high-throughput primer design to enable successful directed sequencing. investigators may specify target regions defined through a rich set of descriptors, such as ensembl accessions and arbitrary genomic coordinates. primer pairs are then selected computationally to produce a minimal amplicon set capable of tiling across the specified target regions. polymerase chain reaction is used in directed sequencing for the discovery of novel polymorphisms. as part of the tiling process, primer pairs are computationally screened to meet the criteria for success with one of two pcr amplification protocols. in the process of improving our sequencing success rate, which currently exceeds 95% for exons, we have discovered novel and accurate computational methods capable of identifying primers that may lead to pcr failures.
 polymerase chain reaction is used in directed sequencing for the discovery of novel polymorphisms. as the first step in pcr directed sequencing, effective pcr primer design is crucial for obtaining high-quality sequence data for target regions. since current computational primer design tools are not fully tuned with stable underlying laboratory protocols, researchers may still be forced to iteratively optimize protocols for failed amplifications after the primers have been ordered. furthermore, potentially identifiable factors which contribute to pcr failures have yet to be elucidated. this inefficient approach to primer design is further intensified in a high-throughput laboratory, where hundreds of genes may be targeted in one experiment.  we have developed a fully integrated computational pcr primer design pipeline that plays a key role in our high-throughput directed sequencing pipeline. investigators may specify target regions defined through a rich set of descriptors, such as ensembl accessions and arbitrary genomic coordinates. primer pairs are then selected computationally to produce a minimal amplicon set capable of tiling across the specified target regions. as part of the tiling process, primer pairs are computationally screened to meet the criteria for success with one of two pcr amplification protocols. in the process of improving our sequencing success rate, which currently exceeds 95% for exons, we have discovered novel and accurate computational methods capable of identifying primers that may lead to pcr failures. we reveal the laboratory protocols and their associated, empirically determined computational parameters, as well as describe the novel computational methods which may benefit others in future primer design research.  the high-throughput pcr primer design pipeline has been very successful in providing the basis for high-quality directed sequencing results and for minimizing costs associated with labor and reprocessing. the modular architecture of the primer design software has made it possible to readily integrate additional primer critique tests based on iterative feedback from the laboratory. as a result, the primer design software, coupled with the laboratory protocols, serves as a powerful tool for low and high-throughput primer design to enable successful directed sequencing. 