Designing and analyzing multiplex PCR primers with openPrimeR

Matthias Döring

2022-04-26

openPrimeR provides functionalities for designing and analyzing multiplex polymerase chain reaction (PCR) primers. In the following, we introduce typical workflows for three application scenarios, namely designing primers, analyzing primers, and comparing primer sets.

Overview of the package

openPrimeR was developed to provide a rational approach for evaluating and designing primers for multiplex PCR such that multiple template sequences are amplified at the same time. The concept of coverage is of critical importance for multiplex PCR as it describes the number of templates that can be amplified with a set of primers. This package was specifically developed to enable researchers to evaluate the coverage of existing sets of primers as well as to design novel primer sets that maximize the coverage with a minimal number of primers. To provide a user-friendly tool, we created a Shiny app, which is available through the openPrimeRui package. The openPrimeR package enables the computation of the most important PCR-related physicochemical properties of primers to gauge whether a set of primers can provide high yields. The following list provides research questions that may be answered using openPrimeR:

Which templates are amplified by the primers?

How well do the primers fulfill desired physicochemical properties?

Among multiple sets of primers, which seems to be the best for a specific task?

What is the smallest set of primers that covers all of the template sequences?

Preliminaries

openPrimeR requires external programs for some features, particularly for computing the physicochemical properties of primers. Please make sure you have the following tools installed on your system such that they are in your system’s path:

If you would like to be able to access the immunoglobulin repository IMGT from the openPrimeR Shiny app, you should additionally fulfill the following dependencies:

openPrimeR will automatically check for all dependencies and inform you about any missing dependencies when the package is attached:

library(openPrimeR)

Note that the tool is still functional if there are missing external programs. However, we recommend that all dependencies are fulfilled to guarantee the best user experience.

Using the Shiny app

If you would like to use the openPrimer shiny application, please install the openPrimeRui package and consider its documentation. In the following, we will solely focus on the openPrimeR package itself and not on the frontend.

Loading data

In order to design primers, we only need to load a set of template sequences and define the target binding regions. To analyze and compare the properties of existing primer sets, we also need to load one or multiple sets of primers. The following table summarizes the possible input data formats for each task:

Task Templates Primers Input file format
Design primers ✓ FASTA, CSV
Analyze primers ✓ ✓ FASTA, CSV
Compare primers ✓ ✓ FASTA, CSV

Loading templates

To load a set of template sequences, we simply input the path to a valid FASTA file and use read_templates(). In the following code snippet, we store the path to a FASTA file that is provided with the package in fasta.file. This file contains the template sequences. In our case, we will load sequences of the human heavy chain immunoglobulin genes:

Having loaded the template sequences successfully, we can investigate the structure of seq.df.simple, which is a Templates object. Since the Templates class is derived from data.frame, you can use it in the same was as any data frame. For example, we can retrieve the header of the first template sequence via

As you can see, the headers of the templates contain several pieces of information that are separated by pipe symbols (|), among others:

To load these annotations into the Templates object, we will now provide additional arguments to read_templates(). We are particularly interested in loading the group information as this information is important for interpreting the results later. In the next code snippet, we use the hdr.structure variable to annotate the meta-information that is provided by the headers of the FASTA file. Moreover, we provide the delim argument to read_templates() to specify that the pipe symbol is used to separated individual fields and supply the id.column argument to identify which field should be used as the identifier in the Templates object:

Since we have specified to load the accession, group, species, and function from the FASTA headers, we can now retrieve these values from seq.df. For example, for the first template we can retrieve the following information:

Note that only the GROUP annotation has an impact on the analysis in terms of visualizing the results later. The other fields can be set arbitrarily and do not have any impact. If there is no metadata that you wish to load, you can simply use the read_templates() call used to declare seq.df.simple. In this case, all templates are considered to belong to a single default group.

Upon loading templates with read_templates(), the primer binding region is set to the first 30 bases for forward primers and to the last 30 bases for reverse primers, where first refers to the 5’ end and last refers to the 3’ end. We can review the target binding regions of forward and reverse primers by accessing seq.df$Allowed_fw or seq.df$Allowed_rev, respectively:

In the following sections, we describe two ways in which the primer binding regions can be defined using assign_binding_regions().

Uniform binding regions

