MCMC-ODPR: Primer design optimization using Markov Chain Monte Carlo sampling
© Kitchen et al.; licensee BioMed Central Ltd. 2012
Received: 26 September 2011
Accepted: 18 October 2012
Published: 5 November 2012
Skip to main content
© Kitchen et al.; licensee BioMed Central Ltd. 2012
Received: 26 September 2011
Accepted: 18 October 2012
Published: 5 November 2012
Next generation sequencing technologies often require numerous primer designs that require good target coverage that can be financially costly. We aimed to develop a system that would implement primer reuse to design degenerate primers that could be designed around SNPs, thus find the fewest necessary primers and the lowest cost whilst maintaining an acceptable coverage and provide a cost effective solution. We have implemented Metropolis-Hastings Markov Chain Monte Carlo for optimizing primer reuse. We call it the Markov Chain Monte Carlo Optimized Degenerate Primer Reuse (MCMC-ODPR) algorithm.
After repeating the program 1020 times to assess the variance, an average of 17.14% fewer primers were found to be necessary using MCMC-ODPR for an equivalent coverage without implementing primer reuse. The algorithm was able to reuse primers up to five times. We compared MCMC-ODPR with single sequence primer design programs Primer3 and Primer-BLAST and achieved a lower primer cost per amplicon base covered of 0.21 and 0.19 and 0.18 primer nucleotides on three separate gene sequences, respectively. With multiple sequences, MCMC-ODPR achieved a lower cost per base covered of 0.19 than programs BatchPrimer3 and PAMPS, which achieved 0.25 and 0.64 primer nucleotides, respectively.
MCMC-ODPR is a useful tool for designing primers at various melting temperatures at good target coverage. By combining degeneracy with optimal primer reuse the user may increase coverage of sequences amplified by the designed primers at significantly lower costs. Our analyses showed that overall MCMC-ODPR outperformed the other primer-design programs in our study in terms of cost per covered base.
Next generation sequencing technology has led to the frequent application of deep sequencing projects, and the use of systems that require a large number of oligonuceotide primers for PCR. The design of a high number of primers is a challenge logistically both in terms of achieving good coverage of target regions and in terms of cost. Although there are a number of primer design programs available, utilizing them for high throughput design can be difficult and financially costly. We aimed to produce a system in which a large number of primers could be designed cost effectively by using the fewest necessary primers, hence the lowest cost, at multiple priming sites where possible whilst maintaining an acceptable level of coverage, and avoiding degeneracy in amplicon targets which overlap in the same regions.
In designing our program we compared our approaches and performance with several other available programs; Primer3 , Batchprimer3 , Primer-BLAST  and PAMPS . The core algorithm for the first three of these programs is that of Primer3. The Primer3 algorithm takes into account the primer size, melting temperature (Tm), GC content, and concentration of monovalent and divalent cations within the PCR reaction mixture, a selection of salt correction formulae and different parameters for simulating the thermodynamics of primer hybridization. Potential primers are then checked by using a mispriming repeat library from the human, rodent or Drosophila genomes, allowing interspersed repeats or other sequence regions to be avoided as primer annealing locations. Primer-BLAST utilizes the Primer3 algorithm and the BLAST local alignment search tool  to ensure only unique primer pairs are selected, thus preventing primers becoming designed around undesired targets such as introns. These two programs output a range of primer pair possibilities for single DNA sequences, but they are not designed for high throughput primer design. The need for high throughput primer design was recognized by You et al. , who produced BatchPrimer3 in which multiple sequences can be input for primer selection, but only one primer pair per input sequence is produced.
Minimizing the cost of a primer design can be achieved by (i) designing degenerate primers able to anneal to a number of related target sequences and (ii) implementing primer reuse utilizing primers that bind to conserved loci that are repeated. Although degeneracy allows for amplification of greater numbers of related sequences, the more degenerate primers are the less specific amplicons will be. Therefore achieving an optimal degree of degeneracy is important to obtaining a suitable trade-off between the number of related sequences amplified and the specificity of these amplicons. A number of variants exist for tackling this problem and achieving a good trade-off for the specificity and sensitivity. The Maximum Coverage Degenerate Primer Design (MCDPD) approach as used in HYDEN  tries to identify a primer of length l and a maximum degeneracy d max that covers a maximum number of sequences, each of length l. The Minimum Degeneracy Degenerate Primer Design (MDDPD) attempts to find a primer of length l and a minimum degeneracy d min that can cover all input sequences with a length equal to or greater than l. The Minimum Primers Degenerate Primer Design (MPDPD) attempts to find the fewest number of primers of length l and a maximum degeneracy of d max for a set of sequences, so that each sequence is covered by at least one primer. Whereas this approach has the constraint that all input sequences have the same length as the primers and may be inadequate in practice, the Multiple Degenerate Primer Design (MDPD) allows the input sequences to have different lengths of greater than l min and attempts to identify primers of length at least l min and degeneracy d max , allowing each sequence to be covered by at least one primer. This was the approach taken by MIPS  and PAMPS . PAMPS is a heuristic, high throughput algorithm which designs degenerate primers through a process of consecutive ad hoc pairwise alignments . This program has been shown to outperform other degenerate primer design systems such as HYDEN and MIPS in terms of computational time.
Algorithms for implementing primer reuse have also been developed: Doi and Imai have described a heuristic algorithm for greedy primer design within multiplex PCR, which attempts to minimize the cost of primers required for multiplex PCR and SNP genotyping . MuPlex is another heuristic algorithm designed for multiplex PCR, which uses a graph based approach to assign the largest number of non-conflicting primers into the fewest ‘cliques’ that can be assigned to multiplex PCR tubes . Lui and Carson have utilized a simulated annealing optimization to maximize primer reuse, which exhaustively searches primer space and aims to converge upon the optimal cost solution . Despite individual cost benefits from either optimization of primer reuse, or automated design of degenerate primers, combining the two techniques is likely to offer additional cost advantages.
Optimizing primer design to make use of degeneracy and multiplexing has been referred to as the Multiple Degenerate Primer Selection Problem (MDPSP), and variants have been shown to be NP-complete. Previous approaches to MDPSP, such as those undertaken by Balla et al  have shown that primer coverage and cost can be improved through approximate (heuristic) greedy algorithms. Jabado et al. provide a heuristic algorithm for degeneracy, Greene SCPrimer . In this method phylogenetic trees are constructed from multiple sequence alignments to identify candidate primers, which are used by a greedy set covering problem (SCP) solving algorithm to determine the minimum set of degenerate primers that may amplify all members of the alignment, so combining degeneracy with primer reuse. Although heuristic approaches generally outperform global optimizations in computation time, the reverse can potentially be true in quality of output. Given that optimization of large multiplexed primer design is not generally time-critical, a global optimization approach seems appropriate.
In order to improve on greedy approaches to MDPSP we present here an algorithm that takes a Markov Chain Monte Carlo (MCMC) approach, which allows sampling through primer parameter space using a probability distribution of acceptance of iterative primer designs. Primers are weighted according to their degree of reuse provided their degeneracy is kept below a user-defined threshold. We have implemented a Metropolis-Hastings algorithm, in which new proposals (e.g. the cost of a primer design) are accepted if they provide a more optimal solution to the current proposal, with the system tending to revert probabilistically, to the current state if the new proposal is more costly. We call the algorithm the Markov Chain Monte Carlo Optimized Degenerate Primer Reuse (MCMC-ODPR) algorithm. We show that the MCMC-ODPR program outperforms Primer3, Primer-BLAST, BatchPrimer3, and PAMPS in terms of cost and in terms of sequence coverage.
The GC content of the primers must be taken into account. Due to the stronger intermolecular forces between the guanine and cytosine nucleotides, Tm is highly dependent on this, which in turn has implications for how well primers hybridize to template sequences and for specificity of amplicons.
Primers should be minimally degenerate with a limited number of degenerate nucleotides per primer.
The forward and reverse pair of primers must not allow amplification of a target that exceeds the desired target length.
The primers should ideally be used multiple times, to keep the financial cost down. However, primer pairs should also be unique, otherwise the specificity of the targets will be reduced.
Primers should not be designed that will hybridize to the target between a forward and reverse primer of an already established amplicon, unless this will lower the cost of the overall design.
Primers should not be designed which are complementary to other primers, i.e. form dimers with other primers.
Primers should not be designed which are likely to self-hybridize and form hairpin structures.
We attempt to minimize cost D by allowing degeneracy in i and by sampling possible designs within design constraints using an MCMC sampler. Greedy approaches to minimizing cost D in a primer set design for multiple target sequences will result in a final cost D which depends in part on which of the target sequences is chosen first, to begin the design, and on which position within each sequence the primer design is begun for that sequence. We sample from uniform prior distributions on both of these parameters. We also allow our sampler to vary the position of individual primers within the design, and the degeneracy of individual bases within each primer within the chosen degeneracy limit, both also with uniform priors. Our sampler allows each of these parameters to vary, initially preferring to sample from the earlier of the four parameters and increasingly preferring to sample from the latter.
Input file parameters for MCMC-ODPR
Genome file for BLAST searching
Max degeneracy (base pairs)
Minumum melting temperature (centigrade)
Maximum melting temperature (centigrade)
Minimum amplicon length (base pairs)
Maximum amplicon length (base pairs)
Initial overlap (base pairs)
Number of optimisations
Maximum gap between sequences (base pairs)
Save interim optimizations?
Output file name
Restart from previous run?
Probability of removing redundant primer pairs
Proportion of iterations to be considered as `early'
Weight greedy methods according to optimization?
Remove non-reusable primers from initial design?
Proportion of failed weight check proposals to accept (heating)
MCMC-ODPR initially goes through three stages of valid primer enumeration. Firstly, it searches all input consensus sequences for all possible primers meeting design constraints, whilst iterating through the input range of melting temperatures. At the second stage primer-dimers are removed and the redundancy of primers and uniqueness of primer pairs is assessed. Primers that have the propensity to self hybridize are subsequently removed. Primers are weighted according to their redundancy, with more highly weighted primers being preferred later in the optimization. In the third stage, for all remaining accepted primers, MCMC-ODPR recursively generates all possible degenerate primers to the input target level of degeneracy and tests whether they match other primers in the data set, to reassess their redundancy. In order to restrict a lack of specificity in the targets of the designed primers, the degree of degeneracy and redundancy must be limited, however. In the case of redundancy a high cut-off value allowing large redundancy values could favor primers being designed from low complexity or microsatellite regions of the target genes. A cut-off of three S or W degenerate bases yielding 23 potential priming sites for degeneracy and a maximum redundancy of 10 are the default values in MCMC-ODPR; these can be changed by the user. After these initial stages, the Metropolis Hastings MCMC optimization begins.
When designing degenerate primers, it is important to maintain a degree of specificity within the primers by limiting the number of degenerate nucleotides per primer, otherwise designed primers will likely hybridize to numerous loci on the target sequence; generating highly unspecific amplicons and reducing the scope for primer reuse and cost effectiveness. Furthermore, if primer Tm is constrained, the possible combinations of degenerate loci that are permitted within a design must also be constrained to avoid generation of degenerate primers with variant Tm. For instance, a primer having one degenerate position allowing an A or a C (IUPAC code M) would generate two primers with different Tm values, which is inappropriate if primers are being designed for a specific Tm. For our purposes, therefore, the degree of possible degeneracy per primer is limited with primers having lower degeneracy, yet higher redundancy being more highly weighted for each proposal. Degenerate oligonucleotides which are permitted in our primer designs are W (bases A and T) and S (bases G and C). Take P to be the set of all possible primers extracted from the candidate gene set. A degeneracy threshold α is implemented and represents the maximum proportion of degenerate nucleotides within the primer. MCMC-ODPR then iteratively generalizes through all possible degenerate variations of the primers within P, whilst not contravening α. This exhaustive set is then searched during the greedy MCMC algorithms described below.
The primer cost optimization itself consists of three different greedy proposal generating methods that are progressively applied more frequently in succession. The first covers all input gene sequences by randomly selecting genes progressively and then picks primers from a random seeded position within the sequence, building amplicon targets from that seed in a periodicity which matches a user defined size range by randomly selecting amplicon sizes within that range to look for subsequent primer positions. This process is repeated, with the subsequent primer set being selected if the global cost is lowered. In this process primer loci may be ‘swapped’, with more reusable primers being accepted over less reusable primers. The second method (per-gene refinement) repeats the seeding process for each gene individually having accepted the global primer arrangement from the first step. The final greedy method adjusts final primer positions by repositioning each primer randomly. Again, if each primer adjustment results in a cost reduction, the adjustment is accepted. To allow the algorithm to converge upon a solution, it is desirable to allow a transition from optimizations that greatly affect a design to more ‘fine tuning’ optimizations. Therefore these methods may be weighted according to how many iterations within the optimization have passed: with the first greedy method being weighted for within the first third of all iterations, the second within the second third of optimizations and the third within the last third. These methods are stochastically chosen, however, hence all three greedy methods can still be chosen randomly by the algorithm at any time, allowing primer set design space to be sampled in different ways.
The primer cost sample space may be visualized as a landscape, with various local cost minima and the global minima. The Metropolis-Hastings algorithm takes a solution proposal and then assesses whether this current proposal is more optimal than previous proposals. If this is the case, the Markov chain is allowed to continue stepping through the landscape. Otherwise the algorithm reverts back to a previous state. An initial random cost, C is proposed. At each iteration the cost associated with the new proposal, C’ is compared to C, with δS = C’-C. Unless δS is negative, the proposal is accepted with a probability proportional to e-δSH; where the heating H is a user input variable (see Table 1). The heating variable is defined as the proportion of primers that have failed their redundancy/weight check proposals that may be accepted, allowing chains to escape local minima and allowing efficient convergence. A variable specifying the proportion of optimizations to consider as ‘young’ is also input to allow convergence to be estimated in mature chains, as an alternative to specifying a fixed number of iterations.
The parameter INITIAL_OVERLAP allows the user to build ‘slack’ into the initial proposal, allowing for a more extensive exploration of the cost landscape before the proposed solution tends towards being trapped in a local optimum. If larger values of INITIAL_OVERLAP are used, the initial proposal becomes more expensive because more primer pairs are used to cover the target sequences, but by initially proposing a more expensive design there is more scope for the algorithm to explore the cost landscape by changing the positions of individual primer pairs, whilst still satisfying the maximum amplicon length constraint. The parameter PROB_REMOVE_REDUNDANT constrains the rate at which this initial slack is removed from the design in cases where two neighboring primer pairs have changed position to the extent that the pairs become nested.
Currently the algorithm only aims to optimize cost, and not coverage. Because of this, as the design begins to approach an optimum cost for a given coverage, it will then likely allow coverage to reduce in order to further minimize the cost by removing one or more primer pairs from the design, for instance, and leaving the corresponding sequence uncovered.
The output of MCMC-ODPR is a single text file, containing an enumeration of all primers designed at each input Tm, the selected degenerate primers, the number of covered bases and the final cost in nucleotides.
The variance in the number of primers obtained from 1–5x redundancy
Frequency of redundancy value when maximum
Number of primers >0
Average number of primers identified
We compared MCMC-ODPR with primer design programs Primer3, Primer-BLAST, BatchPrimer3 and PAMPS. For primer design from single input sequences MCMC-ODPR was compared to Primer3, Primer-BLAST, and for multiple input sequences BatchPrimer3 and PAMPS. Tm = 60 was used for the design of primers from MCMC-ODPR, with 10,000 iterations. All of the programs we compared with MCMC-ODPR produced many more primers than would be required to amplify the target regions, with many amplicon targets occurring in the same area. Consequently, the total primer cost was noted, but the outputs were processed with a script such that non-overlapping primer pairs were selected that represented the maximum coverage of the target gene sequences.
Performance of primer design software for single sequence input
Size of gene (bp)
Number of primers in design
Cost per covered base (nucleotides)
C-repeat binding factor 3 like protein
Primer3 designed 9 to 33 primer pairs per input sequence, from which 4 to 10 primers were selected to provide optimal coverage without overlap, Table 3. MCMC-ODPR achieved a lower cost per unit coverage than Primer 3 in all sequences except Dehydrin 9, where Primer 3 achieved a value of 0.2 compared with MCMC-ODPR, which achieved a value of 0.21. Primer-BLAST designed the number of primers specified in each case (12, 36 and 14), from which an optimal subset 6 to 14 was selected. The resulting cost of primer for the coverage achieved was consistently lower in MCMC-ODPR. Furthermore the amount of target sequence that was covered by the PCR systems designed by MCMC-ODPR was over 95% in all cases, whereas the maximum coverage achieved by Primer3 and Primer-BLAST was 40.1% and 52.7%, respectively.
Performance of primer design software for multiple sequence input
Number of primers designed
Cost per covered base (nucleotides)
Execution time (minutes)*
PAMPS designs degenerate primers according to ad-hoc pairwise alignments. The output of PAMPS is a list of degenerate primers consisting of a range of Tm with no information on which sequence its primers have been generated from. A script was therefore written to match the primers generated by PAMPS to our sequence dataset, keeping the same amplicon length constraints as above and only choosing primers with Tm = 60 degrees centigrade.
PAMPS designed 1205 primers that gave maximal coverage, which was comparable to the number generated by MCMC-ODPR at 1521 primers. PAMPS achieved the lowest coverage of 18.72% of the input sequences with a cost of 25,284 nucleotides. MCMC-ODPR outperformed PAMPS in terms of cost per unit coverage, with this value being calculated to be 0.64 for PAMPS.
MCMC_ODPR execution time with different numbers of sequences input
Number of input sequences
Execution time (minutes)
MCMC-ODPR is a useful addition to the suite of primer design programs available in particular reducing costs through improved solutions to the MDPSP. The bench tests carried out in this study show that MCMC-ODPR consistently produced the lowest cost primer design all cases but two, but even then the slightly higher cost achieved a much larger coverage. Consequently, MCMC-ODPR was the most economically efficient primer design in all cases. The comparisons with other software show that by combining degeneracy with optimal primer reuse the user may increase the coverage of the sequences obtained by the designed primers at significantly lower costs. MCMC-ODPR’s enumeration through an input temperature range will allow the user to observe the trade off between Tm and the suitability of the sequence space for invoking reuse. However, as MCMC-ODPR is an optimization technique that utilizes the Metropolis-Hastings algorithm for optimizing primer reuse, it is intrinsically slow when compared to heuristic algorithms, and may possibly require a number of hours for multiple sequences, especially if a range of Tm has been input.
Project name: MCMC-ODPR
Project home page [permissions have been granted for open (anonymous) access]: http://www2.warwick.ac.uk/fac/sci/lifesci/research/archaeobotany/downloads/MCMC_ODPR
Operating system(s): MacOSX
Programming language: Perl, Java
License: GNU General Public License (GNU-GPL)
Any restrictions to use by non-academics: none.
This work was supported by the NERC (NE/G005974/1).
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.