“Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) is a methodology of differential abundance (DA) analysis for microbial absolute abundance data. ANCOM-BC estimates the unknown sampling fractions, corrects the bias induced by their differences among samples, and identifies taxa that are differentially abundant according to the covariate of interest.”

For more details, please refer to the ANCOM-BC paper.

# Subset to baseline
pseq <- readRDS("C:/Congo/psCongo_V4.rds")
metadata <- read.csv(file="C:/Congo/Congo_metadata_V4.csv", header = TRUE, sep=",")
row.names(metadata) <- metadata$sampleId

range01 <- function(x){(x-min(x))/(max(x)-min(x))}

set.seed(11)

scores <- read.csv(file="C:/Congo/final_maternal_stress_scores_20220131.csv", header = TRUE, sep=",")
names(scores)[names(scores) == "id"] <- "Individual"

scores <- scores[scores$Individual != "G6", ]
scores <- scores[scores$Individual != "G18", ]
scores <- scores[scores$Individual != "G32", ]
scores <- scores[scores$Individual != "G44", ]
scores <- scores[scores$Individual != "G49", ]

scaledscores <- scores
scaledscores$general_trauma <- range01(scores$general_trauma)
scaledscores$sexual_events <- range01(scores$sexual_events)
scaledscores$stress <- range01(scores$stress)
scaledscores$ptsd <- range01(scores$ptsd)
scaledscores$violence <- range01(scores$violence)
scaledscores$depression <- range01(scores$depression)
scaledscores$coping <- range01(scores$coping)
scaledscores$anxiety <- range01(scores$anxiety)
scaledscores$pregnancy <- range01(scores$pregnancy)

scaledscores$Composite <- scaledscores$general_trauma + scaledscores$sexual_events + scaledscores$anxiety + scaledscores$depression + scaledscores$ptsd + scaledscores$stress + scaledscores$violence + scaledscores$pregnancy

Composite_categorical <- arules::discretize(scaledscores[,11], method="cluster", breaks = 2, labels = c("Low", "High"))
table(Composite_categorical)
## Composite_categorical
##  Low High 
##   26   21
scores$Composite_categorical <- as.factor(Composite_categorical) ## appends sample table

df <- merge(metadata, scores[,c(10,11)], by="Individual", all.x=TRUE)
row.names(df) <- df$sampleId
df <- sample_data(df)
pseq <- merge_phyloseq(pseq, df)

pseq_t <- transform_sample_counts(pseq, function(x) x/sum(x)*100)

pseq_t <- subset_samples(pseq_t, Type == "sample")

6 Weeks

Antibiotic usage and Baby Sex are the only covariates.

sixw_data <- subset_samples(pseq_t, Time == "6W")
sixw_data = aggregate_taxa(sixw_data, "Species")

sample_data(sixw_data)[,c("Individual", "X6_BabySex", "Med_Antibiotic")]
##        Individual X6_BabySex Med_Antibiotic
## G1_6W          G1   Msichana              0
## G10_6W        G10   Msichana              0
## G11_6W        G11   Msichana              0
## G12_6W        G12     Kijana              0
## G13_6W        G13   Msichana              0
## G14_6W        G14   Msichana              0
## G15_6W        G15     Kijana              0
## G17_6W        G17     Kijana              0
## G19_6W        G19   Msichana              0
## G2_6W          G2     Kijana              0
## G20_6W        G20   Msichana              0
## G21_6W        G21   Msichana              0
## G22_6W        G22   Msichana              0
## G24_6W        G24     Kijana              0
## G25_6W        G25     Kijana              0
## G26_6W        G26     Kijana              0
## G27_6W        G27     Kijana              0
## G29_6W        G29   Msichana              0
## G3_6W          G3     Kijana              0
## G30_6W        G30   Msichana              0
## G33_6W        G33     Kijana              0
## G34_6W        G34     Kijana              0
## G35_6W        G35     Kijana              0
## G37_6W        G37   Msichana              0
## G38_6W        G38   Msichana              0
## G4_6W          G4     Kijana              0
## G40_6W        G40     Kijana              0
## G43_6W        G43     Kijana              0
## G45_6W        G45   Msichana              0
## G46_6W        G46   Msichana              0
## G47_6W        G47   Msichana              0
## G48_6W        G48   Msichana              0
## G50_6W        G50     Kijana              1
## G51_6W        G51   Msichana              0
## G52_6W        G52   Msichana              0
## G8_6W          G8     Kijana              0
## G9_6W          G9     Kijana              0
outsixw = ancombc(phyloseq = sixw_data, formula = "Composite_categorical + Med_Antibiotic + X6_BabySex", p_adj_method = "holm", zero_cut = 0.9, lib_cut = 0, group = "Composite_categorical", struc_zero = TRUE, neg_lb = TRUE, tol = 1e-5, 
              max_iter = 100, conserve = TRUE, alpha = 0.05, global = TRUE)

