内容

1关于APAlyzer

APAlyzer是一个使用RNA测序数据进行选择性多聚腺苷酸(APA)事件生物信息学分析的工具包。我们的主要方法是比较PolyA_DB数据库(https://exon.apps.wistar.org/PolyA_DB/v3/）(Wang et al. 2017,2018)．目前的版本使用RNA-seq数据来检查3 '非翻译区域(3 ' utr)和内含子中的APA事件。编码区域用于基因表达计算。

2程序安装

APAlyzer应使用BiocManager安装:

如果(!"BiocManager" %in% rownames(installed.packages()) install.packages("BiocManager") BiocManager::install("APAlyzer")

或者，它也可以安装如下:

R CMD安装APAlyzer.tar.gz

此外，用户还可以直接从GitHub安装APAlyzer开发版:

BiocManager::安装(“RJWANGbioinfo / APAlyzer”)

安装后，APAlyzer可用于:

库(APAlyzer)

3.样本数据和PAS参考

3．1RNA-seq BAM文件

该包读取BAM文件以获得不同基因组区域的读覆盖信息。的示例中指定示例BAM文件的路径TBX20BamSubset(Bindreither 2018)包数据。在本例中，BAM文件对应于鼠标RNA-seq数据(映射到mm9)。

suppressMessages(library(" tbx20bam子集"))suppressMessages(library("Rsamtools")) flsall = getBamFileList() flsall
## SRR316184 /home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316184“SRR316185 ##”/home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316185。“SRR316186 ##”/home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316186。SRR316187 /home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316187“## SRR316188 ##”/home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316188“SRR316189 ##”/home/biocbuild/bbs-3.17-bioc/R/library/ tbx20bam子集/extdata/SRR316189.bam”

3.2基因组参考

我们的包装需要基因组中的PAS引用(包括3 ' utr和内含子)。我们已经预先构建了小鼠基因组(mm9)的参考文件，可以从该文件加载extdata：

library(repmis) URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="mm9_REF. "RData”source_data (paste0 (URL,文件,“?生= True”))
##从https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/mm9_REF.RData?raw=True下载数据
下载的数据文件的SHA-1哈希值为:## fa117a77976a6931e4755a3182544df1d1de7eb6
## [1] " reftrraw " "dfIPA" " dflle "

这extdata覆盖3 ' utr APA区域(refUTRraw)、IPA区域(dfIPA)和3 ' -most外显子区域(dfLE)。的refUTRraw为包含6列3'UTR PASs基因组信息的数据帧:

头(refUTRraw, 2)
## gene_symbol铬链第一最后cdsend ## 1 Xkr4 chr1 - 3204561 3189814 3206103 ## 2 Lypla1 chr1 + 4835237 4836816 4835097

dfIPA为内含子PASs的8列数据帧;“上游ss”表示距离IPA最近的5 '或3 '拼接点，“下游ss”表示距离IPA最近的3 '拼接点:

头(dfIPA, 2)
## PASid gene_symbol Chrom链Pos上行ss下行ss ## 5 chr4:32818455:+ Mdn1 chr4 + 32818455 32818305 32818889 ## 10 chr4:32837613:+ Mdn1 chr4 + 32837613 32837452 32837887

dfLE是包含5列3 '最小外显子的数据帧;“LEstart”表示最后3’外显子的开始基因组位置。

头(dfLE, 2)
## gene_symbol铬链LEstart TES ## 1 0610009O20Rik chr18 + 38421534 38425118 ## 2 0610010F05Rik chr11 - 23465136 23462079

除了小鼠mm9之外，我们的套件还有小鼠mm10、人类hg38和人类hg19基因组的预构建版本:

URL = " https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/ "文件= " hg19_REF。RData”source_data (paste0 (URL,文件,“?生= True”))
##从https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/hg19_REF.RData?raw=True下载数据
下载数据文件的## SHA-1哈希值为:## 7ba8e91beba626318e48736692d6442480d1ef06
## [1] "refUTRraw_hg19" "dfIPA_hg19" "dfLE_hg19"

更多的构建前参考可以在参考和测试数据repo中找到:https://github.com/RJWANGbioinfo/PAS_reference_RData_and_testing_data

3．3同时建立3'UTR和内含子PAS参考区

为了量化PAS的相对表达，我们需要为它们构建参考区域，尽管这可以在以前的版本中单独构建。我们还提供了一个新函数REF4PAS从APAlyzer 1.2.1开始立即构建这些区域:

refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19')，] dfIPAraw=dfIPA[which(dfIPA$Chrom=='chr19')，] dfLEraw= dflle [which(dfIPA$Chrom=='chr19')，] PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw) UTRdbraw=PASREF$UTRdbraw dfIPA=PASREF$dfIPA dflle =PASREF$dfIPA dflle

在这种情况下，UTRdbraw为3'UTR APA分析的参考区域，而dfIPA而且dfLE在内含子APA分析中是必需的。

3.4使用GTF文件构建3'UTR PAS和IPA参考

尽管我们强烈建议用户使用从PolyA_DB生成的引用区域。从APAlyzer 1.2.1开始，我们还提供了一个新的函数，可以帮助用户直接从基因注释GTF文件中建立他们的参考，我们希望这可以帮助到尚未被PolyA_DB覆盖的物种:

##建立参考范围3'UTR PASs在鼠标下载。file(url='ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.gtf.gz'， destfile=' mus_muscle . grcm38.99 .gtf.gz') GTFfile=" mus_muscle . grcm38.99 .gtf.gz" PASREFraw=PAS2GEF(GTFfile) refUTRraw=PASREFraw$refUTRraw dfIPAraw=PASREFraw$dfIPA dfLEraw=PASREFraw$ dfIPA dfLEraw= pasref4pas (refUTRraw,dfIPAraw,dfLEraw)

从APAlyzer 1.9.1开始，可以设置两种注释方法AnnoMethod =要注释PAS， ' legacy '或' V2 '(默认)，' legacy '是以前版本中使用的旧方法，而' V2 '是更新和改进的版本。

43'UTRs的APA分析

4．1构建aUTR和cUTR引用

为了计算3'UTR APA相对表达量(RE)，我们首先需要使用定义aUTR和cUTR的参考区域refUTRraw．由于示例数据只包含chr19上的映射信息，我们可以只放大到chr19上的参考区域:

refUTRraw = refUTRraw ((refUTRraw铬美元= =“chr19”),]UTRdbraw = REF3UTR (refUTRraw)

的REF3UTR函数返回每个基因包含aUTR(pPAS到dPAS)和cUTR(cdsend到pPAS)区域的基因组范围:

头(UTRdbraw, 2)
## seqnames ranges strand | gene_symbol gene_name ##    |   ## 10082 chr19 5492231-5495277 - | 1700020D05Rik 1700020D05Rik__aUTR ## ------- ## seqinfo: 21个来自未指定基因组的序列;## seqnames ranges strand | gene_symbol gene_name ##    |   ## 100821 chr19 5495277-5502867 - | 1700020D05Rik 1700020D05Rik__cUTR ## ------- # seqinfo: 21个序列来自一个未指定的基因组;没有seqlengths

4.2相对表达式计算

一旦确定了cUTR和aUTR区域，就可以计算出每个基因的3'UTR APA的REPASEXP_3UTR：

DFUTRraw=PASEXP_3UTR(UTRdbraw, flsall, Strandtype="forward")
##[1]“SRR316184, Strand: forward, finished“##[1]”SRR316185, Strand: forward, finished“##[1]”SRR316186, Strand: forward, finished“##[1]”SRR316187, Strand: forward, finished“##[1]”SRR316188, Strand: forward, finished“##[1]”SRR316189, Strand: forward, finished“##[1]”SRR316189, Strand: forward, finished

的PASEXP_3UTR3UTR需要两个输入:1)aUTR和cUTR参考区域，2)BAM文件路径。除了输入之外，还可以使用定义序列的链Strandtype．详细的使用方法还可以通过命令获取PASEXP_3UTR ?.输出数据帧包含每个基因的reads count(以aUTR或cUTR为单位)、RPKM(以aUTR或cUTR为单位)和RE (log2(aUTR/cUTR)):

头(DFUTRraw, 2)
# # 1 # # gene_symbol SRR316184_areads SRR316184_aRPKM SRR316184_creads 1700020 d05rik 0 0 31 # # 2 1810009 a15rik 0 0 0 # # SRR316184_cRPKM SRR316184_3UTR_RE SRR316185_areads SRR316185_aRPKM # 118.4358 # 1负0 0 # # 2 0.0000南0 0 # # SRR316185_creads SRR316185_cRPKM SRR316185_3UTR_RE SRR316186_areads 20 74.34672负0 # # 1 # # 2 0 0.00000南0 # # SRR316186_aRPKM SRR316186_creads SRR316186_cRPKM SRR316186_3UTR_RE # # 1 0 23 87.90994负无穷到# # 2 0 0 0.00000南# # SRR316187_areads SRR316187_aRPKMSRR316187_creads SRR316187_cRPKM 33 # # 1 0 0 0.00000 88.85001 # # 2 0 0 0 # # SRR316187_3UTR_RE SRR316188_areads SRR316188_aRPKM SRR316188_creads 33 # # 1负0 0 # # 2南0 0 0 # # SRR316188_cRPKM SRR316188_3UTR_RE SRR316189_areads SRR316189_aRPKM # 77.68085 # 1负0 0 # # 2 0.00000南0 0 # # SRR316189_creads SRR316189_cRPKM SRR316189_3UTR_RE # 31.76981 # 1 11负无穷到# # 2 0 0.00000南

5内含子APA分析

5.1构建内含子polyA引用

IPA的分析需要两个基因组区域:IPA区域和3 ' -most外显子。如上所述，小鼠和人类基因组中的这些区域已经预先构建在包中:

#鼠标(mm9): URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="mm9_REF. "RData”source_data (paste0 (URL,文件,“?生= True”))
##从https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/mm9_REF.RData?raw=True下载数据
下载的数据文件的SHA-1哈希值为:## fa117a77976a6931e4755a3182544df1d1de7eb6
## [1] " reftrraw " "dfIPA" " dflle "
#human(hg19): URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="hg19_REF. "RData”source_data (paste0 (URL,文件,“?生= True”))
##从https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/hg19_REF.RData?raw=True下载数据
下载数据文件的## SHA-1哈希值为:## 7ba8e91beba626318e48736692d6442480d1ef06
## [1] "refUTRraw_hg19" "dfIPA_hg19" "dfLE_hg19"

5.2相对表达式计算

与3'UTR APA类似，IPAs的RE可由PASEXP_IPA计算:PASEXP_IPA：

dfIPA=dfIPA[which(dfIPA$Chrom=='chr19')，] dfIPA=dfIPA[which(dfIPA$Chrom=='chr19')，] IPA_OUTraw=PASEXP_IPA(dfIPA, dfIPA, flsall, Strandtype="forward"， nts=1)

注意，作为IPA的一个特定特性，可以使用' nts= '(传递给Rsubread:: featurerts的参数，检查)设置更多线程? Rsubread:: featureCounts以提高计算速度。详细用法可通过命令获取PASEXP_IPA ?．

输出数据帧包含IPA上游(a)、IPA下游(b)和3 ' -most外显子区域(c)的读计数和读密度。IPA的RE计算为log2((a - b)/c)。

头(IPA_OUTraw, 2)
# # 1 # # gene_symbol IPAID PASid 9130011 e15rik chr19:45894071: -9130011 e15rik chr19:45894071: - # # 2 9130011 e15rik chr19:45952090: -9130011 e15rik chr19:45952090: - - - - - - # # SRR316184_IPA_UPreads SRR316184_IPA_UPRPK SRR316184_IPA_DNreads 1 1.445087 0 # # 2 # # 104 12.453598 26 # # SRR316184_IPA_DNRPK SRR316184_LEreads SRR316184_LERPK SRR316184_IPA_RE # # 96 0.0000000 96 59.47955 - -5.363166 # # 2 0.8276828 59.47955 - -2.355049 # # SRR316185_IPA_UPreads SRR316185_IPA_UPRPK SRR316185_IPA_DNreads 0.00000 - 1 # # # # 1 0212114．4892820 ## SRR316185_IPA_DNRPK SRR316185_LEreads SRR316185_LERPK SRR316185_IPA_RE ## 1 1.6611296 96 59.47955 0.000000 ## 2 0.6366791 96 59.47955 -2.102237 ## SRR316186_IPA_UPreads SRR316186_IPA_UPRPK SRR316186_IPA_DNreads ## 1 0 0.00000 0 ## 2 108 12.93258 13 ## SRR316186_IPA_DNRPK SRR316186_LEreads SRR316186_LERPK SRR316186_IPA_RE ## 1 0.0000000 69 42.75093 0.000000 ## 2 0.4138414 69 42.75093 -1.771866 ## SRR316187_IPA_UPreads SRR316187_IPA_UPRPK SRR316187_IPA_DNreads ## 1 0 0.00000 0 ## 2 199 23.82948 22 ## SRR316187_IPA_DNRPK SRR316187_LEreads SRR316187_LERPK SRR316187_IPA_RE ## 1 0.000000 40 24.78315 0.00000000 ## 2 0.700347 40 24.78315 -0.09964814 ## SRR316188_IPA_UPreads SRR316188_IPA_UPRPK SRR316188_IPA_DNreads ## 1 1 1.445087 0 ## 2 207 24.787451 15 ## SRR316188_IPA_DNRPK SRR316188_LEreads SRR316188_LERPK SRR316188_IPA_RE ## 1 0.0000000 55 34.07683 -4.5595631 ## 2 0.4775093 55 34.07683 -0.4872446 ## SRR316189_IPA_UPreads SRR316189_IPA_UPRPK SRR316189_IPA_DNreads ## 1 1 1.445087 0 ## 2 226 27.062627 15 ## SRR316189_IPA_DNRPK SRR316189_LEreads SRR316189_LERPK SRR316189_IPA_RE ## 1 0.0000000 46 28.50062 -4.3017653 ## 2 0.4775093 46 28.50062 -0.1003744

6APA事件的显著性分析

一旦获得了每个样本的读覆盖率信息，就可以比较两个不同组之间的APA调控差异。在本分析中，有两种类型的实验设计:1)无重复;2)重复。根据设计生成样例表:

