内容

#介绍

核糖体足迹，由Jonathan Weissman和Nicholas Ingolia开发1，通过直接量化目前由80S核糖体(核糖体保护片段，RPFs)结合的编码序列来测量翻译。2在真核生物中，rpf的大小约为28-nt，其中核糖体的p位点通常位于reads 5 '端起的13号位置。在细菌中，Allen等人能够更准确地从读取的3 '端识别p位点3.．

核糖体轮廓的示意图表示。

Bioconductor已经提供了几个软件包，包括，riboSeqR，4RiboProfiling5而且ORFik6．这些软件包在分析核糖体足迹数据方面功能强大。的ORFik包也可以使用CageSeq数据来寻找新的转录起始位点。RiboWaltz7是另一个流行的包，它基于R和Bioconductor。

为了帮助研究人员快速评估核糖体分析数据的质量，我们开发了ribosomeProfilingQC包。的ribosomeProfilingQC软件包可以用来方便地做出诊断图，以检查映射质量和帧移。此外，它还可以对核糖体分析数据进行预处理，用于后续的差异分析。我们已经尝试使这个包尽可能的用户友好，唯一需要的输入是一个bam文件的核糖体足迹和RNAseq数据映射到基因组。

请注意，以下所有分析都是基于已知的开放阅读框(ORF)注释。包中提供的样本数据映射到斑马鱼UCSC danRer10组件;与此程序集相关的所有代码都将被删除淡黄色背景突出显示为清晰。

1快速启动

下面是一个使用ribosomeProfilingQC用核糖测序数据的子集。

第一次安装ribosomeProfilingQC以及运行示例所需的其他包。请注意，这里使用的示例数据集来自斑马鱼。要使用来自不同物种或不同组件的数据集运行分析，请安装相应的Bsgenome和TxDb。例如，要分析与mm10对齐的鼠标数据，请安装BSgenome.Mmusculus.UCSC。TxDb.Mmusculus.UCSC.mm10.knownGene。还可以通过函数生成TxDb对象makeTxDbFromGFF从本地GFF文件，或者makeTxDbFromUCSC，makeTxDbFromBiomart,makeTxDbFromEnsembl，从网上资源中GenomicFeatures包中。

library(BiocManager) BiocManager::install(c("ribosomeProfilingQC"， "AnnotationDbi"， "Rsamtools"， "BSgenome.Drerio.UCSC. "danRer10”、“TxDb.Drerio.UCSC.danRer10。refGene”、“motifStack”))

如果安装有问题ribosomeProfilingQC，请先检查你的R版本。的ribosomeProfilingQC包装要求R >= 4.0。

R.version
平台x86_64-pc-linux-gnu ## arch x86_64 ## os linux-gnu ##系统x86_64, linux-gnu ##状态## #主4 ##次0.4 ##年2021 ##月02 ##日15 ## svn rev 80002 ##语言r# #版本。R版本4.0.4 (2021-02-15
##加载库(ribosomeProfilingQC)库(AnnotationDbi)库(Rsamtools)

1.1负载基因组

在本手册中，我们将使用鱼类基因组。

库(BSgenome.Drerio.UCSC. danrer10) ##集合基因组，Drerio是BSgenome.Drerio.UCSC的缩写。danRer10基因组<- Drerio

如果您的程序集是人类hg38，请加载人类库，

库(BSgenome.Hsapiens.UCSC.hg38)基因组<- Hsapiens. hg38

如果您的程序集是鼠标mm10，请加载鼠标库，

库(BSgenome.Mmusculus.UCSC.mm10)基因组<- Mmusculus. mm10

1.2准备注释光盘

这个函数prepareCDS用于为下游分析从TxDb对象。

##对应BSgenome.Drerio.UCSC。txdb <- txdb . drerio . ucsc .danRer10. refgenerefGene ##给它一个简短的名字CDS <- prepareCDS(txdb)

如果您的程序集是Human hg38，请尝试加载库，

txdb <- txdb . hsapiens . ucsc .hg38. knowngeneknownGene ##给它一个简短的名字cd <- prepareCDS(txdb)

如果您的程序集是鼠标mm10，请尝试加载库，