res = outsixw$res
res_global = outsixw$res_global

tab_diff = res$diff_abn
#colnames(tab_diff) = col_name
tab_diff %>% 
  datatable(caption = "Differentially Abundant Taxa from the Primary Result")
da6w <- tab_diff[which(tab_diff$Composite_categoricalHigh=="TRUE"),]

da6w$Composite_categoricalHigh <- "Six Week"

Create look up table to replace species names

look_data <- subset_samples(pseq_t, Time == "6W")
ps_species <- tax_glom(look_data, "Species")
lookuptable <- tax_table(ps_species)
lookuptable <- as.data.frame(lookuptable)

lookuptable$GenusL <- substr(lookuptable$Genus, 1, 1)

lookuptable$SpeciesL <- paste(lookuptable$GenusL, lookuptable$Species, sep=". ")
res <- outsixw$res

res.or_p <- res$q_val

res.or_p <- res.or_p[order(res.or_p$Composite_categoricalHigh, na.last=NA),]
taxa_sig <- row.names(res.or_p[1:15,])
taxa_sig <- unlist(taxa_sig)

ps.taxa.rel.sig <- prune_taxa(taxa_sig, sixw_data)

matrix <- as.matrix(data.frame(otu_table(ps.taxa.rel.sig)))
metadata_sub <- data.frame(sample_data(ps.taxa.rel.sig))
annotation_col = data.frame(
    Composite = as.factor(metadata_sub$Composite_categorical),  
    check.names = FALSE
)
rownames(annotation_col) = rownames(metadata_sub)

annotation_row = data.frame(
    Phylum = as.factor(tax_table(ps.taxa.rel.sig)[, "Phylum"])
)
rownames(annotation_row) <- rownames(matrix)

phylum_col = RColorBrewer::brewer.pal(length(levels(annotation_row$Phylum)), "Paired")
names(phylum_col) = levels(annotation_row$Phylum)
ann_colors = list(
    Composite = c(`High` = "red", `Low` = "blue"),
    Phylum = phylum_col
)

## Fix the presentation of some of the names
rownames(matrix)[rownames(matrix) == "Bacteria_Firmicutes_Bacilli_Staphylococcales_Staphylococcaceae_Staphylococcus_haemolyticus"] <- "haemolyticus"

rownames(matrix)[rownames(matrix) == "Bacteria_Firmicutes_Negativicutes_Veillonellales-Selenomonadales_Veillonellaceae_Veillonella_dispar"] <- "dispar"

col_fun = colorRamp2(c(min(matrix), 1, max(matrix)/2, max(matrix)), c("white", "#c994c7", "#e7298a", "#67001f"))

## use lookup table to replacenames
row.names(matrix) <- lookuptable$SpeciesL[match(row.names(matrix), lookuptable$Species)]

## remove the timepoint reference in the individual names
for ( col in 1:ncol(matrix)){
    colnames(matrix)[col] <-  sub("_.*", "", colnames(matrix)[col])
}

sixweekda <- ComplexHeatmap::pheatmap(matrix, name="Relative\nAbundance\n(%)", 
                         annotation_col = annotation_col, 
                         annotation_row = annotation_row, 
                         annotation_colors = ann_colors,
                         col = col_fun,
                         column_title = "A:  6 Weeks                                                                            \nIndividuals",
                         column_title_gp = gpar(fontsize = 20),
                         row_title = "Taxa",
                         row_title_gp = gpar(fontsize = 20),
                         heatmap_legend_param = list(at = c(min(matrix), 1, max(matrix)/2, max(matrix)),labels = c(min(matrix), "1", round(max(matrix)/2, digits=0), round(max(matrix), digits=0))))

sixweekda

## save 1000 wide by 500 tall
#ps.taxa.rel.sig

