t检验相信大家应该都不陌生。不管是大学里面的数理与统计，还是研究生阶段的生物统计学，里面都会提到t检验。
小编也给大家总结过一些统计学相关的知识
R语言里面也有专门做t检验的函数，t.test。
t检验的应用场景也很多，比如我们经常做的差异表达分析就可以使用t检验来做。当我们手上有很多基因的时候，该如何做t检验会更有效率呢？今天小编就给大家介绍三个批量做t检验的方法。我们这里使用的数据是
这篇文章中用到的m6a甲基化相关的16个基因在TCGA-CHOL（胆管癌）中的表达情况。其实这里我们是做了16次t检验才得到每个基因的p值的。
首先我们把16个m6a基因得表达谱读进来，最后一列为样本类型，也是我们待会做t检验时候的分组依据。具体如何得到这个表达矩阵可以参考
#读取16个m6a甲基化相关基因在CHOL中的表达量
m6a_expr_type=read.table(file="m6a_expr_with_type.txt",header=T,sep="\t",row.names=1)
#获取16个m6a基因的名字，最后一列为样本类型
m6a_sym=names(m6a_expr_type)[1:(ncol(m6a_expr_type)-1)]
方法一、原始一点的方法，for循环
#生成一个空向量来存放计算出的p值
pval=c()
#for循环16次计算每个基因的p值
for(gene in m6a_sym){
#根据type来将样本分成两组
p=t.test(m6a_expr_type[,gene]~m6a_expr_type$type)$p.value
#存放p值
pval=c(pval,p)
}
#输出p值看看
pval
方法二、使用plyr和reshape2
#如果没有安装plyr和reshape2这两个R包，先去掉下面两行的#，运行进行安装
#BiocManager::install("plyr")
#BiocManager::install("reshape2")
#加载plyr和reshape2包
library(plyr)
library(reshape2)
#melt对m6a_expr_type数据格式进行转换
ddply(melt(m6a_expr_type),"variable",
function(x) {
w <- t.test(value~type,data=x)
with(w,data.frame(statistic,p.value))
})
你会发现跟我们使用for循环得到的结果是一致的
方法三、使用rstatix和reshape2
#如果没有安装dplyr，rstatix和reshape2这三个R包，先去掉下面三行的#，运行进行安装
#BiocManager::install("dplyr")
#BiocManager::install("rstatix")
#BiocManager::install("reshape2")
#加载dplyr，rstatix和reshape2这三个R包
library(dplyr)
library(rstatix)
library(reshape2)
result=melt(m6a_expr_type) %>%
group_by(variable) %>%
t_test(value ~ type)
#输出result
result
你会发现跟前面使用for循环和ddply方法得到的结果是一样的
再给大家分享两个小技巧，在计算原始p值的同时，我们还能计算校正之后的p值
#使用fdr方法对原始p值进行校正
result=melt(m6a_expr_type) %>%
group_by(variable) %>%
t_test(value ~ type) %>%
adjust_pvalue(method = "fdr")
你会发现在这张表的最后两列，我们得到了原始的p值和经过FDR方法校正之后的p值
在下面这张图上其实显示的是将p值转换成相应的*（星号），前面我们也给大家介绍过☞【R语言】P值转换成***
其实这里我们可以一次性通过rstatix这个包得到原始p值，FDR校正之后的p值以及转换成对应的***。
#一次性得到原始p值，FDR校正之后的p值以及转换成对应的***
result=melt(m6a_expr_type) %>%
group_by(variable) %>%
t_test(value ~ type) %>%
adjust_pvalue(method = "fdr") %>%
add_significance("p.adj")
#输出result
result
这样我们就可以直接将***画在图上了，具体画图方法可以参考
获取m6a_expr_with_type.txt ☞【R语言】三种批量做T检验的方法
参考资料：
网友评论