内容

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")

或者，它也可以按以下方式安装:

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
/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316184。/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316185。/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316186。/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316187。/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316188。## SRR316189 ## "/home/biocbuild/ bs-3.15-bio /R/library/ tbx20bam子集/extdata/SRR316189.bam"

3.2基因组参考

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

库(repmis) URL = " https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/ "文件= " mm9_REF。RData”source_data (paste0 (URL,文件,“?生= True”))
从:https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/mm9_REF.RData?raw=True下载数据
## fa117a77976a6931e4755a3182544df1d1de7eb6
## [1] "refUTRraw" "dfIPA" " dfl "

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

头(refUTRraw, 2)
## gene_symbol Chrom Strand First Last cdsend ## 1 Xkr4 chr1 - 3204561 3189814 3206103 ## 2 Lypla1 chr1 + 4835237 4836816 4835097

dfIPA是包含8列Intronic PASs的数据帧;'上游'指离IPA最近的5 '或3 '拼接点，'下游'指离IPA最近的3 '拼接点:

头(dfIPA, 2)
## PASid gene_symbol Chrom Strand Pos upstreamSS downstreamSS ## 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 Chrom Strand 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"

更多的预构建参考可以在参考和测试数据回购中找到: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= dfIPA[which(dfIPA$Chrom=='chr19')，] PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw) UTRdbraw=PASREF$UTRdbraw dfIPA=PASREF$dfIPA dfale =PASREF$dfIPA dfale =PASREF$dfIPA dfale =PASREF$dfIPA dfale)

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

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

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

##建立鼠标中3'UTR PASs的参考范围download.file(url='ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.gtf.gz'， destfile='Mus_musculus.GRCm38.99.gtf.gz') GTFfile="Mus_musculus.GRCm38.99.gtf.gz" PASREFraw=PAS2GEF(GTFfile) refUTRraw= PAS2GEF $refUTRraw dfIPAraw=PASREFraw$dfIPA dfLEraw=PASREFraw$ dfIPA PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw)

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

43 ' utr患者APA分析

4．1构建aUTR和cUTR引用

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

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

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

头(UTRdbraw, 2)
##长度为2的GRangesList对象:## $ ' 1700020D05Rik__aUTR ' ## GRanges对象具有1个范围和2个元数据列:## seqnames ranges strand | gene_symbol gene_name ##    |   ## 10082 chr19 5492231-5495277 - | 1700020D05Rik 1700020D05Rik__aUTR ## ------- ## seqinfo: 21个序列来自一个未指定的基因组;no seqlength## ## $ ' 1700020D05Rik__cUTR ' ## GRanges对象具有1个范围和2个元数据列:## 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 =“转发”)
## [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"

的PASEXP_3UTR3UTR需要两个输入:1)aUTR和cUTR参考区域，2)BAM文件的路径。除了输入之外，还可以使用定义序列的链Strandtype．详细的用法也可以通过命令获取PASEXP_3UTR ?.输出数据帧包括每个基因的读取计数(在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/ "文件= " mm9_REF。RData”source_data (paste0 (URL,文件,“?生= True”))
从:https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/mm9_REF.RData?raw=True下载数据
## fa117a77976a6931e4755a3182544df1d1de7eb6
## [1] "refUTRraw" "dfIPA" " dfl "
#人类(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"

5.2相对表达式的计算

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

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

注意，作为IPA的一个特定特性，可以使用' nts= '(传递给Rsubread:: featucounts的参数，检查)设置更多线程? 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))) sampleTable1 .frame
## samplename条件## 1 SRR316184 NT ## 2 SRR316185 NT ## 3 SRR316186 NT ## 4 SRR316187 KD ## 5 SRR316188 KD ## 6 SRR316189 KD
#不复制样本表sampleTable2 = data.frame(samplename = c("SRR316184"，"SRR316187")， condition = c("NT"，"KD")
## 1 SRR316184 NT ## 2 SRR316187 KD

6.1APA在3 ' utr显著调节

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

#使用非重复设计分析KD和NT组之间的3'UTR APA test_3utsing =APAdiff(sampleTable2,DFUTRraw, conKET='NT'， trtKEY='KD'， PAS='3UTR'， CUTreads=0, p_adjust_methods="fdr")

的APAdiff函数需要两个输入:1)定义样本组/条件的样本表，2)读取aUTRs和cUTRs的覆盖率信息，这些信息可以通过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

输出包含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-value<0.05;' NC '是剩下的基因。对于3'UTR大小的变化，' UP '表示3'UTR变长，' DN ' 3'UTR变短。

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

#使用多副本设计分析KD和NT组之间的3'UTR APA test_3utmuti =APAdiff(sampleTable1, DFUTRraw, conKET='NT'， trtKEY='KD'， PAS='3UTR'， CUTreads=0, p_adjust_methods="fdr"， MultiTest='unpaired t-test') head(test_3utmuti,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 =,例如:

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

在方框图中，RED是在' 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)
### violle ggplot(dfplot, aes(x = APA, y = RED)) + geom_violle (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("cumulative fraction")+ geom_vline(xintercept=0, linetype="虚线"，color ="灰色")+ theme_bw() + geom_hline(yintercept=0.5, linetype="虚线"，color ="灰色")

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

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

8.1构建编码区域引用

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

8.2计算表达式

DFGENEraw = GENEXP_CDS (CDSdbraw、flsall Strandtype =“转发”)
## [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"

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

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

##提取一个成对端## bam文件的3质数最对齐，并保存到一个新的bam文件Bamfile='/path/to/inputdir/input。' Outdir='/path/to/ Outdir /' StrandType="forward-reverse" ## "forward-reverse"， or "reverse-forward" or "NONE"

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