#使用复制sampleTable1 = data.frame(samplename = c(names(flsall))， condition = c(rep("NT"，3)，rep("KD"，3)))
## samplename条件## 1 SRR316184 NT ## 2 SRR316185 NT ## 3 SRR316186 NT ## 4 SRR316187 KD ## 5 SRR316188 KD ## 6 SRR316189 KD
= data.frame(samplename = c("SRR316184"，"SRR316187")， condition = c("NT"，"KD")
## SRR316184 NT ## 2 SRR316187 KD

6.13 ' utr显著调节APA

为分析无重复的样本(例中为KD组和NT组)之间的3'UTR APA，sampleTable2使用。这里使用的函数被调用APAdiff(详细信息可通过命令获取APAdiff ?)．它将首先通过样本表来确定它是复制设计还是非复制设计。然后进行APA同情。

test_3utsing =APAdiff(sampleTable2,DFUTRraw, conKET='NT'， trtKEY='KD'， PAS='3UTR'， CUTreads=0, p_adjust_methods="fdr")

的APAdiff函数需要两个输入:1)定义样本分组/条件的样本表，2)读取autr和cutr的覆盖率信息，可以通过PASEXP_3UTR从前面的步骤。组名，即治疗或控制，可以定义为btrtKEY =而且conKET =；所分析的PAS类型应定义为不是=；用于aUTR和cUTR的读截止由CUTreads =默认值为0;定义了p值校正的方法p_adjust_methods =默认值为' fdr '。在非复制设计中，APA模式将在两个样本之间进行比较，输出将显示在一个数据帧中:

