Footsteps on my way !
perl/linux/测序分析

用R实现批量差异分析(t检验和方差分析),自己算P值

对于二代数据的表达差异分析,理论上应该用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 <- “NA”

log2_FC<- “NA”

}else{

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

Pvalue<-y$p.value

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

}

}

# 对p value进行FDR校正

fdr=p.adjust(Pvalue, “BH”)

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

out<-cbind(a,log2_FC,Pvalue,fdr)

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

 

 

##########################我的理解##################################
我感觉代码有问题,经过我改完以后代码如下(我的数据是2个重复,我的测试数据: test_of_yangl.txt ):
最后我觉得这么算的P值不靠谱,不推荐使用自己算P的值

##setwd(“F:/”)
a=read.table("RA24C.txt",header=T,sep="\t") 
a
#预生成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:3])==0&&sd(a[i,4:5])==0){
 Pvalue[i] <-"NA"
 log2_FC[i]<-"NA"
 }else{
 y=t.test(as.numeric(a[i,2:3]),as.numeric(a[i,4:5]))
 Pvalue[i]<-y$p.value
 log2_FC[i]<-log2((mean(as.numeric(a[i,2:3]))+0.001)/(mean(as.numeric(a[i,4:5]))+0.001)) 
 }
}
# 对p value进行FDR校正
fdr[i]=p.adjust(Pvalue[i], "BH") 
# 在原文件后面加入log2FC,p value和FDR,共3列;
out<-cbind(a,log2_FC,Pvalue,fdr) 
write.table(out,file="RA24C.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<-factor(c(rep(“c”,3),rep(“m”,3)))

# 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] <- “NA”

log2_FC[i]<- “NA”

}else{

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

Pvalue[i]<-summary(y)[[1]][,5][1]

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

}

}

# 对p value进行FDR校正

fdr=p.adjust(Pvalue, “BH”)

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

out<-cbind(a,log2_FC,Pvalue,fdr)

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

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

脚本和测试数据我们已经上传了,请点击“阅读原文”跳转到论坛下载。

尊重他人劳动成果,转载请注明出处:Bluesky's blog » 用R实现批量差异分析(t检验和方差分析),自己算P值

分享到:更多 ()

评论 4

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址
  1. #0

    请问作者,方差分析中Pvalue[i]<-summary(y)[[1]][,5][1],改代码中[[1]][,5][1]是什么意思?还望告知,十分感谢

    jiandan4年前 (2017-12-11)回复
  2. #0

    你好,方差分析Pvalue[i]<-summary(y)[[1]][,5][1]中,[[1]][,5][1]怎么理解啊?非常感谢

    abcd3年前 (2018-10-07)回复
  3. #0

    一个错误的代码编写,还给转到到处都是,害人不浅

    • 杨李

      哪儿有问题可以指出来,大家共同改正一下

      杨李2年前 (2019-04-19)回复