datasetsig <- phyloseq2meco(sixw_data)

datasetsig$tidy_dataset()

taxa_sig <- c(taxa_sig, "haemolyticus", "dispar")

t1 <- trans_abund$new(dataset = datasetsig, taxrank = "Species", input_taxaname = taxa_sig)

## use lookup table to replacenames
t1$data_taxanames <- lookuptable$SpeciesL[match(t1$data_taxanames, lookuptable$Species)]
t1$data_abund$Taxonomy <- lookuptable$SpeciesL[match(t1$data_abund$Taxonomy, lookuptable$Species)]


sixweekboxplot <- t1$plot_box(group = "Composite_categorical") + scale_y_continuous(trans='log10') + theme_bw() +theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("blue", "red")) + scale_color_manual(values=c("black", "black")) + guides(fill=guide_legend(title="Composite\nStress")) + ggtitle("\nTaxa") + theme(plot.title = element_text(hjust = 0.5, margin=margin(50,0,20,0))) + theme(text = element_text(size = 20, face="bold")) 

sixweekboxplot

ggsave("C:/Congo/Manuscript/sixweekbox.png", plot=sixweekboxplot, height=1500, width=1500, units="px", dpi=300)

3 Months

Antibiotic and Baby Sex as covariates.

threem_data <- subset_samples(pseq_t, Time == "3M")
threem_data = aggregate_taxa(threem_data, "Species")

sample_data(threem_data)[,c("Individual", "X6_BabySex", "Med_Antibiotic")]
##        Individual X6_BabySex Med_Antibiotic
## G10_3M        G10   Msichana              0
## G11_3M        G11   Msichana              0
## G12_3M        G12     Kijana              0
## G13_3M        G13   Msichana              0
## G14_3M        G14   Msichana              0
## G15_3M        G15     Kijana              0
## G17_3M        G17     Kijana              0
## G19_3M        G19   Msichana              0
## G2_3M          G2     Kijana              0
## G21_3M        G21   Msichana              0
## G22_3M        G22   Msichana              0
## G23_3M        G23   Msichana              0
## G24_3M        G24     Kijana              0
## G25_3M        G25     Kijana              0
## G26_3M        G26     Kijana              0
## G29_3M        G29   Msichana              0
## G3_3M          G3     Kijana              0
## G30_3M        G30   Msichana              0
## G31_3M        G31     Kijana              0
## G33_3M        G33     Kijana              0
## G34_3M        G34     Kijana              0
## G35_3M        G35     Kijana              0
## G36_3M        G36     Kijana              0
## G38_3M        G38   Msichana              0
## G4_3M          G4     Kijana              0
## G40_3M        G40     Kijana              0
## G41_3M        G41     Kijana              0
## G42_3M        G42     Kijana              1
## G45_3M        G45     Kijana              0
## G46_3M        G46   Msichana              0
## G47_3M        G47   Msichana              0
## G48_3M        G48   Msichana              0
## G5_3M          G5     Kijana              0
## G50_3M        G50     Kijana              1
## G52_3M        G52   Msichana              0
## G7_3M          G7     Kijana              0
## G8_3M          G8     Kijana              0
## G9_3M          G9     Kijana              0
outthreem = ancombc(phyloseq = threem_data, formula = "Composite_categorical + Med_Antibiotic + X6_BabySex", 
              p_adj_method = "holm", zero_cut = 0.9, lib_cut = 0, 
              group = "Composite_categorical", struc_zero = TRUE, neg_lb = TRUE, tol = 1e-5, 
              max_iter = 100, conserve = TRUE, alpha = 0.05, global = TRUE)

res = outthreem$res
res_global = outthreem$res_global

