对于二代数据的表达差异分析,理论上应该用reads counts进行计算。这个在我们论坛的专题帖已经有解释:
第14期“基因表达量计算和差异表达分析(下)”【视频】
www.omicshare.com/forum/thread-236-1-12.html
同时,在我们OS-tools已经有基于edgeR软件的差异分析工具。但依然有网友问,如果手头没有reads count数据,而只有RPKM/FPKM值该怎么办?
这个时候,就只能退而求其次,使用t检验或者方差分析。但注意,这两种检验是基于正态分布的检验方法,是不适用于二代数据的,对低丰度基因的检验会产生大量假阳性。不到万不得已不要使用这类方法。
如果非使用不可,可以:
使用以下的R脚本进行批量差异检验(t检验或方差分析);
请将在两组样本中表达量RPKM值均低于1的基因过滤掉(t检验和方差分析在低丰度基因中,假阳性过高,P value不可靠)。
请确定你认真看了上面两点使用建议,再开始看代码。
# t检验的代码如下:
a=read.table(“all_fpkm.txt”,header=T,sep=”\t”)
#预生成2个长度与输入文件行数相同的全为0的向量,将用于存储p value和差异倍数(log2FC)
Pvalue
log2_FC
# 2~4列是处理组1,5~7列是处理组2;
#将使用循环对每一行进行t检验
#如果某一行两组的标准差都等于0,将无法进行t