寻找核糖体暂停位点

引言

❝
核糖体暂停或停滞是指在蛋白质合成过程中，核糖体在 mRNA 上移动的速度减慢或暂时停止。这种现象可以由多种因素引起，以下是一些常见的原因：
❝
mRNA 结构：
二级结构：mRNA 分子的复杂二级结构（如发夹结构）可能会物理阻碍核糖体的移动。
稀有密码子：一些密码子在细胞内对应的 tRNA 较少，导致核糖体在这些密码子上等待合适的 tRNA。
tRNA 供给：
tRNA 短缺：某些 tRNA 在细胞中的浓度较低，导致与这些 tRNA 匹配的密码子处发生停滞。
tRNA 修饰问题：tRNA 的修饰不足可能影响其功能，从而影响翻译效率。
翻译调控：
翻译调控机制：细胞可能通过调控翻译速度来控制蛋白质的合成。例如，通过调节翻译起始或延伸阶段的因子。
信号肽：某些信号肽序列可能会减缓核糖体的移动，以便于协同翻译后修饰或定位。
翻译后修饰和折叠：
共翻译折叠：某些蛋白质在合成过程中需要共翻译折叠，可能导致核糖体停滞以便完成折叠。
翻译后修饰：需要翻译后修饰的蛋白质可能在核糖体上停留更长时间。
环境压力和细胞状态：
营养缺乏：营养缺乏或其他环境压力可能导致翻译速度降低。
应激反应：在细胞应激条件下，核糖体可能会暂停以调节蛋白质合成。
毒素和药物：
抗生素：某些抗生素通过与核糖体结合来抑制蛋白质合成，导致核糖体停滞。
内源性抑制剂：一些内源性分子可以直接与核糖体或翻译因子结合，从而抑制翻译。
核糖体暂停或停滞是细胞调控蛋白质合成的一种重要机制，它可以影响蛋白质的表达水平、折叠和功能。了解这些机制对于研究基因表达调控和开发抗生素等药物具有重要意义。
网上找了很多，发现计算暂停位点的其实很少，目前看到一个 PausePred 来计算暂停得分的（pausing score）的，发表在 RNA 期刊上：

下面是计算原理：

❝
该算法一特定大小的滑动窗口（1000nt）来计算每个位置核糖体的暂停得分，考虑相邻区间的核糖体密度作为背景。我们自己来实现这个算法，下面是自己写的 python 代码示例：

#

4_calculate pausing score according to pausepred methods

#

for key,pos_dict in trans_id_sorted.items():
pos_dict_new = {}

# loop to calculate pausing score for each position
for pos,exp in pos_dict.items():
    pos_n = pos + int(window)
    pos_n_mid = pos + int(window)/2
    pos_n_extend = pos + 1.5*int(window)

    sum_1 = sum([pp for pi,pp in pos_dict.items() if pi >= pos and pi <= pos_n])
    sum_2 = sum([pp for pi,pp in pos_dict.items() if pi >= pos_n_mid and pi <= pos_n_extend])

    if sum_1*sum_2 == 0:
        si = 0
    else:
        si = (int(window)/2)*exp*((sum_1 + sum_2)/(sum_1*sum_2))

    # save in dict
    pos_dict_new[pos] = [exp,si]

# save in dict
trans_id_sorted[key] = pos_dict_new

R 代码示例，测试感觉循环运行速度较慢：
plyr::ldply(seq_along(tid),function(x){
tmp2 <- subset(filtered_df,trans_id == tid[x]) %>%
dplyr::arrange(trans_pos)

  # calculating pausing score
  # window = 100
  # x = 1
  purrr::map_df(1:nrow(tmp2),function(x){
    tmp3 <- tmp2[x,]

    pos <- tmp3$trans_pos

    sum_density_window <- tmp2 %>%
      dplyr::filter(trans_pos >= pos & trans_pos <= pos + 100)
    sum_density_window <- sum(sum_density_window$norm_exp)


    sum_density_window2 <- tmp2 %>%
      dplyr::filter(trans_pos >= pos + 100/2 & trans_pos <= pos + (3/2)*100)
    sum_density_window2 <- sum(sum_density_window2$norm_exp)

    si <- (100/2)*tmp3$norm_exp*((sum_density_window + sum_density_window2)/(sum_density_window*sum_density_window2))

    # add pasuing score
    tmp3$pasuing_score <- si

    return(tmp3)
  }) -> pasuing_df

  return(pasuing_df)
}) -> calculated_df

最后采用 python 来实现这个功能。
安装

install.packages(“devtools”)

devtools::install_github(“junjunlab/RiboProfiler”)

or

remotes::install_github(“junjunlab/RiboProfiler”)

library(RiboProfiler)
测试

前面我们复现的文章可以找到一些有明显暂停位点的基因，作者也展示了几个，我们通过计算暂停得分来找找其他的基因看看。
前面预处理还是一样的：

perform QC analysis

pre_qc_data(longest_trans_file = “longest_info_extend.txt”,
bam_file = bam_file,
out_file = paste(sample_name,”.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

offset <- PsiteOffsetCheck(qc_data = qc_df,read_length = c(28,28))
offset$plot
offset$summary_offset

summary_offset <- offset$summary_offset
summary_offset$Offsets <- -12

sample readLengths Offsets bamFiles bamLegends

1 eIF5Ad 28 -12 eIF5Ad eIF5Ad

2 WT 28 -12 WT WT

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)

function test

norm_df_sub <- norm_df %>%
dplyr::filter(sample %in% c(“WT”,”eIF5Ad”))
❝
get_pausing_site 函数来计算暂停得分，我这里窗口大小设置为 500nt, 得到结果后可以筛选 pausing_score>50 的位置看看：:

test

pause_df <- get_pausing_site(longest_trans_file = “longest_info_extend.txt”,
normed_file = norm_df_sub,
min_counts = 64,
window = 500)

check

head(pause_df,3)

# A tibble: 3 × 5

gene_name trans_id trans_pos norm_exp pausing_score

1 CCC2 YDR270W_mRNA 435 14.4 9.15

2 CCC2 YDR270W_mRNA 441 14.4 9.20

3 CCC2 YDR270W_mRNA 444 14.4 9.25

pause_df_ft <- pause_df %>%
dplyr::filter(pausing_score >= 50)

再结合前面提到的 track_plot 函数选几个基因可视化看看是否真的和结果一样：
track_df <- get_track_df(longest_trans_file = “longest_info_extend.txt”,
select_gene = c(“CCC2″,”RAM2″,”PMA1″,”PAP2”),
normed_file = norm_df_sub)

check

head(track_df,3)

sample gene_name trans_id transpos density type

1 eIF5Ad CCC2 YDR270W_mRNA 1506 14.35714286 ribo

2 eIF5Ad RAM2 YKL019W_mRNA 702 4.59420290 ribo

3 eIF5Ad PMA1 YGL008C_mRNA 898 0.03012654 ribo

track_plot(signal_data = track_df,
gene_anno = “longest_info_extend.txt”)

可以，有点意思。
再看几个：
track_df <- get_track_df(longest_trans_file = “longest_info_extend.txt”,
select_gene = c(“ZWF1″,”IES2″,”ECM38″,”PAC10”),
normed_file = norm_df_sub)

track_plot(signal_data = track_df,
gene_anno = “longest_info_extend.txt”)

结尾

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