tab_diff = res$diff_abn
#colnames(tab_diff) = col_name
tab_diff %>% 
  datatable(caption = "Differentially Abundant Taxa 
            from the Primary Result")
da3m <- tab_diff[which(tab_diff$Composite_categoricalHigh=="TRUE"),]

da3m$Composite_categoricalHigh <- "Three Month"

Create look up table to replace species names

look_data <- subset_samples(pseq_t, Time == "3M")
ps_species <- tax_glom(look_data, "Species")
lookuptable <- tax_table(ps_species)
lookuptable <- as.data.frame(lookuptable)

lookuptable$GenusL <- substr(lookuptable$Genus, 1, 1)

lookuptable$SpeciesL <- paste(lookuptable$GenusL, lookuptable$Species, sep=". ")
res <- outthreem$res

res.or_p <- res$q_val

res.or_p <- res.or_p[order(res.or_p$Composite_categoricalHigh, na.last=NA),]
taxa_sig <- row.names(res.or_p[1:11,])
taxa_sig <- unlist(taxa_sig)
ps.taxa.rel.sig <- prune_taxa(taxa_sig, threem_data)

matrix <- as.matrix(data.frame(otu_table(ps.taxa.rel.sig)))
metadata_sub <- data.frame(sample_data(ps.taxa.rel.sig))
annotation_col = data.frame(
    Composite = as.factor(metadata_sub$Composite_categorical),  
    check.names = FALSE
)
rownames(annotation_col) = rownames(metadata_sub)

annotation_row = data.frame(
    Phylum = as.factor(tax_table(ps.taxa.rel.sig)[, "Phylum"])
)
rownames(annotation_row) <- rownames(matrix)

phylum_col = RColorBrewer::brewer.pal(length(levels(annotation_row$Phylum)), "Paired")
names(phylum_col) = levels(annotation_row$Phylum)
ann_colors = list(
    Composite = c(`High` = "red", `Low` = "blue"),
    Phylum = phylum_col
)

rownames(matrix)[rownames(matrix) == "Bacteria_Firmicutes_Bacilli_Staphylococcales_Staphylococcaceae_Staphylococcus_haemolyticus"] <- "haemolyticus"

rownames(matrix)[rownames(matrix) == "Bacteria_Firmicutes_Negativicutes_Veillonellales-Selenomonadales_Veillonellaceae_Veillonella_dispar"] <- "dispar"


## remove the timepoint reference in the individual names
for ( col in 1:ncol(matrix)){
    colnames(matrix)[col] <-  sub("_.*", "", colnames(matrix)[col])
}

col_fun = colorRamp2(c(min(matrix), 1, max(matrix)/2, max(matrix)), c("white", "#c994c7", "#e7298a", "#67001f"))

## use lookup table to replacenames
row.names(matrix) <- lookuptable$SpeciesL[match(row.names(matrix), lookuptable$Species)]

threemonthda <- ComplexHeatmap::pheatmap(matrix, name="Relative\nAbundance\n(%)", 
                         annotation_col = annotation_col, 
                         annotation_row = annotation_row, 
                         annotation_colors = ann_colors,
                         col = col_fun,
                         column_title = "B:  3 Months                                                                        \nIndividuals",
                         column_title_gp = gpar(fontsize = 20),
                         row_title = "Taxa",
                         row_title_gp = gpar(fontsize = 20),
                         heatmap_legend_param = list(at = c(min(matrix), 1, max(matrix)/2, max(matrix)),labels = c(min(matrix), "1", round(max(matrix)/2, digits=0), round(max(matrix), digits=0))))

threemonthda

## save 1000 wide by 500 tall
#ps.taxa.rel.sig

datasetsig <- phyloseq2meco(threem_data)

datasetsig$tidy_dataset()

taxa_sig <- c(taxa_sig, "haemolyticus", "dispar")

t1 <- trans_abund$new(dataset = datasetsig, taxrank = "Species", input_taxaname = taxa_sig)
#t1$plot_box(group = "Composite_categorical") + scale_y_continuous(trans='log10')

## use lookup table to replacenames
t1$data_taxanames <- lookuptable$SpeciesL[match(t1$data_taxanames, lookuptable$Species)]
t1$data_abund$Taxonomy <- lookuptable$SpeciesL[match(t1$data_abund$Taxonomy, lookuptable$Species)]

threemonthboxplot <- t1$plot_box(group = "Composite_categorical") + scale_y_continuous(trans='log10')+ theme_bw() +theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("blue", "red")) + scale_color_manual(values=c("black", "black")) + guides(fill=guide_legend(title="Composite\nStress")) + ggtitle("\nTaxa") + theme(plot.title = element_text(hjust = 0.5, margin=margin(50,0,20,0))) + theme(text = element_text(size = 20, face="bold")) 

threemonthboxplot

ggsave("C:/Congo/Manuscript/threemonthbox.png", plot=threemonthboxplot, height=1500, width=1500, units="px", dpi=300)

6 Months

Antibiotic usage and Baby Sex are the only covariates.