txdb <- txdb . mmusculus . ucsc .mm10. knowngeneknownGene ##给它一个简短的名字cd <- prepareCDS(txdb)

还可以通过以下方法从gtf文件创建TxDb对象GenomicFeatures: makeTxDbFromGFF函数。要获得GTF文件，您可以从运用或获得在线文件信息通过AnnotationHub．这里我们使用一个准备好的TxDb对象进行测试。

创建一个只包含chr1信息的小TxDb对象。library(GenomicFeatures) txdb <- makeTxDbFromGFF(system. library);file("extdata"， "Danio_rerio.GRCz10.91.chr1.gtf.gz"， package="ribosomeProfilingQC")， organism =" Danio rerio"， chrominfo = seqinfo(Drerio)["chr1"]， taxonomyId = 7955) cd <- prepareCDS(txdb)

1.3.输入

的输入ribosomeProfilingQC是bam文件。准备bam文件，不同于riboSeqR转录组，ribosomeProfilingQC使用bam文件映射到整个基因组。为了获得正确的映射reads，首先尝试通过bowtie2使用以下参数将适配器修剪的序列映射到基因组组装:-local -ma 5 -mp 8,4 -rdg 7,7 -rfg 7,7 -fr -nofw，然后从Ensembl和Repeatmasker注释中过滤映射到rRNA、tRNA、snRNA、snoRNA和misc_RNA的reads。然后，通过tophat2将干净的reads映射到基因组组装，参数如下:-library-type fr-firststrand - transcriptome_index =Transcriptome_data/genome。因为ribo-seq的库类型通常是特定于字符串的，请确保用正确的库类型映射读取。

从ribosomeProfilingQC包bamfilename <- system. library(Rsamtools) ##输入bamFile文件(“extdata”、“RPF.WT.1。， package="ribosomeProfilingQC") ##对于您自己的数据，请设置bamfilename为您的文件路径。例如，你的bam文件位于C:\mydata\a。你要设置bamfilename = "C:\\mydata\\a. "或者你可以通过## setwd("C:\\mydata") ##来改变你的工作目录，然后设置bamfilename = "a.b bam" yieldSize <- 10000000 bamfile <- bamfile (bamfilename, yieldSize = yieldSize)

1.4估计P点

如上图所示，核糖体的P位点位于第13位(如果使用RNase I)，但在不同的实验中，P位点可能会因酶的选择、细胞类型等不同的实验条件而发生移位。的estimatePsite函数可以用来检查P点。的estimatePsite函数将搜索读取过程中出现的开始/停止密码子。的estimatePsite从5 '端开始搜索时，只会使用12、13或14作为最佳P站点候选。

estimatePsite(bamfile, CDS，基因组)
## [1]

已经证明，对于某些酶，如MNase，从3 '端估计P位点效果更好3.．的estimatePsite当从3 '端搜索时，将使用15、16或17作为最佳P站点候选。

estimatePsite(bamfile, CDS, genome, anchor = "3end")
## [1] -16

1.5绘制启动/停止窗口

