干货｜用R实现批量差异分析（t检验和方差分析）

对于二代数据的表达差异分析，理论上应该用reads counts进行计算。这个在我们论坛的专题帖已经有解释：  

第14期“基因表达量计算和差异表达分析（下）”【视频】

www.omicshare.com/forum/thread-236-1-12.html

同时，在我们OS-tools已经有基于edgeR软件的差异分析工具。但依然有网友问，如果手头没有reads count数据，而只有RPKM/FPKM值该怎么办？

这个时候，就只能退而求其次，使用t检验或者方差分析。但注意，这两种检验是基于正态分布的检验方法，是不适用于二代数据的，对低丰度基因的检验会产生大量假阳性。不到万不得已不要使用这类方法。

如果非使用不可，可以：

1. 使用以下的R脚本进行批量差异检验（t检验或方差分析）；

2. 请将在两组样本中表达量RPKM值均低于1的基因过滤掉（t检验和方差分析在低丰度基因中，假阳性过高，P value不可靠）。

请确定你认真看了上面两点使用建议，再开始看代码。

# t检验的代码如下：

a=read.table('all_fpkm.txt',header=T,sep='\t')  

#预生成2个长度与输入文件行数相同的全为0的向量，将用于存储p value和差异倍数（log2FC）

Pvalue<-c(rep(0,nrow(a)))>

log2_FC<-c(rep(0,nrow(a)))>

# 2~4列是处理组1,5~7列是处理组2；

#将使用循环对每一行进行t检验

#如果某一行两组的标准差都等于0，将无法进行t检验，所以这些行使用NA替代

#每一行计算得到p value和log2FC将被加入原文件的后两列；

#计算log2FC时，每组均值加0.001，是为了防止分母为0导致bug；

for(i in 1:nrow(a)){

        if(sd(a[i,2:4])==0&&a[i,5:7]==0){

                Pvalue <->

                log2_FC<->

        }else{

                y=t.test(as.numeric(a[i,2:4]),as.numeric(a[i,5:7]))

                Pvalue<>

                log2_FC<-log2((mean(as.numeric(a[i,2:4]))+0.001)>

        }

}

# 对p value进行FDR校正

fdr=p.adjust(Pvalue, 'BH') 

# 在原文件后面加入log2FC，p value和FDR,共3列；

out<>

write.table(out,file='ttest.out.xls',quote=FALSE,sep='\t',row.names=FALSE)

#方差分析的代码如下，与t检验相比，逻辑相同，只是有些细微的修改。

a=read.table('all_fpkm.txt',header=T,sep='\t')  

#预生成2个长度与输入文件行数相同的全为0的向量，将用于存储p value和差异倍数（log2FC）

Pvalue<-c(rep(0,nrow(a)))>

log2_FC<-c(rep(0,nrow(a)))>

#确定分组信息

type<>

# 2~4列是处理组1,5~7列是处理组2；

#将使用循环对每一行进行方差检验

#两组表达量都是0的基因，不检验；

#每一行计算得到p value和log2FC将被加入原文件的后两列；

#计算log2FC时，每组均值加0.001，是为了防止分母为0导致bug；

for(i in 1:nrow(a)){

        if(sum(a[i,2:4])==0&&sum(a[i,5:7])==0){

                Pvalue[i] <->

                log2_FC[i]<->

        }else{

                y=aov(as.numeric(a[i,2:7])~type)

                Pvalue[i]<>

                log2_FC[i]<-log2((mean(as.numeric(a[i,2:4]))+0.001)>

        }

}

# 对p value进行FDR校正

fdr=p.adjust(Pvalue, 'BH') 

# 在原文件后面加入log2FC，p value和FDR,共3列；

out<>

write.table(out,file='anova.out.xls',quote=FALSE,sep='\t',row.names=FALSE)

另外，如果没有R基础的同学，看这一期视频，学完就基本可以保证正确使用这个脚本：

第12期在线交流 R语言入门——软件简介与实操【视频】

www.omicshare.com/forum/thread-133-1-1.html

备注：理论上如果你理解了这个方法，可以进行任何组间的差异检验。只要把检验方法的部分的相关语句替换就可以了，例如可以替换为秩和检验、相关性检验等；