组间分析—T检验、R语言绘图
2021-04-30 本文已影响0人
fafu生信小蘑菇
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文件
转化之后的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 生信 小蘑菇
