微生物群落的改变-多个群落到一个群落距离计算展示

背景讲解-水稻时间序列

参考刘老师推文：https://mp.weixin.qq.com/s/WUtV6XnzA4a6PNiROjrWEQ

图2. 田间水稻根系微生物组在8-10周趋于稳定。

A-D. 对两个水稻品种分别在两地进行的连续微生物组调查结果相关分析，发现8-10周后群落结构趋于稳定。

E. 所有时间点距离水稻最后取样点的Bray-Curtis距离。发现土壤呈现小幅波动变化，而根系呈现出先快后慢，逐渐趋近的变化规律。

F. 不同水稻品种在两个地点间的距离变化，发现土壤差异稳定，而根系微生物组差异随时间增长而趋于一致。

E-F基于vegan包计算的所有样品两两间Bray-Curtis距离。分别挑选距离终点的距离，和两地间的距离与时间序列上的关系，并采用ggplot2可视化散点图，并添加拟合曲线方便观察变化规律。

这些R包都是常见的包，其中ggClusterNet包提供数据，EasyStat用于距离的差异分析。这个个包均在github上，大家使用install_github安装，如果网络问题安装不上，请在微信搜一搜中搜索：玩转R包，即可看到我写的一篇R包安装指南，尝试过后相信大部分R包都可以安装上。

我对这个部分代码进行了整理

  • 并做了简单的封装，使用起来更加容易，只需要phyloseq对象即可

  • 其次对距离之间差异做了统计分析。

  • 如果有多个处理或者多种情况，可以进行多组同时分析，同时出图，方便比对

可以

#--读入文件
library(phyloseq)
library(vegan)
library(tidyverse)
require(splines)
require(MASS)

library(ggClusterNet)
library(EasyStat)

data(ps)

  • tre1 为指定的一组微生物群落，即选择这组微生物群落，计算这组样本的距离。

  • refer 指定一个参考组，也就是其他全部组均以这个组为参考，计算相对距离。

    这里只能指定一个refer，因为只能有一个参考。

注意：phyloseq对象中map文件分组列名字必须为Group。

map = as.data.frame(sample_data(ps))
head(map)
map$Group
map$Site
tre1 = c("WT","KO","OE")
datab = data.frame(tre1 = tre1)

refer = "WT"

下面一组代码用于计算和出图

plot = list()
for (i in 1:ncol(datab)) {
tre = as.character(datab[[i]])
map = as.data.frame(sample_data(ps))
map$ID = row.names(map)
maps <- dplyr::filter(as.tibble(map),Group%in% tre)
row.names(maps) = maps$ID
ps_sub = ps
sample_data(ps_sub) = as.data.frame(maps);ps_sub

ps_rela = transform_sample_counts(ps_sub , function(x) x / sum(x) );ps_rela

# 提取otu表格，筛选数据
vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
otu= as.data.frame(t(vegan_otu(ps_rela)))

head(otu)
#-添加第一个时间的均值作为标准，求取均值列加入其后
mapsub = as.data.frame(sample_data(ps_rela))
#挑选time为1的样本
ID = mapsub$ID[mapsub$Group == refer]
#挑选样本
head(otu)
ID = as.character(ID)
index = otu[,ID] # 筛选并计算均值
head(index,n = 5)

index_mean= rowMeans(index)
otu$index_mean=index_mean # 均值添加到OTU表
head(otu)
# 计算包括终点均值的所有样品bray距离
bray_curtis = vegan::vegdist(t(otu), method = "bray")
bray_curtis= as.matrix(bray_curtis)

#--------------整理数据-
# 计算各组内部差异程度
# 建立一个存储组内差异的数据框
dat = t(as.data.frame(c("sampleA","sampleB","0.15","group","genosite")))
colnames(dat) = c("sampleA","sampleB","distance","group","type")
rownames(dat) = c("test")

# 每个样品与final对应的距离
# ii = 1
for (ii in tre){
group = rownames(mapsub[mapsub$Group %in% ii,])
# m = 1
for (m in 1:(length(group))) {
x = c(group[m],"index_mean",bray_curtis[group[m],"index_mean"],ii,tre)
dat=rbind(dat,x)
}
}
dat = as.data.frame(dat[-1,], stringsAsFactors=F) # 删除首行框架数据

# dat = dat[dat$distance != 0,] #
# 距离转换为3位数值
dat$distance=round(as.numeric(as.character(dat$distance)), digits=3)
# 分组添加levels，分组有时间顺序
dat$group = factor(dat$group, levels=unique(dat$group))

all = dat
# 开始尝试散点图拟合
all1 = all
all1$group = as.numeric(all1$group)
iris_groups<- group_by(all, group)
all3<- dplyr::summarise(iris_groups, mean(distance))
colnames(all3) = c("group","mean")
datawt = data.frame(ID =all$sampleA,group = all$group,distance = all$distance)

result= aovMcomper (data = datawt, i= 3,method_Mc = "Tukey")
# 提取多重比较结果
result[[1]]
PlotresultBox = aovMuiBoxP(data = datawt, i= 3,sig_show ="abc",result = result[[1]])
#提取图片
p = PlotresultBox[[1]]
p = p + geom_line(aes(x=group, y=mean,group = 1),data = all3) +geom_point(aes(x=group, y=mean),all3,alpha=1, size=4,pch = 21,fill="red")
p <- p + theme_bw() + theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)

