假如现在你手头有如下文件(test.txt),只有俩样品RPKM_A (对照) 和RPKM_B (处理), 值为标准化后的RPKM。
1. 根据foldchange直接筛选
之前在一篇中文文献中见到有人用这种方法,作者自定义差异基因的标准:至少有一组RPKM值大于5,且满足foldchange(差异倍数) > 2,我们可以在LInux中直接可以用awk进行过滤,其实Excel、R中也可以操作,根据个人习惯吧。代码如下:
### 上调基因######## # 提取B组大于等于5,A组等于0的基因。 less test.txt | gawk '{if (($2==0)&&($3>=5)) print $0}' > up.txt # 提取A、B俩组至少有一组大于等于5,且B组值/A组值大于等于2 less test.txt | gawk '{if (($2!=0)&&($3!=0)) print $0}'|gawk '{if (($2>=5)||($3>=5)) print $0}'|sed '1d'|gawk '{if ($3/$2>=2) print $0}' >> up.txt ### 下调基因######### # 提取A组大于等于55,B组等于0的基因 less test.txt | gawk '{if (($2>=5)&&($3==0)) print $0}' > down.txt # 提取A、B俩组至少有一组大于等于5,且A组值/B组值大于等于2 less test.txt | gawk '{if (($2!=0)&&($3!=0)) print $0}'|gawk '{if (($2>=5)||($3>=5)) print $0}'|sed '1d'|gawk '{if ($2/$3>=2) print $0}' >> down.txt
2. edgeR包
这种方法我在提到过,edgeR包可以做无重复的差异分析,不过需要认为指定一个dispersion值(设置BCV值),这样得到的结果比较主观,不同的人就可以有不同的结果。通常如果是实验控制的好的人类数据,那么选择BCV=0.4,比较好的模式生物选择BCV=0.1。参考
代码如下:
library(edgeR) ##跟DESeq2一样,导入数据,预处理(用了cpm函数) exprSet<- read.table(file = "test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE) group_list <- factor(c(rep("Contral",1),rep("Treat",1))) ##设置分组信息,并做TMM标准化 exprSet <- DGEList(counts = exprSet, group = group_list) bcv = 0.1 #设置BCV值 et <- exactTest(exprSet, dispersion=bcv^2) write.csv(topTags(et, n = nrow(exprSet$counts)), 'result.csv', quote = FALSE) #输出主要结果
结果文件如下:
3. Gfold软件
地址:https://zhanglab.tongji.edu.cn/softwares/GFOLD/index.html
Gfold软件应该是做没有生物学重复样本首选的软件,该软件由同济大学开发的,网站 往下拉可以看到该软件的几个功能,其中Example3为鉴定无重复的数据的差异基因。另外,这个软件不支持Windows 版本,是基于linux的一个安装软件,所以我们需要在linux环境下使用
安装GSL
使用Gfold之前必须先安装Gsl,如图下载最新的版本
# 安装最新版的 wget http://mirrors.ocf.berkeley.edu/gnu/gsl/gsl-latest.tar.gz tar -zxv -f gsl-latest.tar.gz .configure --prefix=/home/pub_guest/hekai/soft_ware/GFOLD/gsl-2.6/ make make check(选做) make install
安装Gfold
- 下载安装包
wget https://zhanglab.tongji.edu.cn/softwares/GFOLD/gfold.V1.1.4.tar.gz tar -zxv - f gfold.V1.1.4.tar.gz cd gfold.V1.1.4
- 我们打开目录下的README文件,可以需要执行下面俩句。注意 /your/installed/path/ 是你安装GSL的路径,把下面命令行中的替换为相应的路径即可
# 打开你的bashrc文件,再最后添加 export CXXFLAGS="-g -O3 -I/your/installed/path/include -L/your/installed/path/lib" export LD_LIBRARY_PATH="/your/installed/path/lib:"$LD_LIBRARY_PATH source ~/.bashrc
- 接着输入make, 发现报错了。
- 不慌,文档里也有提醒如果报错 直接输入以下命令行即可
# If it happens, follow step 1 again. If error remains, try the following command: g++ -O3 -Wall -g main.cc -o gfold -lgsl -lgslcblas -I/your/installed/path/include -L/your/installed/path/lib
- 这是我们就会发现目录下多了一个gfold软件,是可执行状态,我们.gfold -h ,可以看到该软件的帮助文档,证明此时已经安装成功了,我们将其添加为环境变量里方便我们使用。
echo 'export PATH=/your/installed/path:$PATH' >>~/.bashrc ###### **/your/installed/path/** 是你安装**gfold.V1.1.4**的路径 source ~/.bashrc
差异分析
- 我们需要准备俩个文件Control.txt和Treat.txt,我们看下处理组Control.txt文件都包含哪些,Gfold输入文件规定必须为5列,前俩列可以分别输入你的基因ID号和Symbol号,俩列内容一样其实也不影响,第三列为原始Counts值,第四列为基因长度,最后一列为标准化的RPKM值,同样对照组Treat.txt文件也按照这样准备。
- 准备完毕后,直接一条命令计算差异。
gfold diff -s1 Control.txt -s2 Treat.txt -o ControlVSTreat.csv
- OK,已经计算完了,我们看下结果文件都有哪些内容。主要一共7列信息,前两列没什么可说,就是gene symbol和gene name,第三列是GFOLD值,相当于log2(Fold Change),该值等于0的基因则记为非差异基因,非0的值才是差异基因,E-FDR是基于重复的Empirical FDR,因此无重复样本的经验FDR均为1。Log2fdc以及后面的RPKM列可以忽略考虑,因为最开始的exon的长度,我们是给定的是一个虚拟的数据。所以真正的确定差异是否显著,主要是看GFOLD值,GFOLD>0,代表case组中高表达,GFOLD<0,代表case组中低表达。后续筛选差异基因如果太多或者太少,我们也可以通过设定GFOLD的阈值来控制,比如设定<-0.3或者大于0.3。
除上述之外几种方法,还有几个比较经典的软件可供选择进行无重复数据的差异分析,比如BMC Bioinformatics发表过的一篇文章中提出了一种新的差异基因分析方法利用贝叶斯方法来推断真实基因表达数的 后验分布。其创新型之一该方法包括了由RNA样品浓度决定的覆盖度参数,之二是真实基因表达量后验分布的比较为寻找差异表达基因提供了一个参照。这种方法针对无重复样本的数据是有一定优势的,这里提供一个链接大家有兴趣的话可以去看该博主的讲解,之后有机会也会尝试一下该软件的使用进行比较。