R语言绘制云雨图

最近学习了用云雨图展示菌群随时间的动态变化，并判断其变化与对照组的差别。

代码来源：

Samara, J., et al. (2022). "Supplementation with a probiotic mixture accelerates gut microbiome maturation and reduces intestinal inflammation in extremely preterm infants." Cell Host Microbe 30(5): 696-711 e695.
    Probiotics are increasingly administered to premature infants to prevent necrotizing enterocolitis and neonatal sepsis. However, their effects on gut microbiome assembly and immunity are poorly understood. Using a randomized intervention trial in extremely premature infants, we tested the effects of a probiotic product containing four strains of Bifidobacterium species autochthonous to the infant gut and one Lacticaseibacillus strain on the compositional and functional trajectory of microbiome. Daily administration of the mixture accelerated the transition into a mature, term-like microbiome with higher stability and species interconnectivity. Besides infant age, Bifidobacterium strains and stool metabolites were the best predictors of microbiome maturation, and structural equation modeling confirmed probiotics as a major determinant for the trajectory of microbiome assembly. Bifidobacterium-driven microbiome maturation was also linked to an anti-inflammatory intestinal immune milieu. This demonstrates that Bifidobacterium strains are ecosystem engineers that lead to an acceleration of microbiome maturation and immunological consequences in extremely premature infants.

输入数据：

metadata.csv

otu_table.csv

代码如下：

library(vegan)
library(ape)
library(plyr)
library(ggridges)
d = read.csv("genus.csv", row.names = 1)#otu_table格式
meta = read.csv("metadata.csv", row.names = 1)#表型数据
vegan_otu = t(d) %>% as.data.frame()#taxa in col
OTU.Bray <- as.matrix(vegdist(vegan_otu, "bray"))
pcoa <- pcoa(OTU.Bray)
pcoa1 <- pcoa$vectors[,1]
pcoa2 <- pcoa$vectors[,2]

meta$pcoa1 = pcoa1
meta$pcoa2 = pcoa2

mu <- ddply(meta, "TP", summarise, grp.mean=mean(pcoa1))#求pcoa1的均值

#排序
sdf <- meta %>%
  arrange(desc(TP))
#绘图
fig <- ggplot(sdf, aes(x = pcoa1, y = factor(TP, levels = c(0,4,3,2,1), labels = c("Control", "T4","T3","T2","T1")), fill = TP)) +
  geom_density_ridges2(rel_min_height=0.01) +
  scale_x_continuous(expand = c(0.01, 0)) +
  scale_y_discrete(expand = c(0.01, 0)) +
  #T1
  annotate("segment", x = mu[2,2], xend = mu[1,2], y = 5.4, yend = 5.4, colour = "black", size=0.5)+ 
  annotate("pointrange", x = mu[2,2], y = 5.4, ymin = 5.4, ymax = 5.4,colour = "red", size = 1)+ 
  #T2
  annotate("segment", x = mu[3,2], xend = mu[1,2], y = 4.6, yend = 4.6, colour = "black", size=0.5)+ 
  annotate("pointrange", x = mu[3,2], y = 4.6, ymin = 4.6, ymax = 4.6,colour = "red", size = 1)+
  #T3
  annotate("segment", x = mu[4,2], xend = mu[1,2], y = 3.5, yend = 3.5, colour = "black", size=0.5)+ 
  annotate("pointrange", x = mu[4,2], y = 3.5, ymin = 3.5, ymax = 3.5,colour = "red", size = 1)+
  #T4
  annotate("segment", x = mu[5,2], xend = mu[1,2], y = 2.5, yend = 2.5, colour = "black", size=0.5)+ 
  annotate("pointrange", x = mu[5,2], y = 2.5, ymin = 2.5, ymax = 2.5,colour = "red", size = 1)+
  #term 6 month
  annotate("segment", x = mu[1,2], xend = mu[1,2], y = 1.5, yend = 7.3, colour = "black", size=0.5)+
  annotate("pointrange", x = mu[1,2], y = 1.5, ymin = 1.5, ymax = 1.5,colour = "red", size = 1)+
  #scale_fill_manual(values = c("mediumpurple3", "plum2", "steelblue4", "steelblue3", "skyblue2"))+
  #scale_color_manual(values = c("mediumpurple3", "plum2", "steelblue4", "steelblue3", "skyblue2"))+
  xlab("Bray-Curtis Dissimilarity (PCoA1)")+
  theme_ridges(font_size = 10, grid = TRUE, center_axis_labels = TRUE) + theme(axis.title.y = element_blank())+
  theme(legend.position = "none")
fig

输出结果：

红点就是pcoa1在组内的均值，可以看到其与对照组的距离，距离越近说明与对照越相似。这个图在菌群的动态分析中挺实用的。