To assign a uniform target binding region to all templates, you can specify positional intervals indicating the binding regions for forward and reverse primers. For forward primers, the interval is specified relative to the template 5’ end, while for reverse primers, the interval is specified relative to the 3’ end. In the following example, we set the binding region of forward primers (fw) to the first 50 template bases and to the last 40 bases for reverse primers (rev):

Note that we have supplied the interval [1,40] to allow binding in the last 40 bases of the templates for reverse primers. This is because the binding region for reverse primers is provided relative to the 3’ end, while the binding region of forward primers is provided relative to the 5’ end. In this way, the reverse binding region can be annotated independent of the length of individual templates.

Let’s verify the different binding regions for forward and reverse primers using the first template sequence:

Individual binding regions

Individual binding regions can be assigned to each template by providing a FASTA file for each primer direction containing the target binding regions for the primers. The FASTA headers of these files should match the headers in the template FASTA file that was provided earlier. In the following example, we use a FASTA file that is provided with the package to define the individual binding regions for the forward primers only. In this case, the FASTA file specifies the leaders of the human heavy chain immunoglobulin sequences that have been loaded in seq.df:

The binding regions for forward primers may now be different for each template. For example, the binding regions for the following two templates occur at different positions in the templates:

Note, that since we did not supply individual binding regions for reverse primers, their binding regions were not adjusted:

Loading and writing settings

Before we can start an analysis, we need to define the analysis settings. openPrimeR supplies predefined XML files specifying default settings for different applications:

In our case, we select the high-stringency primer design conditions for Taq polymerase and load the DesignSettings object with read_settings()):

You can use str(settings) to explore the structure of settings:

Slot Getter/Setter Purpose
Input_Constraints constraints() Desired values for primer properties
Input_Constraint_Boundaries constraintLimits() Limits for relaxing constraints during primer design
Coverage_Constraints cvg_constraints() Constraints for estimating the coverage
PCR_conditions PCR() Experimental PCR conditions
constraint_settings conOptions() Settings for evaluating the constraints

Since the settings object contains all the relevant information for the analysis of primers, you should review and possibly customize the settings before starting an analysis. Particularly make sure that the PCR conditions specified in the settings agree with your experimental conditions. For example, the PCR sodium ion concentration can be accessed via PCR(settings)$Na_concentration. Moreover, the coverage constraints should typically contain one constraint (e.g. coverage_model, primer_efficiency, or annealing_DeltaG). Moreover, for designing primers, you should select a coverage constraint that provides highly specific coverage calls. For the purpose of this vignette, however, we will not apply any coverage constraints and instead stringently limit the number of mismatches between primers and templates:

You should also make sure that the requirements physicochemical constraints for high-quality primers fulfill your expectations. For example, if we do not want to filter designed primers using the GC clamp criterion, we can remove the requirement for a GC clamp in the following way:

We may also want to design only primers of a specific length. To generate only primers of length 25, we could specify this via

For designing primers we may also want to prevent any mismatch binding. To implement this, we simply set

For more possible customizations, please refer to the documentation of the settings class, which can be accessed via ?DesignSettings.

Having customized the settings, we can store the modified settings to disk in the following way:

The next time we want to perform an analysis, we can simply load the stored settings using the read_settings function using out.file as an argument.

Designing primers

Since the loaded template sequences contain only the variable region of the human immunoglobulins and not the constant region, we will restrict ourselves to designing only forward primers. To design primers, we need only a single function, namely design_primers(), which requires the settings specified in design.settings and the templates and binding regions provided in template.df. By setting mode.directionality to ‘fw’, we specify that we want to design forward primers only. The other possible choices are rev for designing only reverse primers and both for designing both forward and reverse primers. Moreover, we subset template.df to the first two template sequences to limit the runtime of the design procedure for this example:

# Design forward primers for the first two templates
optimal.primers <- design_primers(template.df[1:2,], mode.directionality = "fw",
                                  settings = design.settings)

optimal.primers is a list in which the optimal primers are stored in optimal.primers$opti and the corresponding filtered primers are stored in optimal.primers$filtered. The optimal primer sets for all evaluated melting temperatures are stored in optimal.primers$all_results.