ThreemostBamfile = " test.3most。bam" IPA_OUTraw=PASEXP_IPA(dfIPA, dfl, 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参考区域

库(repmis) URL = " https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/ "文件= " 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, dfl, flsall, Strandtype="invert"， nts=1)

10.5APA在3 ' utr显著调节

############# 3utr APA ################# sampleTable = data.frame(samplename = c('Heart_rep1'， 'Heart_rep2'， 'Heart_rep3'， 'Heart_rep4'， ' tests_rep1 '， ' tests_rep2 '， ' tests_rep3 '， ' tests_rep4 ')， condition = c(rep("Heart"，4)， rep(" tests "，4)) test_3UTRAPA=APAdiff(sampleTable, ut_apa_out, conKET='Heart'， trtKEY=' tests '， PAS=' 3utr '， CUTreads=5)
表(test_3UTRAPA APAreg美元)
## ## dn nc up ## 490 719 80
APAVolcano(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")

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

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

12会话信息

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

sessionInfo ()
## R版本4.2.0 RC (22-04-19 r82224) ##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.4 LTS ## ##矩阵产品:default ## BLAS: /home/biocbuild/bbs-3.15-bio /R/lib/libRblas. ##因此## LAPACK: /home/biocbuild/bbs-3.15-bio /R/lib/libRlapack。因此## ## 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] base ## ##其他附加的包:# [7] repmis_0.5 TBX20BamSubset_1.31.0 Rsamtools_2.12.0 ## [10] Biostrings_2.64.0 XVector_0.36.0 genomics icranges_1 .48.0 ## [13] GenomeInfoDb_1.32.0 IRanges_2.30.0 S4Vectors_0.34.0 ## [16] BiocGenerics_0.42.0 APAlyzer_1.10.0 BiocStyle_2.24.0 ## ##通过命名空间加载(并没有附加):# [9] splines_4.2.0 R.methodsS3_1.8.1 ## [11] cachem_1.0.6 geneplotter_1.74.0 ## [13] jsonlite_1.8.0 annotate_1.74.0 ## [15] dbplyr_2.1.1 png_0.1-7 ## [17] R.oo_1.24.0 BiocManager_1.30.17 ## [19] compiler_4.2.0 httr_1.4.2 ## [21] assertthat_0.2.1 Matrix_1.4-1 ## [23] fastmap_1.1.0 cli_3.3.0 ## [25] htmltools_0.5.2 prettyunits_1.1.1 ## [27] tools_4.2.0gtable_0.3.0 # # [29] glue_1.6.2 GenomeInfoDbData_1.2.8 # # [31] dplyr_1.0.8 rappdirs_0.3.3 # # [33] Rcpp_1.0.8.3 jquerylib_0.1.4 # # [35] vctrs_0.4.1 rtracklayer_1.56.0 # # [37] xfun_0.30 stringr_1.4.0 # # [39] lifecycle_1.0.1 restfulr_0.0.13 # # [41] xml_3.99 - 0.9 zlibbioc_1.42.0 # # [43] scales_1.2.0 BSgenome_1.64.0 # # [45] VariantAnnotation_1.42.0 hms_1.1.1 # # [47] MatrixGenerics_1.8.0 parallel_4.2.0 # # [49] SummarizedExperiment_1.26.0 RColorBrewer_1.1-3 # # [51] yaml_2.3.5 curl_4.3.2 # # [53]## [55] biomaRt_2.52.0 stringi_1.7.6 ## [57] RSQLite_2.2.12 highr_0.9 ## [59] genefilter_1.78.0 BiocIO_1.6.0 ## [61] filelock_1.0.2 BiocParallel_1.30.0 ## [63] Rsubread_2.10.0 rlang_1.0.2 ## [65] pkgconfig_2.0.3 matrixStats_0.62.0 ## [67] bitops_1.0-7 evaluate_0.15 ## [69] lattice_0.20-45 purrr_0.3.4 ## [71] labeling_0.4.2基因组校正s_1.32.0 ## [73] bit_4.0.4 tidyselect_1.1.2 ## [75] plyr_1.8.7 magrittr_2.0.3 ## [79] R6_2.5.1magick_2.7.3 # # [81] generics_0.1.2 DelayedArray_0.22.0 # # [83] DBI_1.1.2 withr_2.5.0 # # [85] pillar_1.7.0 survival_3.3-1 # # [87] KEGGREST_1.36.0 rcurl_1.98 - 1.6 # # [89] tibble_3.1.6 crayon_1.5.1 # # [91] HybridMTest_1.40.0 utf8_1.2.2 # # [93] BiocFileCache_2.4.0 rmarkdown_2.14 # # [95] progress_1.2.2 locfit_1.5 - 9.5 # # [97] grid_4.2.0 data.table_1.14.2 # # [99] blob_1.2.3 digest_0.6.29 # # [101] xtable_1.8-4 R.cache_0.15.0 # # [103] tidyr_1.2.0 R.utils_2.11.0 # # [105] munsell_0.5.0 bslib_0.3.1

13确认

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

参考文献

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

Bindreither, d . 2018。tbx20bam子集:来自"Tbx20"实验的Bam文件的子集．

王锐，R. Nambiar, D.郑，田斌。2017。PolyA_DB 3编目多重基因组深度测序鉴定的裂解和多聚腺苷酸化位点。核酸的研究46 (D1): D315-D319。https://doi.org/10.1093/nar/gkx1000．

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