头(test_3UTRsing, 2)
## gene_symbol RED pvalue p_adj APAreg ## 3 1810055G02Rik -0.1195816 1.0000000 1.0000000 NC ## 4 2700081O15Rik 0.3862662 0.3292784 0.8007581 NC
表(test_3UTRsing APAreg美元)
## ## dn nc up ## 31 246 8 .

输出包含4列:“基因符号”描述基因信息;RED为两组间的相对表达差异;“pvalue”是基于Fisher精确检验的统计显著性;' p_adj '为FDR调节的pvalue， ' APAreg '为基因中3'UTR APA调节模式。我们在“APAreg”中定义了3种类型，“UP”表示治疗组的aUTR丰度(在本例中为“KD”)比对照组(在本例中为“NT”)至少高5%，且“pvalue”<0.05;“DN”表示处理组的aUTR丰度比对照组低5%，p值<0.05;' NC '是剩下的基因。对于3'UTR大小变化，' UP '表示3'UTR变长，' DN ' 3'UTR变短。

对于重复设计，我们使用t检验进行显著性分析。然而，其他基于负二项数据分布的工具，如
DEXSeq(安德斯、雷耶斯和W. 2012)也可以用。

使用多重重复设计分析KD和NT组之间的3'UTR APA test_3UTRmuti=APAdiff(sampleTable1, DFUTRraw, conKET='NT'， trtKEY='KD'， PAS='3UTR'， CUTreads=0, p_adjust_methods="fdr"， MultiTest='unpaired t-test') head(test_3UTRmuti,2)
## gene_symbol RED pvalue p_adj APAreg ## 3 1810055G02Rik -0.7631689 0.08590526 0.4961374 NC ## 4 2700081O15Rik 0.1872116 0.07691742 0.4961374 NC
表(test_3UTRmuti APAreg美元)
## ## dn nc up ## 24 226

