鉴定差异翻译效率基因之 Riborex

1、介绍

今天介绍一个 Ribo-seq 中用于快速鉴定差异翻译效率基因的软件：Riborex。文章在 2017 年发表在 Bioinformatics 上：

标题：

Riborex: fast and flexible identification of differential translation from Ribo-seq data^[1]

Riborex 基于在已有的差异分析软件 DESeq2、edgeR、Voom 的框架下，根据负二项模型（negative binomial distribution）和广义线性模型（GLM）来寻找翻译效率差异的基因。作者对不同的翻译效率计算软件和 Riborex 里 3 种不同的方法鉴定的基因进行 overlap，都有大部分的重叠：

此外对不同数据集的 Ribo-seq 数据使用不同软件分析，相对于之前的几款软件，所耗费的时间也是最快的：

不同的数据集之间，Riborex 的准确率相对于其它软件也有着良好的效果，此外 Xtail 软件也表现不错：

总之，Riborex 是作为一个开源的 R 包实现的，它的准确性与现有的最精确的方法不相上下，但是速度要快得多。

2、安装

包的信息及软件存放在 githup 网址上：

smithlabcode/riborex：^[2] https://github.com/smithlabcode/riborex

作者在 3 个月前还更新了 README 文件，此外在 vignetts 文件夹里有个 manual 使用帮助文档，可以参考使用。

该软件依赖 DESeq2、edgeR 和 fdrtool 3 个包，提前装好，作者强烈建议使用 conda 安装：

$ conda install -c bioconda r-riborex

R 里安装：

# 1、安装包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
install.packages('fdrtool')
# 从githup安装
library(devtools)
options(unzip='internal')
devtools::install_github('smithlabcode/riborex')

# 2、加载R包
library(riborex)

3、使用

输入文件需要 RNA-seq 和 Ribo-seq 的定量 count 表达矩阵，可以使用 HT-seq 或 featureCount 等软件得到，接下来使用软件来鉴定差异翻译效率的基因：

# 加载测试数据
data(riborexdata)
# 获取rna-seq和ribo-seq count
RNACntTable <- rna
RiboCntTable <- ribo
# 查看数据
head(RNACntTable,3)
                   BN_336 BN_337 BN_338 BN_339
ENSRNOG00000000017      7     11      4      4
ENSRNOG00000000024   2467   2478   3258   2316
ENSRNOG00000000033    206    282    330    244
head(RiboCntTable,3)
                   BN_341 BN_342 BN_343 BN_344
ENSRNOG00000000017     15      5     10      2
ENSRNOG00000000024   5206   5921   2864   1985
ENSRNOG00000000033     30     30     23     13
# 分组
rnaCond <- c("control", "control", "treated", "treated")
riboCond <- c("control", "control", "treated", "treated")
# 差异分析，使用DESeq2
res.deseq2 <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, engine = "DESeq2")
DESeq2 mode selected
combining design matrix
applying DESeq2 to modified design matrix
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

# 查看结果
head(res.deseq2,3)
log2 fold change (MLE): EXTRA1 treated vs control
Wald test p-value: EXTRA1 treated vs control
DataFrame with 3 rows and 6 columns
                     baseMean log2FoldChange     lfcSE      stat    pvalue      padj
                    <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
ENSRNOG00000000017    7.69493       1.533227  1.566415  0.978813  0.327672  0.736753
ENSRNOG00000000024 3544.23807      -0.102151  0.202054 -0.505564  0.613163  0.901926
ENSRNOG00000000033  111.43945       0.265252  0.503397  0.526924  0.598246  0.896407

检查 p 值的分布：

• 良好的 p 值分布图应该是在 0 侧有个最高的峰，往后面则是均匀分布，具体可参考：什么，你算出的 P-value 看上去像齐天大圣变的庙？：

# 检查p值分布
hist(res.deseq2$pvalue, main = 'DESeq2 unadjusted p-values', xlab='Unadjusted p-values',
     col = '#F4C7AB')

统计差异结果：

# 统计差异结果
summary(res.deseq2)
out of 13916 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 1107, 8%
LFC < 0 (down)     : 1217, 8.7%
outliers [1]       : 0, 0%
low counts [2]     : 540, 3.9%
(mean count < 6)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

保存差异结果：

# 保存结果
write.table(res.deseq2, "riborex_res_deseq2.txt", quote=FALSE,sep = 't')

使用 edgeR 来差异分析：

# 使用edgeR来差异分析
res.edgeR <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "edgeR")
edgeR mode selected
combining design matrix
applying edgeR to modified design matrix

# 查看结果
head(res.edgeR$table,3)
                         logFC    logCPM         LR    PValue       FDR
ENSRNOG00000000017  1.36290149 -2.456968 1.80649827 0.1789289 0.4627345
ENSRNOG00000000024 -0.30127172  6.212404 2.13670654 0.1438103 0.4037250
ENSRNOG00000000033  0.07178854  1.313235 0.02631769 0.8711269 0.9636097

# 使用edgeRD来对RNA-seq和Ribo-seq分别评估离散度
res.edgeRD <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "edgeRD")

使用 limma 包的 Voom 来进行差异分析：

# 使用limma包的Voom来进行差异分析：
res.voom <- riborex(RNACntTable, RiboCntTable, rnaCond, riboCond, "Voom")
Voom mode selected
combining design matrix
applying Voom to modified design matrix

