Molecular Cell 文章结果复现（终）

引言

❝
给前面文章复现的结果收个尾，差不多文章里大部分的生信图都复现的差不多了。学有所得，总之是花费了很多时间和精力的。之前的可以参看下面推文链接：
Molecular Cell 文章部分结果复现
Molecular Cell 文章部分结果复现(续)
计算 Ribosome 距离二肽/三肽基序距离的密度
安装

install.packages(“devtools”)

devtools::install_github(“junjunlab/RiboProfiler”)

or

remotes::install_github(“junjunlab/RiboProfiler”)

library(RiboProfiler)
富集 motif 的序列可视化

occupancy_file = c(“peptide_motif/WT_tripeptide_occupancy.txt”,
“peptide_motif/eIF5Ad_tripeptide_occupancy.txt”)
sampe_name = c(“WT”,”eIF5Ad”)

triAmino_motif_plot(occupancy_file = occupancy_file,
sample_name = c(“WT”,”eIF5Ad”),
top_motif = 29)

metagene

整体的分布情况：

使用 metagene_plot 函数可视化，先看看距离 start codon 的：
library(ggplot2)

non_hs <- norm_df %>%
dplyr::filter(sample %in% c(“WT”,”eIF5Ad”))

hs <- norm_df %>%
dplyr::filter(sample %in% c(“WTHS”,”eIF5AdHS”))

metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,100),
mode = “nt”)

高盐条件处理的：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,100),
mode = “nt”)

可以把 mode 改为 codon 单位：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,500),
mode = “codon”)

设置 collapse = T 合并展示：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-25,500),
mode = “codon”,
collapse = T)

展示离 stop codon 的分布：
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-500,25),
mode = “codon”,
collapse = T,
type = “sp”)

文章展示如下，右侧为特定位置放大展示：

我们绘图看看：
pnon_hs <- metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-100,50),
mode = “nt”,
collapse = T,
type = “sp”)

p_hs <- metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-100,50),
mode = “nt”,
collapse = T,
type = “sp”)

library(patchwork)

pnon_hs / p_hs

然后把特定范围放大展示：

close up view

library(ggforce)

pnon_hs_close_up <-
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = non_hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-100,50),
mode = “nt”,
collapse = T,
type = “sp”) +
facet_zoom(xlim = c(-5,50),ylim = c(0,1),zoom.size = 0.5)

p_hs_close_up <-
metagene_plot(longest_trans_file = “longest_info_extend.txt”,
normed_file = hs,
min_counts = 64,
min_cds_length = 600,
relative_distance = c(-100,50),
mode = “nt”,
collapse = T,
type = “sp”) +
facet_zoom(xlim = c(-5,50),ylim = c(0,1),zoom.size = 0.5)

pnon_hs_close_up / p_hs_close_up

以上就是文献主要的图，文章的大多数图就基本复现完成了。
track plot

❝
前面提到 codon_track_plot 函数是以 codon 为单位在 CDS 区域绘制核糖体占据密度的峰图的。有时候你可能需要查看是否在 UTR 区域是否存在 翻译 情况，我们可以直接传给 track_plot 函数即可。

首先使用 get_track_df 函数获取对应基因的绘图数据：
track_df <- get_track_df(longest_trans_file = “longest_info_extend.txt”,
select_gene = c(“THO1″,”PCL6″,”CHC1”),
normed_file = norm_df)

check

head(track_df,3)

sample gene_name trans_id transpos density type

1 eIF5Ad CHC1 YGL206C_mRNA 2369 3.359513 ribo

2 eIF5Ad PCL6 YER059W_mRNA 498 7.784284 ribo

3 eIF5Ad CHC1 YGL206C_mRNA 4410 3.359513 ribo

然后就可以进行可视化了：
track_plot(signal_data = track_df,
gene_anno = “longest_info_extend.txt”)

下面是文章选择基因的峰图：

结尾

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