组间分析—T检验、R语言绘图

作者: fafu生信小蘑菇 | 来源:发表于2021-04-30 23:30 被阅读0次

组间分析—T检验、R语言绘图
R语言ggstatsplot包做T检验
R语言T检验
R语言t检验
【R画图学习21.5】单因素方差分析评估条件
生信入门学习笔记day4@2021.06.28
R语言-配对t检验
看看不同分析工具的单样本t检验（R、SPSS、Python、Ex
2021-06-18 两组间差异的非参数检验之Wilcox秩和检
宏基因组、扩增子组间差异检验，终于有人讲清楚了！！！

T检验

在做微生物群落分析时经常要寻找不同群落之间的差异，我们会经常关心两个或者多个群落之间是否显著不同。就两个微生物群落间的比较，T检验是最常用的方法之一。

本次推文主要内容：

一、什么是T检验

二、T检验如何与群落微生物的alpha多样性结合

三、如何使用R语言进行T检验以及绘制箱线图

一、什么是T检验

t检验，又称student t检验（Student's t test），主要用于样本含量较小（例如n < 30），总体标准差σ未知的正态分布。用于统计量服从正态分布，但方差未知的情况，从而比较两个平均数的差异是否显著。

前提条件：要求样本服从正态分布或近似正态分布，不然需要利用一些数据转化（如取对数、开根号等）试图将其转化为服从正态分布的数据。如果还是不满足正态分布，只能利用非参数检验方法。

这里就简单了解一下

二、T检验如何与群落微生物的alpha多样性结合

Alpha多样性是指一个特定区域或生态系统内的多样性，是反映丰富度和均匀度的综合指标。Alpha多样性主要与两个因素有关：一是种类数目，即丰富度；二是多样性，群落中个体分配上的均匀性。

群落丰富度（Community richness）指数主要包括Chao指数和Ace指数。
 Chao和Ace越大，说明群落中含有的OTU数目越多，群落的丰富度越大。
群落多样性（Community diversity）指数，包括Shannon指数和Simpson指数。
Shannon值越大，说明群落多样性越高，但是Simpson 指数值较大，这说明群落多样性较低

我们可以T检验来比较两个微生物群落中alpha多样性指数，来反映两个群落在物种丰富度和均匀度方面的差异，来推断两个微生物群落的差异大小。

三、如何使用R语言进行T检验以及绘制箱线图

这里准备了一组示例文件

“alpha.txt”为16S细菌群落测序所获得的部分alpha多样性指数数据，包含3列信息：Sample为样本名称；Chao指数、Ace指数、Shannon指数和Simpson指数等四个alpha多样性指数。

“group.txt”为样本分组信息，Sample为样本名称；group为样本的分组信息。

alpha.txt文件

alpha.txt

group.txt文件

group.txt

接下来，我们使用R语言运行T检验，来观察样本之间的各alpha多样性指数是否存在显著不同

1.样本数据的预处理

首先，我们在excel中查看一下

补充一点知识：

reshape2R包主要有两个主要的功能：melt和cast

melt：将wide-format数据“熔化”成long-format数据；

cast：获取long-format数据“重铸”成wide-format数据。

#reshape2包可以实现在宽格式（wide-format）和长格式（long-format）之间转换数据。
library(reshape2)
 
#读入文件，合并分组信息，数据重排
#stringsAsFactors = FALSE表示在读取数据时，遇到字符串之后，不将其转化为factors，仍然保留为字符串
alpha <- read.table('alpha.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE)
group  <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE)
alpha1 <- melt(merge(alpha, group, by = 'sample'), id = c('sample', 'group'))

alpha.txt文件

alpha

group.txt文件

group

转化之后的alpha文件

#这里我们查看 group 1 和 group 2 的 alpha多样性指数是否存在显著差异
#选择出需要比较的分组1和2
richness_ACE <- subset(alpha, variable == 'ACE' & group %in% c('1', '2'))
richness_ACE$group <- factor(richness_ACE$group)

richness_Chao1 <- subset(alpha, variable == 'Chao1' & group %in% c('1', '2'))
richness_Chao1$group <- factor(richness_Chao1$group)

richness_Simpson <- subset(alpha, variable == 'Simpson' & group %in% c('1', '2'))
richness_Simpson$group <- factor(richness_Simpson$group)