在复制设计中，附加参数MultiTest可用于设置统计方法。在当前版本中，有3种方法可以设置:“unpaired t-test”(默认)，“paired t-test”，“ANOVA”。运行之后APAdiff，会生成几列:‘RED’是两组平均相对表达式的差值;' pvalue '是通过定义的方法计算的p值MultiTest．“APAreg”是预测的APA调控方向;' UP '定义为' RED ' >0且' pvalue ' <0.05;而' DN '是相反的;' NC '是剩下的基因。

6.2在内含子中显著调控APA

IPA比较与APA使用的3'UTR相似APAdiff(1)使用IPA表达式作为输入，(2)需要将' PAS= '定义为' IPA '，(3)对每个IPA进行分析。值得注意的是，IPA的调控方向与3'UTR APA的调控方向相似。这意味着“UP”被定义为IPA的上调(RED > 0);“DN”是反义词;' NC '是剩下的基因。

无重复的KD组和NT组之间的IPA分析如下:

test_IPAsing=APAdiff(sampleTable2, IPA_OUTraw, conKET='NT'， trtKEY='KD'， PAS='IPA'， CUTreads=0, p_adjust_methods="fdr") head(test_IPAsing,2)
## gene_symbol PASid RED pvalue p_adj APAreg ## 2 9130011E15Rik chr19:45952090:- 2.255401 1.790010e-12 3.615821e-10 UP ## 5 Ablim1 chr19:57151655:- 0.626637 5.392007e-01 1.000000e+00 NC

