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:
# Define the FASTA primer file to load
primer.location <- system.file("extdata", "IMGT_data", "primers", "IGHV",
"Ippolito2012.fasta", package = "openPrimeR")
# Load the primers
primer.df <- read_primers(primer.location, fw.id = "_fw")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 tcaacacaacggttcccagttaEvaluation 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:
# Allow off-target coverage
conOptions(settings)$allowed_other_binding_ratio <- c("max" = 1.0)
# Evaluate all constraints found in 'settings'
constraint.df <- check_constraints(primer.df, template.df,
settings, active.constraints = names(constraints(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():
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))
#> [1] "Coverage (n=8): 99.35%; Subset Coverage (n=3): 98.71%"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:
# View the number of terminal GCs for primers failing the GC constraint
constraint.df$gc_clamp_fw[!constraint.df$EVAL_gc_clamp]
#> [1] 0 0 0
# View the desired number of terminal GCs:
constraints(settings)$gc_clamp
#> min max
#> 1 3If 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:
filtered.df <- filter_primers(constraint.df, template.df, settings,
active.constraints = c("gc_clamp", "gc_ratio"))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():
# Define the path to the output file
my.file <- tempfile(fileext = ".pdf")
# Store a PDF report for 'constraint.df' in 'my.file'
create_report(constraint.df, template.df, my.file,
settings, sample.name = "My analysis")Note that this function requires pandoc as well as LaTeX such that rmarkdown can create the PDF report.