The primer design function can be customized in many ways. The init.algo argument can be used to specify how the initial set of primers is generated. By default, this is done by extracting substrings from the template sequences (‘naive’). A tree-based initialization strategy producing degenerate primers can be activated by setting init.algo to ‘tree’, which is favorable for related template sequences.

Another important argument is required.cvg, which defines the desired coverage ratio of the templates. By default, this is set to 1 indicating that 100% of templates should be covered. If the desired coverage cannot be reached with a primer set fulfilling all properties specified via the settings argument, the constraints are relaxed in order to reach the target coverage. This behavior can be deactivated by setting required.cvg to 0.

The opti.algo argument specifies the algorithm that is used for optimizing the primer sets. By default, this is a greedy algorithm (‘Greedy’), but we can also select an integer linear program (‘ILP’). The ILP ensures that designed primer sets are minimal, but this comes at the cost of a worst-case exponential runtime. Hence, the ILP is recommended for small template sets. For larger sets of templates, the default (‘Greedy’) should be used, which provides an approximate solution in polynomial runtime.

To conclude the primer design task, let’s store the designed primers as a FASTA file:

out.file <- tempfile("my_primers", fileext = ".fasta")
write_primers(optimal.primers$opti, out.file)

Analyzing primers

In the following, we will analyze the physicochemical properties of an existing primer set. Since the set of primers we have designed earlier is quite small, we will load a new set of primers first.

Loading primers

Similarly to templates, primers can also be loaded from a FASTA file. Similarly to ensuring that the information in the header is loaded correctly via read_templates(), you should make sure that the primer directionalities are correctly annotated in the FASTA file. This means that the header of every primer should contain a keyword that uniquely identifies the primer to either be a forward or reverse primer. The FASTA file we are going to load contains the forward primers that were designed for the variable region of human IGH genes by Ippolito et al.. Since all primers are forward primers, the headers of all primers in the FASTA file are annotated with the ‘_fw’ keyword. Therefore, we set the fw.id argument for read_primers() accordingly:

We can view the identifiers and sequences of the primers in the following way:

Evaluation of biochemical constraints

To determine which biochemical constraints are fulfilled by the loaded primers, we can use check_constraints(), which requires only the set of primers, a set of templates, and a settings object. Since we want to determine the coverage of the primers along the full template sequence, we will set the allowed ratio of off-coverage events to 100% by modifying the constraint settings. We then call check_constraints() in order to determine all constraints defined in settings:

Note that we could have also computed only a subset of the constraints specified via settings, if we had passed the active.constraints argument to check_constraints(). The constraint.df variable provides a Primers data frame containing the results from evaluating all constraints provided via settings (e.g. primer coverage, GC ratio, melting temperature). We can retrieve the values of individual constraints by accessing the corresponding columns in constraint.df. In the following, we will explore the coverage of templates in more detail.

Primer coverage

The number of templates that are covered by each primer is annotated in the primer_coverage column:

To identify the primers that cover individual templates, we can annotate the coverage information in the template data frame using update_template_cvg():

Now, the primer_coverage column is also available in template.df such that we can identify the number of primers that cover the first five templates in the set:

The overall ratio of covered template sequences can be retrieved via get_cvg_ratio():

The output indicates that the primer set is expected to amplify 99.35% of templates according to the currently specified conditions for primer coverage (at most 3 mismatches). We can learn more about the coverage by computing further statistics. For example, we may be interested in identifying which groups of templates are expected to be amplified and which are not. For this purpose, we can use get_cvg_stats(), which provides information on the number of covered template sequences according to their group annotation:

Group Coverage
Total 154 of 155 (99.35%)
IGHV1 25 of 25 (100%)
IGHV2 21 of 21 (100%)
IGHV3 55 of 56 (98.21%)
IGHV4 47 of 47 (100%)
IGHV5 3 of 3 (100%)
IGHV6 2 of 2 (100%)
IGHV7 1 of 1 (100%)

In this case, we see that all but a single template of IGHV3 does not seem to be covered by the primer set. A visualization of the coverage can be obtained via

This plot shows three bars:

  • Identity Coverage: The ratio of templates that are covered by primers that are perfectly complementary to the templates
  • Expected Coverage: The expected ratio of covered templates under the current conditions for computing the coverage
  • Available Templates: The number of available template sequences per group of templates