使用复制数据分析KD组和NT组之间的IPA如下:

test_IPAmuti=APAdiff(sampleTable1, IPA_OUTraw, conKET='NT'， trtKEY='KD'， PAS='IPA'， CUTreads=0, p_adjust_methods="fdr"， MultiTest='unpaired t-test') head(test_IPAmuti,2)
## gene_symbol PASid RED pvalue p_adj APAreg ## 2 9130011E15Rik chr19:45952090:- 1.84729507 0.001294432 0.2808917 UP ## 5 Ablim1 chr19:57151655:- -0.09466715 0.900027947 0.9480877 NC

7分析结果的可视化

从APAlyzer 1.2.1开始，我们提供了两个新的函数APAVolcano而且APABox供用户使用火山图和箱形图绘制他们的RED结果。在火山图中，用户还可以标记使用最多的基因顶级=或一组特定的基因使用markergenes =，例如:

apolcano (test_3utsing, PAS='3UTR'， Pcol = "pvalue"， top=5, main='3UTR APA')

在箱形图中，红色位于' UP '、' DN '和' NC '基因上:

APABox(test_3UTRsing, xlab = "APAreg"， ylab = "RED"， plot_title = NULL)

除了火山图和箱形图外，APA比较结果还可以用箱形图、小提琴图或CDF曲线绘制。对于前面的3'UTR APA和IPA比较输出，首先需要构建绘图数据帧:

