Molecular Cell 文章部分结果复现(续)

引言

❝
codon occupancy(密码子占有率) 指的是核糖体在不同密码子上停留的多少。反映了核糖体在特定密码子上的相对停留时间。高 occupancy 表示核糖体在该密码子上停留时间较长，低 occupancy 则表示停留时间较短。与翻译的延伸密切相关，比如核糖体在一些特定位置密码子或者二肽/三肽基序会发生暂停。
❝
根据已有的研究文献，核糖体在某些特定密码子或密码子序列上会出现显著的暂停。这些暂停位点通常与特定的生物学功能或翻译调控机制相关。以下是一些常见的暂停位点类型：
稀有密码子： 例如，在大肠杆菌中，AGG、AGA（精氨酸）和 CUA（亮氨酸）等稀有密码子。 在哺乳动物细胞中，CUA（亮氨酸）和 CGA（精氨酸）等密码子。
连续碱基相同的密码子： 多个脯氨酸密码子（如 CCC-CCC-CCC） 多个甘氨酸密码子（如 GGG-GGG-GGG） 这些序列可能导致核糖体滑移或暂停
带电氨基酸序列： 连续的碱性氨基酸（如赖氨酸 AAA、AAG） 连续的酸性氨基酸（如谷氨酸 GAA、GAG）
特定的三肽或四肽序列： PPP（脯氨酸-脯氨酸-脯氨酸） PPG（脯氨酸-脯氨酸-甘氨酸） PP(X)P，其中 X 可以是任何氨基酸
依赖 eIF5A 的位点： 主要是含有脯氨酸和甘氨酸的序列 例如，PPG、PPP、PGG 等
与蛋白质折叠相关的位点： 通常位于结构域边界 可能包含疏水氨基酸簇
依赖 RACK1 的位点： 某些含有精氨酸的序列
与 mRNA 二级结构相关的位点： 通常位于强烈的 mRNA 二级结构附近
终止密码子附近： 特别是在某些读穿现象（readthrough）中观察到 起始密码子附近： 通常与翻译起始的调控相关 内部重新起始位点： 某些病毒和细胞内 mRNA 的特征
依赖 Shine-Dalgarno 序列的位点（原核生物）： 通常与翻译暂停和重新起始相关 D-氨基酸引起的暂停： 在某些特殊情况下，如细菌细胞壁合成中
与翻译后修饰相关的位点： 例如，某些需要共翻译修饰的位点
应激条件下特异的暂停位点： 如热休克或氧化应激条件下
物种特异的暂停位点： 不同物种可能有其特有的暂停模式 这些暂停位点的确切序列和影响程度可能因物种、细胞类型和生理条件而异。此外，许多研究发现，暂停并不总是由单个密码子决定，而是由更长的序列上下文共同影响。因此，在解释 codon occupancy 数据时，需要考虑更广泛的序列环境和生物学背景。
前面复现的文章也有相关的结果（二肽/三肽基序附近的位置）：

计算三肽的暂停得分（筛选规则：过滤出现频率小于 30/100 次的基序）：

方法部分：

❝
我们自己去做这样一个计算，看看是否能得到文章类似的结果。不求一模一样，毕竟这种科研性的文章能复现百分之七八十就已经很不错了。
安装

install.packages(“devtools”)

devtools::install_github(“junjunlab/RiboProfiler”)

or

remotes::install_github(“junjunlab/RiboProfiler”)

library(RiboProfiler)
codon occupancy

首先从 sam/bam 提取数据：
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”,
seq_type = “singleEnd”)

load qc data

qc_df <- load_qc_data()

check

head(qc_df,3)

length framest relst framesp relsp feature trans_pos trans_id counts sample group

1 27 0 2862 0 -270 3 2913 YLR249W_mRNA 30 eIF5Ad NA

2 28 0 1443 0 -1569 3 1494 YDR270W_mRNA 1 eIF5Ad NA

3 24 0 123 0 -252 3 174 YKL192C_mRNA 1 eIF5Ad NA

然后看一下 psite 的 offset，我们看 28nt 的 3nt 周期性较好，就选用这个片段来分析：
offset <- PsiteOffsetCheck(qc_data = qc_df,read_length = c(28,28))
offset$plot

offset 信息：
offset$summary_offset

sample readLengths Offsets bamFiles bamLegends

1 eIF5Ad 28 -12 eIF5Ad eIF5Ad

2 WT 28 -12 WT WT

使用 adjust_offset 函数矫正 offset：
offset_df <- adjust_offset(offset_df = offset$summary_offset,
qc_df = qc_df)
对每个位置的 counts 进行标准化：
norm_df <- get_nomalized_counts(longest_trans_file = “longest_info_extend.txt”,
qc_df = offset_df)

head(norm_df,3)

length framest relst framesp relsp feature trans_pos trans_id counts sample group norm_exp

1 27 0 2870 0 -270 3 2921 YLR249W_mRNA 30 eIF5Ad NA 0.4015918

2 28 0 1455 0 -1569 3 1506 YDR270W_mRNA 1 eIF5Ad NA 6.9953596

3 28 0 651 0 -309 3 702 YKL019W_mRNA 1 eIF5Ad NA 2.1564626

