其他
Ribo-seq 质控软件：ribosomeProfilingQC
杨花落尽子规啼
1、介绍
今天来介绍一款对 Ribo-seq 进行质量控制和绘图的 R 包：ribosomeProfilingQC。
通过对核糖体保护的片段（ribosome-protected fragments, RPFs）进行高通量测序，然后进行定量，即可获得 mRNA 的翻译水平，在真核生物中，一般核糖体保护的片段大小在 28nt 左右：
mRNA 翻译过程大概分为：1、翻译起始。2、翻译延伸。3、翻译终止。3 个过程，核糖体的 p 位 通常位于 5 末端 的 13 位：
2、使用
使用前先安装需要的 R 包，R 版本需要 4.0 以上：
BiocManager::install(c("ribosomeProfilingQC",
"AnnotationDbi", "Rsamtools",
"BSgenome.Drerio.UCSC.danRer10",
"TxDb.Drerio.UCSC.danRer10.refGene",
"motifStack"))
library(ribosomeProfilingQC)
library(AnnotationDbi)
library(Rsamtools)
2.1：加载基因组
# If your assembly is Human hg38 please load the human library,
library(BSgenome.Hsapiens.UCSC.hg38)
genome <- Hsapiens
# If your assembly is Mouse mm10 please load the mouse library,
library(BSgenome.Mmusculus.UCSC.mm10)
genome <- Mmusculus
# 测试数据需要的基因组
library(BSgenome.Drerio.UCSC.danRer10)
## set genome, Drerio is a shortname for BSgenome.Drerio.UCSC.danRer10
genome <- Drerio
2.2：准备 CDS 注释文件
# If your assembly is Human hg38 please try to load the library,
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene ## give it a short name
CDS <- prepareCDS(txdb)
# If your assembly is Mouse mm10 please try to load the library,
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene ## give it a short name
CDS <- prepareCDS(txdb)
# 测试数据需要的
## which is corresponding to BSgenome.Drerio.UCSC.danRer10
library(TxDb.Drerio.UCSC.danRer10.refGene)
txdb <- TxDb.Drerio.UCSC.danRer10.refGene ## give it a short name
CDS <- prepareCDS(txdb)
你可以通过 GenomicFeatures::makeTxDbFromGFF 函数从自己本地 gtf 文件创建 TxDb 对象：
library(GenomicFeatures)
txdb <- makeTxDbFromGFF(gtffile, format="gtf")
2.3：输入文件
需要比对到 基因组 的 bam 文件，这里使用包自带的测试文件：
读取 bam 文件：
# 加载测试数据
bamfilename <- system.file("extdata", "RPF.WT.1.bam",
package="ribosomeProfilingQC")
# 设置每次读取数里
yieldSize <- 10000000
# 读取bam文件
bamfile <- BamFile(bamfilename, yieldSize = yieldSize)
2.4：评估 P 位点
核糖体的 P 位点位于 13 位(如果使用 RNase I)，但在不同的实验中，由于酶的选择和细胞类型等实验条件的不同，P 位点可能会发生移动。estimatePsite 函数可以从 5'末端寻找 reads 上的 start/stop codons 来推测 P 位点，12、13、14 作为最好的 P 候选位点：
estimatePsite(bamfile, CDS, genome)
[1] 13
已经证明，对于某些酶，如 MNase，从 3 端估计 P 位点效果更好。当从 3 端搜索时，estimatePsite 将使用 15、16 或 17 作为最佳 P 位点候选：
estimatePsite(bamfile, CDS, genome, anchor = "3end")
[1] -16
2.5：绘制起始密码子/终止密码子的窗口
eadsEndPlot 可以帮助确定真正最好的 P site。在下面的例子中，开始密码子在读码的 5’端位置 9，reads 的 3’端位置 19。这意味着有大量的核糖体停靠在翻译起始位置，大部分读取长度为 28nt：
readsEndPlot(bamfile, CDS, toStartCodon=TRUE)
readsEndPlot(bamfile, CDS, toStartCodon=TRUE, fiveEnd=FALSE)
以下的分布模式意味着大量的基因处于活跃的表达中：
2.6：获取 P 位点坐标
getPsiteCoordinates 函数用于读取所有 P 位点坐标。理想情况下，最好的位点应该是 13。为了检验数据质量，我们设置 bestpsite = 13：
pc <- getPsiteCoordinates(bamfile, bestpsite = 13)
2.7：片段长度分布
理想的核糖体保护片段长度为 27 ~ 29 nt。要检查 reads 长度分布，使用以下函数：
readsLen <- summaryReadsLength(pc)
可以筛选特定的长度做下游分析：
## for this QC demo, we will only use reads length of 28-29 nt.
pc.sub <- pc[pc$qwidth %in% c(28, 29)]
2.8：正义链和反义链绘图
大部分 Ribo-seq 的建库是 链特异性 的，所以大部分会比对到正义链上：
strandPlot(pc.sub, CDS)
2.9：基因位置分布
对于核糖体足迹，大多数读取应该比对到 CDS 区域。readsDistribution 函数将显示不同基因组元件中的 P 位点:CDS、5' UTR、3' UTR、其他类型外显子、内含子、启动子、下游或基因间区域：
pc.sub <- readsDistribution(pc.sub, txdb, las=2)
2.10：5’UTR/CDS/3’UTR 的元基因图
5' UTR、CDS 和 3' UTR 区域的 reads 分布可以用元基因图表示：
cvgs.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")
metaPlot(cvgs.utr5, cvgs.CDS, cvgs.utr3, sample=1)
2.11：阅读框
函数 assignReadingFrame 用于设置位于已知注释 cds 内的 P 位点的读框。plotDistance2Codon 函数可以用来绘制翻译起始或终止位点的阅读框使用情况。函数 plotframeddensity 可以用来绘制不同阅读框的比例。这些图可以帮助你反复检查 P 位点的位置是否正确。如果它是正确的，大多数读取应该分配到 frame0 ：
pc.sub <- assignReadingFrame(pc.sub, CDS)
plotDistance2Codon(pc.sub)
plotFrameDensity(pc.sub)
用原始数据来查看：
pc <- assignReadingFrame(pc, CDS)
plotFrameDensity(pc)
函数 plotTranscript 可用于查看特定转录本的阅读框分布：
plotTranscript(pc.sub, c("ENSDART00000161781", "ENSDART00000166968",
"ENSDART00000040204", "ENSDART00000124837"))
2.12：ORFscore vs coverageRate
ORFscore 可以用来量化 rpf 对给定 CDS 的 frame 0 的偏倚分布。整个 CDS 的覆盖率（Coverage rate）可以帮助研究人员检查整个 CDS 的 RPFs 分布。覆盖率是通过测量>= 1 in-frame 的 CDS 的比例来确定的。如果覆盖率约为 1，则整个 CDS 被活性核糖体覆盖：
cvg <- frameCounts(pc.sub, coverageRate=TRUE)
ORFscore <- getORFscore(pc.sub)
## Following code will plot the ORFscores vs coverage.
plot(cvg[names(ORFscore)], ORFscore,
xlab="coverage ORF1", ylab="ORF score",
type="p", pch=16, cex=.5, xlim=c(0, 1))
3、数据不好的示例
4、下游分析
4.1：计数（Count for RPFs）
library(ribosomeProfilingQC)
library(AnnotationDbi)
path <- system.file("extdata", package="ribosomeProfilingQC")
# 加载bam文件
RPFs <- dir(path, "RPF.*?\\.[12].bam$", full.names=TRUE)
# 加载注释文件
gtf <- file.path(path, "Danio_rerio.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 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(cnts$annotation[rownames(cnts$RPFs), ], cnts$RPFs),
# "counts.csv")
要获得 GTF 文件，可以从 ensemble 下载它，或者通过 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
gtf
[1] "ftp://ftp.ensembl.org/pub/release-103/gtf/homo_sapiens/Homo_sapiens.GRCh38.103.gtf.gz"
## for mouse mm10
mm10 <- query(ah, c("Ensembl", "GRCm38", "gtf"))
mm10 <- mm10[length(mm10)]
gtf <- mcols(mm10)$sourceurl
gtf
[1] "ftp://ftp.ensembl.org/pub/release-102/gtf/mus_musculus/Mus_musculus.GRCm38.102.gtf.gz"
4.2：RPFs 的差异分析
# install edgeR by BiocManager::install("edgeR")
library(edgeR)
# sample groups: KD:knockdown; CTL:Control
gp <- c("KD", "KD", "CTL", "CTL")
y <- DGEList(counts = cnts$RPFs, group = gp)
y <- calcNormFactors(y)
design <- model.matrix(~0+gp)
colnames(design) <- sub("gp", "", colnames(design))
design
CTL KD
1 0 1
2 0 1
3 1 0
4 1 0
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$gp
[1] "contr.treatment"
y <- estimateDisp(y, design)
## To perform quasi-likelihood F-tests:
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit)
# set n=nrow(qlf) to pull all results.
topTags(qlf, n=3)
Coefficient: KD
logFC logCPM F PValue FDR
ENSDARG00000103054 -11.16762 8.682141 86767.19 6.128825e-16 2.261536e-13
ENSDARG00000074275 -10.96103 8.689056 63046.41 1.931312e-15 3.563270e-13
ENSDARG00000043247 -11.66550 8.621404 54103.55 3.346689e-15 4.032524e-13
## To perform likelihood ratio tests:
fit <- glmFit(y, design)
lrt <- glmLRT(fit)
# set n=nrow(lrt) to pull all results.
topTags(lrt, n=3)
Coefficient: KD
logFC logCPM LR PValue FDR
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 0
4.3：RPFs 和 RNA-seq 一起计数
# 数据路径
path <- system.file("extdata", package="ribosomeProfilingQC")
# ribo data
RPFs <- dir(path, "RPF.*?\\.[12].bam$", full.names=TRUE)
# rna data
RNAs <- dir(path, "mRNA.*?\\.[12].bam$", full.names=TRUE)
# gtf data
gtf <- file.path(path, "Danio_rerio.GRCz10.91.chr1.gtf.gz")
## make sure that the order of the genes listed in the bam files for RPFs and RNAseq data is the same.
cnts <- countReads(RPFs, RNAs, gtf, level="tx")
# 查看结果
head(cnts$RPFs,3)
RPF.KD1.1.bam RPF.KD1.2.bam RPF.WT.1.bam RPF.WT.2.bam
ENSDART00000000087 8 5 24 4
ENSDART00000002886 41 25 10 8
ENSDART00000003317 14 3 2 0
# mRNA
head(cnts$mRNA,3)
mRNA.KD1.1.bam mRNA.KD1.2.bam mRNA.WT.1.bam mRNA.WT.2.bam
ENSDART00000164359 2 3 5 2
ENSDART00000157701 0 0 2 0
ENSDART00000158290 0 0 0 2
# 注释
head(cnts$annotation,3)
GeneID Chr Start
1 ENSDART00000164359 chr1;chr1;chr1;chr1;chr1;chr1 6408;6892;9558;10081;11550;11751
2 ENSDART00000157701 chr1;chr1;chr1 6642;6892;9558
3 ENSDART00000158290 chr1;chr1 6642;6892
End Strand Length
1 6760;6955;9694;10191;11625;12027 -;-;-;-;-;- 1018
2 6760;6955;9919 -;-;- 545
3 6760;7335 -;- 563
# 保存数据
rn <- cnts$annotation$GeneID
write.csv(cbind(cnts$annotation,
cnts$RPFs[match(rn, rownames(cnts$RPFs)), ],
cnts$mRNA[match(rn, rownames(cnts$mRNA)), ]),
"counts.csv")
4.4：计算翻译效率（Translational Efficiency （TE））
# 计算fpkm
fpkm <- getFPKM(cnts)
# 计算TE
TE <- translationalEfficiency(fpkm)
head(TE$TE,3)
RPF.KD1.1.bam RPF.KD1.2.bam RPF.WT.1.bam RPF.WT.2.bam
ENSDART00000000087 1.254362e-01 0.1056103 0.1209225 0.06922225
ENSDART00000002886 1.016359e+04 2734.3877431 1887.9600718 679.32052417
ENSDART00000003317 1.008050e+03 631.8404789 186.5107979 1.00000000
4.5：计算差异翻译效率
使用 limma 包：
library(limma)
# sample groups: KD:knockdown; CTL:Control
gp <- c("KD", "KD", "CTL", "CTL")
# 取log2
TE.log2 <- log2(TE$TE + 1)
# 分组矩阵
design <- model.matrix(~0+gp)
colnames(design) <- sub("gp", "", colnames(design))
fit <- lmFit(TE.log2, design)
fit2 <- eBayes(fit)
# set number=nrow(fit2) to pull all results
topTable(fit2, number=3)
CTL KD AveExpr F P.Value adj.P.Val
ENSDART00000138070 23.095739 22.15413 22.624936 15378.78 6.130721e-06 0.000514134
ENSDART00000152424 7.818822 1.00000 4.409411 13393.30 7.314603e-06 0.000514134
ENSDART00000020327 13.376763 14.42394 13.900350 12287.78 8.165402e-06 0.000514134
使用 edgeR 包：
library(edgeR)
gp <- c("KD", "KD", "CTL", "CTL")
design <- model.matrix(~0+gp)
colnames(design) <- sub("gp", "", colnames(design))
y <- lapply(cnts[c("RPFs", "mRNA")], function(.cnt){
.y <- DGEList(counts = .cnt,
samples = data.frame(gp=gp))
.y <- calcNormFactors(.y)
.y <- estimateDisp(.y, design)
.y
})
## To perform likelihood ratio tests:
fit <- lapply(y, glmFit, design = design)
lrt <- lapply(fit, glmLRT)
# difference between CTL and KD for RPFs
topTags(lrt[["RPFs"]], n=3)
Coefficient: KD
logFC logCPM LR PValue FDR
ENSDART00000160650 -18.74358 9.514299 11432.193 0 0
ENSDART00000110824 -15.31215 8.284378 7025.735 0 0
ENSDART00000075230 -14.92803 8.151413 4590.713 0 0
# difference between CTL and KD for mRNA
topTags(lrt[["mRNA"]], n=3)
Coefficient: KD
logFC logCPM LR PValue FDR
ENSDART00000157701 -15.79362 8.460275 19653.85 0 0
ENSDART00000158290 -15.79362 8.460789 19653.85 0 0
ENSDART00000161139 -15.79362 8.307693 19653.85 0 0
results <- lapply(lrt, topTags, p.value = 0.05, n = nrow(lrt[[1]]))
RPFs_unique <- setdiff(rownames(results[["RPFs"]]),
rownames(results[["mRNA"]]))
mRNA_zero <- rownames(cnts[["mRNA"]])[rowSums(cnts[["mRNA"]]) == 0]
RPFs_unique <- RPFs_unique[!RPFs_unique %in% mRNA_zero]
## top 3 events for translation level
head(results[["RPFs"]][RPFs_unique, ], n=3)
Coefficient: KD
logFC logCPM LR PValue FDR
ENSDART00000139003 -12.92114 7.834459 1747.696 0 0
ENSDART00000127426 -12.69363 8.589542 14080.615 0 0
ENSDART00000147683 -11.89208 8.254962 6028.894 0 0
剩下的其它功能大家自行去探索吧，R 包在 Bioconductor 上。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦, 数据和代码已上传至 QQ群，欢迎加入下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
欢迎小伙伴留言评论！
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀CIRCexplorer3: 对 circRNA 进行相对定量
◀circRNA-seq：CIRCexplorer2 使用指南（二）
◀circRNA-seq：CIRCexplorer2 使用指南（一）
◀...