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.
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:
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:
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.
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.
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|
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
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
# Show the binding region of the first template for forward primers seq.df$Allowed_fw #>  "atggactggacctggagcatccttttcttg" # Show the corresponding interval in the templates c(seq.df$Allowed_Start_fw, seq.df$Allowed_End_fw) #>  1 30 # Show the binding region of the first template for reverse primers seq.df$Allowed_rev #>  "cgacacggccgtgtattactgtgcgagaga" # Show the corresponding interval in the templates c(seq.df$Allowed_Start_rev, seq.df$Allowed_End_rev) #>  324 353
In the following sections, we describe two ways in which the primer binding regions can be defined using
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 (
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 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
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:
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
You can use
str(settings) to explore the structure of
||Desired values for primer properties|
||Limits for relaxing constraints during primer design|
||Constraints for estimating the coverage|
||Experimental PCR conditions|
||Settings for evaluating the constraints|
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.
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
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.
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:
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
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.
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:
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.
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
We can view the identifiers and sequences of the primers in the following way:
print(primer.df[, c("ID", "Forward")]) #> ID Forward #> 2 Ippolito2012|VH1|1_fw caggtccagctkgtrcagtctgg #> 1 Ippolito2012|VH157|2_fw caggtgcagctggtgsartctgg #> 3 Ippolito2012|VH2|3_fw cagrtcaccttgaaggagtctg #> 5 Ippolito2012|VH3|4_fw gaggtgcagctgktggagwcy #> 7 Ippolito2012|VH4|5_fw caggtgcagctgcaggagtcsg #> 6 Ippolito2012|VH4-DP63|6_fw caggtgcagctacagcagtggg #> 8 Ippolito2012|VH6|7_fw caggtacagctgcagcagtca #> 4 Ippolito2012|VH3N|8_fw tcaacacaacggttcccagtta
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
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
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.
The number of templates that are covered by each primer is annotated in the
To identify the primers that cover individual templates, we can annotate the coverage information in the template data frame using
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
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:
|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:
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.
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
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[]. 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
original.cvg <- as.character(get_cvg_ratio(constraint.df, template.df, as.char = TRUE)) subset.cvg <- as.character(get_cvg_ratio(my.primer.subset, template.df, as.char = TRUE)) print(paste0("Coverage (n=", nrow(constraint.df), "): ", original.cvg, "; Subset Coverage (n=", nrow(my.primer.subset), "): ", subset.cvg)) #>  "Coverage (n=8): 99.35%; Subset Coverage (n=3): 98.71%"
The binding regions of the primers in the templates can be visualized with
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.
We can determine which primers passed the physicochemical constraints that were supplied to the
check_constraints function through visual inspection via
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
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.
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
Note that this function requires
pandoc as well as
LaTeX such that
rmarkdown can create the PDF report.
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:
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)
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:
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:
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:
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.
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