提取转录本的 CDS 序列：
FetchSequence(gene_file = “longest_info.txt”,
genome_file = “./Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa”,
output_file = “sac_cds.fa”,
type = “cds”,
coding_type = “NT”)
然后计算 codon occupancy：
codon_occupancy(qc_file = norm_df,
longest_trans_file = “longest_info_extend.txt”,
min_counts = 32,
cds_fasta_file = “sac_cds.fa”)

list.files(“codon_occupancy/”)
[1] “eIF5Ad_codon_occupancy.txt” “eIF5Ad_codon_pos_exp.txt”
[3] “WT_codon_occupancy.txt” “WT_codon_pos_exp.txt”
绘图，其中 SP 代表 stop codon：
list.files(path = “codon_occupancy/”,pattern = “occupancy.txt”,full.names = T)

library(ggh4x)

codon_occupancy_plot(codon_occupancy_file = c(“codon_occupancy/eIF5Ad_codon_occupancy.txt”,
“codon_occupancy/WT_codon_occupancy.txt”),
sample_name = c(“eIF5Ad”,”WT”),
compare_var = c(“eIF5Ad”,”WT”),
codon_type = “amino”)

你可以使用 rm_stopcodon 去除：
codon_occupancy_plot(codon_occupancy_file = c(“codon_occupancy/eIF5Ad_codon_occupancy.txt”,
“codon_occupancy/WT_codon_occupancy.txt”),
sample_name = c(“eIF5Ad”,”WT”),
compare_var = c(“eIF5Ad”,”WT”),
codon_type = “amino”,
rm_stopcodon = T)

用 codon 展示：
codon_occupancy_plot(codon_occupancy_file = c(“codon_occupancy/eIF5Ad_codon_occupancy.txt”,
“codon_occupancy/WT_codon_occupancy.txt”),
sample_name = c(“eIF5Ad”,”WT”),
compare_var = c(“eIF5Ad”,”WT”),
codon_type = “codon”)

分面展示：
codon_occupancy_plot(codon_occupancy_file = c(“codon_occupancy/eIF5Ad_codon_occupancy.txt”,
“codon_occupancy/WT_codon_occupancy.txt”),
sample_name = c(“eIF5Ad”,”WT”),
compare_var = c(“eIF5Ad”,”WT”),
codon_type = “codon”,
facet = T)

tri amino acids pause score

❝
前面提到 三肽 上的核糖体密度也和核糖体的延伸相关。但是文章里并没有说如何选取三肽的。因为基因的氨基酸序列如果是 3 的倍数，则刚好可以三个三个取，但是如果不是 3 的倍数，那么如何取呢，还是就是从头往后三个三个取，直到最后一个后者两个氨基酸直接不要。暂时我是按这种方法算的。
我们看看氨基酸序列长度是否是 3 的倍数的数量：
amino = pyfastx.Fastx(“sac_cds_aa.fa”)
ct3 = 0
non_ct3 = 0
for name,seq in amino:
if len(seq)%3 == 0:
ct3 += 1
else:
non_ct3 += 1

print(ct3)
print(non_ct3)

❝
可以看到有 2/3 的基因的氨基酸序列并不是 3 的倍数，所以这个 tri-peptide 的选取还是有点说法的。文章不说，我们也不清楚。
计算 tripeptide pause score：
peptide_motif_score(amino_file = “sac_cds_aa.fa”,
codon_exp_file = c(“codon_occupancy/WT_codon_pos_exp.txt”,
“codon_occupancy/eIF5Ad_codon_pos_exp.txt”),
sample_name = c(“WT”,”eIF5Ad”),
occurrence_threshold = 50)
绘图：
triAmino_scater_plot(occupancy_file = c(“peptide_motif/WT_tripeptide_occupancy.txt”,
“peptide_motif/eIF5Ad_tripeptide_occupancy.txt”),
sampe_name = c(“WT”,”eIF5Ad”),
mark_motif = c(“DDP”,”DPP”,”PPP”,”PDK”,”PDG”,
“GPP”,”PDP”,”DNP”,”DGP”,”SPP”))

和文章还是比较类似的。
codon track plot

文章里的 track plot 也可以画一下：
library(patchwork)

codon_track_plot(codon_exp_file = c(“codon_occupancy/eIF5Ad_codon_pos_exp.txt”,
“codon_occupancy/WT_codon_pos_exp.txt”),
longest_trans_file = “longest_info_extend.txt”,
sample_name = c(“eIF5Ad”,”WT”),
select_gene = c(“THO1″,”PCL6″,”CHC1”))

文章里标注了 THO1 这个基因暂停位点的三肽基序为 DGP，我们拿自己的数据验证一下：

read.delim(“codon_occupancy/eIF5Ad_codon_pos_exp.txt”,header = F) %>%
dplyr::filter(V1 == “YER063W_mRNA”) %>%
dplyr::arrange(dplyr::desc(V3)) %>%
head()

V1 V2 V3

1 YER063W_mRNA 100 94.07746

2 YER063W_mRNA 175 40.25282

3 YER063W_mRNA 99 34.23803

4 YER063W_mRNA 118 13.88028

5 YER063W_mRNA 101 11.56690

6 YER063W_mRNA 49 11.10423

可以看到 density 最高的是 codon pos 100 的位置，我们看看氨基酸序列 100 位置：

❝
实际上是 101 的位置，还是有些出入的，再一个和文章分析方法也有区别，我们使用的是 5‘end 去做分析的，文章使用的 3’end。怎么说呢，有想法的欢迎交流讨论。

结尾

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