sixm_data <- subset_samples(pseq_t, Time == "6M")
sixm_data = aggregate_taxa(sixm_data, "Species")

sample_data(sixm_data)[,c("Individual", "X6_BabySex", "Med_Antibiotic")]
##        Individual X6_BabySex Med_Antibiotic
## G1_6M          G1   Msichana              0
## G10_6M        G10   Msichana              0
## G11_6M        G11   Msichana              1
## G13_6M        G13   Msichana              1
## G14_6M        G14   Msichana              0
## G16_6M        G16   Msichana              1
## G17_6M        G17     Kijana              0
## G2_6M          G2     Kijana              0
## G24_6M        G24     Kijana              0
## G26_6M        G26     Kijana              0
## G27_6M        G27     Kijana              0
## G28_6M        G28   Msichana              0
## G3_6M          G3     Kijana              1
## G30_6M        G30   Msichana              0
## G31_6M        G31     Kijana              0
## G35_6M        G35     Kijana              0
## G36_6M        G36     Kijana              0
## G39_6M        G39     Kijana              1
## G40_6M        G40     Kijana              0
## G42_6M        G42     Kijana              1
## G43_6M        G43     Kijana              0
## G46_6M        G46   Msichana              0
## G48_6M        G48   Msichana              0
## G50_6M        G50     Kijana              1
## G52_6M        G52   Msichana              0
## G7_6M          G7     Kijana              0
## G9_6M          G9     Kijana              0
outsixm = ancombc(phyloseq = sixm_data, formula = "Composite_categorical + Med_Antibiotic + X6_BabySex", 
              p_adj_method = "holm", zero_cut = 0.9, lib_cut = 0, 
              group = "Composite_categorical", struc_zero = TRUE, neg_lb = TRUE, tol = 1e-5, 
              max_iter = 100, conserve = TRUE, alpha = 0.05, global = TRUE)

res = outsixm$res
res_global = outsixm$res_global

