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

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循环得到的结果是一致的

方法三、使用rstatixreshape2

#如果没有安装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值

获取m6a_expr_with_type.txt ☞【R语言】三种批量做T检验的方法

参考资料:

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

☞统计学知识大梳理

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

☞R入门教程——cookbook for R

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

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

【R语言】P值转换成***

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,657评论 6 505
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,889评论 3 394
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,057评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,509评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,562评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,443评论 1 302
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,251评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,129评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,561评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,779评论 3 335
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,902评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,621评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,220评论 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,838评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,971评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,025评论 2 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,843评论 2 354

推荐阅读更多精彩内容