分析差异翻译基因时,为了保证结果的正确性,除了使用Xtail
软件外,也用riborex
分析了一遍,该软件封装了常规转录组常用的三种差异分析框架DESeq2
、edgeR
、Voom
。另外,该软件还封装了fdrtool
可以用来矫正pvalue
。具体的代码示例可以参考软件文档,这里就不演示了,主要说说使用软件时需要注意的事项。
详细的使用说明书:
https://github.com/smithlabcode/riborex/blob/master/vignettes/riborex.pdf
DESeq2
做为常规转录组最常用的差异表达基因分析的框架,确实是一个不可多得的选择。同时,也对翻译组差异翻译基因的鉴定有着很大的影响,很多翻译组的差异分析软件也是基于这些框架而来。因此,这里说的注意事项是选择DESeq2
做为riborex
的分析引擎时,需要注意的地方。
用DESeq2
做差异分析时,最后提取差异结果时,默认设计矩阵中前面的分组做为参考,也可以提供contrast
参数,指定如何提取结果,确定哪一个分组做为参考,弄不好得到的结果可能跟预期相反。例如,分组变量名为group
,含有case
、ctrl
两个水平,那么最后提取差异时,可以用两种方式指定:
res <- results(dds, contrast=c("group", "case", "ctrl"))
# 或者下面的方式
res <- results(dds, name="group_case_vs_ctrl")
要明白放在后面的组名如ctrl
是做为对照组的,如果case
放在后面就会导致上下调基因反过来。虽然不是什么大问题,但有可能引起误导而浪费不必要的时间。
那么,riborex
使用DESeq2
做为引擎分析差异基因时,最后也要根据contrast
提取差异结果。例如,软件文档示例数据的分组情况如下:
rnaCond <- c("control", "control", "treated", "treated")
riboCond <- c("control", "control", "treated", "treated")
虽然软件也提供了contrast
参数用于指定如何提取,但好像该参数并那么好使。默认情况下,contrast
无需指定,与DESeq2
默认情况一样,前面的分组做为参考,即control
为对照。但是,类似文档中只有Multi-factor experiment
的示例格式指定这个参数:
然而,不是Multi-factor experiment
,该怎么办呢?比如样本分组情况如下:
rnaCond <- c("case", "case", "ctrl", "ctrl")
riboCond <- c("case", "case", "ctrl", "ctrl")
通常习惯以对照组为参考,看处理组如何变化,默认的结果是以case
做为参考,这显然有点反人类。如果想得到以ctrl
为参考的结果,该怎么办呢?
因为软件有contrast
参数,所以第一时间想到用这个参数提取想要的结果,可是看了一圈发现不知道如何提供参数,也是无奈。当然,此路不通还有别的选择,比如可以将样本顺序颠倒一下就可以了,这样默认情况就是想要的结果。
没办法,咱还是想搞清楚怎么提供参数,就去看了看软件的源码,关于contrast
参数的代码如下:
这段代码的else
部分代码写法有点怪异,提供的contrast
第一个元素要在combinedCond
列名里面,而咱们根本不知道combinedCond
内容是什么。既然倒腾了,那就倒腾明白吧,可以通过下面的方式提供参数:
contrast <- c('cond', 'case', 'ctrl')
res <- riborex(rna, ribo, rnaCond, riboCond, "DESeq2", contrast = contrast)
contrast
第一个元素得是cond
,已经由代码的逻辑固定了。
事实上,riborex
只是封装了DESeq2
,咱们也可以自行构建设计矩阵,直接使用DESeq2
来做差异翻译基因分析。
protocol + condition + condition:protocol
protocol
为RNAseq
和riboseq
各样本的分组向量,condition
为样本所属的数据类型分组向量,condition:protocol
为前面两个的交互效应。
riborex
与xtail
的结果重合度挺好,而xtail
的结果更为保守,想要多一些的结果可以选择前者。