xtail | 差异翻译基因分析

  基因翻译，即遗传密码从转录本到蛋白质的过程，也是功能基因能够发挥自身作用前的必经之路。基因的翻译效率收到很多因素的调控，这其中最直接的是基因本身的转录丰度，正常情况下，翻译速度与转录丰度成正比。所以，与普通rnaseq寻找差异表达基因不同，寻找差异翻译基因，目的是找到受功能调控而翻译速度变化的基因，首先要排除转录丰度变化引起翻译速度变化的基因。因此，差异翻译基因分析时，会同时利用普通的rna-seq与类似ribo-seq(专门用于定量正在翻译的基因的翻译丰度) 定量结果。

  翻译基因的定量核心目的是捕获ribosome protected fragment (RPF)，如下面的流程示意图所示：

  对于像ribo-seq这样的技术，由于其定量的特殊性，会有一些相应质控方式来反映数据的特征与质量。比如，三碱基周期，基因翻译过程中核糖体移动时每隔三个碱基 (一个密码子) 停顿一次，因此正确的RPF信号在CDS上应该呈现三碱基周期，即高-低-低模式如下图C：

  就像普通转录组测序一样，翻译组测序最主要的目的也是寻找差异基因。理解差异翻译之前，先说一下翻译效率(TE = RPF/mRNA)，即翻译丰度占总转录丰度的比值。很多差异软件从TE这个角度出发，分析不同条件下TE的变化(Δlog2R)来寻找差异翻译基因，比如今天要说的软件xtail也有这方面的考虑。此外，该软件还别出心裁的考虑了另外一个角度，即Δlog2FC，这个出发点假设可以理解为，如果基因的翻译丰度变化是由于实验条件导致而不是转录丰度，那么，不同条件下翻译丰度会有变化，而转录丰度没有变化。综合来看，(Δlog2R)和Δlog2FC两种方式都是为了寻找差异翻译基因，理论上应该殊途同归，但实际的做法是取两者中相对更保守的结果作为最终的结果。故而，与其他软件相比，xtail得到的结果相对保守一点。

  本来想安装releases版本的xtail，结果一直安装不上，最后选择了github上的开发版。看着参考文档的示例代码，觉得xtail使用本身不难，实际上却很快就遇到了问题：

mrna_cnt <- read.table('mrna_cnt.txt', header=T, row.names=1, seq='\t')
head(mrna_cnt, 3)
                     case1 case2 ctrl1 ctrl2
ENSMUSG00000000001.4  3968  4546  5738  6049
ENSMUSG00000000003.15    2     1     10     13
ENSMUSG00000000028.15 2203  3112  2836  3390

ribo_cnt <- read.table('ribo_cnt.txt', header=T, row.names=1, seq='\t')
head(ribo_cnt, 3)
                         case1    case2    ctrl1    ctrl2
ENSMUSG00000048696.11      139       45       50      188
ENSMUSG00000028830.14        2       89      459       31
ENSMUSG00000020530.14        5      310      868      418

condition <- c("case","case","ctrl","ctrl")
res_xtail <- xtail(mrna_cnt , ribo_cnt, condition, baseLevel="ctrl", bins=1000, threads=2)
Calculating the library size factors
1. Estimate the log2 fold change in mrna
converting counts to integer mode
2. Estimate the log2 fold change in rpf
3. Estimate the difference between two log2 fold changes
4. Estimate the log2 ratio in first condition
converting counts to integer mode
5. Estimate the log2 ratio in second condition
converting counts to integer mode
6. Estimate the difference between two log2 ratios
Error: logical subscript contains NAs

  这个问题第一次使用软件时，大概率都会遇到，比如网上就有人遇到过，也给出了一种解决方法。

  不可否认，提高minMeanCount可以解决问题，但提高到多少不好说，得看两个表达谱的数据情况。所以，简单一句话，用这个方法不直接，因为阈值设高会丢弃很多基因，设置低了还是会出错。

  那么，到底是什么引发错误的呢？确认问题才能更好地解决问题，于是乎倒腾了一番，最终确定到了问题所在：

# fold change better than ratio
fc_result <- res$pvalue_v1 > res$pvalue_v2
log2FC_determine_num <- sum(fc_result)
res[fc_result,"log2FC_TE_final"] <- res[fc_result,"log2FC_TE_v1"]
res[fc_result,"pvalue_final"] <- res[fc_result,"pvalue_v1"]
# ratio better than fold change
ratio_result <- res$pvalue_v1 <= res$pvalue_v2
log2R_determine_num <- sum(ratio_result)
res[ratio_result,"log2FC_TE_final"] <- res[ratio_result,"log2FC_TE_v2"]
res[ratio_result,"pvalue_final"] <- res[ratio_result,"pvalue_v2"]

  在pvalue_v1与pvalue_V2比较时，因为有些pvalue的值是NA，从而造成log2FC_TE_final赋值这一行代码出错。解决了这个问题，便可以正常运行并得到结果：

restab <- resultsTable(res_xtail)
head(restab, 3)
DataFrame with 6 rows and 8 columns
                      log2FC_TE_v1   pvalue_v1 ctrl_log2TE.1 log2FC_TE_v2
                         <numeric>   <numeric>     <numeric>    <numeric>
ENSMUSG00000000001.4     -1.442377 3.21985e-01    -1.5386621    -1.446969
ENSMUSG00000000028.15    -0.822741 6.82438e-01     0.0632993    -0.718147
ENSMUSG00000000037.17     0.698955 3.99256e-01     0.6652733     0.620557
                      pvalue_v2 log2FC_TE_final pvalue_final pvalue.adjust
                      <numeric>       <numeric>    <numeric>     <numeric>
ENSMUSG00000000001.4   0.388114       -1.452377     0.392985      0.999861
ENSMUSG00000000028.15  0.640160       -0.822741     0.792438      0.998961
ENSMUSG00000000037.17  0.710053        0.602557     0.700053      0.998961

  结果看上去好像没什么问题，但仔细一看与文档中resultsTable默认提取的结果相比，好像多出一列ctrl_log2TE.1，看着列名后缀盲猜应该是合并数据时有相同的列名。不过，这个应该仅是保留最终结果时的瑕疵没有影响，可以直接忽略。

  xtail返回结果是一个对象，如果想要完整的结果数据框，可以直接用res_xtail$resultsTable命令获取。现在再回过头看前面的问题，在不改代码的前提下，最好的方式应该是先过滤数据，保留在各条件下都有表达值的基因。