CORNAS：一种快速简单鉴定无重复转录组差异基因的方法

还记得上次文章的最后提到CORNAS这种方法吗？最近刚好在Github上看到了这个项目，就花了点时间看了下文档感觉操作也比较简单，这里记录一下使用过程，大家共同学习一下。

介绍

该软件17年发表在BMC上，是一种快速贝叶斯方法，它可以根据样品的覆盖参数来估计真实基因计数的后验分布，从而提高计算差异表达基因的准确性，更多的原理大家可以去看文章或者中文解读。

image

下载

使用软件之前，我们最好先看下该软件Github中的Readme文件，看一下我们需要下载哪些东西并且怎么使用。下图第一段中内容可以得知CORNAS是用R语言写用来快速计算无重复转录组数据差异基因，我们可以直接命令行模式中运行R脚本，也可以在R语言编辑器(常用Rstudio)里运行，这个包包含了俩个主要文件：CORNAS.R和cornas.config。软件测试基于R版本3.0，运行基因列表数据也是基于mapping的原始数据，无需标准化。这里我们先将这个文件压缩包下载下来。

image-20200929205449914

使用方法

1. 方法一 (Shell中运行):

    Rscript CORNAS.R <config> <datatable>
   

2. 方法二(Rstudio中运行)

    source("/path/to/CORNAS.R")
    cornas("/path/to/config" , "/path/to/datatable")
   

根据方法二我们可以看出CORNAS软件其实是作者写的由一些函数构成的R脚本，并非是个R包，我们就先加载CONRNAS.R脚本中的函数，然后通过其中的cornas函数运行即可，核心参数为俩个文件：config配置文件和datatable数据文件。所以我们只要按照规定的文件格式将自己的文件整理一致就可以轻松使用了，还是蛮简单的。

CORNAS.R

图片.png

贴心的是作者在文档的example_run目录里准备了这俩个示例文件，我们在使用前可以先运行一下了解运行过程，之后再替换文件内容即可。 OK，我们先分别看下这几个文件的内容格式，然后运行一遍。

图片.png

1. datatable数据文件

   文件有四列信息，第一列可以是基因/转录本的ID，第二列为对应的长度信息，第三列和第四列分别为没有重复的俩样本。注意文件不能包含标题行，而且第二列的长度信息可以不需要，文件分隔以Tab形式分隔。

   
   图片.png

2. contig配置文件

   这个软件比较特别的地方就在于此，也是它的核心部分，我们需要简单配置一个文件，运行自己的数据时候只需要改动我圈出的部分即可，上面对应的数字代表所在的列数，而下面的每个样品的覆盖度。

   关于覆盖度的概念，作者提到(The sequencing coverage is the number of total reads (observed counts) divided by the actual amount of fragments present in the PCR mix)，意思即为样品存在的实际片段数量除以观察到的总数值，这里他提供了俩种选择：1. 如果你提前已经计算了样品A和B的覆盖度值，那就直接填写在配置文件中 ； 2. 如果没有这俩个样品的覆盖度也OK，这里程序会根据你datatable文件中每个样品的总reads量计算这个覆盖度的值。

除此之外，作者另外还有俩个可选的参数，我们可以增添到配置文件下方 (其实建议默认就好)：

Alpha：如果你想改变alpha值从默认的99%。降低该参数的值可以增加鉴定差异表达基因的灵敏度，数量会变多。

Fold change：差异倍数，默认是1.5，降低该值也可以增加差异表达基因的数量。

image

运行

方式一

Linux中以命令行的模式运行

Rscript CORNAS.R cornas.config.test4 test4_kidneyliver_example.tab >cornas_test4_example2.out

方式二

Rstudio中运行

source("/path/to/CORNAS.R")
cornasExample1 <- cornas("cornas.config.test4" , "test4_kidneyliver_example.tab")

以上俩种方法根据个人运行环境自由选择，代码都很简单一行搞定，运行速度也非常快大概10几秒就可以了，我们现在直接看下结果文件都包含了哪些东东，也就是对应示例文件中的cornas_test4_example2.out文件。

image-20200929220429044

结果文件就是上图一排**************号下面的内容，上面内容就是你当时设置的参数及一些重要的值，我们重点看下面这些，从左到右有十列内容，重点其实为前4列的内容。

1. Gene_Name = 输入文件的基因/转录本ID
2. DEG_call = 通过设置阈值鉴定该基因是否为差异基因: 是 (Y) or 否 (N).
3. Express_higher = 如果 DEG_call == Y, 显示哪个样品的表达量最高 (A or B), 否则为 "-".
4. Fold_difference = 高表达的样品与另一样品之间的差异倍数fold difference/change.、
5. A_O-count = 输入文件中样品A的表达量
6. B_O-count = 输入文件中样品B的表达量
7. A_T-lower = 样本A中真实计数后验分布的下界。
8. A_T-upper = 样本A中真实计数后验分布的上界。
9. B_T-lower = 样本B中真实计数后验分布的下界.
10. B_T-upper = 样本B中真实计数后验分布的上界。