论文信息
- 论文题目: Deterministic colonization arises early during the transition of soil bacteria to the phyllosphere and is shaped by plant-microbe interactions
- DOI: 10.1186/s40168-025-02090-1
- 复现目标: Figure 4 微生物富集结果
数据来源
- 原始测序数据: NCBI SRA 数据库 BioProject: PRJNA1174406
- 分析脚本和元数据: Figshare 平台公开数据包
数据加载与预处理
# 加载必要的包
library(vegan)
library(ggplot2)
library(RColorBrewer)
library(tidyr)
library(dplyr)
library(ggpubr)
# 加载数据
MeganData <- read.table(file="./ForBarchart.txt", header=T, sep="\t")
# 数据预处理
PlotLevel <- "Order"
MeganDataPlot <- MeganData %>% count(Sample, Order)
MeganDataPlot$ra <- MeganDataPlot$n / 30
MeganDataPlot$ra[which(MeganDataPlot$Sample == "PB2PB")] <- MeganDataPlot$n[which(MeganDataPlot$Sample == "PB2PB")] / 31
# 设置颜色
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(10)
names(mycolors) <- unique(MeganDataPlot$Order)
Taxscale <- scale_fill_manual(name = "Order",values = mycolors)
# 设置因子水平
MeganDataPlot$Sample <- factor(MeganDataPlot$Sample, levels=c("NG2NG", "NG2PB", "NG2Col", "PB2PB", "PB2NG", "PB2Col"))
print("数据加载成功")
Sample Order n ra
1 NG2Col Burkholderiales 1 0.03333333
2 NG2Col Corynebacteriales 4 0.13333333
3 NG2Col Hyphomicrobiales 3 0.10000000
4 NG2Col Micrococcales 8 0.26666667
5 NG2Col Pseudomonadales 1 0.03333333
6 NG2Col Sphingomonadales 11 0.36666667
CrossInoc 柱状图
# 绘制柱状图
ggplot(data=MeganDataPlot) + geom_bar(aes(x=Sample, y=ra, fill=Order), stat="identity") +
theme_classic(base_size=18) + ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust=1)) +
Taxscale
CrossInoc 气泡图
# 气泡图数据准备
MeganDataPlot <- MeganData %>% count(Sample, Genus)
MeganDataPlot$Sample <- factor(MeganDataPlot$Sample, levels=c("NG2NG", "NG2PB", "NG2Col", "PB2PB", "PB2NG", "PB2Col"))
MeganDataPlot$Origin <- substr(MeganDataPlot$Sample, start=1, stop = 2)
MeganDataPlot_NG <- MeganDataPlot[which(MeganDataPlot$Origin == "NG"),]
MeganDataPlot_PB <- MeganDataPlot[which(MeganDataPlot$Origin == "PB"),]
# NG 组气泡图
NG_bubble <- ggplot(MeganDataPlot_NG, aes(x =Sample, y = Genus)) +
geom_point(aes(size=n),shape=21, fill="black") +
theme_classic() +
scale_size_area(max_size = max(MeganDataPlot_NG$n))
# PB 组气泡图
PB_bubble <- ggplot(MeganDataPlot_PB, aes(x =Sample, y = Genus)) +
geom_point(aes(size=n),shape=21, fill="black") +
theme_classic() +
scale_size_area(max_size = max(MeganDataPlot_PB$n))
# 组合图
ggarrange(NG_bubble, PB_bubble + rremove("ylab"), widths = c(2, 2), common.legend = T)
统计分析
Origin 组间差异分析
# 数据准备
MeganData$Origin <- substr(MeganData$Sample, start=1, stop = 2)
MeganDataStat <- MeganData %>% count(Origin, Genus)
# 卡方检验
pvals_genera_originNGvPB <- c()
for(i in 1:length(unique(MeganDataStat$Genus))){
MeganDataStat_test <- MeganDataStat[which(MeganDataStat$Genus == unique(MeganDataStat$Genus)[i]),]
MeganDataStat_test$rest <- 90 - MeganDataStat_test$n
conttab <- rbind(MeganDataStat_test$n, MeganDataStat_test$rest)
test <- chisq.test(conttab)
pvals_genera_originNGvPB <- c(pvals_genera_originNGvPB,test$p.value)
}
# 结果汇总
pvals_genera_originNGvPB <- data.frame(pvals_genera_originNGvPB, unique(MeganDataStat$Genus))
pvals_genera_originNGvPB$padj <- p.adjust(pvals_genera_originNGvPB$pvals_genera_originNGvPB, method="fdr")
print("Origin 组间差异分析结果：")
head(pvals_genera_originNGvPB)
pvals_genera_originNGvPB unique.MeganDataStat.Genus. padj
1 1.640690e-01 2.331507e-01
2 1.759367e-20 Acidovorax 5.278101e-20
3 1.243805e-19 Brevibacillus 3.052975e-19
4 4.560565e-01 Curtobacterium 5.597058e-01
5 1.000000e+00 Frondihabitans 1.000000e+00
6 5.448859e-18 Janthinobacterium 1.225993e-17
NG 组内差异分析
pvals_genera_NG <- c()
for(i in 1:length(unique(MeganDataPlot_NG$Genus))){
MeganDataPlot_NG_test <- MeganDataPlot_NG[which(MeganDataPlot_NG$Genus == unique(MeganDataPlot_NG$Genus)[i]),]
MeganDataPlot_NG_test$rest <- 30 - MeganDataPlot_NG_test$n
conttab <- rbind(MeganDataPlot_NG_test$n, MeganDataPlot_NG_test$rest)
test <- chisq.test(conttab)
pvals_genera_NG <- c(pvals_genera_NG,test$p.value)
}
pvals_genera_NG <- data.frame(pvals_genera_NG, unique(MeganDataPlot_NG$Genus))
pvals_genera_NG$padj <- p.adjust(pvals_genera_NG$pvals_genera_NG, method="fdr")
print("NG 组内差异分析结果：")
pvals_genera_NG unique.MeganDataPlot_NG.Genus. padj
1 6.047729e-01 6.719698e-01
2 5.381482e-01 Curtobacterium 6.719698e-01
3 7.697974e-01 Janthinobacterium 7.697974e-01
4 5.852511e-01 Leifsonia 6.719698e-01
5 3.186355e-07 Methylobacterium 7.080790e-07
6 7.697974e-01 Microbacterium 7.697974e-01
PB 组内差异分析
pvals_genera_PB <- c()
for(i in 1:length(unique(MeganDataPlot_PB$Genus))){
MeganDataPlot_PB_test <- MeganDataPlot_PB[which(MeganDataPlot_PB$Genus == unique(MeganDataPlot_PB$Genus)[i]),]
MeganDataPlot_PB_test$rest <- 30 - MeganDataPlot_PB_test$n
conttab <- rbind(MeganDataPlot_PB_test$n, MeganDataPlot_PB_test$rest)
test <- chisq.test(conttab)
pvals_genera_PB <- c(pvals_genera_PB,test$p.value)
}
pvals_genera_PB <- data.frame(pvals_genera_PB, unique(MeganDataPlot_PB$Genus))
pvals_genera_PB$padj <- p.adjust(pvals_genera_PB$pvals_genera_PB, method="fdr")
print("PB 组内差异分析结果：")
pvals_genera_PB unique.MeganDataPlot_PB.Genus. padj
1 6.376282e-01 7.970352e-01
2 3.186355e-07 Brevundimonas 7.080790e-07
3 2.004595e-01 Clavibacter 3.340992e-01
4 3.186355e-07 Herbiconiux 7.080790e-07
5 2.065286e-06 Leifsonia 4.130572e-06
6 1.000000e+00 Methylobacterium 1.000000e+00
复现结论
本次复现成功生成了论文 Figure 4 相关的所有图表，验证了原作者结果的可复现性。主要发现包括：
- CrossInoc 柱状图展示了不同处理组的微生物相对丰度分布
- 气泡图揭示了属水平的微生物群落结构差异
- 统计分析证实了 Origin 组间和组内的显著差异