plot[[i]] = p

}

如果有多个图，这里使用ggarrange组图，只有一个就无所谓了。

library(ggpubr)
p = ggarrange(plotlist = plot, common.legend = TRUE, legend="right", ncol = 1)
p

这里我对Group列多弄了几个分组，来演示一下多个分析和出图的用法。

map = as.data.frame(sample_data(ps))
head(map)
map$Group = paste(map$Group,c(1,2),sep = "")
sample_data(ps)= map

tre1 = c("WT1","KO1","OE1")
tre2 = c("WT1","KO2","OE2")

datab = data.frame(tre1 = tre1,tre2 = tre2)

refer = "WT1"

下面一组代码用于计算和出图

plot = list()
for (i in 1:ncol(datab)) {
tre = as.character(datab[[i]])
map = as.data.frame(sample_data(ps))
map$ID = row.names(map)
maps <- dplyr::filter(as.tibble(map),Group%in% tre)
row.names(maps) = maps$ID
ps_sub = ps
sample_data(ps_sub) = as.data.frame(maps);ps_sub

ps_rela = transform_sample_counts(ps_sub , function(x) x / sum(x) );ps_rela

# 提取otu表格，筛选数据
vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
otu= as.data.frame(t(vegan_otu(ps_rela)))

head(otu)
#-添加第一个时间的均值作为标准，求取均值列加入其后
mapsub = as.data.frame(sample_data(ps_rela))
#挑选time为1的样本
ID = mapsub$ID[mapsub$Group == refer]
#挑选样本
head(otu)
ID = as.character(ID)
index = otu[,ID] # 筛选并计算均值
head(index,n = 5)

index_mean= rowMeans(index)
otu$index_mean=index_mean # 均值添加到OTU表
head(otu)
# 计算包括终点均值的所有样品bray距离
bray_curtis = vegan::vegdist(t(otu), method = "bray")
bray_curtis= as.matrix(bray_curtis)

#--------------整理数据-
# 计算各组内部差异程度
# 建立一个存储组内差异的数据框
dat = t(as.data.frame(c("sampleA","sampleB","0.15","group","genosite")))
colnames(dat) = c("sampleA","sampleB","distance","group","type")
rownames(dat) = c("test")

# 每个样品与final对应的距离
# ii = 1
for (ii in tre){
group = rownames(mapsub[mapsub$Group %in% ii,])
# m = 1
for (m in 1:(length(group))) {
x = c(group[m],"index_mean",bray_curtis[group[m],"index_mean"],ii,tre)
dat=rbind(dat,x)
}
}
dat = as.data.frame(dat[-1,], stringsAsFactors=F) # 删除首行框架数据

# dat = dat[dat$distance != 0,] #
# 距离转换为3位数值
dat$distance=round(as.numeric(as.character(dat$distance)), digits=3)
# 分组添加levels，分组有时间顺序
dat$group = factor(dat$group, levels=unique(dat$group))

all = dat
# 开始尝试散点图拟合
all1 = all
all1$group = as.numeric(all1$group)
iris_groups<- group_by(all, group)
all3<- dplyr::summarise(iris_groups, mean(distance))
colnames(all3) = c("group","mean")
datawt = data.frame(ID =all$sampleA,group = all$group,distance = all$distance)

result= aovMcomper (data = datawt, i= 3,method_Mc = "Tukey")
# 提取多重比较结果
result[[1]]
PlotresultBox = aovMuiBoxP(data = datawt, i= 3,sig_show ="abc",result = result[[1]])
#提取图片
p = PlotresultBox[[1]]
p = p + geom_line(aes(x=group, y=mean,group = 1),data = all3) +geom_point(aes(x=group, y=mean),all3,alpha=1, size=4,pch = 21,fill="#FF7F00")
p <- p + theme_bw() + theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)

plot[[i]] = p

}

如果有多个图，这里使用ggarrange组图，只有一个就无所谓了。

library(ggpubr)
p = ggarrange(plotlist = plot, common.legend = TRUE, legend="right", ncol = 2)
p

加主编微信 加入群聊

关于微生信生物 你想要的都在这里

微生信生物

赞 (0)

相关推荐