For the loaded set of primers, the identity coverage is nearly as high as the expected coverage, which indicates that the analyzed primer set should have a high amplification fidelity. Note that for IGHV5, the identity coverage is 0%, while the expected coverage is 100%. This just shows that there is no primer that is fully complementary to any of the IGHV5 templates, however, according to the coverage conditions there are primers that are expected to cover all of the IGHV5 templates via mismatch binding.

Optimal primer subsets

Typically one wants to avoid large primer sets for reasons of cost and priming efficiency. We provide a method for reducing the size of an existing primer set by computing optimal subsets with respect to the coverage of templates. Using this approach, it is possible to limit the size of a primer set without sacrificing any coverage. We can compute all optimal subsets of a primer set for which the coverage has already been annotated using subset_primer_set():

primer.subsets is a list whose i-th entry contains the primer set of size i. For example, the optimal subset of size 3 could be accessed via primer.subsets[[3]]. We can easily determine the most suitable subset by plotting the coverage of all optimal subsets:

The plot visualizes two types of information. The line graph shows the total percentage of covered templates, while the stacked bars indicate the coverage of individual primers. Since some of the primers in the set cover the same templates, the cumulative coverage of the stacked bars exceeds 100% for subsets of size 3 or larger, although the total coverage already seems to be saturated for the subset of size 3. It seems that subsets with more than 3 primers offer only additional redundant coverage of the templates. Hence, we might decide to select the primer subset of size 3 as it seems to achieve the same coverage as the full primer set:

We can verify that the coverages of the full primer set and the subset of size 3 seem to match using get_cvg_ratio():

Binding regions

The binding regions of the primers in the templates can be visualized with plot_primer_binding_regions():

The x-axis of the plot shows the binding positions of the primers relative to the target binding regions. Position -1 indicates the end of the target binding region and we can see that all primers bind beyond the target region. Since we have specified the leader of the immunoglobulins as the target region, this just means that all primers bind at the beginning of the exon, which ensures that the full antibody sequence is recovered by the primers. Typical primer sets for amplifying immunoglobulins will target the exon, that is, the region directly following the leader. To investigate the individual positions of primer binding for individual templates, we can use plot_primer(). In this case, we just provide the first primer and the first ten templates to the plotting function in order to restrict the dimension of the plot:

The plot shows primers that cover a template as arrows above the corresponding templates. If a template is covered it is shown as a black line and otherwise it is shown in grey.

Constraint evaluation

We can determine which primers passed the physicochemical constraints that were supplied to the check_constraints function through visual inspection via plot_constraints():

The plot shows which constraints are passed (blue) and which constraints are failed (red) for every primer. For example, looking at the column indicating the GC clamp constraint, we see that only Ippolito2012|VH3|4_fw and Ippolito2012|VH6|7_fw and Ippolito2012|VH3N|8_fw did not fulfill our requirements for the GC clamp. This is because both primers have no GC clamp, although we required a GC clamp with a length between 1 and 3 in the settings:

If you are wondering why the specificity of all evaluated primer seems so low, this can be explained by the binding regions of the primers. As we have seen earlier, all primers bind outside the target binding region. Therefore, the specificity of every primer is 0% and the constraint on the specificity of the primers is never fulfilled.

We can not only can visualize whether a primer passed a specific constraint, but also view the distribution of values corresponding to a specific property. For example, let us take a look at the number of terminal GCs in the primers by creating a histogram:

The y-axis shows the number of primers exhibiting a certain number of GCs at their 3’ ends. The dashed lines indicate the desired extent of the GC clamp according to the settings object that we passed to plot_constraint().

Filtering primers

All our previous evaluations were performed without requiring the primers to actually fulfill any of the constraints we postulated. We can filter primers according to a set of biochemical constraints such that only primers fulfilling all requirements are retained. For example, if we want to select only primers fulfilling the requirements for GC clamp and melting temperature range, we could obtain the filtered data set in the following way:

Now, we could perform further analyses on this data set. For example, using get_cvg_ratio(), we could determine that the percentage of templates that are covered by the filtered primer set is only 13.55% since only 1 primer remains after filtering.

Report generation

We can create a PDF report that summarizes the analysis of a set of primers by passing an analyzed primer set with its corresponding templates and settings to create_report():

Note that this function requires pandoc as well as LaTeX such that rmarkdown can create the PDF report.

