Loading...
墨滴

fafu生信小蘑菇

2021/04/30  阅读:17  主题:橙心

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

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

在做微生物群落分析时经常要寻找不同群落之间的差异,我们会经常关心两个或者多个群落之间是否显著不同。就两个微生物群落间的比较,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
alpha.txt

group.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
alpha

group.txt文件

group
group

转化之后的alpha文件

转化之后的alpha文件
转化之后的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检验
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图
QQ图

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

T检验

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

独立样本的t检验

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

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

  2. 该函数默认两组数据间相互独立(默认参数paired = FALSE),执行独立样本的t检验。

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

  4. **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检验

从图上我们可以看到独立样本的 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
箱线图-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
箱线图-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 生信 小蘑菇

fafu生信小蘑菇

2021/04/30  阅读:17  主题:橙心

作者介绍

fafu生信小蘑菇