test_3UTRmuti APA =“3 'utr”test_IPAmuti APA美元=“音标”dfplot = rbind (test_3UTRmuti [c(“红”,APA)], test_IPAmuti [c(“红”,APA)])

制作小提琴图和CDF曲线ggplot2：

库(ggplot2)
###小提琴ggplot(dfplot, aes(x = APA, y = RED)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2)+ theme_bw() + geom_hline(yintercept=0, linetype="虚线"，color ="红色")

###CDF ggplot(dfplot, aes(x = RED, color = APA)) + stat_ecdf(geom =" step") + ylab("累积分数")+ geom_vline(xintercept=0, linetype="虚线"，color ="灰色")+ theme_bw() + geom_hline(yintercept=0.5, linetype="虚线"，color ="灰色")

8利用编码区进行基因表达分析

APA经常参与基因表达调控。为了比较不同样本中的基因表达与APA，我们的包提供了一个简单的函数，使用映射到编码序列的RNA-seq读取来评估表达变化。

8.1构建编码区域参考

suppressMessages(library("GenomicFeatures")) suppressMessages(library("org.Mm.eg.db")) extpath = system。file("extdata"， "mm9.chr19.refGene.R.DB"， package="APAlyzer") txdb=loadDb(extpath, packageName='GenomicFeatures') IDDB = org.Mm.eg.db CDSdbraw=REFCDS(txdb,IDDB)
## 'select()'返回键和列之间的many:1映射

8.2表达式计算

DFGENEraw=GENEXP_CDS(CDSdbraw, flsall, Strandtype="forward")
##[1]“SRR316184, Strand: forward, finished“##[1]”SRR316185, Strand: forward, finished“##[1]”SRR316186, Strand: forward, finished“##[1]”SRR316187, Strand: forward, finished“##[1]”SRR316188, Strand: forward, finished“##[1]”SRR316189, Strand: forward, finished“##[1]”SRR316189, Strand: forward, finished

9从对端bam文件中提取3'most对齐

从APAlyzer 1.5.5开始，我们提供了一个名为ThreeMostPairBam供用户从配对端的bam文件中提取三个质数对齐，并将其保存到新的bam文件中。例如:

##提取一个成对端## bam文件的3 '大部分对齐，并保存到一个新的bam文件Bamfile='/path/to/inputdir/input。bam' Outdir='/path/to/ Outdir /' StrandType="forward-reverse" ## "forward-reverse"，或"reverse-forward"或"NONE" ThreeMostPairBam (BamfilePath=Bamfile, OutDirPath=Outdir, StrandType='forward-reverse')

输出bamfile(一个名为*.3most的bam文件。Bam位于Outdir)可作为输入文件PASEXP_3UTR．将此bam文件用于PASEXP_IPA，用户需要设置SeqType到“ThreeMostPairEnd”，例如:

ThreemostBamfile = " test.3most。bam" IPA_OUTraw=PASEXP_IPA(dfIPA, dfLE, ThreemostBamfile, Strandtype="forward"， SeqType="ThreeMostPairEnd")

10完整分析示例:小鼠睾丸与心脏的APA分析

10.1关于这个数据集

为了使用我们的软件包提供APA分析的完整教程，我们现在已经准备了一个测试数据集，通过向下采样小鼠心脏(GSM900193)和睾丸(GSM900199)的RNA-Seq数据:

样品标识 SRRID 样品名称 下采样读数
GSM900199 SRR453175 Heart_Rep1 500万年
GSM900199 SRR453174 Heart_Rep2 500万年
GSM900199 SRR453173 Heart_Rep3 500万年
GSM900199 SRR453172 Heart_Rep4 500万年
GSM900193 SRR453143 Testis_rep1 500万年
GSM900193 SRR453142 Testis_rep2 500万年
GSM900193 SRR453141 Testis_rep3 500万年
GSM900193 SRR453140 Testis_rep4 500万年

