ribosomeProfilingQC指南
2020年8月6日
包
ribosomeProfilingQC 1.0.6
内容
#介绍
核糖体的碳足迹,由乔纳森·斯曼和尼古拉斯Ingolia1、措施翻译直接量化的编码序列目前受80年代核糖体(ribosome-protected片段,rpf)2。在真核生物中,围绕28-nt rpf的大小,核糖体的p区通常在13从5 '端读取位置。在细菌中,艾伦等人能够更准确地确定p区3 '末端的读取3。
核糖体剖析的示意图表示。
在Bioconductor已经有几个包可用,包括,riboSeqR4,RiboProfiling5和ORFik6。这些包是强大的核糖体的碳足迹的分析数据。的ORFik包还可以用于转录起始点的寻找新的使用CageSeq数据。RiboWaltz7是另一个基于R和Bioconductor流行包。
帮助研究者快速评估核糖体分析数据的质量,我们开发了ribosomeProfilingQC包。的ribosomeProfilingQC包块可用于轻松地做出诊断检查映射质量和转移。此外,它可以为随后的微分分析预处理核糖体分析数据。我们试图让这个包尽可能友好的和唯一的输入是bam文件核糖体的足迹和RNAseq数据映射到基因组。
后请注意,所有的分析都是基于已知的开放阅读框(ORF)注释。包中提供的示例数据映射到斑马鱼UCSC danRer10组装;本届大会将相关的所有代码以淡黄色背景突出显示为清晰。
1快速启动
这是一个例子ribosomeProfilingQC与ribo-seq数据的一个子集。
第一次安装ribosomeProfilingQC和其他包需要运行示例。请注意,这里使用的示例数据集来自斑马鱼。运行分析与数据集从不同种类或不同的装配,请安装相应的Bsgenome和TxDb。例如,mm10分析鼠标数据一致,请安装BSgenome.Mmusculus.UCSC。mm10, TxDb.Mmusculus.UCSC.mm10.knownGene。你也可以通过函数生成TxDb对象makeTxDbFromGFF从本地文件人造石铺地面,或makeTxDbFromUCSC,makeTxDbFromBiomart,makeTxDbFromEnsembl,从网络资源GenomicFeatures包中。
库(BiocManager) BiocManager:安装(c (“ribosomeProfilingQC”、“AnnotationDbi”,“Rsamtools”、“BSgenome.Drerio.UCSC。danRer10”、“TxDb.Drerio.UCSC.danRer10。refGene”、“motifStack”))如果你在安装有困难ribosomeProfilingQC,请检查你的版本。的ribosomeProfilingQC包需要R > = 4.0。
R.versionx86_64-pc-linux-gnu # # # # _ # #平台拱x86_64 os linux-gnu # # # #系统x86_64, linux-gnu # # # #状态主要4 # #小0.2 # # # # 2020年06月22 # # # #天svn转速78730 R # # # #语言版本。字符串R版本4.0.2(2020-06-22)# #昵称再次起飞# #装载库库(ribosomeProfilingQC)图书馆(AnnotationDbi)图书馆(Rsamtools)1。1负载基因组
在本手册中,我们将使用鱼基因组。
图书馆(BSgenome.Drerio.UCSC.danRer10) # #设置基因组,Drerio是BSgenome.Drerio.UCSC的缩写。< - Drerio danRer10基因组如果你的装配是人类hg38请加载库,
< - Hsapiens库(BSgenome.Hsapiens.UCSC.hg38)基因组如果你的装配是鼠标mm10请加载鼠标库,
< - Mmusculus库(BSgenome.Mmusculus.UCSC.mm10)基因组1。2准备annotaiton cd
这个函数prepareCDS准备的信息用于下游分析从TxDb对象。
# # BSgenome.Drerio.UCSC对应。(TxDb.Drerio.UCSC.danRer10.refGene) txdb < - TxDb.Drerio.UCSC.danRer10 danRer10图书馆。refGene # #给它一个短名称cd < - prepareCDS (txdb)如果你的装配是人类hg38请尝试加载库,
库(TxDb.Hsapiens.UCSC.hg38.knownGene) txdb < - TxDb.Hsapiens.UCSC.hg38。knownGene # #给它一个短名称cd < - prepareCDS (txdb)如果你的装配是鼠标mm10请尝试加载库,
库(TxDb.Mmusculus.UCSC.mm10.knownGene) txdb < - TxDb.Mmusculus.UCSC.mm10。knownGene # #给它一个短名称cd < - prepareCDS (txdb)您还可以创建一个TxDb gtf文件对象GenomicFeatures: makeTxDbFromGFF函数。让GTF文件,您可以下载它运用或者在线文件信息通过AnnotationHub。这里我们使用一个准备TxDb对象进行测试。
# #创建一个小TxDb对象只包含chr1信息。库(GenomicFeatures) txdb < makeTxDbFromGFF(系统。文件(“extdata”、“Danio_rerio.GRCz10.91.chr1.gtf。广州”,包= " ribosomeProfilingQC”),生物= "鲐鱼类”,chrominfo = seqinfo (Drerio) (“chr1”), taxonomyId = 7955) cd < - prepareCDS (txdb)1。3输入
的输入ribosomeProfilingQC是bam文件。准备bam文件,不同riboSeqR包,问读取映射到转录组,ribosomeProfilingQC使用bam文件映射到整个基因组。正确映射读取,第一次尝试将适配器修剪序列映射到基因组大会通过bowtie2以下参数:——马5 mp 8 4读数7,7 -rfg 7, 7 - fr -nofw然后鱼片读取映射到rRNA, tRNA,核内小rna, snoRNA misc_RNA运用和Repeatmasker注释。之后,清洁读取映射到基因组大会通过tophat2以下参数:库类型fr-firststrand -transcriptome-index = Transcriptome_data /基因组。因为图书馆类型的ribo-seq重要strand-specific,请确保映射读取正确的库类型。
库(Rsamtools) # #输入的bamFile ribosomeProfilingQC包bamfilename < -系统。文件(“extdata”、“RPF.WT.1。bam”,包= " ribosomeProfilingQC”) # #为你自己的数据,请bamfilename设置为你的文件路径。# #例如,bam文件位于C: \ mydata \。# #你想要设置的bam bamfilename = " C: \ \ mydata \ \。bam“# #或由# # setwd你可以改变你的工作目录(C: \ \ mydata) # #,然后设置bamfilename = "。bam“yieldSize < - 10000000 bamfile < bamfile (bamfilename, yieldSize = yieldSize)1。4估计P网站
在上面的图,P的核糖体是13(如果使用核糖核酸酶I)的位置。然而,在不同的实验中,P网站可能是由于转向各种实验条件如酶和细胞类型的选择。的estimatePsite函数可以用来检查网站。的estimatePsite函数将搜索启动/停止密码子,发生在读取。的estimatePsite只能使用12、13或14最佳P网站从5 '末端候选人搜索时。
estimatePsite (bamfile、cd、基因组)# # 13 [1]它已经表明,某些酶,如MNase,估计3 '末端的P网站更好的工作3。的estimatePsite最好将使用15,16或17 P网站从3 '末端候选人搜索时。
estimatePsite (bamfile、cd、基因组,锚=“3”)# # -16年[1]1。5启动/停止的窗户
的readsEndPlot函数将图5 '末端或3 '端读取从cd启动/停止位置。当分配阅读框没有区别的读取,如果你设置最佳P网站13日或10或16(从5 '末端)。的readsEndPlot可以帮助用户确定真正的最佳Psite。在下面的示例中,起始密码子是丰富的位置9从5 '末端的读取和位置19 3 '端读取。这意味着有很多核糖体停靠的翻译起始位置和大部分的读取长度28元。
核糖体停靠TSS
readsEndPlot (bamfile、cd、toStartCodon = TRUE)
readsEndPlot (bamfile、cd、toStartCodon = TRUE, fiveEnd = FALSE)
如果你看到后分布,这意味着大量的基因活跃表达。
积极的表达
如果你看到一个警告或错误消息抱怨染色体序列的不一致,请核实你和正确使用TxDb对象基因组组装。如果这个警告消息是染色体的补丁你不感兴趣,你可以忽略的警告消息。
1。6阅读所有P站点坐标
的getPsiteCoordinates函数是用于读取所有P站点坐标。理想情况下,bestpsite应该13。要测试的数据质量,我们设置bestpsite = 13。
pc < - getPsiteCoordinates (bamfile bestpsite = 13)1。7片段大小分布
Ribosome-protected碎片应该27 29-nt长。检查片段大小分布,使用以下功能:
readsLen < - summaryReadsLength (pc)
1.7.1上过滤器的读取片段大小
过滤读取长度为下游分析,使用下面的脚本:
# # QC演示,我们将只使用读取长度28 - 29日nt.电脑。子< -个人电脑(pc qwidth美元% % c (28、29)]1。8感觉/反义链阴谋
大部分的阅读应该映射到有义链因为strand-specific ribo-seq库。
strandPlot (pc。子,CDS)
1。9基因组元素分布
核糖体的碳足迹,大部分的阅读应该映射到cd。的readsDistribution函数将显示P站点位置在不同基因组元素:cd、5 'utr 3 'utr,其他类型外显子,内含子、启动子下游,或基因间区域。下游高百分比显示,有很高比例的替代聚腺苷酸化网站从注释数据使用。很大一部分intronic地区显示intron-retaining成绩单的可能性。
电脑。子< - readsDistribution (pc。子,txdb las = 2)
1.10Metagene分析5 'utr / cd / 3 'utr情节
metagene图可以显示读取分布5 'utr, cd和3 'utr地区。
汇报工作。utr5 < - coverageDepth (rpf = bamfilename gtf = txdb地区=“utr5”)总经理汇报。cd < - coverageDepth (rpf = bamfilename gtf = txdb地区=“cd”)总经理汇报。utr3 < - coverageDepth (rpf = bamfilename gtf = txdb地区=“utr3”) metaPlot(汇报工作。utr5,汇报工作。cd、汇报工作。utr3、样品= 1)
1.11阅读框
函数assignReadingFrame用于设置P网站的阅读框位于已知的带注释的cd。的plotDistance2Codon函数可以用来绘制在转录起始和终止网站阅读框使用。函数plotFrameDensity可用于每一帧中所有的rpf崩溃。这些情节可以帮助你仔细检查,如果p区位置是否正确。如果它是正确的,应该分配给frame0大部分读取。
电脑。子< - assignReadingFrame (pc。子,CDS)plotDistance2Codon(pc.sub)
plotFrameDensity (pc.sub)
确定有多少生读取映射和P网站在0号帧中。
pc < - assignReadingFrame(电脑、cd) plotFrameDensity (pc)
函数plotTranscript可以用来查看阅读框分布记录。
plotTranscript (pc。子,c (“ENSDART00000161781”、“ENSDART00000166968”,“ENSDART00000040204”,“ENSDART00000124837”))
1.12ORFscore vs coverageRate
ORFscore2可以用来量化rpf的偏差分布对第一帧的cd。整个cd的覆盖率可以帮助研究人员检查rpf分布在整个cd。保险是由测量的比例在坐标系CDS头寸> = 1读取。如果报道是关于1,整个光盘是由活跃的核糖体。
总经理汇报< - frameCounts (pc。子,coverageRate = TRUE) ORFscore < - getORFscore (pc.sub) # #代码将情节ORFscores vs报道。# #试试通过删除“#”。#地块(总经理汇报[名称(ORFscore)], ORFscore, # xlab =“ORF1报道”,ylab =“羊痘疮分数”,#类型=“p”, pch = 16, cex =。5,xlim = c (0,1))2坏的情况下
在这里,我们显示核糖体的碳足迹数据质量差数据和用于下游分析不应该睡觉。
bamfilename < -系统。文件(“extdata”、“RPF.chr1.bad。bam”,包= " ribosomeProfilingQC ") yieldSize < - 10000000 bamfile < - bamfile (bamfilename, yieldSize = yieldSize)电脑< - getPsiteCoordinates (bamfile bestpsite = 13)电脑。子< -个人电脑(pc qwidth美元% % c(27 28 29)] # #在这个例子中,大部分的读取映射到反义链# #这可能表明有一些问题映射步骤strandPlot (pc。子,CDS)
# #在这个做法,大部分读取映射到# # inter-genic区域而不是cd,这可能表明,核糖体保护# #片段没有被恰当地隔离/选择电脑。子< - readsDistribution (pc。子,txdb las = 2)
# # P选择适当的地点也很重要。# #如果我们分配错了P的站点位置帧映射# #可能受到影响。pc < - getPsiteCoordinates (bamfile 12)电脑。子< -个人电脑(pc qwidth美元% % c(27 28 29)电脑。子< - assignReadingFrame (pc。子,CDS)plotDistance2Codon(pc.sub)
plotFrameDensity (pc.sub)
3准备下游分析
3所示。1rpf只
3.1.1数rpf
与RNAseq下游分析包括微分分析、比较,等等。函数frameCounts将生成一个为每个记录或基因数向量,可用于鉴别分析。countReads可用于计数ribo-seq多个文件。
库(ribosomeProfilingQC)库(AnnotationDbi) < -道路系统。文件(“extdata”,包=“ribosomeProfilingQC”) RPF < - dir(路径,“RPF。* ? \ \ [12]。bam美元”,full.names = TRUE) gtf < -文件。路径(路径,“Danio_rerio.GRCz10.91.chr1.gtf.gz”)碳纳米管< - countReads (rpf gtf = gtf水平=“基因”,bestpsite = 13, readsLen = c(28、29))头(碳纳米管rpf美元)# # RPF.KD1.1。bam RPF.KD1.2。bam RPF.WT.1。bam RPF.WT.2。bam # # ENSDARG00000000086 8 5 24日4 # # ENSDARG00000000606 23 10 0 0 # # ENSDARG00000001470 1 0 4 4 # # ENSDARG00000002285 5 5 0 1 # # ENSDARG00000002377 116 47 143 64 # # ENSDARG00000002634 4 3 0 0# #保存碳纳米管,请尝试以下代码通过删除“#”# write.csv (cbind(碳纳米管美元注释[rownames(碳纳米管rpf美元),),碳纳米管rpf美元),#“counts.csv”)让GTF文件,您可以下载它运用或者在线文件信息通过AnnotationHub。
BiocManager::安装(“AnnotationHub”)图书馆(AnnotationHub)啊= AnnotationHub() # #为人类hg38 hg38 < -查询(啊,c(“运用”、“GRCh38 gtf)) hg38 < - hg38[长度(hg38)] gtf < - mcols (hg38)美元sourceurl # #鼠标mm10 mm10 < -查询(啊,c(“运用”、“GRCm38 gtf)) mm10 < - mm10[长度(mm10)] gtf < - mcols mm10 sourceurl美元3.1.2只有rpf微分分析
库(磨边机)# #安装由BiocManager刨边机::安装(磨边机)医生< - c (“KD”、“KD”,“CTL”、“CTL”) # #样品组:KD:击倒;CTL:控制y < - DGEList(数量=碳管rpf美元,组= gp) y < - calcNormFactors (y)设计< - model.matrix (0 ~ + gp) colnames(设计)< -子(“全科医生”、“colnames(设计)y <——estimateDisp (y,设计)# #执行quasi-likelihood野生:符合< - glmQLFit (y,设计)qlf < - glmQLFTest(适合)topTags (qlf, n = 3) #组n = nrow (qlf)把所有的结果。# #系数:KD # # logFC logCPM F PValue罗斯福# # ENSDARG00000103054 -11.16762 8.682141 86767.21 6.128821 2.261535 e-16 e-13 # # ENSDARG00000074275 e15汽油1.931310 -10.96103 8.689056 63046.43 3.563268 e-13 # # ENSDARG00000043247 e-13 4.032526 -11.66550 8.621404 54103.55 3.346689 e15汽油# #进行似然比测试:符合< - glmFit (y,设计)轻轨交通< - glmLRT(适合)topTags(轻轨交通,n = 3) #集n = nrow(轻轨)把所有的结果。# # # #系数:KD logFC logCPM LR PValue罗斯福# # ENSDARG00000027355 -18.73631 9.551459 12085.481 0 0 # # ENSDARG00000053222 -14.92366 8.172169 6682.324 0 0 # # ENSDARG00000037748 -14.36672 7.796084 1603.511 0 03.1.3可变剪接、翻译起始和聚腺苷酸化
覆盖< - coverageDepth (rpf [grepl (“KD1 | WT”, rpf)], gtf = txdb水平=“基因”,地区=“扩展功能”)group1 < - c (“RPF.KD1.1”、“RPF.KD1.2”) group2 < - c (“RPF.WT。1”、“RPF.WT.2”) # #子集数据,仅供示例运行覆盖< -拉普(覆盖率,函数(.ele){# #不运行这个真实数据.ele美元覆盖< -拉普(。避署报道,美元[',i = seq.int (50)) .ele农庄组织< - .ele农庄美元[seq.int (50)] .ele}) se < - spliceEvent(覆盖、group1 group2)表(se美元类型)# # # # aSE # # 135plotSpliceEvent (se, se特性[1],美元保险,group1, group2)
3所示。2rpf和RNA-seq
3.2.1之上通过计数
3.2.1.1rpf和RNA-seq计数
的countReads功能可以用来估算ribo-seq和RNA-seq数据的多个文件。
路径< -系统。文件(“extdata”,包=“ribosomeProfilingQC”) RPF < - dir(路径,“RPF。* ? \ \ [12]。bam美元”,full.names = TRUE) rna < - dir(路径,“mRNA。* ? \ \ [12]。bam美元”,full.names = TRUE) gtf < -文件。路径(路径,“Danio_rerio.GRCz10.91.chr1.gtf.gz”)# #确保基因的顺序列出的bam文件rpf # #和RNAseq数据是相同的。碳纳米管< countReads (rpf rna, gtf水平= tx) # #拯救碳纳米管,请尝试以下代码通过删除“#”注释GeneID #美元得到rn # < -奈米碳管写。csv (cbind(#碳纳米碳管注释美元,美元rpf[匹配(rn, rownames(碳纳米管rpf)美元),),#碳纳米管美元mRNA(匹配(rn, rownames(碳纳米管mRNA)美元),]),#“counts.csv”)3.2.1.2转化效率(TE)
核糖体的绝对水平入住率强烈与编码和非编码RNA水平相关记录。介绍了转化效率8显示的相关性。TE的比例是规范化核糖体足迹丰富mRNA密度。一个常见的归一化法是使用片段每千碱基记录每百万映射读取(FPKM)。
fpkm < - getFPKM(碳纳米管)TE < - translationalEfficiency (fpkm)3.2.1.3微分分析TE
我们假设log2转化转化效率计算的比率rpf mrna有线性相关与真正的转化效率。然后,我们使用limma包测试微分转化效率。
库(limma)医生< - c (“KD”、“KD”,“CTL”、“CTL”) # #样品组:KD:击倒;CTL:控制TE。log2 < - log2 (TE TE + 1美元)#地块(TE。log2 [1], TE。log2 [3], # xlab = colnames (TE.log2) [1], ylab = colnames (TE.log2)[3], #主要=“转化效率”,pch = 16, cex = 5)设计< - model.matrix (0 ~ + gp) colnames(设计)< -子(“全科医生”、“colnames(设计)<——lmFit (TE。log2、设计)fit2 < - ebay(适合)topTable (fit2数量= 3)# #设置数量= nrow (fit2)把所有的结果# # CTL KD AveExpr F P。值# # ENSDART00000138070 e-06 # # ENSDART00000152424 23.095739 22.15413 22.624936 15378.78 6.130721 7.818822 1.00000 4.409411 13393.30 7.314603 e-06 # # ENSDART00000020327 13.376763 14.42394 13.900350 12287.78 8.165402 e-06 # # adj.P。Val 0.000514134 # # # # ENSDART00000138070 ENSDART00000152424 0.000514134 # # ENSDART00000020327 0.0005141343.2.2通过覆盖
3.2.2.1最大的n - m转化效率
如果我们把相关mrna或rpf转化效率水平
计算中所有重要的记录,我们将发现TE规范化。它显示了一个更高的价值低表达了成绩单和低价值在高度表达了成绩单。
plotTE (TE、样品= 2,xaxis =“mRNA log2 = TRUE, pch = 16, cex = 5)
# plotTE (TE、样品= 2,xaxis =“rpf log2 = TRUE, pch = 16, cex = 5)这个问题可以通过计算来解决最大值(TE max)在最高度ribosome-occupied 90元窗口内的功能8。请注意,TE马克斯的归一化法是不再FPKM。
汇报< - coverageDepth (rpf rna, txdb) TE90 < - translationalEfficiency(总经理汇报,窗口= 90,normByLibSize = TRUE) plotTE (TE90、样品= 2,xaxis =“mRNA log2 = TRUE, pch = 16, cex = 5)
plotTE (TE90示例= 2,xaxis =“爱国阵线”,log2 = TRUE, pch = 16, cex = 5)
上面的例子是TE90 cd地区。下面的代码展示了如何计算3 TE90 'utr地区。
汇报工作。utr3 < - coverageDepth (rpf, rna、txdb、地区=“utr3”) TE90。utr3 < - translationalEfficiency(汇报工作。utr3,窗口= 90)# #代码将情节TE90 3 'utr地区。# #试试通过删除“#”。# plotTE (TE90。utr3、样品= 2,xaxis = mRNA, log2 = TRUE, pch = 16, cex = 5) # plotTE (TE90。utr3、样品= 2,xaxis =“爱国阵线”,log2 = TRUE, pch = 16, cex = 5)3.2.2.2核糖体发布评分(RRS)
RRS的比率计算rpf RNA-seq读)(规范化的cd的rpf 3 'utr。因为很难定义cd地区非编码rna, RRS功能的非编码rna不能计算ribosomeReleaseScore。
RRS < - ribosomeReleaseScore (TE90 TE90。utr3 log2 = TRUE) # #代码将比较RSS 2样品。# #试试通过删除“#”。#地块(RRS [1]——[3], # xlab =“KD1 log2 RRS转换”,# ylab =“WT1 log2转换RRS”) # #代码将显示RSS TE90。# #试试通过删除“#”。#地块(RRS [1], log2 (TE90 TE美元[rownames (RRS), 1]), # xlab =“KD1 log2转换RSS”, # ylab =“KD1 log2转换TE”)3.2.2.3Metagene分析情节
情节metagene覆盖cd、5 'utr和3 'utr。你会注意到卢旺达爱国阵线的覆盖率数据更丰富的cd与相应的RNAseq数据
汇报工作。utr5 < - coverageDepth (rpf, rna、txdb、地区=“utr5”) # #样不同数量绘制metagene分析# #对不同样品# metaPlot(汇报工作。utr5总经理汇报,汇报工作。utr3示例= 2,xaxis = rpf) metaPlot(汇报工作。utr5总经理汇报,汇报工作。utr3示例= 2,xaxis =“信使rna”)
4片段长度组织相似性得分(棉花)1
牙线可以用来比较读取长度的分布等背景的基因。基因簇可以提取gtf /人造石铺地面文件下载运用。
# #文档:https://useast.ensembl.org/Help/Faq?id = 468 gtf < -进口(“Danio_rerio.GRCz10.91.gtf.gz”)gtf文件也可以通过下载AnnotationHub
BiocManager::安装(“AnnotationHub”)图书馆(AnnotationHub)啊= AnnotationHub() # #为人类hg38 hg38 < -查询(啊,c(“运用”、“GRCh38 gtf)) hg38 < - hg38[[长度(hg38)]] # #鼠标mm10 mm10 < -查询(啊,c(“运用”、“GRCm38 gtf)) mm10 < - mm10[[长度(mm10)]] # #因为TxDb.Mmusculus.UCSC.mm10基因id。knownGene和# # TxDb.Hsapiens.UCSC.hg38。entriz_id knownGene # #, gene_id entriz_id mm10或hg38需要改变。库(ChIPpeakAnno)库(org.Mm.eg.db) mm10 gene_id <美元——ChIPpeakAnno: xget (mm10 gene_id美元,org.Mm.egENSEMBL2EG)库(org.Hg.eg.db) hg38 gene_id <美元——ChIPpeakAnno: xget (hg38 gene_id美元,org.Mm.egENSEMBL2EG)gtf < - gtf [! is.na (gtf gene_id美元)]gtf < - gtf (gtf gene_id美元!= " ")# #蛋白质编码< - gtf gene_id美元(美元gtf transcript_biotype % % c (“IG_C_gene”、“IG_D_gene”,“IG_J_gene”、“IG_LV_gene”、“IG_M_gene”、“IG_V_gene”、“IG_Z_gene”、“nonsense_mediated_decay”、“nontranslating_CDS”、“non_stop_decay”、“protein_coding”、“TR_C_gene”、“TR_D_gene”、“TR_gene”、“TR_J_gene”、“TR_V_gene”)] # #线粒体基因水户< - gtf gene_id美元[grepl (“^ mt \ \”, gtf gene_name美元)| gtf transcript_biotype美元% % c (“Mt_tRNA”、“Mt_rRNA”)] # #长非编码lincRNA < - gtf gene_id美元(美元gtf transcript_biotype % % c (“3 prime_overlapping_ncrna”、“lincRNA”、“ncrna_host”、“non_coding”)] # #短非编码sncRNA < - gtf gene_id美元(美元gtf transcript_biotype % % c (“microrna的”、“miRNA_pseudogene”、“misc_RNA”、“misc_RNA_pseudogene”、“Mt_rRNA”、“Mt_tRNA”、“Mt_tRNA_pseudogene”、“ncRNA”、“pre_miRNA”、“RNase_MRP_RNA”、“RNase_P_RNA”、“核糖体rna”、“rRNA_pseudogene”、“scRNA_pseudogene”、“snlRNA”、“snoRNA”、“snRNA_pseudogene”、“SRP_RNA”、“tmRNA”、“tRNA”、“tRNA_pseudogene”、“核糖酶”,“scaRNA”、“sRNA”)] # #假基因假基因< - gtf gene_id美元(美元gtf transcript_biotype % % c (“disrupted_domain”、“IG_C_pseudogene”、“IG_J_pseudogene”、“IG_pseudogene”、“IG_V_pseudogene”、“processed_pseudogene”、“假基因”、“transcribed_processed_pseudogene”、“transcribed_unprocessed_pseudogene”、“translated_processed_pseudogene”、“translated_unprocessed_pseudogene”、“TR_J_pseudogene”、“TR_V_pseudogene”、“unitary_pseudogene”、“unprocessed_pseudogene”)] danrer10。注释< -列表(蛋白质=独特(蛋白质),水户=独特(水),lincRNA =独特(lincRNA) sncRNA =独特(sncRNA),假基因=独特(假基因))这里我们加载pre-saved chr1注释的示例代码。
danrer10。注释< - readRDS(系统。文件(“extdata”、“danrer10.annotations。rds”,包= " ribosomeProfilingQC "))库(ribosomeProfilingQC)库(GenomicFeatures) # # cd注释txdb < - makeTxDbFromGFF(系统做好准备。文件(“extdata”、“Danio_rerio.GRCz10.91.chr1.gtf。广州”,包= " ribosomeProfilingQC”),生物= "鲐鱼类”,chrominfo = seqinfo (Drerio) (“chr1”), taxonomyId = 7955) cd < - prepareCDS (txdb)库(Rsamtools) # #输入的bamFile ribosomeProfilingQC包bamfilename < -系统。文件(“extdata”、“RPF.WT.1。bam”,包= " ribosomeProfilingQC”) # #为你自己的数据,请bamfilename设置为你的文件路径。yieldSize < - 10000000 bamfile < bamfile (bamfilename yieldSize = yieldSize)电脑< - getPsiteCoordinates (bamfile bestpsite = 13) readsLengths < - 20:34 fl < -牙线(pc, ref = danrer10。注释蛋白质美元、cd = cd, readLengths = readsLengths水平=“基因”,画= FALSE) fl.max < - t (fl [c (1, which.max (fl cooks.distance美元)),as.character (readsLengths)]) matplot (fl。max =“l”型x = readsLengths xlab =“片段长度”,ylab =“阅读的分数”,上校= c(“灰色”、“红色”)、lwd = 2, lty = 1)传说(“topright”,传说= c (“ref”、“选择基因”),坳= c(“灰色”、“红色”)、lwd = 2, lty = 1, cex = 5)
引用
1。Ingolia: T。et al。核糖体分析揭示了普遍的翻译之外的带注释的蛋白质编码基因。细胞的报道8日,1365 - 1379 (2014)。
2。Bazzini, A。et al。识别小orf使用核糖体足迹和脊椎动物进化的保护。在EMBO杂志33岁的981 - 993 (2014)。
3所示。穆罕默德,F。,Green, R. & Buskirk, A. R. A systematically-revised ribosome profiling method for bacteria reveals pauses at single-codon resolution.Elife8日,e42591 (2019)。
4所示。钟,b . Y。et al。使用duplex-specific在核糖核酸酶分析和用户友好的软件包ribo-seq数据分析。核糖核酸21日,1731 - 1745 (2015)。
5。Popa,。et al。RiboProfiling: bioconductor包标准ribo-seq管道处理。F1000Research5,(2016)。
6。人类uORFome Tjeldnes, h·阿特拉斯和监管组织。卑尔根大学(2018)。
7所示。Lauria F。et al。RiboWaltz:优化核糖体核糖体p区定位的分析数据。PLoS计算生物学14日,e1006169 (2018)。
8。Ingolia: T。,Ghaemmaghami, S., Newman, J. R. & Weissman, J. S. Genome-wide analysis in vivo of translation with nucleotide resolution using ribosome profiling.科学324年,218 - 223 (2009)。