1安装

此包需要在本地安装Primer3和BLASTn。TAPseq使用Primer3 v.2.5.0和blastn v.2.6.0开发和测试。强烈建议使用Primer3 >= 2.5.0!早期版本需要一个primer3_config目录，在调用与Primer3交互的函数时需要指定该目录。源代码和安装说明可在以下地址找到:

Primer3:https://github.com/primer3-org/primer3/releases
BLASTn:https://www.ncbi.nlm.nih.gov/books/NBK279690

请先安装这些工具，并将它们添加到您的PATH中。如果你不想将这些工具添加到你的“全局”PATH中，你可以将以下代码添加到你的. rprofile文件中。这将在您启动新会话时将工具添加到R中的PATH中。

Sys。setenv(PATH = paste(Sys.getenv("PATH")， "/full/ PATH /to/primer3-x.x x. "x/src"， sep = ":"))setenv(PATH = paste(Sys.getenv("PATH")， "/full/ PATH /to/blast+/ncbi-blast-x.x x。X +/bin"， sep = ":"))

或者，您可以在调用TAPseq函数(TAPseqInput()， designPrimers()， checkPrimers())时指定第三方软件的路径作为参数。

2记录序列

让我们从为目标基因组的引物设计创建转录本序列开始。首先，我们试图推断可能的多聚腺苷酸(polyA)位点为选定的靶基因。

库(TAPseq)库(GenomicRanges)库(BiocParallel) #染色体11基因组区域数据的基因注释(“chr11_genes”)#转换为每个基因包含注释外显子的GRangesList。为了节省时间，我们只使用了chr11基因的一小部分。使用完整的集合来做一个更现实的例子target_genes <- split(chr11_genes, f = chr11_genes$gene_name) target_genes <- target_genes[18:27] # K562 Drop-seq read data(这只是R包中的一个小示例文件)dropseq_bam <- system. zip . zip文件(“extdata”、“chr11_k562_dropseq。bam"， package = "TAPseq") #注册后端并行化寄存器(multicoream (workers = 5)) #从Drop-seq数据中推测polyA站点polyA_sites <- inferPolyASites(target_genes, bam = dropseq_bam, polyA_downstream = 50, wdsize = 100, min_cvrg = 1, parallel = TRUE)

这里我们假设所有推断的polyA位点都是正确的。在实际应用中，建议手动检查polyA站点预测并删除任何明显的假阳性。最简单的方法是将polyA站点导出到.bed。然后可以将该文件与目标基因注释和.bam文件一起加载到基因组浏览器(例如IGV)中，以直观地检查polyA位点预测。

library(rtracklayer) export(polyA_sites, con = "path/to/polyA_sites. "Bed "， format = " Bed ") export(unlist(target_genes)， con = "path/to/target_genes. "， format = " Gtf ")

一旦我们对我们的polyA预测感到满意，我们就用它们来截断目标基因转录本。所有重叠任何polyA位点的靶基因转录本在polyA位点被截断，因为我们假设它们标记了转录本3 '端。默认情况下，最下游的polyA站点被用来截断一个转录本。如果每个基因有多个转录本重叠polyA位点，则从所有重叠转录本生成合并的截断转录本模型。如果给定基因的转录本与任何polyA位点都没有重叠，则从该基因的所有转录本建立一致的转录本模型。

# truncate_txs <- truncateTxsPolyA(target_genes, polyA_sites = polyA_sites, parallel = TRUE)

为了完成这一部分，我们需要提取被截断的转录本的序列。这里我们使用Bioconductor的BSgenome来导出转录序列，但是基因组序列也可以从fasta文件导入到DNAStringSet对象中。

库(BSgenome) #人类基因组BSgenome对象(需要从Bioconductor安装)hg38 <- getBSgenome("BSgenome. hsapiens . ucsc .hg38") #获取所有截断转录本的序列txs_seqs <- getTxsSeq(truncated_txs, genome = hg38)

3.设计外部底漆

Primer3使用boulder-IO记录进行输入和输出(参见:http://primer3.org/manual.html)．TAPseq实现了TsIO和TsIOList对象，用于在R的S4类系统中存储TAP-seq引物设计的Primer3输入和输出。在引物设计过程中，它们作为引物3的用户界面。

首先，我们根据获得的转录本序列创建Primer3输入，用于外部正向引物设计。所有PCR反应都使用了反向引物，这样Primer3就可以选择适合它工作的正向引物。默认情况下，该功能使用10x Beads-oligo-dT和正确的底漆，并据此选择最佳、最低和最高熔化温度。如果使用另一种协议(例如Drop-seq)，则可以对这些参数进行调整(TAPseqInput ?)．

——TAPseqInput(txs_seqs, target_annot = truncated_txs, product_size_range = c(350, 500)， primer_num_return = 5)

我们现在可以使用Primer3设计引物，并将它们添加到对象中。

#为每个目标基因设计5个外部引物outer_primers <- designPrimers(outer_primers)

输出TsIO对象包含设计引物和预期扩增子:

#获取特定基因的引物和pcr产品tapseq_primers(outer_primers$HBE1) pcr_products(outer_primers$HBE1) #这些也可以用于所有基因tapseq_primers(outer_primers) pcr_products(outer_primers)

我们现在已经成功设计了每个靶基因5个外引物。要为每个基因选择最好的引物，我们只需选择惩罚得分最低的引物。然而，我们现在将使用BLAST来尝试估计每个引物的潜在脱靶引物。我们将利用这一点来为每个基因选择最好的引物，即脱靶最少的引物。

4爆炸引物

创建一个包含所有注释转录本和染色体序列的BLAST数据库需要几分钟(大约2Gb的空闲存储空间)。对于这个例子，我们可以生成一个有限的BLAST数据库，只包含所有靶基因的转录本和11号染色体的序列。

# 11号染色体序列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)

在实际应用中，可以生成一个包含整个基因组和转录组的数据库。数据库可以保存在一个位置，用于多个引物设计，而不必每次都重新构建它。

库(BSgenome) #人类基因组BSgenome对象(需要从Bioconductor安装)hg38 <- getBSgenome("BSgenome. hspiens . ucsc .hg38") #下载和导入gencode hg38注释url <- "ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz" annot <- import(url，# create blast数据库blastdb <- file.path(tempdir()， "blastdb") createBLASTDb(genome = hg38, annot = tx_exons, blastdb = blastdb)

一旦我们有了BLAST数据库，我们就可以用它来估计所有引物的外显子、内含子和基因间脱靶的数量。请注意，外显子和内含子脱靶的数量是在基因水平上提供的，即一个外显子脱靶意味着一个引物可能与另一个基因的外显子区结合。如果脱靶结合位点与两个基因的外显子重叠，它将被计算两次，因为它可以与两个基因的转录本结合。

#现在我们可以爆炸我们的外部引物对创建数据库outer_primers <- blastPrimers(outer_primers, blastdb = blastdb, max_mismatch = 0, min_aligned = 0.75) #引物现在包含估计脱靶tapseq_primers(outer_primers$IFITM1)

为了最终确定我们的外部引物，我们想要选择每个目标基因的最佳引物，即具有最少的外显子、内含子和基因间脱靶(按此顺序)的引物。

#根据潜在的off-targets数量选择每个目标基因的最佳引物best_outer_primers <- pickPrimers(outer_primers, n = 1, by = "off_targets") #每个对象现在只包含最佳引物tapseq_primers(best_outer_primers$IFITM1)

5内底漆设计

为了设计嵌套内引物，我们只需要在较小的产品尺寸范围内重复相同的过程。

#创建新的TsIO对象内部引物,注意不同产品尺寸inner_primers < - TAPseqInput (txs_seqs, target_annot = truncated_txs product_size_range = c(150、300),primer_num_return = 5) #设计内部引物inner_primers < - designPrimers (inner_primers) #高炉内部引物inner_primers < - blastPrimers (inner_primers, blastdb = blastdb max_mismatch = 0, min_aligned = 0.75) #选择最佳引物/目标基因best_inner_primers < - pickPrimers (inner_primers n = 1,By = "off_targets")

完成了!我们成功地为我们的目的基因面板设计了TAP-seq外、内引物集。

6多路复用的兼容性

作为额外的步骤，我们可以验证所设计的引物集是否兼容PCR多路复用。为此，我们使用Primer3的" check_primers "功能:

#使用checkPrimers为每个可能的引物对运行Primer3的"check_primers"功能outer_comp <- checkPrimers(best_outer_primers) inner_comp <- checkPrimers(best_inner_primers)

例如，我们现在可以绘制出每一对的估计互补性分数。我们突出了得分高于47的配对，这在引物设计过程中被Primer3认为是“关键”的(参见Primer3了解更多信息)。

库(dplyr)库(ggplot2) #内外互补分数合并到一个data.frame comp < - bind_rows(外= outer_comp内= inner_comp .id =“套”)#添加变量对任何complemetarity得分高于47 comp < - comp % > %变异(high_compl = if_else (primer_pair_compl_any_th > 47 | primer_pair_compl_end_th > 47岁真的=“高”,假=“ok”))% > %变异(high_compl =因素(high_compl水平= c(“好”,“高”)))#情节互补分数ggplot (comp,aes (x = primer_pair_compl_any_th, y = primer_pair_compl_end_th)) + facet_wrap(~集,ncol = 2) + geom_point (aes(颜色= high_compl),α= 0.25)+ scale_color_manual(值= c(“黑色”,“红色”),删除= FALSE) + geom_hline (aes (yintercept = 47),颜色=“darkgray线型=“冲”)+ geom_vline (aes (xintercept = 47),颜色=“darkgray线型=“冲”)+实验室(title =“互补分数TAP-seq底漆组合”,颜色=“互补性”)+ theme_bw ()

互补得分在47以上的引物对均未发现，可用于多重pcr。

7出口引物

为了完成工作流程，我们可以将设计好的引物导出到data.frames中，而data.frames可以被写入.csv文件，方便引物集的存储。

# create data.frames for outer和inner primer outer_primers_df <- primerDataFrame(best_outer_primers) inner_primers_df <- primerDataFrame(best_inner_primers) #生成的data.frames包含所有相关的primer数据头(outer_primers_df)

此外，我们可以用引物结合位点的基因组坐标创建BED轨迹，以便在基因组浏览器中查看。

#创建床跟踪外部和内部引物和自定义颜色outer_primers_track < - createPrimerTrack (best_outer_primers,颜色=“steelblue3”)inner_primers_track < - createPrimerTrack (best_inner_primers,颜色=“goldenrod1”)#输出data.frames包含线在床上格式为每个引物(outer_primers_track)负责人(inner_primers_track) #出口跟踪请全部文件(“写入标准输出,替换文件名)exportPrimerTrack (outer_primers_track,exportPrimerTrack(inner_primers_track, con = "")

或者，可以轻松地为外部和内部引物创建引物轨道，并将其写入一个.bed文件，以便两个引物集可以在相同的轨道中查看，但具有不同的颜色。

exportPrimerTrack(createPrimerTrack(best_outer_primers, color = "steelblue3")， createPrimerTrack(best_inner_primers, color = "goldenrod1")， con = "path/to/primer .bed")

8会话信息

本插图中的所有输出都是在以下条件下产生的:

sessionInfo() #> R version 4.2.1(2022-06-23) #>平台:x86_64-pc-linux-gnu (64-bit) #>运行在:Ubuntu 20.04.5 LTS #> #>矩阵产品:默认#> BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。所以#> LAPACK: /home/biocbuild/bbs-3.16-bio /R/lib/libRlapack。so #> #> locale: #> [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# > [3] LC_TIME=en_GB LC_COLLATE= c# > [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 #> [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# > [9] LC_ADDRESS=C LC_TELEPHONE= c# > [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION= c# > #>附加基础包:#> [1]stats4 stats graphics grDevices utils datasets methods #>[8]基础#> #>其他附加包:#> [11] IRanges_2.32.0 S4Vectors_0.36.0 #> [13] BiocGenerics_0.44.0 TAPseq_1.10.0 #> [15] BiocStyle_2.26.0 #> #>通过命名空间加载(且未附加):#> [1] bitops_1.0-7 matrixstats_0.0.2 #> b[5] progress_1.2.2 httr_1.4.4 #> [7] tools_4.2.1 bslib_0.4.0 #> [9] utf8_1.2.2 R6_2.5.1 #> [11] colorspace_2.0-3 DBI_1.1.3 #> [13] withr_2.5.0 tidyselect_1.2.0 #> [15] prettyunits_1.1.1 bit_4.0.4 #> [17] curl_4.3.3 compiler_4.2.1 #> [19] cli_3.4.1 Biobase_2.58.0 #> [21] xml2_1.3.3 DelayedArray_0.24.0 #> [23] labeling_0.4.2 bookdown_0.29 #> [25] sass_0.4.2 scales_1.2.1 #> [27] rappdirs_0.3.3 string_1 .4.1 #> [29]> [33] htmltools_0.5.3 MatrixGenerics_1.10.0 #> [31] rmarkdown_2.17 pkgconfig_2.0.3 #> [35] highr_0.9 dbplyr_2.2.1 #> [37] fastmap_1.1.0 rlang_1.0.6 #> [39] RSQLite_2.2.18 farver_2.1.1 #> [41] jquerylib_0.1.4 BiocIO_1.8.0 #> [43] generics_0.1.3 jsonlite_1.8.3 #> [45] RCurl_1.98-1.9 magrittr_2.0.3 #> [47] GenomeInfoDbData_1.2.9 Matrix_1.5-1 #> [49] Rcpp_1.0.9 munsell_0.5.0 #> b[51] fansi_1.0.3 lifecycle_1.0.3 #> [53] stringi_1.7.8 yaml_2.3.6 #> [55]summary - edexperiment_1 .28.0 zlibbioc_1.44.0 #> [57] biocfilecache_2.0 grid_4.2.1 #> [61] crayon_1.5.2 lattice_0.20-45 #> [63] GenomicFeatures_1.50.0 hms_1.1.2 #> [65] KEGGREST_1.38.0 magick_2.7.3 #> [67] knitr_1.40 pillar_1.8.1 #> [69] rjson_0.2.21 codetools_0.2-18 #> [71] biomaRt_2.54.0 XML_3.99-0.12 #> [73] BiocManager_1.30.19 png_0.1-7 #> [79] vctrs_0.5.0 gtable_0.3.1 #> [79] purrr_0.3.5 tidyr_1.2.1 #> [81] assertthat_0.2.1 #cachem_1.0.6 #> [83] xfun_0.34 restfulr_0.0.15 #> [85] tibble_3.1.8 genome alignments_1.34.0 #> [87] AnnotationDbi_1.60.0 memoise_2.0.1 #> [89] ellipsis_0.3.2