10.2下载bam文件

download_testbam () flsall < - dir (getwd(),“本”)flsall < -paste0 (getwd(),“/”,flsall)名称(flsall) < -gsub(‘本’,”迪尔(getwd(),“本”))

10.3构建PAS参考区域

library(repmis) URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="mm9_REF. "RData”source_data (paste0 (URL,文件,“?生= True”))PASREF=REF4PAS(refUTRraw,dfIPA,dfLE) UTRdbraw=PASREF$UTRdbraw dfIPA=PASREF$dfIPA dfLE=PASREF$dfLE

10.43'UTR APA和IPA相对表达量的计算

UTR_APA_OUT=PASEXP_3UTR(UTRdbraw, flsall, Strandtype="invert") IPA_OUT=PASEXP_IPA(dfIPA, dflle, flsall, Strandtype="invert"， nts=1)

10.53 ' utr显著调节APA

############# 3utr APA ################# sampleTable = data.frame(samplename = c('Heart_rep1'， 'Heart_rep2'， 'Heart_rep3'， 'Heart_rep4'， 'Testis_rep1'， 'Testis_rep2'， 'Testis_rep3'， 'Testis_rep4')， condition = c(rep("Heart"，4)， rep("Testis"，4))) test_3UTRAPA=APAdiff(sampleTable,UTR_APA_OUT, conKET='Heart'， trtKEY='Testis'， PAS=' 3utr '， CUTreads=5)
表(test_3UTRAPA APAreg美元)
## ## dn nc up ## 490 719 80
apvolcano (test_3UTRAPA, PAS='3UTR'， Pcol = "pvalue"， plot_title='3UTR APA')

APABox(test_3UTRAPA, xlab = "APAreg"， ylab = "RED"， plot_title = NULL)

10.6内含子显著调控APA

############# 异丙醇  ################# test_IPA = APAdiff (sampleTable IPA_OUT conKET =‘心’,trtKEY =“睾丸”,不等于“音标”,CUTreads = 5)
表(test_IPA APAreg美元)
## ## dn nc up ## 125 924 312
apvolcano (test_IPA, PAS='IPA'， Pcol = "pvalue"， plot_title='IPA')

APABox(test_IPA, xlab = "APAreg"， ylab = "RED"， plot_title = NULL)

11常见问题

如何生成BAM文件列表进行分析?

包含BAM文件名和文件路径的BAM文件列表。假设所有BAM文件都存储在bamdir中，那么BAM文件列表可以通过以下方式获得:

Flsall =dir(bamdir，".bam") Flsall =paste0(bamdir, Flsall) names(Flsall)=dir(bamdir，".bam")

为什么我得到错误消息时，我试图获得txdb使用makeTxDbFromUCSC?

您可以尝试升级您的Bioconductor，或使用GTF加载基因组注释，或使用' . r.b b '文件加载预构建基因组注释，例如mm9. refgene . r.b b。

12会话信息

会话信息记录了当前文档生成过程中使用的所有包版本。

