引言

❝
ribo-seq 数据的质量对于后续的分析非常重要。一般在进行后续深入分析的时候都需要提前评估数据的质量是否符合常规 ribo-seq 的特征。分享一篇浙大毛原辉老师的文章。

目前下面几种 ribo-seq 的测序技术，不同技术应用所研究的问题有所不一样：

关于检测翻译效率的变化水平，一些由于核糖体在 CDS 上发生暂停而导致核糖体密度的升高，所带来的偏差有待进一步评估：

小开放阅读框的检测（smORFs）也非常依赖于 riboseq 数据的质量：

核糖体保护片段的大小收到 RNA 酶消化时间的影响，消化时间不够会导致过长的长度，此外在有些特殊情况下（氨基酸缺乏）也会产生较短的片段：

被保护片段应该大部分处于 CDS 区域内：

一些 5UTR 区域的可能是翻译起始，3UTR 的则可能是遇到终止密码子导致通读或者重新进行翻译起始，不过这种概率比较低：

被保护片段的对框（in-frame）情况也是一个重要的指标：：

对框的计算收到 5‘end reads 位置计算的影响：

不同建库的方法也和数据质量关系非常相关：

Dual-ligation 方法提高了 5’end 的准确性：

毛老师开发了一个新的方法大大提高了对框的比例：

安装

install.packages(“devtools”)

devtools::install_github(“junjunlab/RiboProfiler”)

or

remotes::install_github(“junjunlab/RiboProfiler”)

library(RiboProfiler)
5‘end 和 3’end 的区别

❝
前面我们复现的 molecular cell 文章使用的就是 3‘end assignment 去分析的，我们拿的是 5’end assignment 去做所有分析的，我们可以看到只有 28nt 长度的片段对框还不错：

然后我们增加几行代码去做这样一个 3‘end assignment 这个事情：

check assign type

if assignType == “end5”
# relative distance
rel2st = exact_pos – start_codon_pos
rel2sp = exact_pos – (stop_codon_pos + 1)
else
# relative distance
rel2st = exact_pos + read_length – 1 – start_codon_pos
rel2sp = exact_pos + read_length – 1 – (stop_codon_pos + 1)
end

assign frame

frame_st = abs(rel2st)%3
frame_sp = abs(rel2sp)%3
指定一下 assignType 参数：
bam_file = list.files(“../3.map-data”,pattern = “.genome.bam$”,full.names = TRUE)
bam_file

1] “../3.map-data/eIF5Ad.sorted.genome.bam” “../3.map-data/eIF5AdHS.sorted.genome.bam”

[3] “../3.map-data/WT.sorted.genome.bam” “../3.map-data/WTHS.sorted.genome.bam”

sample_name <- sapply(strsplit(list.files(“../3.map-data/”,pattern = “.sorted.genome.bam$”),
split = “.sorted.genome.bam”), “[“,1)
sample_name

[1] “eIF5Ad” “eIF5AdHS” “WT” “WTHS”

perform QC analysis

pre_qc_data(longest_trans_file = “longest_info_extend.txt”,
bam_file = bam_file[c(1,3)],
out_file = paste(sample_name[c(1,3)],”.qc.txt”,sep = “”),
mapping_type = “genome”,
assignType = “end3”,
seq_type = “singleEnd”)

load qc data

qc_df <- load_qc_data()
然后检测一下对框情况：

length with frame

qc_plot(qc_data = qc_df,type = “length_frame”,
facet_wrap_list = list(ncol = 2))

❝
可以看到对框情况改善了很多。同时把对框信息可以展示在 aggregation plot 里：
offset <- PsiteOffsetCheck(qc_data = qc_df,
read_length = c(25,30),
relative_distance = c(-20,50))
offset$plot
offset$summary_offset

summary_offset <- offset$summary_offset
summary_offset$Offsets <- c(“15, 15, 15, 15, 15, 15”)

offset_df <- adjust_offset(offset_df = summary_offset,
qc_df = qc_df)

norm_df <- get_nomalized_counts(longest_trans_file = “longest_info_extend.txt”,
qc_df = offset_df)

这里发现其实用峰最高的地方来评估不同片段长度的 offset 感觉存在一些问题，那个最高的峰前面其实还有比较高的峰，所以是不是最前面才算是 offset 呢？
aggragation plot
画个正常的看看：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = norm_df,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,500),
mode = “nt”)

添加 frame 对框的信息：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = norm_df,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,500),
frame = T,
mode = “nt”)

转换为 codon 单位：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = norm_df,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,500),
frame = T,
mode = “codon”)

结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。