TAP-seq primer design workflow

Installation

This package requires local installations of Primer3 and BLASTn. TAPseq was developed and tested using Primer3 v.2.5.0 and blastn v.2.6.0. It’s strongly suggested to use Primer3 >= 2.5.0! Earlier versions require a primer3_config directory, which needs to be specified whenever calling functions interacting with Primer3. Source code and installation instructions can be found under:

Please install these tools first and add them to your PATH. If you don’t want to add the tools to your “global” PATH, you can add the following code to your .Rprofile file. This should add the tools to your PATH in R whenever you start a new session.

Sys.setenv(PATH = paste(Sys.getenv("PATH"), "/full/path/to/primer3-x.x.x/src", sep = ":"))

Sys.setenv(PATH = paste(Sys.getenv("PATH"), "/full/path/to/blast+/ncbi-blast-x.x.x+/bin",
sep = ":"))

Alternatively you can specify the paths to 3rd party software as arguments when calling TAPseq functions (TAPseqInput(), designPrimers(), checkPrimers()).

Transcript sequences

Let’s start by creating transcript sequences for primer design for a target gene panel. First we try to infer likely polyadenylation (polyA) sites for the selected target genes.

library(TAPseq)
library(GenomicRanges)
library(BiocParallel)

# gene annotations for chromosome 11 genomic region
data("chr11_genes")

# convert to GRangesList containing annotated exons per gene. for the sake of time we only use
# a small subset of the chr11 genes. use the full set for a more realistic example
pcr_products(outer_primers$HBE1) # these can also be accessed for all genes tapseq_primers(outer_primers) pcr_products(outer_primers) We have now successfully designed 5 outer primers per target gene. To select the best primer per gene we could just pick the primer with the lowest penalty score. However, we will now use BLAST to try to estimate potential off-target priming for every primer. We will use this to then select the best primer per gene, i.e. the primer with the fewest off-targets. BLAST primers Creating a BLAST database containing all annotated transcripts and chromosome sequences takes a couple of minutes (and ~2Gb of free storage). For this example, we can generate a limited BLAST database, containing only transcripts of all target genes and the sequence of chromosome 11. # chromosome 11 sequence chr11_genome <- DNAStringSet(getSeq(hg38, "chr11")) names(chr11_genome) <- "chr11" # create blast database blastdb <- file.path(tempdir(), "blastdb") createBLASTDb(genome = chr11_genome, annot = unlist(target_genes), blastdb = blastdb) In a real-world application one would generate a database containing the entire genome and transcriptome. The database can be saved in a location and used for multiple primer designs whithout having to rebuild it everytime. library(BSgenome) # human genome BSgenome object (needs to be istalled from Bioconductor) hg38 <- getBSgenome("BSgenome.Hsapiens.UCSC.hg38") # download and import gencode hg38 annotations url <- "ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz" annot <- import(url, format = "gtf") # extract exon annotations for protein-coding genes to build transcripts tx_exons <- annot[annot$type == "exon" & annotgene_type == "protein_coding"] # create blast database blastdb <- file.path(tempdir(), "blastdb") createBLASTDb(genome = hg38, annot = tx_exons, blastdb = blastdb) Once we have our BLAST database, we can use it to estimate the number of exonic, intronic and intergenic off-targets for all primers. Note, that the number of exonic and intronic off-targets is provided on the level of genes, i.e. one exonic off-target means that a given primer might bind in exonic region(s) of one other gene. If the off-target binding site overlaps with exons of two genes, it will be counted twice, as it could bind to transcripts of two genes. # now we can blast our outer primers against the create database outer_primers <- blastPrimers(outer_primers, blastdb = blastdb, max_mismatch = 0, min_aligned = 0.75) # the primers now contain the number of estimated off-targets tapseq_primers(outer_primersIFITM1)

To finalize our set of outer primers, we want to choose the best primer per target gene, i.e. the one with the fewest exonic, intronic and intergenic off-target hits (in that order).

# select best primer per target gene based on the number of potential off-targets
best_outer_primers <- pickPrimers(outer_primers, n = 1, by = "off_targets")

# each object now only contains the best primer
tapseq_primers(best_outer_primers\$IFITM1)

To design nested inner primers we simply need to repeat the same procedure with a smaller product size range.

# create new TsIO objects for inner primers, note the different product size
inner_primers <- TAPseqInput(txs_seqs, target_annot = truncated_txs,
product_size_range = c(150, 300), primer_num_return = 5)

inner_primers <- designPrimers(inner_primers)

# blast inner primers
inner_primers <- blastPrimers(inner_primers, blastdb = blastdb, max_mismatch = 0,
min_aligned = 0.75)

# pick best primer per target gene
best_inner_primers <- pickPrimers(inner_primers, n = 1, by = "off_targets")

Done! We succesfully designed TAP-seq outer and inner primer sets for our target gene panel.

Multiplex compatibility

As an additional step, we can verify if the designed primer sets are compatible for PCR multiplexing. For that we use Primer3’s “check_primers” functionality:

# use checkPrimers to run Primer3's "check_primers" functionality for every possible primer pair
outer_comp <- checkPrimers(best_outer_primers)
inner_comp <- checkPrimers(best_inner_primers)

We can now for instance plot the estimated complementarity scores for every pair. We highlight pairs with a score higher than 47, which is considered “critical” by Primer3 during primer design (see Primer3 for more information).

library(dplyr)
library(ggplot2)

# merge outer and inner complementarity scores into one data.frame
comp <- bind_rows(outer = outer_comp, inner = inner_comp, .id = "set")

# add variable for pairs with any complemetarity score higher than 47
comp <- comp %>%
mutate(high_compl = if_else(primer_pair_compl_any_th > 47 | primer_pair_compl_end_th > 47,
true = "high", false = "ok")) %>%
mutate(high_compl = factor(high_compl, levels = c("ok", "high")))

# plot complementarity scores
ggplot(comp, aes(x = primer_pair_compl_any_th, y = primer_pair_compl_end_th)) +
facet_wrap(~set, ncol = 2) +
geom_point(aes(color = high_compl), alpha = 0.25) +
scale_color_manual(values = c("black", "red"), drop = FALSE) +
geom_hline(aes(yintercept = 47), colour = "darkgray", linetype = "dashed") +
geom_vline(aes(xintercept = 47), colour = "darkgray", linetype = "dashed") +
labs(title = "Complementarity scores TAP-seq primer combinations",
color = "Complementarity") +
theme_bw()

No primer pairs with complementarity scores above 47 were found, so they should be ok to use in multiplex PCRs.

Export primers

To finish the workflow, we can export the designed primers to data.frames, which can be written to .csv files for easy storage of primer sets.

# create data.frames for outer and inner primers
outer_primers_df <- primerDataFrame(best_outer_primers)
inner_primers_df <- primerDataFrame(best_inner_primers)

# the resulting data.frames contain all relevant primer data
head(outer_primers_df)

Furthermore we can create BED tracks with the genomic coordinates of the primer binding sites for viewing in a genome browser.

# create BED tracks for outer and inner primers with custom colors
outer_primers_track <- createPrimerTrack(best_outer_primers, color = "steelblue3")
inner_primers_track <- createPrimerTrack(best_inner_primers, color = "goldenrod1")

# the output data.frames contain lines in BED format for every primer
exportPrimerTrack(inner_primers_track, con = "")
exportPrimerTrack(createPrimerTrack(best_outer_primers, color = "steelblue3"),
con = "path/to/primers.bed")