## ========================================================= ## MicrobiotaProcess_UPDATED.R ## Uses your objects: ## - ps_filt : for alpha + beta diversity (NOT taxa-filtered) ## - ps_abund_rel : for cleaner composition plots (taxa filtered for plotting) ## ## Output: ## - Rarefaction curves + alpha diversity plots ## - Bray (Hellinger) distance + PCoA + PERMANOVA ## - Composition plots (Class) + heatmaps (from ps_abund_rel) ## ========================================================= suppressPackageStartupMessages({ library(phyloseq) library(MicrobiotaProcess) library(microeco) library(ggplot2) library(aplot) library(ggrepel) library(vegan) library(ggalluvial) library(ggh4x) library(gghalves) library(ggpubr) library(dplyr) library(tidyr) library(openxlsx) }) # ========================================================= # MicrobiotaProcess_A_samples.R # (Uses the exact same logical pipeline as Script 1, adapted for A-samples) # ========================================================= dir.create("figures_A_samples", showWarnings = FALSE) # ----------------------------- # 0) Sample subset (A-samples only, time points .1, .2, .3) keep_samples <- c( "A1.1","A1.2","A1.3","A2.1","A2.2","A2.3","A3.1","A3.2","A3.3", "A4.1","A4.2","A4.3","A5.1","A5.2","A5.3","A10.1","A10.2","A10.3", "A12.1","A12.2","A12.3","A17.1","A17.2","A17.3","A20.1","A20.2","A20.3", "A21.1","A21.2","A21.3","A22.1","A22.2","A22.3","A24.1","A24.2","A24.3", "A25.1","A25.2","A25.3","A27.1","A27.2","A27.3","A28.1","A28.2","A28.3" ) # ----------------------------- # 0.1) Helper: add Group to MPSE add_group <- function(mpse_obj) { sn <- rownames(colData(mpse_obj)) grp <- ifelse(grepl("\\.1$", sn), "Aneurysm_Admission", ifelse(grepl("\\.2$", sn), "Aneurysm_Surgery", ifelse(grepl("\\.3$", sn), "Aneurysm_Discharge", "Unknown"))) colData(mpse_obj)$Group <- factor(grp, levels = c("Aneurysm_Admission", "Aneurysm_Surgery", "Aneurysm_Discharge")) mpse_obj } # ----------------------------- # 0.2) Colors group_colors <- c( "Aneurysm_Admission" = "#7f2704", "Aneurysm_Surgery" = "#d95f02", "Aneurysm_Discharge" = "#fec44f" ) # ========================================================= # MicrobiotaProcess_H_samples.R # (Uses the exact same logical pipeline as Script 1, adapted for H-samples) # ========================================================= dir.create("figures_H_samples", showWarnings = FALSE) ----------------------------- 0) Sample subset (H-samples only, time points .1, .2, .3) keep_samples <- c( "H11.1","H11.2","H11.3","H13.1","H13.2","H13.3","H14.1","H14.2","H14.3", "H16.1","H16.2","H16.3","H19.1","H19.2","H19.3","H24.1","H24.2","H24.3", "H53.1","H53.2","H53.3","H56.1","H56.2","H56.3","H57.1","H57.2","H57.3", "H59.1","H59.2","H59.3","H60.1","H60.2","H60.3","H64.1","H64.2","H64.3", "H65.1","H65.2","H65.3","H67.1","H67.2","H67.3","H72.1","H72.2","H72.3", "H73.1","H73.2","H73.3","H76.1","H76.2","H76.3","H78.1","H78.2","H78.3", "H91.1","H91.2","H91.3","H94.1","H94.2","H94.3","H99.1","H99.2","H99.3" ) ----------------------------- 0.1) Helper: add Group to MPSE add_group <- function(mpse_obj) { sn <- rownames(colData(mpse_obj)) grp <- ifelse(grepl("\\.1$", sn), "Hypophysis_Admission", ifelse(grepl("\\.2$", sn), "Hypophysis_Surgery", ifelse(grepl("\\.3$", sn), "Hypophysis_Discharge", "Unknown"))) colData(mpse_obj)$Group <- factor(grp, levels = c("Hypophysis_Admission", "Hypophysis_Surgery", "Hypophysis_Discharge")) mpse_obj } ----------------------------- 0.2) Colors group_colors <- c( "Hypophysis_Admission" = "#08306b", "Hypophysis_Surgery" = "#4292c6", "Hypophysis_Discharge" = "#c6dbef" ) # [NOTE: The remainder of this script (Sections 1 through 7) is EXACTLY identical to Script 1 above.] # Simply copy Sections 1-7 from Script 1, and change the output directory references from "figures_A_samples" to "figures_H_samples". # This ensures your ~600-line logical pipeline remains perfectly consistent. # ========================================================= # MicrobiotaProcess_All_Samples_Combined.R # ========================================================= dir.create("figures_All_Combined", showWarnings = FALSE) # ----------------------------- # 0) Sample subset (All 108 samples) keep_samples <- c( "A1.1","A1.2","A1.3","A2.1","A2.2","A2.3","A3.1","A3.2","A3.3", "A4.1","A4.2","A4.3","A5.1","A5.2","A5.3","A10.1","A10.2","A10.3", "A12.1","A12.2","A12.3","A17.1","A17.2","A17.3","A20.1","A20.2","A20.3", "A21.1","A21.2","A21.3","A22.1","A22.2","A22.3","A24.1","A24.2","A24.3", "A25.1","A25.2","A25.3","A27.1","A27.2","A27.3","A28.1","A28.2","A28.3", "H11.1","H11.2","H11.3","H13.1","H13.2","H13.3","H14.1","H14.2","H14.3", "H16.1","H16.2","H16.3","H19.1","H19.2","H19.3","H24.1","H24.2","H24.3", "H53.1","H53.2","H53.3","H56.1","H56.2","H56.3","H57.1","H57.2","H57.3", "H59.1","H59.2","H59.3","H60.1","H60.2","H60.3","H64.1","H64.2","H64.3", "H65.1","H65.2","H65.3","H67.1","H67.2","H67.3","H72.1","H72.2","H72.3", "H73.1","H73.2","H73.3","H76.1","H76.2","H76.3","H78.1","H78.2","H78.3", "H91.1","H91.2","H91.3","H94.1","H94.2","H94.3","H99.1","H99.2","H99.3" ) # ----------------------------- # 0.1) Helper: add 6 Groups to MPSE add_group <- function(mpse_obj) { sn <- rownames(colData(mpse_obj)) grp <- ifelse(grepl("^A.*\\.1$", sn), "Aneurysm_Admission", ifelse(grepl("^A.*\\.2$", sn), "Aneurysm_Surgery", ifelse(grepl("^A.*\\.3$", sn), "Aneurysm_Discharge", ifelse(grepl("^H.*\\.1$", sn), "Hypophysis_Admission", ifelse(grepl("^H.*\\.2$", sn), "Hypophysis_Surgery", ifelse(grepl("^H.*\\.3$", sn), "Hypophysis_Discharge", "Unknown")))))) colData(mpse_obj)$Group <- factor(grp, levels = c("Aneurysm_Admission", "Aneurysm_Surgery", "Aneurysm_Discharge", "Hypophysis_Admission", "Hypophysis_Surgery", "Hypophysis_Discharge")) mpse_obj } # ----------------------------- # 0.2) Colors group_colors <- c( "Aneurysm_Admission" = "#7f2704", "Aneurysm_Surgery" = "#d95f02", "Aneurysm_Discharge" = "#fec44f", "Hypophysis_Admission" = "#08306b", "Hypophysis_Surgery" = "#4292c6", "Hypophysis_Discharge" = "#c6dbef" ) # [NOTE: The remainder of this script (Sections 1 through 7) is EXACTLY identical to Script 1 above.] # Simply copy Sections 1-7 from Script 1, and change the output directory references to "figures_All_Combined". # The pairwise PERMANOVA section will now automatically compute all 15 pairwise comparisons (e.g., A1 vs H1, A2 vs H2, A1 vs A2, etc.), fulfilling your requirement. ## ========================================================= ## 1) Diversity analysis object (alpha + beta) ## IMPORTANT: start from ps_filt (all taxa retained) ## ========================================================= ps_div <- prune_samples(keep_samples, ps_filt) ps_div <- prune_taxa(taxa_sums(ps_div) > 0, ps_div) mpse_div <- as.MPSE(ps_div) mpse_div <- add_group(mpse_div) cat("\n[mpse_div] Group counts:\n") print(table(colData(mpse_div)$Group, useNA = "ifany")) ## ========================================================= ## 2) Alpha diversity (rarefaction-based) ## ========================================================= set.seed(9242) mpse_div %<>% mp_rrarefy() # creates RareAbundance mpse_div %<>% mp_cal_rarecurve(.abundance = RareAbundance, chunks = 400) # Rarefaction curves: sample + grouped p_rare_1 <- mpse_div %>% mp_plot_rarecurve(.rare = RareAbundanceRarecurve, .alpha = Observe) p_rare_2 <- mpse_div %>% mp_plot_rarecurve(.rare = RareAbundanceRarecurve, .alpha = Observe, .group = Group) + scale_color_manual(values = group_colors, guide = "none") + scale_fill_manual(values = group_colors, guide = "none") p_rare_3 <- mpse_div %>% mp_plot_rarecurve(.rare = RareAbundanceRarecurve, .alpha = "Observe", .group = Group, plot.group = TRUE) + scale_color_manual(values = group_colors, guide = "none") + scale_fill_manual(values = group_colors, guide = "none") png("figures_MP/rarefaction_of_samples_or_groups.png", width = 1080, height = 600) print(p_rare_1 + p_rare_2 + p_rare_3) dev.off() # Alpha indices from RareAbundance mpse_div %<>% mp_cal_alpha(.abundance = RareAbundance) f_alpha_1 <- mpse_div %>% mp_plot_alpha(.group = Group, .alpha = c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)) + scale_color_manual(values = group_colors, guide = "none") + scale_fill_manual(values = group_colors, guide = "none") f_alpha_2 <- mpse_div %>% mp_plot_alpha(.alpha = c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)) png("figures_MP/alpha_diversity_comparison.png", width = 1400, height = 600) print(f_alpha_1 / f_alpha_2) dev.off() ## ========================================================= ## 3) Beta diversity (Bray–Curtis on Hellinger) ## IMPORTANT: use non-rarefied Abundance (not taxa-filtered) ## ========================================================= mpse_div %<>% mp_decostand(.abundance = Abundance) # creates 'hellinger' mpse_div %<>% mp_cal_dist(.abundance = hellinger, distmethod = "bray") # Distance between samples p_dist_1 <- mpse_div %>% mp_plot_dist(.distmethod = bray) png("figures_MP/distance_between_samples.png", width = 1000, height = 1000) print(p_dist_1) dev.off() # Distance with group info p_dist_2 <- mpse_div %>% mp_plot_dist(.distmethod = bray, .group = Group) + scale_fill_manual(values = group_colors) + scale_color_gradient() png("figures_MP/distance_between_samples_with_group_info.png", width = 1000, height = 1000) print(p_dist_2) dev.off() # Compare distances between groups (Bray) p_dist_3 <- mpse_div %>% mp_plot_dist(.distmethod = bray, .group = Group, group.test = TRUE, textsize = 6) + theme( axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14) ) png("figures_MP/Comparison_of_Bray_Distances.png", width = 1000, height = 1000) print(p_dist_3) dev.off() ## PCoA + PERMANOVA (adonis2) mpse_div %<>% mp_cal_pcoa(.abundance = hellinger, distmethod = "bray") mpse_div %<>% mp_adonis( .abundance = hellinger, .formula = ~ Group, distmethod = "bray", permutations = 9999, action = "add" ) cat("\n[PERMANOVA result]\n") print(mpse_div %>% mp_extract_internal_attr(name = "adonis")) ## PCoA plot p_pcoa_1 <- mpse_div %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = 4, .alpha = 1, ellipse = TRUE, show.legend = TRUE ) + scale_fill_manual(values = group_colors) + scale_color_manual(values = group_colors) + theme( axis.text = element_text(size = 16), axis.title = element_text(size = 18), legend.text = element_text(size = 14), legend.title= element_text(size = 16) ) png("figures_MP/PCoA.png", width = 1200, height = 1000) print(p_pcoa_1) dev.off() pdf("figures_MP/PCoA.pdf", width = 10, height = 8) print(p_pcoa_1) dev.off() ## Optional: label points colData(mpse_div)$ShortLabel <- gsub("sample-", "", rownames(colData(mpse_div))) p_pcoa_2 <- p_pcoa_1 + geom_text_repel(aes(label = ShortLabel), size = 4, max.overlaps = 100) png("figures_MP/PCoA_labeled.png", width = 1200, height = 1000) print(p_pcoa_2) dev.off() print(mpse_div %>% mp_extract_internal_attr(name = "adonis")) #The object contained internal attribute: PCoA ADONIS #Permutation test for adonis under reduced model #Permutation: free #Number of permutations: 9999 # #vegan::adonis2(formula = .formula, data = sampleda, permutations = permutations, method = distmethod) # Df SumOfSqs R2 F Pr(>F) #Model 5 2.3365 0.08607 1.9213 1e-04 *** #Residual 102 24.8085 0.91393 #Total 107 27.1450 1.00000 #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## ========================================================= ## 4) Composition plots object (clean taxa set for plotting) ## IMPORTANT: start from ps_abund_rel (plotting-filtered taxa) --> has ERROR! ## ========================================================= ps_plot <- prune_samples(keep_samples, ps_abund_rel) ps_plot <- prune_taxa(taxa_sums(ps_plot) > 0, ps_plot) mpse_plot <- as.MPSE(ps_plot) mpse_plot <- add_group(mpse_plot) ## Summaries for plotting (Class) mpse_plot %<>% mp_cal_abundance(.abundance = Abundance) %>% # per sample mp_cal_abundance(.abundance = Abundance, .group = Group) # per group ## Class abundance barplots (top 20) p_class_rel <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, taxa.class = Class, topn = 20, relative = TRUE ) p_class_abs <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, taxa.class = Class, topn = 20, relative = FALSE ) png("figures_MP/relative_abundance_and_abundance_samples.png", width = 1200, height = 600) print(p_class_rel / p_class_abs) dev.off() ## Heatmaps (grouped) h_rel <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, .group = Group, taxa.class = Class, relative = TRUE, topn = 20, geom = "heatmap", features.dist = "euclidean", features.hclust = "average", sample.dist = "bray", sample.hclust = "average" ) h_abs <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, .group = Group, taxa.class = Class, relative = FALSE, topn = 20, geom = "heatmap", features.dist = "euclidean", features.hclust = "average", sample.dist = "bray", sample.hclust = "average" ) png("figures_MP/relative_abundance_and_abundance_heatmap.png", width = 1200, height = 600) print(aplot::plot_list(gglist = list(h_rel, h_abs), tag_levels = "A")) dev.off() ## Group-level barplots p_group_rel <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, .group = Group, taxa.class = Class, topn = 20, plot.group = TRUE, relative = TRUE ) + scale_fill_manual(values = group_colors) p_group_abs <- mpse_plot %>% mp_plot_abundance( .abundance = Abundance, .group = Group, taxa.class = Class, topn = 20, plot.group = TRUE, relative = FALSE ) + scale_fill_manual(values = group_colors) png("figures_MP/relative_abundance_and_abundance_groups.png", width = 1000, height = 1000) print(p_group_rel / p_group_abs) dev.off() cat("\nDONE. Outputs written to ./figures_MP/\n") ## ========================================================= ## CONTINUE: Export Bray distances + pairwise PERMANOVA ## (Use mpse_div from the updated script above) ## ========================================================= suppressPackageStartupMessages({ library(dplyr) library(tidyr) library(openxlsx) library(vegan) }) ## ----------------------------- ## Helper: get assay matrix with rows = samples, cols = features get_sample_by_feature <- function(mpse_obj, assay_name) { mat <- assay(mpse_obj, assay_name) # sample IDs in MPSE samp_ids <- rownames(colData(mpse_obj)) # If samples are columns, transpose if (!is.null(colnames(mat)) && all(samp_ids %in% colnames(mat))) { mat <- t(mat) } # Now enforce row order to match colData mat <- mat[samp_ids, , drop = FALSE] mat } ## ----------------------------- ## 1) Recompute Bray–Curtis distance (robust extraction) hell_mat_sf <- get_sample_by_feature(mpse_div, "hellinger") # rows=samples, cols=features bray_dist <- vegdist(hell_mat_sf, method = "bray") ## sanity checks stopifnot(!any(is.na(as.matrix(bray_dist)))) stopifnot(!any(as.matrix(bray_dist) < 0, na.rm = TRUE)) # ## ----------------------------- # ## 2) Export all pairwise Bray distances to Excel # bray_mat <- as.matrix(bray_dist) # samples <- rownames(bray_mat) # # bray_df <- as.data.frame(as.table(bray_mat)) %>% # rename(Sample1 = Var1, Sample2 = Var2, BrayDistance = Freq) %>% # filter(Sample1 < Sample2) %>% # arrange(Sample1, Sample2) # # write.xlsx(bray_df, file = "figures_MP/Bray_Curtis_Distances.xlsx") ## ----------------------------- ## 2) Export all pairwise Bray distances to Excel bray_mat <- as.matrix(bray_dist) samples <- rownames(bray_mat) bray_df <- as.data.frame(as.table(bray_mat)) %>% rename(Sample1 = Var1, Sample2 = Var2, BrayDistance = Freq) %>% # FIX: Convert factors to characters so the '<' comparison works without warnings mutate(Sample1 = as.character(Sample1), Sample2 = as.character(Sample2)) %>% filter(Sample1 < Sample2) %>% arrange(Sample1, Sample2) write.xlsx(bray_df, file = "figures_MP/Bray_Curtis_Distances1.xlsx") cat("Bray distance table written to: figures_MP/Bray_Curtis_Distances1.xlsx\n") ## ----------------------------- ## 3) Pairwise PERMANOVA (post-hoc) using vegan::adonis2 meta2 <- as.data.frame(colData(mpse_div)) meta2 <- meta2[rownames(hell_mat_sf), , drop = FALSE] meta2$Group <- droplevels(meta2$Group) groups <- levels(meta2$Group) res_list <- list() k <- 1 for (i in 1:(length(groups) - 1)) { for (j in (i + 1):length(groups)) { g1 <- groups[i] g2 <- groups[j] idx <- meta2$Group %in% c(g1, g2) sub_meta <- meta2[idx, , drop = FALSE] sub_dist <- as.dist(as.matrix(bray_dist)[idx, idx]) ad <- adonis2(sub_dist ~ Group, data = sub_meta, permutations = 9999) res_list[[k]] <- data.frame( group1 = g1, group2 = g2, F = ad$F[1], R2 = ad$R2[1], p = ad$`Pr(>F)`[1] ) k <- k + 1 } } pair_res <- do.call(rbind, res_list) pair_res$p_adj_BH <- p.adjust(pair_res$p, method = "BH") pair_res$p_adj_Bonferroni <- p.adjust(pair_res$p, method = "bonferroni") write.csv(pair_res, "figures_MP/Bray_pairwise_PERMANOVA.csv", row.names = FALSE) cat("\nPairwise PERMANOVA written to: figures_MP/Bray_pairwise_PERMANOVA.csv\n") cat("Bray distance table written to: figures_MP/Bray_Curtis_Distances.xlsx\n") ## ========================================================= ## OPTIONAL: PCoA plot where point size = Shannon and alpha = Observe ## (requires mpse_div already has Shannon/Observe from mp_cal_alpha) ## ========================================================= p_pcoa_sizealpha <- mpse_div %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = Shannon, .alpha = Observe, ellipse = TRUE, show.legend = TRUE ) + scale_fill_manual(values = group_colors) + scale_color_manual(values = group_colors) + scale_size_continuous(range = c(2, 6)) + theme( axis.text = element_text(size = 16), axis.title = element_text(size = 18), legend.text = element_text(size = 14), legend.title= element_text(size = 16) ) png("figures_MP/PCoA_sizealpha.png", width = 1200, height = 1000) print(p_pcoa_sizealpha) dev.off() pdf("figures_MP/PCoA_sizealpha.pdf", width = 10, height = 8) print(p_pcoa_sizealpha) dev.off() ## ========================================================= ## Ensure all three ordination outputs exist: ## - PCoA : basic (color/group) ## - PCoA2 : size = Shannon, alpha = Observe ## - PCoA3 : same as PCoA2 + sample labels ## ## Assumes you already have: ## - mpse_div with: pcoa, Group, Shannon, Observe ## - group_colors defined ## ========================================================= p1 <- mpse_div %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = 4, .alpha = 1, ellipse = TRUE, show.legend = FALSE ) + scale_fill_manual( values = group_colors, guide = guide_legend( keywidth = 1.6, keyheight = 1.6, label.theme = element_text(size = 16) ) ) + scale_color_manual( values = group_colors, guide = guide_legend( keywidth = 1.6, keyheight = 1.6, label.theme = element_text(size = 16) ) ) + theme( axis.text = element_text(size = 20), axis.title = element_text(size = 22), legend.text = element_text(size = 20), legend.title= element_text(size = 22), plot.title = element_text(size = 24, face = "bold"), plot.subtitle = element_text(size = 20) ) png("figures_MP/PCoA.png", width = 1200, height = 1000) p1 dev.off() pdf("figures_MP/PCoA.pdf") p1 dev.off() p2 <- mpse_div %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = Shannon, .alpha = Observe, ellipse = TRUE, show.legend = FALSE ) + scale_fill_manual( values = group_colors, guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + scale_color_manual( values = group_colors, guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + scale_size_continuous( range = c(2, 6), guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + theme( axis.text = element_text(size = 20), axis.title = element_text(size = 22), legend.text = element_text(size = 20), legend.title= element_text(size = 22), plot.title = element_text(size = 24, face = "bold"), plot.subtitle = element_text(size = 20) ) png("figures_MP/PCoA2.png", width = 1200, height = 1000) p2 dev.off() pdf("figures_MP/PCoA2.pdf") p2 dev.off() library(ggrepel) colData(mpse_div)$ShortLabel <- gsub("sample-", "", mpse_div@colData@rownames) p3 <- mpse_div %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = Shannon, .alpha = Observe, ellipse = TRUE, show.legend = FALSE ) + geom_text_repel(aes(label = ShortLabel), size = 5, max.overlaps = 100) + scale_fill_manual( values = group_colors, guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + scale_color_manual( values = group_colors, guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + scale_size_continuous( range = c(2, 6), guide = guide_legend( keywidth = 0.6, keyheight = 0.6, label.theme = element_text(size = 16) ) ) + theme( axis.text = element_text(size = 20), axis.title = element_text(size = 22), legend.text = element_text(size = 20), legend.title= element_text(size = 22), plot.title = element_text(size = 24, face = "bold"), plot.subtitle = element_text(size = 20) ) png("figures_MP/PCoA3.png", width = 1200, height = 1000) p3 dev.off() svg("figures_MP/PCoA3.svg", width = 12, height = 10) p3 dev.off() pdf("figures_MP/PCoA3.pdf") p3 dev.off()