richness_Shannon <- subset(alpha, variable == 'Shannon' & group %in% c('1', '2'))
richness_Shannon$group <- factor(richness_Shannon$group)
#用view查看数据集，也可以用head
View(richness_ACE)
View(richness_Chao1)
View(richness_Shannon)
View(richness_Simpson)

多样性指数

可以在上图中看到四列：

Sample：样本名称；

group：两组的分组名称，且分组列已转化为因子类型

variable：alpha多样性指数ACE

value:ACE指数的数值。

2.正态性假设检验

前面提到过， T检验的前提是数据必须符合正态分布模型。

因此在执行t检验之前我们需要去验证数据分布是否具有正态性。若数据不符合正态分布，则t检验将无法适用于该数据（此时可以考虑转化数据，或者使用非参数的检验方法）。

验证数据是否符合正态分布的方法很多，比如直方图、QQ图以及Shapiro-Wilk检验等，以下展示 QQ图和Shapiro-Wilk检验。

Shapiro-Wilk检验——球性检验

该方法类似于使用线性回归的方法，是检验其回归曲线的残差，来验证数据的分布是否正态性分布。R语言中可用函数shapiro.test()执行Shapiro-Wilk检验,，若结果p值大于0.05，则数据分布符合正态性。

##Shapiro-Wilk 检验，当且仅当  两者 p 值都大于 0.05 时表明数据符合正态分布
shapiro_ACE <- tapply(richness_ACE$value, richness_ACE$group, shapiro.test)
shapiro_ACE
shapiro_Chao1 <- tapply(richness_Chao1$value, richness_Chao1$group, shapiro.test)
shapiro_Chao1
shapiro_Simpson <- tapply(richness_Simpson$value, richness_Simpson$group, shapiro.test)
shapiro_Simpson
shapiro_Shannon <- tapply(richness_Shannon$value, richness_Shannon$group, shapiro.test)
shapiro_Shannon
shapiro_ACE$'1'$p.value
shapiro_ACE$'2'$p.value
shapiro_Chao1$'1'$p.value
shapiro_Chao1$'2'$p.value
shapiro_Simpson$'1'$p.value
shapiro_Simpson$'2'$p.value
shapiro_Shannon$'1'$p.value
shapiro_Shannon$'2'$p.value

Shapiro-Wilk检验

从上图可以看到，四种alpha多样性指数都大于0.05，所以这些数据分布符合正态性，可以进行下一步的T检验

QQ图检验

使用car包中的qqPlot()，绘制QQ图查看数值分布。在图中横坐标是标准的正态分布值，纵坐标是我们自己数据的值。如果所有的点都离直线很近，落在置信区间内（图中虚线部分，默认展示95%置信区间），表明数据符合正态性。

##qq 图验证数据的正态性
library(car) #需要用到car包的qqplot
par(mfrow=c(2,2)) #设置2乘2的格式