的readsEndPlot函数将绘制从cd的开始/停止位置移位的5 '末端或3 '末端读取。如果你将最佳P点设置为13或10或16(从5 '端开始)，那么在为大多数读取分配阅读帧时没有区别。的readsEndPlot可以帮助用户确定真正的最佳Psite。在下面的例子中，起始密码子从读取的5 '末端的-9位置富集，从读取的3 '末端的19位置富集。这意味着有很多核糖体停靠在转译起始位置，大多数reads长度为28 nt。

核糖体对接TSS

readsEndPlot(bamfile, CDS, toStartCodon=TRUE)

readsEndPlot(bamfile, cd, toStartCodon=TRUE, fiveEnd=FALSE)

如果你看到下面的分布，说明有很多基因在活跃表达。

积极的表达

如果您看到一条关于染色体序列不一致的警告或错误消息，请验证您正在使用具有正确基因组组装的TxDb对象。如果此警告消息针对的是您不感兴趣的补丁染色体，则可以忽略警告消息。

1.6读取所有P点坐标

的getPsiteCoordinates函数用于读取所有P点坐标。理想情况下，最佳位置应该是13。为了测试数据质量，我们设置bestpsite = 13。

pc <- getPsiteCoordinates(bamfile, bestpsite = 13)

1.7碎片大小分布

核糖体保护的片段理想情况下应该是27 - 29 nt长。查看分片大小分布，使用以下函数:

readsLen <- summaryReadsLength(pc)

1.7.1上根据片段大小过滤读取

使用下面的脚本根据读取的长度进行下游分析:

##对于这个QC演示，我们将只使用28-29 nt. pc的读取长度。Sub <- pc[pc$qwidth %in% c(28,29)]

1.8正义/反义链图

因为ribo-seq库是特定于链的，所以大多数读取应该映射到感知链。

strandPlot (pc。子,CDS)

1.9基因组元素分布

对于核糖体足迹，大多数读取应该映射到CDS区域。的readsDistribution函数将显示P位点在不同基因组元件中的位置:CDS, 5'UTR, 3'UTR，其他类型外显子，内含子，启动子，下游或基因间区。较高的下游百分比表明注释数据中有较高比例的替代多聚腺苷酸位点使用。内含子区域的高比例表明可能存在保留内含子的转录本。

电脑。sub <- readsDistribution(pc。子，txdb, las=2)

1.105'UTR /CDS/ 3'UTR的元基因分析图

元基因图可以显示5'UTR、CDS和3'UTR区域的reads分布。

汇报工作。utr5 <- coverageDepth(RPFs = bamfilename, gtf = txdb, region="utr5") cvgs。CDS <- coverageDepth(RPFs = bamfilename, gtf = txdb, region=" CDS ") cvgs。utr3 <- coverageDepth(RPFs = bamfilename, gtf = txdb, region="utr3")utr5,汇报工作。cd、汇报工作。utr3、样品= 1)

1.11阅读框

函数assignReadingFrame用于为位于已知带注释的cd中的P个站点设置读取帧。的plotDistance2Codon函数可用于绘制转录起始或终止位点的阅读框使用情况。函数plotFrameDensity可以用来折叠每一帧中的所有rpf。这些图可以帮助你再次检查p点位置是否正确。如果是正确的，大多数读取应该分配到frame0。

电脑。sub <- assignReadingFrame(pc。子,CDS)plotDistance2Codon(pc.sub)

plotFrameDensity (pc.sub)

来确定有多少原始读取映射到帧0中的P个位置。

- assignReadingFrame(pc, CDS) plotFrameDensity(pc)

函数plotTranscript可用于查看给定文本的阅读帧分布。

plotTranscript (pc。sub, c("ENSDART00000161781"， "ENSDART00000166968"， "ENSDART00000040204"， "ENSDART00000124837"))

1.12ORFscore vs coverageRate

ORFscore2可以用来量化rpf对给定CDS的第一帧的偏向分布。整个CDS的覆盖率可以帮助研究人员检查整个CDS的RPFs分布。覆盖率是通过测量帧内CDS位置>= 1读取的比例来确定的。如果覆盖率约为1，则整个CDS都被活性核糖体覆盖。

cvg <- frameCounts(pc。sub, coverageRate=TRUE) ORFscore <- getORFscore(pc.sub) ##下面的代码将绘制ORFscore与覆盖率的关系。试着去掉“#”。#plot(cvg[名称(ORFscore)]， ORFscore， # xlab="覆盖率ORF1"， ylab="ORF评分"，# type="p"， pch=16, cex=。5, xlim=c(0,1))

2坏的情况下

在这里，我们展示了核糖体足迹数据是低质量的数据，不应该用于下游分析。

