【R语言】三种批量做T检验的方法

生信交流平台

发布于 2022-09-21 17:43:12

1.5K0

发布于 2022-09-21 17:43:12

文章被收录于专栏：用户7627119的专栏用户7627119的专栏

t检验相信大家应该都不陌生。不管是大学里面的数理与统计，还是研究生阶段的生物统计学，里面都会提到t检验。

小编也给大家总结过一些统计学相关的知识

☞统计学中数据分析方法汇总！

☞统计学知识大梳理

☞100个统计学 & R语言学习资源网站

R语言里面也有专门做t检验的函数，t.test。

☞R入门教程——cookbook for R

☞R语言入门-工欲善其事必先利其器

t检验的应用场景也很多，比如我们经常做的差异表达分析就可以使用t检验来做。当我们手上有很多基因的时候，该如何做t检验会更有效率呢？今天小编就给大家介绍三个批量做t检验的方法。我们这里使用的数据是

☞m6a甲基化相关基因boxplot并显示p值

这篇文章中用到的m6a甲基化相关的16个基因在TCGA-CHOL（胆管癌）中的表达情况。其实这里我们是做了16次t检验才得到每个基因的p值的。

首先我们把16个m6a基因得表达谱读进来，最后一列为样本类型，也是我们待会做t检验时候的分组依据。具体如何得到这个表达矩阵可以参考

☞m6a甲基化相关基因boxplot并显示p值

#读取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甲基化相关基因boxplot并显示p值