Comparing primer sets

To compare existing primer sets, it is necessary to precompute all constraints of interest for each of the primer sets, which can be done with check_constraints() as we have demonstrated earlier. Instead of evaluating multiple primer sets, we will simply load pre-evaluated sets of primers and template sets that were stored as CSV files. For example, we could have stored our previous evaluation results in the following way in order to load these data later:

primer.xml <- tempfile("my_primers", fileext =".csv")
write_primers(constraint.df, primer.xml, "CSV")
template.xml <- tempfile("my_templates", fileext = ".csv")
write_templates(constraint.df, template.xml, "CSV")

For the following example, we will simply load primers and templates that are shipped with the openPrimeR package:

# Define the primer sets we want to load
sel.sets <- c("Glas1999", "Rubinstein1998", "Cardona1995", "Persson1991", "Ippolito2012", "Scheid2011")
# List all available IGH primer sets
primer.files <- list.files(path = system.file("extdata", "IMGT_data", "comparison", 
                           "primer_sets", "IGH", package = "openPrimeR"),
                pattern = "*\\.csv", full.names = TRUE)
# Load all available primer sets
primer.data <- read_primers(primer.files)
# Select only the sets defined via 'sel.sets'
sel.idx <- which(names(primer.data) %in% sel.sets)
primer.data <- primer.data[sel.idx]
# Provide a set of templates for every primer set
template.files <- rep(system.file("extdata", "IMGT_data", "comparison", "templates", 
                              "IGH_templates.csv", package = "openPrimeR"), 
                              length(primer.data))
template.data <- read_templates(template.files)

Both, primer.data and template.data are lists containing primer and template data frames, respectively. Note that these lists should always have the same lengths, that is, every primer set should have an associated template set and vice versa. Having loaded the primer and template data sets, we can plot an overview of the constraints that are fulfilled by the primers in each set:

plot_constraint_fulfillment(primer.data, settings, plot.p.vals = FALSE)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

In this plot, each facet corresponds to a primer set and each physicochemical constraint is shown as a colored bar whose height indicates the percentage of primers fulfilling the constraint. An intuitive interpretation of the plot is that sets with high-quality primers should have many high bars. For example, we can quickly see that the primers from Persson et al. (1991) do not have any primers fulfilling our requirements for the ratio of GC, which should be between 1 and 3 according to the loaded settings in order to ensure similar binding behaviors of the primers.

The distribution of each evaluated constraint can be investigated in more detail. For example, to investigate the influence of the GC ratios on the melting temperatures of the primers in each set, we can create the following boxplot:

plot_constraint(primer.data, settings, active.constraints = c("gc_ratio", "melting_temp_range"))
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

In the boxplot, each dot corresponds to a single primer and the boxes show the 1st, 2nd, and 3rd quartiles, from bottom to top. Since we have provided the current constraint settings to the plotting function, the desired ranges for GC ratios and melting temperatures are indicated as horizontal dashed lines in the plot. The plot shows that the GC ratios from the primers in the set from Cardona et al. (1995) are all in the desired range and have a small spread, while the GC ratios of other primer sets have much larger spreads, for example those from the set of Glas et al. (1999). The visualization reveals the association between GC ratios and melting temperatures. The primer sets from Cardona, Persson, and Rubinstein all have small deviations in their GC ratios effecting small deviations in their melting temperatures, while the other primer sets exhibit high variances for both properties.

Remember that the plot we generated using plot_constraint_fulfillment() revealed that most of the loaded primer sets failed the specificity constraint. By plotting the binding regions for each of the primer sets through plot_primer_binding_regions() we can find out why this is the case:

plot_primer_binding_regions(primer.data, template.data)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
#> `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

The plot reveals that only the primers from Scheid et al. mostly bind in the target region, while the other primer sets all bind outside the target region. Moreover, since we have evaluated the primers allowing for off-target binding, we find that the forward primers from Glas et al. do not seem to target the 5’ region of the templates but are rather spread along the length of the templates, a property that would be detrimental when trying to amplify complete antibody cDNA.

Want to learn more?

For brevity’s sake, we did not provide explanations and examples for all possible uses of the package. If you would like to learn more about how you can use openPrimeR, please consider using our interactive tutorial. The tutorial is provided as a Shiny app, which can be started via runTutorial().