#读取 OTU 丰度表
otu <- read.delim('otu_table.txt', row.names = 1, stringsAsFactors = FALSE)
otu <- data.frame(t(otu))
#读取样本分组
group <- read.delim('group.txt', stringsAsFactors = FALSE)
###############观察和比较药物处理前后两组的菌群结构是否存在明显的不同
##使用 vegan 包的置换多元方差分析（PERMANOVA），比较用药前后两组中肠道菌群结构是否存在不同
library(vegan)
## Warning: package 'vegan' was built under R version 4.1.0
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.0
## Loading required package: lattice
## This is vegan 2.5-7
#执行 PERMANOVA，详情 ?adonis
#样本间的成对距离选择使用物种多样性中最常使用的 Bray-curtis 距离，并通过 999 次置换估计 p 值
adonis_drug <- adonis(otu~drug, group, distance = 'bray', permutations = 999)
adonis_drug
##
## Call:
## adonis(formula = otu ~ drug, data = group, permutations = 999, distance = "bray")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## drug 1 0.03599 0.035994 0.88132 0.03852 0.693
## Residuals 22 0.89851 0.040841 0.96148
## Total 23 0.93450 1.00000
##使用 vegan 包执行主坐标分析（PCoA），观察前后两组中肠道菌群结构的差异
#library(vegan)
#为了与上文 PERMANOVA 对应起来，这里同样计算了样本间成对的 Bray-curtis 距离，详情 ?vegdist
bray_dis <- vegdist(otu, method = 'bray')
#PCoA 排序，详情 ?cmdscale
pcoa <- cmdscale(bray_dis, k = (nrow(otu) - 1), eig = TRUE)
#提取前两轴的贡献度
pcoa_exp <- pcoa$eig/sum(pcoa$eig)
pcoa1 <- paste('PCoA axis1 :', round(100*pcoa_exp[1], 2), '%')
pcoa2 <- paste('PCoA axis2 :', round(100*pcoa_exp[2], 2), '%')
#提取前两轴的的坐标，并添加样本的分组信息
site <- data.frame(pcoa$point)[1:2]
site$sample <- rownames(site)
site <- merge(site, group, by = 'sample')
names(site)[2:3] <- c('pcoa1', 'pcoa2')
#ggplot2 绘制二维平面图展示 PCoA 结果
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.0
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 4.1.0
p <- ggplot(data = site) +
geom_point(aes(x = pcoa1, y = pcoa2, color = drug), size = 2) + #绘制两轴的样本点
geom_text_repel(aes(x = pcoa1, y = pcoa2, label = sample, color = drug), size = 2.5,
box.padding = unit(0.3, 'lines'), show.legend = FALSE) + #添加样本标签
scale_color_manual(limits = c('Before', 'After'), values = c('#D27FB2', '#764697')) + #两组的颜色
theme(panel.grid = element_blank(), panel.background = element_blank(),
axis.line = element_line(color = 'black'), legend.key = element_blank()) +
labs(x = pcoa1, y = pcoa2, color = '')
p
#在图中添加两组的质心（分别计算两组样本在 PCoA 两轴坐标的均值），并以中心大点显示在图中
group_average <- aggregate(cbind(pcoa1, pcoa2)~drug, data = site, FUN = mean)
p1 <- p +
stat_ellipse(aes(x = pcoa1, y = pcoa2, color = drug), level = 0.95, linetype = 2, show.legend = FALSE) + #绘制两组的 95% 置信椭圆
geom_point(data = group_average, aes(x = pcoa1, y = pcoa2, color = drug), size = 5, show.legend = FALSE) #绘制中央质心大点
#在图中添加上文的 PERMANOVA 的结果（P 值）
p1 <- p1 +
annotate('text', label = 'PERMANOVA', x = 0.18, y = 0.15, size = 3) +
annotate('text', label = sprintf('italic(P) == %.3f', adonis_drug$aov.tab[1,6]), x = 0.18, y = 0.13, size = 3, parse = TRUE)
p1
###############比较受试者个体水平的肠道菌群结构差异
#执行 PERMANOVA 比较受试者个体水平的肠道菌群结构差异
#样本间的成对距离选择使用物种多样性中最常使用的 Bray-curtis 距离，并通过 999 次置换估计 p 值
adonis_individual <- adonis(otu~individual, group, distance = 'bray', permutations = 999)
adonis_individual
##
## Call:
## adonis(formula = otu ~ individual, data = group, permutations = 999, distance = "bray")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## individual 1 0.11899 0.118987 3.2099 0.12733 0.001 ***
## Residuals 22 0.81551 0.037069 0.87267
## Total 23 0.93450 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#继续在上述 PCoA 结果图中以虚线连接来自相同个体的配对样本
site_before <- subset(site, drug == 'Before')[c('individual', 'pcoa1', 'pcoa2')]
site_after <- subset(site, drug == 'After')[c('individual', 'pcoa1', 'pcoa2')]
site_ab <- merge(site_before, site_after, by = 'individual')
p2 <- p +
geom_segment(data = site_ab, aes(x = pcoa1.x, y = pcoa2.x, xend = pcoa1.y, yend = pcoa2.y),
color = 'gray', size = 0.3, linetype = 2)
#在图中添加上文的 PERMANOVA 结果（P 值）
p2 <- p2 +
annotate('text', label = 'PERMANOVA', x = 0.13, y = 0.15, size = 3) +
annotate('text', label = sprintf('italic(P) == %.3f', adonis_individual$aov.tab[1,6]), x = 0.13, y = 0.13, size = 3, parse = TRUE)
p2
##备注：
#刚才是先画的散点，后画的虚线，故虚线图层遮挡在散点的上方
#若想让虚线位于散点的下方，作图时先画虚线，后画散点即可，如下所示
ggplot(data = site) +
geom_segment(data = site_ab, aes(x = pcoa1.x, y = pcoa2.x, xend = pcoa1.y, yend = pcoa2.y),
color = 'gray', size = 0.3, linetype = 2) + #以虚线连接配对样本
geom_point(aes(x = pcoa1, y = pcoa2, color = drug), size = 2) + #绘制两轴的样本点
geom_text_repel(aes(x = pcoa1, y = pcoa2, label = sample, color = drug), size = 2.5,
box.padding = unit(0.3, 'lines'), show.legend = FALSE) + #添加样本标签
scale_color_manual(limits = c('Before', 'After'), values = c('#D27FB2', '#764697')) + #两组的颜色
theme(panel.grid = element_blank(), panel.background = element_blank(),
axis.line = element_line(color = 'black'), legend.key = element_blank()) +
labs(x = pcoa1, y = pcoa2, color = '') +
annotate('text', label = 'PERMANOVA', x = 0.13, y = 0.15, size = 3) + #添加 PERMANOVA 的 P 值信息
annotate('text', label = sprintf('italic(P) == %.3f', adonis_individual$aov.tab[1,6]), x = 0.13, y = 0.13, size = 3, parse = TRUE)
#ref https://github.com/lyao222lll/sheng-xin-xiao-bai-yu-2021/tree/main/2021