Bamfilename <- system。文件(“extdata”、“RPF.chr1.bad。bam"， package="ribosomeProfilingQC") yieldSize <- 10000000 bamfile <- bamfile (bamfilename, yieldSize = yieldSize) pc <- getPsiteCoordinates(bamfile, bestpsite = 13) pc。sub <- pc[pc$qwidth %in% c(27,28,29)] ##在本例中，大多数读取都映射到反义strand ##，这可能表明在映射步骤strandPlot(pc. c.)中存在一些问题。子,CDS)

在本例中，大多数读取都映射到基因间区域##，而不是CDS，这可能表明核糖体保护的##片段没有被正确分离/选择。sub <- readsDistribution(pc。子，txdb, las=2)

##选择合适的P站点也很关键。如果我们分配错误的P点位置，帧映射可能会受到影响。pc <- getPsiteCoordinates(bamfile, 12)Sub <- pc[pc$qwidth %in% c(27,28,29)] pc。sub <- assignReadingFrame(pc。子,CDS)plotDistance2Codon(pc.sub)

plotFrameDensity (pc.sub)

3.为下游分析做准备

3.1rpf只

3.1.1rpf计数

下游分析包括差分分析、与RNAseq的比较等。函数frameCounts将为每个转录本或基因生成计数向量，可用于差异分析。countReads可用于计数多个文件的核糖-seq。

library(ribosomeProfilingQC) library(AnnotationDbi) path <- system. library(ribosomeProfilingQC)file("extdata"， package="ribosomeProfilingQC") RPFs <- dir(path， "RPF.*?\\.[12]. txt ")bam$"， full.names=TRUE)path(path， " danio_rreio . grcz10.91 .chr1.gtf.gz") cnts <- countReads(RPFs, gtf=gtf, level="gene"， bestpsite=13, readsLen=c(28,29)) head(cnts$RPFs)
# # RPF.KD1.1。bam RPF.KD1.2。bam RPF.WT.1。bam RPF.WT.2。bam ## ENSDARG00000000086 85 24 4 ## ENSDARG00000000606 23 10 00 ## ENSDARG00000001470 10 4 4 ## ENSDARG00000002285 5 01 ## ENSDARG00000002377 116 47 143 64 ## ENSDARG00000002634 43 00
##保存碳纳米管，请尝试删除'#' # write.csv(cbind(cnts$annotation[rownames(cnts$RPFs)，]， cnts$RPFs)， # "counts.csv")

要获得GTF文件，您可以从运用或获得在线文件信息通过AnnotationHub．

BiocManager::install("AnnotationHub") library(AnnotationHub) ah = AnnotationHub() ## for human hg38 hg38 <- query(ah, c("Ensembl"， "GRCh38"， "gtf")) hg38 <- hg38[length(hg38)] gtf <- mcols(hg38)$sourceurl ## for mouse mm10 mm10 <- query(ah, c("Ensembl"， "GRCm38"， "gtf")) mm10 <- mm10[length(mm10)] gtf <- mcols(mm10)$sourceurl

3.1.2仅对rpf进行微分分析

库(edgeR) ##安装edgeR by BiocManager::install("edgeR") gp <- c("KD"， "KD"， "CTL"， "CTL") ##样本组:KD:击倒;CTL:Control y <- DGEList(counts = cnts$RPFs, group = gp) y <- calcNormFactors(y) design <- model.matrix(~0+gp) colnames(design) <- sub("gp"， ""， colnames(design)) y <- estimateDisp(y, design) ##执行准似然f检验:fit <- glmQLFit(y, design) qlf <- glmQLFTest(fit) topTags(qlf, n=3) # set n=nrow(qlf)来拉出所有结果。
##系数:KD ## logFC logCPM F PValue FDR ## ENSDARG00000103054 -11.16762 8.682141 86767.21 6.128821e-16 2.261535e-13 ## ENSDARG00000074275 -10.96103 8.689056 63046.43 1.931310e-15 3.563268e-13 ## ENSDARG00000043247 -11.66550 8.621404 54103.55 3.346689e-15 4.032526e-13
要执行似然比测试:fit <- glmFit(y, design) lrt <- glmLRT(fit) topTags(lrt, n=3) # set n=nrow(lrt)以拉出所有结果。
##系数:KD ## logFC logCPM LR PValue FDR ## ENSDARG00000027355 -18.73631 9.551459 12085.481 00 ## ENSDARG00000053222 -14.92366 8.172169 6682.324 00 ## ENSDARG00000037748 -14.36672 7.796084 1603.511 00

3.1.3选择性剪接，翻译起始和多聚腺苷酸化

coverage <- coverageDepth(RPFs[grepl("KD1|WT"， RPFs)]， gtf=txdb, level="gene"， region="feature with extension") group1 <- c("RPF.KD1.1"， "RPF.KD1.2") group2 <- c("RPF.WT. wd1.2 ")1"， "RPF.WT.2") ##子集数据，对于样本只运行覆盖<- lapply(覆盖，函数(.ele){##不运行这个真实的数据。ele$coverage <- lapply(. ele)。ele$coverage， ' ['， i=seq.int(50))) .ele$granges <- .ele$granges[seq.int(50)] .ele}) se <- spliceEvent(coverage, group1, group2)表(se$type)
## ## aSE ## 129
plotSpliceEvent(se, se$feature[1]， coverage, group1, group2)

3.2RPFs和RNA-seq

3.2.1之上通过计数

3.2.1.1计算RPFs和RNA-seq

的countReads函数可以用来计数多个文件的核糖-seq和RNA-seq数据。

路径<- system. Path。file("extdata"， package="ribosomeProfilingQC") RPFs <- dir(path， "RPF.*?\\.[12]. txt ")bam美元”,full.names = TRUE) rna < - dir(路径,“mRNA。* ? \ \[12]。bam$"， full.names=TRUE)路径(路径,“Danio_rerio.GRCz10.91.chr1.gtf.gz”)
确保rpff# #和RNAseq数据的bam文件中列出的基因顺序是相同的。cnts <- countReads(RPFs, rna, gtf, level="tx") ##要保存碳纳米管，请尝试删除'#' # rn <- cnts$annotation$GeneID # write.csv(cbind(cnts$annotation， # cnts$RPFs[match(rn, rownames(cnts$RPFs)))，]， # cnts$mRNA[match(rn, rownames(cnts$mRNA))，])， # " numbers .csv")

3.2.1.2翻译效率(TE)

核糖体占用的绝对水平与编码和非编码转录本的RNA水平密切相关。介绍了平移效率8为了显示相关性。TE是归一化核糖体足迹丰度与mRNA密度的比值。一种常用的归一化方法是使用每千碱基的片段数每百万的转录本映射读数(FPKM)。

fpkm <- getFPKM(cnts) TE <- translationalEfficiency(fpkm)

3.2.1.3TE的微分分析

我们假设，通过rpf与mrna的比值计算的log2转换翻译效率与实际翻译效率呈线性相关。然后我们使用limma包来测试微分平移效率。

库(limma) gp <- c("KD"， "KD"， "CTL"， "CTL") ##样本组:KD:敲除;CTL:控制TE。log2 <- log2(TE$TE + 1) #plot(TE。log2[， 1]， TE。log2[， 3]， # xlab=colnames(TE.log2)[1]， ylab=colnames(TE.log2)[3]， # main="翻译效率"，pch=16, cex=.5) design <- model.matrix(~0+gp) colnames(设计)<- sub("gp"， ""， colnames(设计))fit <- lmFit(TE. log2)[1]， ylab=colnames(TE.log2)[3])log2, design) fit2 <- eBayes(fit) topTable(fit2, number=3) ## set number=nrow(fit2)提取所有结果
## CTL KD AveExpr F P.Value ## ens达特00000138070 23.095739 22.15413 22.624936 15378.78 6.130721e-06 ## ens达特00000152424 7.818822 1.00000 4.409411 13393.30 7.314603e-06 ## ens达特00000020327 13.376763 14.42394 13.900350 12287.78 8.165402e-06 # adj.P.Val ## ens达特00000138070 0.000514134 ## ens达特00000152424 0.000514134 ## ens达特000020327 0.000514134

3.2.2通过覆盖

3.2.2.1最大N-mer平动效率

如果我们绘制mrna或rpf水平与翻译效率的相关性
通过抄本内的所有计数计算，我们会发现TE没有很好地标准化。它在低表达的转录本中表现出较高的值，而在高表达的转录本中表现出较低的值。

plot (TE, sample=2, xaxis="mRNA"， log2=TRUE, pch=16, cex=.5)

#plot (TE, sample=2, xaxis="RPFs"， log2=TRUE, pch=16, cex=.5)

这个问题可以通过计算一个特征中核糖体占用最多的90 nt窗口的最大值(TE max)来解决8．请注意，TE max的归一化方法不再是FPKM了。

cvgs <- coverageDepth(RPFs, rna, txdb) TE90 <- translationalEfficiency(cvgs, window = 90, normByLibSize=TRUE) plot (TE90, sample=2, xaxis="mRNA"， log2=TRUE, pch=16, cex=.5)

plot (TE90, sample=2, xaxis="RPFs"， log2=TRUE, pch=16, cex=.5)

上面的例子是针对CDS地区的TE90。下面的代码演示如何计算3'UTR区域的TE90。

汇报工作。utr3 <- coverageDepth(RPFs, rna, txdb, region="utr3") TE90。utr3 <- translationalEfficiency(cvgs. utr3)下面的代码将为3'UTR区域绘制TE90。试着去掉“#”。# plotTE (TE90。utr3，样本=2,xaxis="mRNA"， log2=TRUE, pch=16, cex=.5) #plot (TE90。utr3, sample=2, xaxis="RPFs"， log2=TRUE, pch=16, cex=.5)

3.2.2.2核糖体释放评分(RRS)

RRS计算为CDS中的rpf(由RNA-seq reads归一化)与3'UTR中的rpf之比。由于非编码rna的CDS区域难以定义，因此无法通过Function计算非编码rna的RRSribosomeReleaseScore．

RRS <- ribosomereleascore (TE90, TE90。下面的代码将比较2个样本的RSS。试着去掉“#”。#plot(RRS[， 1]， RRS[， 3]， # xlab="log2转换后的KD1 RRS "， # ylab="log2转换后的WT1 RRS ") ##下面的代码将显示沿着TE90的RSS。试着去掉“#”。#plot(RRS[， 1]， log2(TE90$TE[rownames(RRS)， 1])， # xlab="log2转换后的KD1的RSS "， # ylab="log2转换后的KD1的TE ")

3.2.2.3元基因分析图

绘制CDS、5'UTR和3'UTR的超基因覆盖率。您将注意到，与相应的RNAseq数据相比，CDS中RPF数据的覆盖范围要丰富得多

汇报工作。utr5 <- coverageDepth(RPFs, rna, txdb, region="utr5") ##设置不同数量的样本来绘制宏基因分析##针对不同的样本#metaPlot(cvgs。ut5 cvgs cvgs。utr3, sample=2, xaxis = "RPFs")ut5 cvgs cvgs。utr3, sample=2, xaxis = "mRNA")

4片段长度组织相似度评分(FLOSS)1

FLOSS可用于比较reads长度的分布与背景(如基因簇)。基因簇可以从下载的gtf/gff文件中提取运用．

##文档:https://useast.ensembl.org/Help/Faq?id=468 gtf <- import("Danio_rerio.GRCz10.91.gtf.gz")

gtf文件也可以通过AnnotationHub

BiocManager::install("AnnotationHub") library(AnnotationHub) ah = AnnotationHub() ## for human hg38 hg38 <- query(ah, c("Ensembl"， "GRCh38"， "gtf")) hg38 <- hg38[[length(hg38)]] ## for mouse mm10 mm10 <- query(ah, c("Ensembl"， "GRCm38"， "gtf")) mm10 <- mm10[[length(mm10)]]] ##因为基因id在TxDb.Mmusculus.UCSC.mm10中。## TxDb.Hsapiens.UCSC.hg38。knownGene ##为entriz_id, mm10或hg38中的gene_id需要修改为entriz_id。library(ChIPpeakAnno) library(org.Mm.eg.db) mm10$gene_id <- ChIPpeakAnno::xget(mm10$gene_id, org.Mm.egENSEMBL2EG) library(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 %in% c("IG_C_gene"， "IG_J_gene"， "IG_LV_gene"， "IG_M_gene"， "IG_V_gene"， "IG_Z_gene"， "nonsense_mediated_decay"， "nontranslating_CDS"， " nonstop_decay "， "protein_coding"， "TR_C_gene"， "TR_D_gene"， "TR_J_gene"， "TR_V_gene")] ##线粒体基因mito <- gtf$gene_id[grepl("^mt\\-"， gtf$gene_name) | gtf$transcript_biotype %in% c("Mt_tRNA"，##长非编码lincRNA <- gtf$gene_id[gtf$transcript_biotype %in% c("3prime_overlapping_ncrna"， "lincRNA"， "ncrna_host"， "non_coding")] ##短非编码sncRNA <- gtf$gene_id[gtf$transcript_biotype %in% c("miRNA"， "miRNA_pseudogene"， "misc_RNA"， "misc_RNA_pseudogene"， "Mt_rRNA"， "Mt_tRNA"， "Mt_tRNA_pseudogene"， "ncRNA"， "pre_miRNA"， "RNase_MRP_RNA"， "RNase_MRP_RNA"， " rna_pseudogene "， "snlRNA"， "snoRNA"， "snRNA_pseudogene"， "SRP_RNA"， "tmRNA"， "tRNA"，" trna_假基因"，"核酶"，"scaRNA"， "sRNA")] ##假基因<- gtf$gene_id[gtf$transcript_biotype %in% c("disrupted_domain"， "IG_C_pseudogene"， "IG_J_pseudogene"， "IG_V_pseudogene"， "经处理的假基因"，"经处理的假基因"，"经处理的假基因"，"经处理的假基因"，"TR_J_pseudogene"， "unitary_pseudogene"， "unprocessed_pseudogene")] danrer10。注释<- list(protein=unique(protein)， mito=unique(mito)， lincRNA=unique(lincRNA)， sncRNA=unique(sncRNA)， pseudogene=unique(pseudogene))

在这里，我们为示例代码加载预先保存的chr1注释。

danrer10。注释<- readRDS(system。文件(“extdata”、“danrer10.annotations。rds"， package = "ribosomeProfilingQC"))
library(ribosomeProfilingQC) library(GenomicFeatures) ## prepare CDS注释txdb <- makeTxDbFromGFF(system. txt)file("extdata"， "Danio_rerio.GRCz10.91.chr1.gtf.gz"， package="ribosomeProfilingQC")， organism =" Danio rerio"， chrominfo = seqinfo(Drerio)["chr1"]， taxonomyId = 7955) CDS <- prepareCDS(txdb) library(Rsamtools) ##输入bamFile从ribosomeProfilingQC包bamfilename <-系统。文件(“extdata”、“RPF.WT.1。， package="ribosomeProfilingQC") ##对于您自己的数据，请设置bamfilename为您的文件路径。yieldSize <- 10000000 bamfile <- bamfile (bamfilename, yieldSize = yieldSize) pc <- getPsiteCoordinates(bamfile, bestpsite = 13) readsLengths <- 20:34 fl <- FLOSS(pc, ref = danrer10。注释$protein, CDS = CDS, readlength =readsLengths, level="gene"， draw = FALSE) fl.max <- t(fl[c(1, where .max(fl$cooks.distance))， as.character(readsLengths)]) matplot(fl. character)max, type =" l"， x=readsLengths, xlab="Fragment Length"， ylab="Fraction of Reads"， col = c("gray"， "red")， lwd = 2, lty = 1) legend("topright"， legend = c("ref"， "selected gene")， col = c("gray"， "red")， lwd = 2, lty = 1, cex=.5)

参考文献

1.北卡罗来纳州，英格利亚。et al。核糖体分析揭示了注释蛋白编码基因之外的普遍翻译。细胞的报道8日,1365 - 1379(2014)。

2.巴齐尼，A。et al。利用核糖体足迹和进化保护识别脊椎动物的小孔。EMBO期刊33岁的981 - 993(2014)。

3.Mohammad, F.， Green, R. & Buskirk, A. R.一种系统修订的细菌核糖体分析方法揭示了单密码子分辨率的暂停。Elife8日,e42591(2019)。

4.钟炳义et al。在核糖体分析中使用双工特异性核酸酶和用于核糖测序数据分析的用户友好的软件包。核糖核酸21日,1731 - 1745(2015)。

5.Popa,。et al。RiboProfiling:用于标准ribo-seq管道处理的生物导体包。F1000Research5，(2016)。

6.人类uORFome及其跨组织调节的图谱。(卑尔根大学，2018)。

7.Lauria F。et al。RiboWaltz:核糖体分析数据中核糖体p位点定位的优化。PLoS计算生物学14，e1006169(2018)。

8.Ingolia, N. T.， Ghaemmaghami, S.， Newman, J. R. & Weissman, J. S.利用核糖体谱分析核苷酸解析的体内翻译全基因组分析。科学324年,218 - 223(2009)。