qqPlot(lm(value~group, data = richness_ACE), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
qqPlot(lm(value~group, data = richness_Chao1), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
qqPlot(lm(value~group, data = richness_Simpson ), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
qqPlot(lm(value~group, data = richness_Shannon), simulate = TRUE, main = 'QQ Plot', labels = FALSE)

QQ图

从上图可以看出，所有数据都在置信区间内，表明数据符合正态性。

T检验

经过上面正态性检验表明我们的数据分布通过了正态假设检验，可以执行t检验。

独立样本的t检验

我们的数据样本是不同物种的群落微生物，所以样本间是相互独立的，所以选用独立样本t检验。

R语言t检验的函数是t.test()

该函数默认两组数据间相互独立（默认参数paired = FALSE），执行独立样本的t检验。

在R语言中的t检验默认假定方差不相等（默认参数var.equal = FALSE），并使用Welsh的修正自由度；可以通过参数“var.equal= TRUE”假定方差相等，并使用合并方差估计。

t.test()默认的备择假设是双侧的（默认参数alternative = 'two.sided'），即执行双侧检验；可分别使用“alternative= 'less'”或“alternative = 'greater'”执行单侧t检验。

我们可以通过doBy包中summaryBy函数来分别计算各组中的均值以及标准差

#通过doBy包中summaryBy函数来分别计算各组中的均值以及标准差
library(doBy)   
ACE1 <- summaryBy(value~group, richness_ACE, FUN = c(mean, sd))
Chao11 <- summaryBy(value~group, richness_Chao1, FUN = c(mean, sd))
Simpson1 <- summaryBy(value~group, richness_Simpson, FUN = c(mean, sd))
Shannon1 <- summaryBy(value~group, richness_Shannon, FUN = c(mean, sd))

标准差

由于各组方差不相等，所以我们执行假设方差不相等的双侧检验

##独立样本的 t 检验
t_test_ACE <- t.test(value~group, richness_ACE, paired = FALSE, alternative = 'two.sided')
t_test_ACE
t_test_Chao1 <- t.test(value~group, richness_Chao1, paired = FALSE, alternative = 'two.sided')
t_test_Chao1
t_test_Simpson <- t.test(value~group, richness_Simpson, paired = FALSE, alternative = 'two.sided')
t_test_Simpson
t_test_Shannon <- t.test(value~group, richness_Shannon, paired = FALSE, alternative = 'two.sided')
t_test_Shannon

t_test_ACE$p.value
t_test_Chao1$p.value
t_test_Simpson$p.value
t_test_Shannon$p.value

独立样本T检验

从图上我们可以看到独立样本的 t 检验的结果，

ACE的p值小于0.05，group1和group2的ACE指数间存在显著不同。

chao1的p值小于0.05，group1和group2的chao1指数间存在极显著不同。

Simpson的p值大于0.05，group1和group2的Simpson指数间没有差异。

Shannon的p值大于0.05，group1和group2的Shannon指数间没有差异。

可视化展示

可以使用箱线图和小提琴图来将两组差异进行可视化展示，这样使用箱线图

#boxplot() 箱线图示例
par(mfrow=c(2,2))
boxplot(value~group, data = richness_ACE, col = c('blue', 'orange'), 
        ylab = 'ACE', xlab = 'Group', main = 't-test: p-value < 0.05')
boxplot(value~group, data = richness_Chao1, col = c('blue', 'orange'), 
        ylab = 'Chao1', xlab = 'Group', main = 't-test: p-value < 0.01')
boxplot(value~group, data = richness_Simpson, col = c('blue', 'orange'), 
        ylab = 'Simpson', xlab = 'Group', main = 't-test: p-value > 0.05')
boxplot(value~group, data =richness_Shannon, col = c('blue', 'orange'), 
        ylab = 'Shannon', xlab = 'Group', main = 't-test: p-value >0.05')

箱线图-boxplot

推荐使用ggplot2 来绘制

library(patchwork)#加载拼图包
library(ggplot2)    #ggplot2 作图

p1<-ggplot(richness_ACE, aes(x = group, y = value, fill = group)) + 
  geom_boxplot(outlier.size = 0.7) +
  labs(x = 'Group', y = 'ACE', title = 't-test: p-value < 0.05')+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) 


p2<-ggplot(richness_Chao1, aes(x = group, y = value, fill = group)) + 
  geom_boxplot(outlier.size = 0.7) +
  labs(x = 'Group', y = 'Chao1', title = 't-test: p-value < 0.01')+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) 

p3<-ggplot(richness_Simpson, aes(x = group, y = value, fill = group)) + 
  geom_boxplot(outlier.size = 0.7) +
  labs(x = 'Group', y = 'Simpson', title = 't-test: p-value > 0.05')+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) 

p4<-ggplot(richness_Shannon, aes(x = group, y = value, fill = group)) + 
  geom_boxplot(outlier.size = 0.7) +
  labs(x = 'Group', y = 'Shannon', title = 't-test: p-value > 0.05')+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) 

p1+p2+p3+p4

箱线图-ggplot2

当然我们也可以绘制柱形图来表示差异

#先分别计算各组中的均值以及标准差，展示为均值 ± 标准差
library(doBy)   #使用其中的 summaryBy() 以方便按分组计算
library(ggplot2)    #ggplot2 作图

ACE1 <- summaryBy(value~group, richness_ACE, FUN = c(mean, sd))
Chao11 <- summaryBy(value~group, richness_Chao1, FUN = c(mean, sd))
Simpson1 <- summaryBy(value~group, richness_Simpson, FUN = c(mean, sd))
Shannon1 <- summaryBy(value~group, richness_Shannon, FUN = c(mean, sd))