# 查看差异结果
head(res.voom,3)
                         logFC    AveExpr          t   P.Value adj.P.Val         B
ENSRNOG00000000017  1.39674189 -2.8437969  1.2847048 0.2268559 0.5243467 -4.928058
ENSRNOG00000000024 -0.30047073  6.0075571 -1.8737990 0.0893910 0.2991861 -5.682138
ENSRNOG00000000033  0.07661457  0.7070877  0.1687028 0.8692764 0.9540355 -6.320171

4、矫正 p 值

对于我们的 p 值分布如果不是那么理想，很有可能鉴定出来的差异基因数量会很少，作者还提供了矫正 p 值的解决办法，让 p 值分布尽可能 “correct” ：

# 操作基本一致
# 载入数据
RNACntTable.corrected <- rna.null
RiboCntTable.corrected <- ribo.null
# 分组
rnaCond.corrected <- c(rep('T0', 3), rep('T24',3))
riboCond.corrected <- rnaCond.corrected
# 差异分析
res.deseq2.corrected <- riborex(RNACntTable.corrected, RiboCntTable.corrected,
                                rnaCond.corrected, riboCond.corrected)
DESeq2 mode selected
combining design matrix
applying DESeq2 to modified design matrix
estimating size factors
estimating dispersions
gene-wise dispersion estimates
...

# 检查p值分布
hist(res.deseq2.corrected$pvalue, main = 'DESeq2 unadjusted p-values', xlab='Unadjusted p-values',
     col = '#F4C7AB')

可以看到这个 p 值分布是向右倾斜逐渐升高的，p < 0.05 基因就会很少，接下来用 fdrtool 进行矫正：

# 矫正
results.corrected <- correctNullDistribution(res.deseq2.corrected)
Step 1... determine cutoff point
Step 2... estimate parameters of null distribution and eta0
Step 3... compute p-values and estimate empirical PDF/CDF
Step 4... compute q-values and local fdr
# 检查矫正后的p值分布
hist(results.corrected$pvalue, main = 'DESeq2 unadjusted p-values after correction',
     xlab='Corrected unadjusted p-values',col = '#F4C7AB')

5、多因子实验设计的差异分析

此外，此软件还支持多个实验条件设计的差异翻译效率的分析：

# 多因子实验设计
rna <- RNACntTable[,c(1,2,3,4,1,2,3,4)]
ribo <- RiboCntTable[,c(1,2,3,4,1,2,3,4)]
# 分组
rnaCond <- data.frame(factor1=(c("control1", "control1", "treated1", "treated1",
                                 "control1", "control1", "control1", "control1")),
                      factor2=(c("control2", "control2", "control2", "control2",
                                 "control2", "control2", "treated2", "treated2")))
riboCond <- data.frame(factor1=(c("control1", "control1", "treated1", "treated1",
                                  "control1", "control1", "control1", "control1")),
                       factor2=(c("control2", "control2", "control2", "control2",
                                  "control2", "control2", "treated2", "treated2")))

进行差异分析需要提供指定的比较分组：

# 进行差异分析需要提供指定的比较分组
contrast = c("factor2", "control2", "treated2")
# 差异分析
res.deseq2 <- riborex(rna, ribo, rnaCond, riboCond, "DESeq2", contrast = contrast)
# 查看统计结果
summary(res.deseq2)
out of 13916 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 1887, 14%
LFC < 0 (down)     : 1987, 14%
outliers [1]       : 0, 0%
low counts [2]     : 270, 1.9%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

# edgeR和edgeRD类似
res.edgeR <- riborex(rna, ribo, rnaCond, riboCond, "edgeR", contrast = contrast)
res.edgeRD <- riborex(rna, ribo, rnaCond, riboCond, "edgeRD", contrast = contrast)

6、绘制火山图

画个火山图对上下调基因可视化一下：

# 转换为data.frame
dat <- as.data.frame(res.deseq2)
# 定义上下调基因
dat$type <- ifelse(abs(dat$log2FoldChange)>1 & dat$pvalue < 0.05,
                   ifelse(dat$log2FoldChange >1 & dat$pvalue < 0.05,'UP','DOWN'),'nonSig')
# 查看上下调基因数量
table(dat$type)
DOWN nonSig     UP
   983  12011    922

# 绘图
library(ggplot2)
library(ggprism)

ggplot(dat,aes(x = log2FoldChange,y = -log10(pvalue))) +
  geom_point(aes(color = type,shape = type),alpha = .5,size = 3) +
  theme_prism() + ylim(0,80) + xlim(-6,6) +
  theme(legend.position = c(0.9,0.9)) +
  scale_color_manual(values = c('UP'= '#F54748','nonSig'='#F8A488','DOWN'='#5AA897')) +
  geom_hline(yintercept = -log10(0.05),linetype = 'dashed',size = 1) +
  geom_vline(xintercept = c(-1,1),linetype = 'dashed',size = 1)

7、实验没有 replicate 能分析吗？

对于实验没有 replicate 是不能做统计学检验分析的，但是 edgeR 是可以对 RNA-seq 在没有生物学重复的情况下进行差异分析，对于差异翻译效率分析的话，我也不知道 Riborex 是否能对于没有重复的条件下进行差异分析。

于是昨天在 githup 上问了作者这个问题，今天中午便很快得到了答复！

很不幸的是，不可以哟。

参考资料