鉴定差异翻译效率基因之 deltaTE 下篇

1、回顾

上期大概介绍了一下文章内容，以及 Ribo-seq 数据分析的基本流程，那么如何评估一个良好的 Ribo-seq 的测序结果呢，文章使用了 ribotish 这个软件对比对好的 bam 文件进行质控：

A.片段长度分布：

此文章图.^[1]

Beyond Read-Counts:Ribo-seq Data Analysis toUnderstand the Functions of the Transcriptome.^[2]

一个良好的 ribo-seq 测序数据的片段分布应该在 29nt 左右之间，因为核糖体保护的 RNA 片段长度就在 28-30nt 左右，当然由于物种的不同也会有一点的偏差。

B.3nt 周期性

核糖体翻译的过程是从 ATG 起始密码子开始翻译，并以 3 个碱基为一个密码子，在 RNA 上 3 个碱基的周期移动，Ribo-seq 是细胞内部翻译状态的一个快照，富集到的所有片段应该也是 3 的倍数，所以比对到的序列密度分布也应该是周期性的变化。此外，RPFs 主要分布在 CDS 区。

图中柱形图有 3 个 frame，分别代表所有密码子的第 1、2、3 个碱基，大部分是落在 frame 0 的位置。

2、实战

在确定好我们的 Ribo-seq 数据的质量较好后，定量然后就可以拿定量结果做差异翻译效率分析来寻找不同条件下的翻译被调控的基因了。

输入数据：

1、RNA-seq 和 Ribo-seq 的原始 count 矩阵：

2、样本分组表：

2.1、安装

首先需要我们在 linux 里安装 DESeq2：

$ R
R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> BiocManager::install("DESeq2")

安装好后输入 q() 退出 R 环境。

$ git clone https://github.com/SGDDNB/translational_regulation.git

下载完后进入文件夹：

$ cd translational_regulation && ls
DTEG.R  README.md  goi_viz.R  sample_data

2.2、使用

查看 DTEG.R 用法：

$ Rscript DTEG.R arg1 arg2 arg3 arg4 arg5 arg6
where command arguments arg1-6 are as follows:
Argument 1 (arg1): Ribo-seq count matrix file path
Argument 2 (arg2): RNA-seq count matrix file path
Argument 3 (arg3): Sample information file path
Argument 4 (arg4): Batch effect covariate: yes=1,or no=0
Argument 5 (arg5): Save Rdata file as a record for future use (optional, Default = 1)
Argument 6 (arg6): Verbose mode (optional,Default = 0)

# Example:
$ Rscript DTEG.R ./ribo_counts.txt ./rna_counts.txt
./sample_info.txt 1

# 进入示例数据文件夹
$ cd sample_data/ && ls
ribo_counts.txt  rna_counts.txt  sample_info.txt

可以看到该文件夹下已经有需要的 3 个数据了，我们一步法搞定分析：

$ Rscript DTEG.R sample_data/ribo_counts.txt sample_data/rna_counts.txt sample_data/sample_info.txt 1
$ ls -l
total 20
-rwxrwxrwx 1 root root 10212 Apr 10 11:25 DTEG.R
-rwxrwxrwx 1 root root  3497 Apr 10 11:25 README.md
drwxrwxrwx 1 root root  4096 Jul  9 20:26 Results
-rwxrwxrwx 1 root root   805 Apr 10 11:25 goi_viz.R
drwxrwxrwx 1 root root  4096 Apr 10 11:25 sample_data
# 多出来一个Results的文件夹
$ tree Results/
Results/
├── Result_figures.pdf
├── fold_changes
│   ├── deltaRNA.txt
│   ├── deltaRibo.txt
│   └── deltaTE.txt
└── gene_lists
    ├── DTEGs.txt
    ├── DTG.txt
    ├── buffered.txt
    ├── exclusive.txt
    ├── forwarded.txt
    └── intensified.txt

结果很快就出来了，我们可以看到：

a、 fold_changes 文件夹下有 3 个文件，分别为 RNA-seq、Ribo-seq 和 deltaTE（差异翻译效率）的差异结果。我放在一起查看：

b、 gene_lists 文件夹下则是分的 6 类差异结果数据，里面都是筛选出来的 gene id，和你的输入文件的行名一致。

c、 Result_figures.pdf 是绘图结果：

2.3、单个基因绘图

作者还提供了 Run goi_vis.R 脚本对单个基因进行在不同条件下的倍数变化的绘图，类似于上图的 D-G 图：

用法：

$ Rscript goi_viz.R arg1 arg2 arg3 arg4
where command arguments arg1-6 are as follows:
Argument 1 (arg1): Ribo-seq fold change file path.
Argument 2 (arg2): RNA-seq fold change file path.
Argument 3 (arg3): TE fold change file path.
Argument 4 (arg4): ENSEMBL gene ID

我们尝试一下画，挑一个 RPFs 显著上调的 ENSG00000134198 基因：

$ Rscript goi_viz.R Results/fold_changes/deltaRibo.txt Results/fold_changes/deltaRNA.txt Results/fold_changes/deltaTE.txt ENSG00000134198
null device
          1

$ ls -l
total 28
-rwxrwxrwx 1 root root 10212 Apr 10 11:25 DTEG.R
-rwxrwxrwx 1 root root  4436 Jul  9 21:07 ENSG00000134198.pdf
-rwxrwxrwx 1 root root  3497 Apr 10 11:25 README.md
drwxrwxrwx 1 root root  4096 Jul  9 20:42 Results
-rwxrwxrwx 1 root root   805 Apr 10 11:25 goi_viz.R
drwxrwxrwx 1 root root  4096 Apr 10 11:25 sample_data

输出一个 ENSG00000134198.pdf 文件：

可以看到 RPFs 和 RNA 水平都是上升的，但 TE 不显著，应该属于 DTG 类型。

最后放一张作者的分类示意图：

作者还提供了一步一步分析的代码，感兴趣的可以去看文章内容，这里就不放了。

参考资料