sessionInfo ()
## R正在开发中(不稳定)(2022-10-25 r83175) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 22.04.1 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas。so ## LAPACK: /usr/lib/x86_64-linux-gnu/ LAPACK /liblapack.so.3.10.0 ## ## 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_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:[4] Biostrings_2.67.0 XVector_0.39.0 genomicranges_1.1.1.0 ## [16] GenomeInfoDb_1.35.0 IRanges_2.33.0 S4Vectors_0.37.0 ## [16] BiocGenerics_0.45.0 APAlyzer_1.13.0 BiocStyle_2.27.0 ## ##通过命名空间加载(并且没有附加):# # # # [1] DBI_1.1.3 bitops_1.0-7 [3] biomaRt_2.55.0 rlang_1.0.6 # # [5] magrittr_2.0.3 matrixStats_0.62.0 # # [7] compiler_4.3.0 RSQLite_2.2.18 # # [9] png_0.1-7 vctrs_0.5.0 # # [11] stringr_1.4.1 pkgconfig_2.0.3 # # [13] crayon_1.5.2 Rsubread_2.13.0 # # [15] fastmap_1.1.0 magick_2.7.3 # # [17] dbplyr_2.2.1 ellipsis_0.3.2 # # [19] labeling_0.4.2 utf8_1.2.2 # # [21] rmarkdown_2.17 purrr_0.3.5 # # [23] bit_4.0.4 xfun_0.34 # # [25] zlibbioc_1.45.0 cachem_1.0.6 # # [27] jsonlite_1.8.3 progress_1.2.2 # # [29]blob_1.2.3 highr_0.9 # # [31] DelayedArray_0.25.0 BiocParallel_1.33.0 # # [33] HybridMTest_1.43.0 parallel_4.3.0 # # [35] prettyunits_1.1.1 VariantAnnotation_1.45.0 # # [37] R6_2.5.1 bslib_0.4.0 # # [39] stringi_1.7.8 RColorBrewer_1.1-3 # # [41] rtracklayer_1.59.0 genefilter_1.81.0 # # [43] jquerylib_0.1.4 Rcpp_1.0.9 # # [45] bookdown_0.29 assertthat_0.2.1 # # [47] SummarizedExperiment_1.29.0 R.utils_2.12.1 # # [49] R.cache_0.16.0 Matrix_1.5-1 # # [51] splines_4.3.0 tidyselect_1.2.0 # # [53] yaml_2.3.6codetools_0.2-18 # # [55] curl_4.3.3 plyr_1.8.7 # # [57] lattice_0.20-45 tibble_3.1.8 # # [59] withr_2.5.0 KEGGREST_1.39.0 # # [61] evaluate_0.17 survival_3.4-0 # # [63] BiocFileCache_2.7.0 xml2_1.3.3 # # [65] pillar_1.8.1 BiocManager_1.30.19 # # [67] filelock_1.0.2 MatrixGenerics_1.11.0 # # [69] generics_0.1.3 rcurl_1.98 - 1.9 # # [71] hms_1.1.2 munsell_0.5.0 # # [73] scales_1.2.1 xtable_1.8-4 # # [75] glue_1.6.2 tools_4.3.0 # # [77] BiocIO_1.9.0 data.table_1.14.4 # # [79] BSgenome_1.67.0 annotate_1.77.0 # #[81] locfit_1. 1.1 - 2.21 genome alignments_1.35.0 ## [83] XML_3.99-0.12 grid_4.3.0 ## [85] tidyr_1.2.1 colorspace_2.0-3 ## [87] GenomeInfoDbData_1.2.9 restfulr_0.0.15 ## [89] cli_3.4.1 rappdirs_0.3.3 ## [91] fansi_1.0.3 dplyr_1.0.10 ## [93] gtable_0.3.1 r.d osss3_1 .8.2 ## [95] DESeq2_1.39.0 sass_0.4.2 ## [97] digest_0.6.30 ggrepel_0.9.1 ## [99] farver_2.1.1 geneplotter_1.77.0 ## [101] rjson_0.2.21 R.oo_1.25.0 ## [103] memoise_2.0.1 htmltools_0.5.3 ## [105] lifecycle_1.0.3 httr_1.4.4 ## [107]bit64_4.0.5

13确认

我们感谢滨天实验室成员的有益讨论和包测试。

参考文献

安德斯，S.， A.雷耶斯，Huber W. 2012。从Rna-Seq数据中检测外显子的差异使用基因组研究22日:4025。https://doi.org/10.1101/gr.133744.111．

宾德莱瑟，D. 2018。tbx20bam子集:来自“Tbx20”实验的Bam文件子集．

王锐，R. Nambiar，郑迪，田波。2017。“PolyA_DB 3在多个基因组中通过深度测序鉴定的切割和多聚腺苷酸位点目录。”核酸研究46 (d1): d315-d319。https://doi.org/10.1093/nar/gkx1000．

王锐，郑德栋，田波。2018。哺乳动物基因中保守的卵裂和多聚腺苷酸事件概要。基因组研究28(10): 1427-41。https://doi.org/10.1101/gr.237826.118．