library(patchwork) #加载拼图包
p1<-ggplot(ACE1, aes(group, value.mean, fill = group))+ 
  geom_col(width = 0.4, show.legend = FALSE)+ 
  geom_errorbar(aes(ymin = value.mean - value.sd, ymax = value.mean + value.sd), width = 0.15, size = 0.5)+ 
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) +
  labs(x = 'Group', y = 'ACE', title = 't-test: p-value < 0.05')+
  scale_y_continuous(expand = c(0,0),limits = c(0,2000))

p2<-ggplot(Chao11, aes(group, value.mean, fill = group))+ 
  geom_col(width = 0.4, show.legend = FALSE)+ 
  geom_errorbar(aes(ymin = value.mean - value.sd, ymax = value.mean + value.sd), width = 0.15, size = 0.5)+ 
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) +
  labs(x = 'Group', y = 'Chao1', title = 't-test: p-value < 0.01')+
  scale_y_continuous(expand = c(0,0),limits = c(0,1800))

p3<-ggplot(Simpson1, aes(group, value.mean, fill = group))+ 
  geom_col(width = 0.4, show.legend = FALSE)+ 
  geom_errorbar(aes(ymin = value.mean - value.sd, ymax = value.mean + value.sd), width = 0.15, size = 0.5)+ 
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) +
  labs(x = 'Group', y = 'Simpson', title = 't-test: p-value > 0.05')+
  scale_y_continuous(expand = c(0,0),limits = c(0,1.2))

p4<-ggplot(Shannon1, aes(group, value.mean, fill = group))+ 
  geom_col(width = 0.4, show.legend = FALSE)+ 
  geom_errorbar(aes(ymin = value.mean - value.sd, ymax = value.mean + value.sd), width = 0.15, size = 0.5)+ 
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        plot.title = element_text(hjust = 0.5)) +
  labs(x = 'Group', y = 'Shannon', title = 't-test: p-value > 0.05')+
  scale_y_continuous(expand = c(0,0),limits = c(0,9))

p1+p2+p3+p4 #拼图

柱形图

好了，就两个微生物群落间的比较，T检验的有关知识就到这里结束了，欢迎期待下一期！！！

谢谢你的阅读
如有不足之处，请联系小蘑菇！
欢迎关注微信公众号：fafu 生信小蘑菇

组间分析—T检验、R语言绘图
T检验在做微生物群落分析时经常要寻找不同群落之间的差异，我们会经常关心两个或者多个群落之间是否显著不同。就两个微...
R语言ggstatsplot包做T检验
R语言用ggstatsplot包做方差分析和绘图R语言ggstatsplot包做卡方检验单样本均值比较 1、点图...
R语言T检验
我们比较的对象是南方和非南方各州，因变量为监禁的概率。一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的...
R语言t检验
本文首发于公众号：医学和生信笔记医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R...
【R画图学习21.5】单因素方差分析评估条件
对于两组数据间的差异分析，最常见的方法就是使用T检验。当拓展到多组（三组及以上）时，使用T检验逐一两两比较的方法无...
生信入门学习笔记day4@2021.06.28
R语言基础 R 语言是一种主要用于统计分析、绘图、数据挖掘的数学编程语言。R language: The R Pr...
R语言-配对t检验
准备data demo 方法一：把配对的数据分别保存在两个数值型向量里（确保数值的配对位置）方法二：把配对的数据...
看看不同分析工具的单样本t检验（R、SPSS、Python、Ex
本章来介绍不同分析工具对于单样本t检验的操作方式和参数解读：一、看看用R语言怎么做检验？前提在题目中已经标出，数...
2021-06-18 两组间差异的非参数检验之Wilcox秩和检
在进行两组数据间的差异分析时，我们通常会想到使用 t检验 [#rd]。但若数据不满足执行t检验的参数假设（例如数据...
宏基因组、扩增子组间差异检验，终于有人讲清楚了！！！
什么是组间差异检验？就是组间的差异分析以及显著性检验，应用统计学上的假设检验方法，检验组间是否有差异及其差异程度。...