tab_diff = res$diff_abn
#colnames(tab_diff) = col_name
tab_diff %>% 
  datatable(caption = "Differentially Abundant Taxa 
            from the Primary Result")
da6m <- tab_diff[which(tab_diff$Composite_categoricalHigh=="TRUE"),]

da6m$Composite_categoricalHigh <- "Six Month"

Create look up table to replace species names

look_data <- subset_samples(pseq_t, Time == "6M")
ps_species <- tax_glom(look_data, "Species")
lookuptable <- tax_table(ps_species)
lookuptable <- as.data.frame(lookuptable)

lookuptable$GenusL <- substr(lookuptable$Genus, 1, 1)

lookuptable$SpeciesL <- paste(lookuptable$GenusL, lookuptable$Species, sep=". ")
res <- outsixm$res

res.or_p <- res$q_val

res.or_p <- res.or_p[order(res.or_p$Composite_categoricalHigh, na.last=NA),]
taxa_sig <- row.names(res.or_p[1:20,])
taxa_sig <- unlist(taxa_sig)
ps.taxa.rel.sig <- prune_taxa(taxa_sig, sixm_data)

matrix <- as.matrix(data.frame(otu_table(ps.taxa.rel.sig)))
metadata_sub <- data.frame(sample_data(ps.taxa.rel.sig))
annotation_col = data.frame(
    Composite = as.factor(metadata_sub$Composite_categorical),  
    check.names = FALSE
)
rownames(annotation_col) = rownames(metadata_sub)

annotation_row = data.frame(
    Phylum = as.factor(tax_table(ps.taxa.rel.sig)[, "Phylum"])
)
rownames(annotation_row) <- rownames(matrix)

phylum_col = RColorBrewer::brewer.pal(length(levels(annotation_row$Phylum)), "Paired")
names(phylum_col) = levels(annotation_row$Phylum)
ann_colors = list(
    Composite = c(`High` = "red", `Low` = "blue"),
    Phylum = phylum_col
)

rownames(matrix)[rownames(matrix) == "Bacteria_Firmicutes_Clostridia_Lachnospirales_Lachnospiraceae_[Ruminococcus] torques group_faecis"] <- "faecis"


## remove the timepoint reference in the individual names
for ( col in 1:ncol(matrix)){
    colnames(matrix)[col] <-  sub("_.*", "", colnames(matrix)[col])
}

col_fun = colorRamp2(c(min(matrix), 1, max(matrix)/2, max(matrix)), c("white", "#c994c7", "#e7298a", "#67001f"))

## use lookup table to replacenames
row.names(matrix) <- lookuptable$SpeciesL[match(row.names(matrix), lookuptable$Species)]

rownames(matrix)[rownames(matrix) == "[. faecis"] <- "R. torques"
rownames(matrix)[rownames(matrix) == "[. innocuum"] <- "C. innocuum"


sixmonthda <- ComplexHeatmap::pheatmap(matrix, name="Relative\nAbundance\n(%)", 
                         annotation_col = annotation_col, 
                         annotation_row = annotation_row, 
                         annotation_colors = ann_colors,
                         col = col_fun,
                         column_title = "C:  6 Months                                                                             \nIndividuals",
                         column_title_gp = gpar(fontsize = 20),
                         row_title = "Taxa",
                         row_title_gp = gpar(fontsize = 20),
                         heatmap_legend_param = list(at = c(min(matrix), 1, max(matrix)/2, max(matrix)),labels = c(min(matrix), "1", round(max(matrix)/2, digits=0), round(max(matrix), digits=0))))

sixmonthda

datasetsig <- phyloseq2meco(sixm_data)

datasetsig$tidy_dataset()

taxa_sig <- c(taxa_sig, "faecis", "[Ruminococcus] torques group_faecis")
taxa_sig <- c(taxa_sig, "R. torques")
taxa_sig <- c(taxa_sig, "C. innocuum")

t1 <- trans_abund$new(dataset = datasetsig, taxrank = "Species", input_taxaname = taxa_sig)
#t1$plot_box(group = "Composite_categorical") + scale_y_continuous(trans='log10')

## use lookup table to replacenames
t1$data_taxanames <- lookuptable$SpeciesL[match(t1$data_taxanames, lookuptable$Species)]
t1$data_abund$Taxonomy <- lookuptable$SpeciesL[match(t1$data_abund$Taxonomy, lookuptable$Species)]

t1$data_taxanames[t1$data_taxanames == "[. innocuum"] <- "C. Innocuum"
t1$data_abund$Taxonomy[t1$data_abund$Taxonomy == "[. innocuum"] <- "C. Innocuum"

t1$data_taxanames[t1$data_taxanames == "[. faecis"] <- "R. torques"
t1$data_abund$Taxonomy[t1$data_abund$Taxonomy == "[. faecis"] <- "R. torques"

sixmonthboxplot <- t1$plot_box(group = "Composite_categorical") + scale_y_continuous(trans='log10')+ theme_bw() +theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_manual(values=c("blue", "red")) + scale_color_manual(values=c("black", "black")) + guides(fill=guide_legend(title="Composite\nStress")) + ggtitle("\nTaxa") + theme(plot.title = element_text(hjust = 0.5, margin=margin(50,0,20,0))) + theme(text = element_text(size = 20, face="bold")) 

sixmonthboxplot

ggsave("C:/Congo/Manuscript/sixmonthbox.png", plot=sixmonthboxplot, height=1500, width=1500, units="px", dpi=300)
## because of an issue with the package complexpheatmap, you need to manually save the heatmap figures the first time, and then load them here for the multifigure. The single panes are always dynamically made each time this markdown document is run.

# plot1 <- readPNG("C:/Congo/Manuscript/sixweekda.png") # saved as 1000 by 500
# plot2 <- readPNG("C:/Congo/Manuscript/threemonthda.png") # saved as 1000 by 500
# plot3 <- readPNG("C:/Congo/Manuscript/sixmonthda.png") # saved as 1000 by 500
# twoplot <- grid.arrange(rasterGrob(plot1),rasterGrob(plot2),rasterGrob(plot3),ncol=1) 

## save 1000 x 1500

## This is the multipane boxplot figure

# plot_grid(sixweekboxplot + theme(legend.background = element_blank(), plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0))) +
#   theme(legend.position = c(0.92, 0.78)) + ggtitle(""), threemonthboxplot + theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0))) +
#   theme(legend.position="none") + ggtitle(""), sixmonthboxplot + theme(plot.title = element_text(hjust = 0.5, margin=margin(0,0,0,0))) +
#   theme(legend.position="none") + ggtitle(""), labels = c('A: 6 Weeks', 'B: 3 Months', 'C: 6 Months'), label_size = 20, ncol=1)

## save 1317 x 1705

Session information

sessionInfo()
R version 4.1.3 (2022-03-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] cowplot_1.1.1     gridExtra_2.3     png_0.1-7         circlize_0.4.15  
 [5] file2meco_0.4.0   microeco_0.12.1   DT_0.24           ANCOMBC_1.4.0    
 [9] qwraps2_0.5.2     magrittr_2.0.3    microbiome_1.16.0 phyloseq_1.38.0  
[13] forcats_0.5.2     stringr_1.4.1     dplyr_1.0.10      purrr_0.3.4      
[17] readr_2.1.2       tidyr_1.2.0       tibble_3.1.8      ggplot2_3.3.6    
[21] tidyverse_1.3.2  

loaded via a namespace (and not attached):
  [1] readxl_1.4.1           backports_1.4.1        systemfonts_1.0.4     
  [4] plyr_1.8.7             igraph_1.3.4           splines_4.1.3         
  [7] crosstalk_1.2.0        GenomeInfoDb_1.30.1    digest_0.6.29         
 [10] foreach_1.5.2          htmltools_0.5.3        arules_1.7-4          
 [13] fansi_1.0.3            googlesheets4_1.0.1    cluster_2.1.4         
 [16] doParallel_1.0.17      tzdb_0.3.0             ComplexHeatmap_2.10.0 
 [19] Biostrings_2.61.2      modelr_0.1.9           matrixStats_0.62.0    
 [22] colorspace_2.0-3       rvest_1.0.3            textshaping_0.3.6     
 [25] haven_2.5.1            rbibutils_2.2.9        xfun_0.32             
 [28] crayon_1.5.1           RCurl_1.98-1.8         jsonlite_1.8.0        
 [31] survival_3.4-0         iterators_1.0.14       ape_5.6-2             
 [34] glue_1.6.2             gtable_0.3.1           gargle_1.2.0          
 [37] zlibbioc_1.39.0        XVector_0.33.0         GetoptLong_1.0.5      
 [40] Rhdf5lib_1.15.2        shape_1.4.6            BiocGenerics_0.40.0   
 [43] scales_1.2.1           DBI_1.1.3              Rcpp_1.0.9            
 [46] clue_0.3-61            stats4_4.1.3           htmlwidgets_1.5.4     
 [49] httr_1.4.4             RColorBrewer_1.1-3     ellipsis_0.3.2        
 [52] pkgconfig_2.0.3        farver_2.1.1           sass_0.4.2            
 [55] dbplyr_2.2.1           utf8_1.2.2             tidyselect_1.1.2      
 [58] labeling_0.4.2         rlang_1.0.6            reshape2_1.4.4        
 [61] munsell_0.5.0          cellranger_1.1.0       tools_4.1.3           
 [64] cachem_1.0.6           cli_3.3.0              generics_0.1.3        
 [67] ade4_1.7-19            broom_1.0.1            evaluate_0.16         
 [70] biomformat_1.22.0      fastmap_1.1.0          yaml_2.3.5            
 [73] ragg_1.2.2             knitr_1.40             fs_1.5.2              
 [76] nlme_3.1-159           xml2_1.3.3             compiler_4.1.3        
 [79] rstudioapi_0.14        reprex_2.0.2           bslib_0.4.0           
 [82] stringi_1.7.8          highr_0.9              lattice_0.20-45       
 [85] Matrix_1.4-1           nloptr_2.0.3           vegan_2.6-2           
 [88] permute_0.9-7          multtest_2.49.0        vctrs_0.4.1           
 [91] pillar_1.8.1           lifecycle_1.0.3        rhdf5filters_1.5.0    
 [94] Rdpack_2.4             jquerylib_0.1.4        GlobalOptions_0.1.2   
 [97] data.table_1.14.2      bitops_1.0-7           R6_2.5.1              
[100] IRanges_2.27.2         codetools_0.2-18       MASS_7.3-58.1         
[103] assertthat_0.2.1       rhdf5_2.37.4           rjson_0.2.21          
[106] withr_2.5.0            S4Vectors_0.31.5       GenomeInfoDbData_1.2.7
[109] mgcv_1.8-40            parallel_4.1.3         hms_1.1.2             
[112] rmarkdown_2.16         googledrive_2.0.0      Rtsne_0.16            
[115] Biobase_2.53.0         lubridate_1.8.0