Queen-killing Melissococcus plutonius study analysis
Brendan Daisley
September 04, 2026
1. Project overview
This RMarkdown document contains the analysis associated with the Case of the Bad Queen Cells project, which investigates an unusual disease outbreak affecting developing honey bee (Apis mellifera) queens in commercial queen-rearing colonies. The project combines pathogen screening, microbiome profiling, microbial isolation, comparative genomics, and experimental infection assays to identify factors associated with queen developmental failure and characterize the microbial and genomic features underlying the disease phenotype.
- Related article: Daisley et al. (2026). Queen-killing
variants of Melissococcus plutonius exhibit phage-mediated genome
diversification and are associated with black queen cell virus
co-infection. In draft
- Related GitHub repository: https://github.com/bdaisley/BadQueenCells
The R code below was used to generate the following figures
presented in the related article:
2. Run auxiliary scripts
Run the bash scripts below to automatically download FASTQ files from SRA, quality filter, and denoise sequences to generate amplicon sequencing variants. Note: This step is for comprehension and can be skipped for the purpose of this tutorial since the counts table and denoised sequences are already available in this repository.
33 #-------------------------------------------------------
34 # V3V4-16S rRNA microbiota profiling workflow
35 #-------------------------------------------------------
36 #bash scripts/V3V4_01_download.sh
37 #bash scripts/V3V4_02_cutadapt.sh
38 #bash scripts/V3V4_03_fastqc.sh
39 #bash scripts/V3V4_04_denoising.sh
40 #Rscript scripts/V3V4_05_classification.R
41
42 #-------------------------------------------------------
43 # Shotgun metagenomic sequencing workflow
44 #-------------------------------------------------------
45
46 #bash scripts/shotgun_01_download.sh
47 #bash scripts/shotgun_02_bowtie2.sh
48 #bash scripts/shotgun_03_fastp.sh
49 #bash scripts/shotgun_04_kraken2.sh
50 #bash scripts/shotgun_05_metacerberus.sh
51 #bash scripts/shotgun_06_eggNOG_mapper2.sh
52
53 #-------------------------------------------------------
54 # Longread genome assembly workflow
55 #-------------------------------------------------------
56 #bash scripts/longread_01_download.sh
57 #bash scripts/longread_02_genome_assembly_pipeline.sh
58 #bash scripts/longread_03_prodigal.sh
59 #bash scripts/longread_04_metacerberus.sh
60 #python scripts/longread_04_metacerberus_hydra.py
61 #bash scripts/longread_05_MOBsuite.sh
62 #python scripts/longread_06_VIBRANT.py
63 #bash scripts/longread_06_VIBRANT.sh
64 3. Load required packages
65 library(dplyr)
66 library(ggplot2)
67 library(stringr)
68 library(readr)
69 library(reshape2)
70 library(ggdendro)
71 library(cowplot)
72 library(dendsort)
73 library(ggiraph)
74 library(ggordiplots)
75 library(igraph)
76 library(NetCoMi)
77 library(phyloseq)
78 library(metagMisc)
79 library(Maaslin2)
80 library(ggrepel)
81 library(rfPermute)
82 library(randomForest)
83 library(pROC)
84 library(caret)
85 library(glmnet)
86 library(readxl)
87
88 theme_set(theme_bw())
89
90 sessionInfo()## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/Toronto
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] readxl_1.4.3 glmnet_4.1-8 Matrix_1.7-1
## [4] caret_6.0-94 pROC_1.18.5 randomForest_4.7-1.1
## [7] rfPermute_2.5.4 ggrepel_0.9.6 Maaslin2_1.18.0
## [10] metagMisc_0.5.0 phyloseq_1.48.0 NetCoMi_1.2.0
## [13] SpiecEasi_1.1.3 igraph_2.0.3 ggordiplots_0.4.3
## [16] glue_1.7.0 vegan_2.6-10 lattice_0.22-5
## [19] permute_0.9-7 ggiraph_0.9.2 dendsort_0.3.4
## [22] cowplot_1.1.3 ggdendro_0.2.0 reshape2_1.4.4
## [25] readr_2.1.5 stringr_1.5.1 ggplot2_4.0.0
## [28] dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.4 matrixStats_1.3.0
## [3] DirichletMultinomial_1.46.0 lubridate_1.9.3
## [5] httr_1.4.7 RColorBrewer_1.1-3
## [7] doParallel_1.0.17 dynamicTreeCut_1.63-1
## [9] tools_4.4.2 backports_1.5.0
## [11] utf8_1.2.4 R6_2.5.1
## [13] lazyeval_0.2.2 mgcv_1.9-1
## [15] rhdf5filters_1.16.0 withr_3.0.1
## [17] gridExtra_2.3 preprocessCore_1.66.0
## [19] fdrtool_1.2.18 qgraph_1.9.8
## [21] WGCNA_1.73 cli_3.6.3
## [23] Biobase_2.64.0 biglm_0.9-3
## [25] sass_0.4.9 robustbase_0.99-3
## [27] mvtnorm_1.2-6 S7_0.2.0
## [29] pbapply_1.7-2 pbivnorm_0.6.0
## [31] systemfonts_1.3.1 yulab.utils_0.2.1
## [33] SPRING_1.0.4 foreign_0.8-86
## [35] scater_1.32.1 parallelly_1.38.0
## [37] decontam_1.24.0 filematrix_1.3
## [39] huge_1.3.5 VGAM_1.1-13
## [41] rstudioapi_0.16.0 impute_1.78.0
## [43] RSQLite_2.3.9 generics_0.1.3
## [45] shape_1.4.6.1 gtools_3.9.5
## [47] GO.db_3.19.1 biomformat_1.32.0
## [49] ggbeeswarm_0.7.2 fansi_1.0.6
## [51] S4Vectors_0.42.1 DECIPHER_3.0.0
## [53] abind_1.4-8 lifecycle_1.0.4
## [55] yaml_2.3.10 SummarizedExperiment_1.34.0
## [57] recipes_1.1.0 rhdf5_2.48.0
## [59] SparseArray_1.4.8 grid_4.4.2
## [61] lavaan_0.6-19 blob_1.2.4
## [63] crayon_1.5.3 beachmat_2.20.0
## [65] KEGGREST_1.44.1 pillar_1.9.0
## [67] knitr_1.48 optparse_1.7.5
## [69] GenomicRanges_1.56.2 pulsar_0.3.11
## [71] corpcor_1.6.10 future.apply_1.11.2
## [73] codetools_0.2-19 fontLiberation_0.1.0
## [75] data.table_1.15.4 MultiAssayExperiment_1.30.3
## [77] vctrs_0.6.5 png_0.1-8
## [79] treeio_1.31.0 Rdpack_2.6.1
## [81] cellranger_1.1.0 gtable_0.3.6
## [83] cachem_1.1.0 gower_1.0.1
## [85] xfun_0.51 prodlim_2024.06.25
## [87] rbibutils_2.3 S4Arrays_1.4.1
## [89] pcaPP_2.0-5 survival_3.7-0
## [91] timeDate_4041.110 SingleCellExperiment_1.26.0
## [93] iterators_1.0.14 hardhat_1.4.0
## [95] lava_1.8.0 bluster_1.14.0
## [97] ipred_0.9-15 nlme_3.1-165
## [99] bit64_4.0.5 fontquiver_0.2.1
## [101] GenomeInfoDb_1.40.1 bslib_0.8.0
## [103] irlba_2.3.5.1 vipor_0.4.7
## [105] rpart_4.1.23 mixedCCA_1.6.2
## [107] colorspace_2.1-1 BiocGenerics_0.50.0
## [109] DBI_1.2.3 Hmisc_5.2-0
## [111] nnet_7.3-19 ade4_1.7-22
## [113] mnormt_2.1.1 tidyselect_1.2.1
## [115] bit_4.0.5 compiler_4.4.2
## [117] htmlTable_2.4.3 BiocNeighbors_1.22.0
## [119] fontBitstreamVera_0.1.1 DelayedArray_0.30.1
## [121] checkmate_2.3.2 scales_1.4.0
## [123] DEoptimR_1.1-3 psych_2.4.12
## [125] quadprog_1.5-8 rappdirs_0.3.3
## [127] digest_0.6.36 rmarkdown_2.27
## [129] XVector_0.44.0 htmltools_0.5.8.1
## [131] pkgconfig_2.0.3 jpeg_0.1-10
## [133] base64enc_0.1-3 sparseMatrixStats_1.16.0
## [135] orca_1.1-3 MatrixGenerics_1.16.0
## [137] fastmap_1.2.0 rlang_1.1.6
## [139] htmlwidgets_1.6.4 UCSC.utils_1.0.0
## [141] DelayedMatrixStats_1.26.0 farver_2.1.2
## [143] jquerylib_0.1.4 jsonlite_1.8.8
## [145] BiocParallel_1.38.0 ModelMetrics_1.2.2.2
## [147] BiocSingular_1.20.0 magrittr_2.0.3
## [149] Formula_1.2-5 scuttle_1.14.0
## [151] GenomeInfoDbData_1.2.12 Rhdf5lib_1.26.0
## [153] Rcpp_1.1.0 ape_5.8
## [155] viridis_0.6.5 gdtools_0.4.4
## [157] stringi_1.8.4 rootSolve_1.8.2.4
## [159] zlibbioc_1.50.0 MASS_7.3-61
## [161] plyr_1.8.9 listenv_0.9.1
## [163] parallel_4.4.2 doSNOW_1.0.20
## [165] Biostrings_2.72.1 splines_4.4.2
## [167] multtest_2.60.0 hms_1.1.3
## [169] fastcluster_1.2.6 stats4_4.4.2
## [171] ScaledMatrix_1.12.0 evaluate_0.24.0
## [173] tzdb_0.4.0 foreach_1.5.2
## [175] getopt_1.20.4 tidyr_1.3.1
## [177] purrr_1.0.2 future_1.34.0
## [179] rsvd_1.0.5 tidytree_0.4.6
## [181] class_7.3-22 glasso_1.11
## [183] viridisLite_0.4.2 snow_0.4-4
## [185] tibble_3.2.1 memoise_2.0.1
## [187] beeswarm_0.4.0 AnnotationDbi_1.66.0
## [189] IRanges_2.38.1 cluster_2.1.6
## [191] corrplot_0.94 TreeSummarizedExperiment_2.12.0
## [193] timechange_0.3.0 globals_0.16.3
## [195] mia_1.12.04. Pathogen screening
Pathogen screening identifies M. plutonius as the primary
disease-associated signal
Targeted pathogen screening was used
to identify candidate etiological agents associated with diseased queen
pupae. These analyses establish M. plutonius as the strongest
disease-associated signal while also revealing its relationship with
BQCV and other potential explanatory variables.
91 #=================================================================================
92 #setwd("/home/brendan/Documents/scripts/brendan/000___GITHUB/2026_09_09___Case_of_the_bad_queen_cell")
93 dir.create("figures", showWarnings=FALSE)
94 #=================================================================================
95
96 pathogen.df <- readr::read_tsv("data/qPCR_pathogen_data.txt") %>% #filter(!is.na(BQCV)) %>%
97 #Set variables of interest as factors
98 mutate_at(vars(disease_status, breeder, builder), as.factor) %>%
99 #Calculate scaled values for downstream analyses
100 mutate(Mplut_scaled = scale(Mplut, center=TRUE, scale=TRUE)[,1]) %>%
101 mutate(Plarv_scaled = scale(Plarv, center=TRUE, scale=TRUE)[,1]) %>%
102 mutate(Crithidia_Lotmaria_scaled = scale(Crithidia_Lotmaria, center=TRUE, scale=TRUE)[,1]) %>%
103 mutate(Napis_scaled = scale(Napis, center=TRUE, scale=TRUE)[,1]) %>%
104 mutate(Ncera_scaled = scale(Ncera, center=TRUE, scale=TRUE)[,1]) %>%
105 mutate(ABPV_scaled = scale(ABPV, center=TRUE, scale=TRUE)[,1]) %>%
106 mutate(BQCV_scaled = scale(BQCV, center=TRUE, scale=TRUE)[,1]) %>%
107 mutate(DWVA_scaled = scale(DWVA, center=TRUE, scale=TRUE)[,1]) %>%
108 mutate(IAPV_scaled = scale(IAPV, center=TRUE, scale=TRUE)[,1]) %>%
109 mutate(KBV_scaled = scale(KBV, center=TRUE, scale=TRUE)[,1]) %>%
110 mutate(LSV1_scaled = scale(LSV1, center=TRUE, scale=TRUE)[,1]) %>%
111 mutate(SBV_scaled = scale(SBV, center=TRUE, scale=TRUE)[,1]) %>%
112 mutate(DWVB_scaled = scale(DWVB, center=TRUE, scale=TRUE)[,1])
113
114 pathogen.df.melt <- melt(pathogen.df, id.vars=c("sample_id", "disease_status"))
115
116 #---------------------------------------------
117 # 1) Pathogen panel
118 #---------------------------------------------
119 pathogen.df.melt.all <- pathogen.df.melt %>%
120 #Filter to only pathogen abundance values
121 filter(grepl("ABPV$|BQCV$|DWVA$|IAPV$|KBV$|LSV1$|SBV$|DWVB$|Mplut$|Plarv$|Napis$|Ncera$|Crithidia_Lotmaria$",
122 variable)) %>%
123 mutate(variable=factor(variable, levels=c("ABPV",
124 "BQCV",
125 "DWVA",
126 "IAPV",
127 "KBV",
128 "LSV1",
129 "SBV",
130 "DWVB",
131 "Mplut",
132 "Plarv",
133 "Napis",
134 "Ncera",
135 "Crithidia_Lotmaria"))) %>%
136 mutate(disease_status=factor(disease_status,
137 levels=c("a_Healthy","b_Diseased"))) %>%
138 mutate(value=as.numeric(value)) %>% mutate(var1 = variable)
139
140 ### Note: Filter out IAPV and Crithidia_Lotmaria for statistical comparisons as these pathogens were not detected in any samples
141
142 df.anova <- pathogen.df.melt.all %>%
143 filter(!grepl("IAPV|Crithidia_Lotmaria", variable)) %>%
144 droplevels()
145
146 #Check normality
147 fit <- lm(value ~ var1 * disease_status, data = df.anova)
148 qq_plot <- ggplotify::as.ggplot(~{
149 qqnorm(residuals(fit))
150 qqline(residuals(fit))
151 })
152 qq_hist <- ggplotify::as.ggplot(~{hist(residuals(fit), breaks = 30, main = "Residual Distribution", xlab = "Residuals")})
153 cowplot::plot_grid(qq_plot, qq_hist)
154 # Two-way ANOVA
155 anova.res <- df.anova %>%
156 rstatix::anova_test(
157 dv = value,
158 between = c(var1, disease_status),
159 type = 2)
160
161 print(anova.res)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 var1 10 308 148.161 2.42e-111 * 0.828
## 2 disease_status 1 308 58.873 2.23e-13 * 0.160
## 3 var1:disease_status 10 308 14.927 1.11e-21 * 0.326162 #Posthoc comparisons with BH (stratified comparisons)
163 pwc1.all <- df.anova %>%
164 group_by(var1) %>%
165 rstatix::emmeans_test(value ~ disease_status,
166 p.adjust.method = "BH",
167 detailed = TRUE ) %>%
168 ungroup() %>%
169 rstatix::adjust_pvalue(p.col = "p",
170 output.col = "p.adj",
171 method = "BH") %>%
172 #rstatix::add_significance("p.adj") %>%
173 rstatix::add_xy_position() %>%
174 #Adjust for plotting on graph below
175 mutate(xmin = (as.numeric(factor(var1, levels = levels(df.anova$var1))))-0.25, xmax = xmin+0.5) %>%
176 #Format to scientific formula for p-values
177 mutate(across(c(p, p.adj), ~ sprintf("%.2e", .x)))
178
179 pwc1.all## # A tibble: 11 × 19
## var1 term .y. group1 group2 null.value estimate se df conf.low
## <fct> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ABPV disease_s… value a_Hea… b_Dis… 0 0.00625 0.301 308 -0.585
## 2 BQCV disease_s… value a_Hea… b_Dis… 0 -3.11 0.301 308 -3.70
## 3 DWVA disease_s… value a_Hea… b_Dis… 0 -0.103 0.301 308 -0.695
## 4 KBV disease_s… value a_Hea… b_Dis… 0 -0.162 0.301 308 -0.753
## 5 LSV1 disease_s… value a_Hea… b_Dis… 0 -0.153 0.301 308 -0.745
## 6 SBV disease_s… value a_Hea… b_Dis… 0 -0.347 0.301 308 -0.939
## 7 DWVB disease_s… value a_Hea… b_Dis… 0 -0.415 0.301 308 -1.01
## 8 Mplut disease_s… value a_Hea… b_Dis… 0 -2.96 0.301 308 -3.55
## 9 Plarv disease_s… value a_Hea… b_Dis… 0 -0.135 0.301 308 -0.727
## 10 Napis disease_s… value a_Hea… b_Dis… 0 -0.228 0.301 308 -0.820
## 11 Ncera disease_s… value a_Hea… b_Dis… 0 -0.0514 0.301 308 -0.643
## # ℹ 9 more variables: conf.high <dbl>, statistic <dbl>, p <chr>, p.adj <chr>,
## # p.adj.signif <chr>, y.position <dbl>, groups <named list>, xmin <dbl>,
## # xmax <dbl>Figure 1B
180 fig1B <- ggplot(df.anova, aes(x=var1, y=value)) +
181 geom_bar(aes(fill=disease_status), stat="summary", fun="mean", width=0.75, alpha=1,
182 position = position_dodge(width = 0.75)) +
183 geom_errorbar(data= . %>% dplyr::group_by(var1, disease_status) %>%
184 dplyr::summarise(mean = mean(value, na.rm = TRUE),
185 se = sd(value, na.rm = TRUE) / sqrt(n()), .groups = 'drop'),
186 aes(x = var1, ymin = mean - se, ymax = mean + se, group = disease_status),
187 width = 0.2, linewidth=0.5, position = position_dodge(width = 0.75), inherit.aes = FALSE) +
188 geom_point(aes(fill=disease_status),
189 size=0.5, alpha=0.5,
190 position = position_jitterdodge(jitter.width=0.3, jitter.height = 0, dodge.width = 0.75)) +
191 ggpubr::stat_pvalue_manual(pwc1.all, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
192 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
193 scale_fill_manual(values=c("a_Healthy" = "#bbf0bb", "b_Diseased" = "#edc2c2")) +
194 theme_classic(base_size=10) +
195 labs(y="Log10 genome equivalents/bee") +
196 ylim(c(0,12)) +
197 labs(caption = rstatix::get_pwc_label(pwc1.all)) +
198 annotate("text", x = Inf, y = 2, label = "LOD", hjust = 0.2,
199 vjust = 0.5, size = 3, colour="grey50") +
200 geom_vline(xintercept=c(7.5,9.5), linetype="solid", linewidth=0.1, colour="black") +
201 geom_hline(yintercept=1.99, linetype="dotted", linewidth=0.2, colour="black") +
202 theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1),
203 axis.title.x = element_blank(),
204 legend.position="none",
205 panel.grid=element_blank())
206
207 fig1B
Figure 1C
208 #---------------------------------------------
209 # Pearson correlation plot
210 #---------------------------------------------
211 fig1C <-ggplot(pathogen.df, aes(x=BQCV_scaled, y=Mplut_scaled)) +
212 geom_jitter(aes(colour=disease_status, fill=disease_status), size=3, width=0.1, height=0.1, stroke=1, shape=21) +
213 smplot2::sm_statCorr(color = 'grey50', corr_method = 'pearson', linetype="dashed", label_x = 0, label_y = 2.5, text_size=3) +
214 scale_fill_manual(values=c("a_Healthy" = scales::alpha("#bbf0bb", 0.5), "b_Diseased" = scales::alpha("#edc2c2", 0.3))) +
215 scale_colour_manual(values=c("a_Healthy" = "#7a9e7a", "b_Diseased" = "#ad6a6a")) + theme_classic(base_size=10) +
216 xlab("BCQV") + ylab("M. plutonius") +
217 theme(panel.grid=element_blank(), legend.position="none")
218 fig1C
219 print(cor.test(pathogen.df$BQCV, pathogen.df$Mplut, method = "pearson"))##
## Pearson's product-moment correlation
##
## data: pathogen.df$BQCV and pathogen.df$Mplut
## t = 5.4191, df = 28, p-value = 8.83e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4785684 0.8552857
## sample estimates:
## cor
## 0.7154793Machine learning models
220 #---------------------------------------------
221 # ROC curve
222 #---------------------------------------------
223 roc.df <- pathogen.df %>%
224 mutate(disease_status=ifelse(disease_status=="a_Healthy", "Healthy", "Diseased")) %>%
225 mutate(disease_status = factor(disease_status))
226
227 set.seed(3)
228
229 # Define number of folds
230 k_folds <- 10
231
232 # Create k-fold cross-validation folds
233 folds <- createFolds(roc.df$disease_status, k = k_folds, list = TRUE, returnTrain = FALSE)
234
235 # Initialize vectors to store predicted probabilities and true labels
236 rf_probs_all <- rep(NA, nrow(roc.df))
237 lr_probs_all <- rep(NA, nrow(roc.df))
238 true_labels <- roc.df$disease_status
239
240 # Cross-validation loop
241 for (i in 1:k_folds) {
242 test_idx <- folds[[i]]
243 train_idx <- setdiff(seq_len(nrow(roc.df)), test_idx)
244
245 train_data <- roc.df[train_idx, ]
246 test_data <- roc.df[test_idx, ]
247
248 # ---- Random Forest ----
249 rf_model <- randomForest(disease_status ~ BQCV_scaled + ABPV_scaled + DWVA_scaled + KBV_scaled +
250 LSV1_scaled + SBV_scaled + DWVB_scaled + Napis_scaled + Ncera_scaled +
251 Mplut_scaled + Plarv_scaled + breeder + builder,
252 data = train_data, ntree = 1000)
253 rf_probs_all[test_idx] <- predict(rf_model, newdata = test_data, type = "prob")[, 2]
254
255 # ---- Lasso Logistic Regression ----
256 X_train <- model.matrix(disease_status ~ BQCV_scaled + ABPV_scaled + DWVA_scaled + KBV_scaled +
257 LSV1_scaled + SBV_scaled + DWVB_scaled + Napis_scaled + Ncera_scaled +
258 Mplut_scaled + Plarv_scaled + breeder + builder, data = train_data)[, -1]
259 X_test <- model.matrix(disease_status ~ BQCV_scaled + ABPV_scaled + DWVA_scaled + KBV_scaled +
260 LSV1_scaled + SBV_scaled + DWVB_scaled + Napis_scaled + Ncera_scaled +
261 Mplut_scaled + Plarv_scaled + breeder + builder, data = test_data)[, -1]
262 y_train <- train_data$disease_status
263
264 lasso_model <- cv.glmnet(X_train, y_train, family = "binomial", alpha = 1, grouped=FALSE)
265 best_lambda <- lasso_model$lambda.min
266
267 lr_probs_all[test_idx] <- predict(lasso_model, newx = X_test, s = best_lambda, type = "response")[, 1]
268 }
269
270 # ---- Evaluate ROC and AUC ----
271 rf_roc <- roc(true_labels, rf_probs_all)
272 lr_roc <- roc(true_labels, lr_probs_all)
273 rf_roc##
## Call:
## roc.default(response = true_labels, predictor = rf_probs_all)
##
## Data: rf_probs_all in 17 controls (true_labels Diseased) < 13 cases (true_labels Healthy).
## Area under the curve: 0.9819274 # Extract ROC data for ggplot
275 rf_df <- data.frame(FPR = 1 - rf_roc$specificities, TPR = rf_roc$sensitivities, Model = "Random Forest")
276 lr_df <- data.frame(FPR = 1 - lr_roc$specificities, TPR = lr_roc$sensitivities, Model = "Lasso Logistic Regression")
277 rf_roc$specificities## [1] 0.00000000 0.05882353 0.11764706 0.17647059 0.23529412 0.29411765
## [7] 0.35294118 0.41176471 0.47058824 0.52941176 0.58823529 0.64705882
## [13] 0.70588235 0.76470588 0.82352941 0.82352941 0.88235294 0.94117647
## [19] 0.94117647 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [25] 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
## [31] 1.00000000278 # Combine for ggplot
279 roc_df <- rbind(rf_df, lr_df) %>% arrange(TPR, .by_group = TRUE) %>% arrange(FPR, .by_group = TRUE)
280
281 auc_rf <- auc(rf_roc)
282 auc_lr <- auc(lr_roc)
283 roc_df## FPR TPR Model
## 1 0.00000000 0.00000000 Random Forest
## 2 0.00000000 0.00000000 Lasso Logistic Regression
## 3 0.00000000 0.07692308 Random Forest
## 4 0.00000000 0.07692308 Lasso Logistic Regression
## 5 0.00000000 0.15384615 Random Forest
## 6 0.00000000 0.15384615 Lasso Logistic Regression
## 7 0.00000000 0.23076923 Random Forest
## 8 0.00000000 0.23076923 Lasso Logistic Regression
## 9 0.00000000 0.30769231 Random Forest
## 10 0.00000000 0.30769231 Lasso Logistic Regression
## 11 0.00000000 0.38461538 Random Forest
## 12 0.00000000 0.38461538 Lasso Logistic Regression
## 13 0.00000000 0.46153846 Random Forest
## 14 0.00000000 0.46153846 Lasso Logistic Regression
## 15 0.00000000 0.53846154 Random Forest
## 16 0.00000000 0.53846154 Lasso Logistic Regression
## 17 0.00000000 0.61538462 Random Forest
## 18 0.00000000 0.61538462 Lasso Logistic Regression
## 19 0.00000000 0.69230769 Random Forest
## 20 0.00000000 0.69230769 Lasso Logistic Regression
## 21 0.00000000 0.76923077 Random Forest
## 22 0.00000000 0.76923077 Lasso Logistic Regression
## 23 0.00000000 0.84615385 Random Forest
## 24 0.00000000 0.84615385 Lasso Logistic Regression
## 25 0.05882353 0.84615385 Random Forest
## 26 0.05882353 0.84615385 Lasso Logistic Regression
## 27 0.05882353 0.92307692 Random Forest
## 28 0.11764706 0.84615385 Lasso Logistic Regression
## 29 0.11764706 0.92307692 Random Forest
## 30 0.11764706 0.92307692 Lasso Logistic Regression
## 31 0.17647059 0.92307692 Random Forest
## 32 0.17647059 0.92307692 Lasso Logistic Regression
## 33 0.17647059 1.00000000 Random Forest
## 34 0.23529412 0.92307692 Lasso Logistic Regression
## 35 0.23529412 1.00000000 Random Forest
## 36 0.29411765 0.92307692 Lasso Logistic Regression
## 37 0.29411765 1.00000000 Random Forest
## 38 0.35294118 0.92307692 Lasso Logistic Regression
## 39 0.35294118 1.00000000 Random Forest
## 40 0.41176471 0.92307692 Lasso Logistic Regression
## 41 0.41176471 1.00000000 Random Forest
## 42 0.41176471 1.00000000 Lasso Logistic Regression
## 43 0.47058824 1.00000000 Random Forest
## 44 0.47058824 1.00000000 Lasso Logistic Regression
## 45 0.52941176 1.00000000 Random Forest
## 46 0.52941176 1.00000000 Lasso Logistic Regression
## 47 0.58823529 1.00000000 Random Forest
## 48 0.58823529 1.00000000 Lasso Logistic Regression
## 49 0.64705882 1.00000000 Random Forest
## 50 0.64705882 1.00000000 Lasso Logistic Regression
## 51 0.70588235 1.00000000 Random Forest
## 52 0.70588235 1.00000000 Lasso Logistic Regression
## 53 0.76470588 1.00000000 Random Forest
## 54 0.76470588 1.00000000 Lasso Logistic Regression
## 55 0.82352941 1.00000000 Random Forest
## 56 0.82352941 1.00000000 Lasso Logistic Regression
## 57 0.88235294 1.00000000 Random Forest
## 58 0.88235294 1.00000000 Lasso Logistic Regression
## 59 0.94117647 1.00000000 Random Forest
## 60 0.94117647 1.00000000 Lasso Logistic Regression
## 61 1.00000000 1.00000000 Random Forest
## 62 1.00000000 1.00000000 Lasso Logistic RegressionFigure 1D
284 fig1D <- ggplot(roc_df %>% mutate(FPR = ifelse(Model == "Lasso Logistic Regression", FPR+0.008, FPR)) %>%
285 mutate(TPR = ifelse(Model == "Lasso Logistic Regression", TPR+0.008, TPR)) %>%
286 mutate(Model = factor(Model, levels=unique(c("Random Forest", "Lasso Logistic Regression")))),
287 aes(x = FPR, y = TPR, color = Model)) +
288 geom_line(size = 1.2, alpha=1) +
289 theme_minimal(base_size=10) +
290 ylab("TPR") + xlab("FPR-1") +
291 scale_color_manual(values = c("Random Forest" = "steelblue",
292 "Lasso Logistic Regression" = "gold"),
293 labels = c(paste0("Random Forest\n(AUC = ", round(auc_rf, 3), ")"),
294 paste0("Lasso Logistic\n(AUC = ", round(auc_lr, 3), ")"))) +
295 geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
296 #theme_classic() +
297 theme(
298 panel.border = element_rect(color="black", fill=NA),
299 panel.grid = element_blank(),
300 legend.position = c(0.95, 0.10), # x, y coordinates (0 = left/bottom, 1 = right/top)
301 legend.justification = c("right", "bottom"),
302 legend.text=element_text(size=7),
303 legend.title=element_text(size=8),
304 legend.background = element_rect(fill = "white", color = "gray80", linewidth = 0.3),
305 legend.box.background = element_rect(color = "black")
306 )
307
308 fig1D
Predictor importance
309 #---------------------------------------------
310 # Variable importance model
311 #---------------------------------------------
312 set.seed(3)
313
314 rf_model <- randomForest(disease_status ~ BQCV_scaled + ABPV_scaled + DWVA_scaled + KBV_scaled +
315 LSV1_scaled + SBV_scaled + DWVB_scaled + Napis_scaled + Ncera_scaled +
316 Mplut_scaled + Plarv_scaled + breeder + builder, data = train_data, ntree = 1000)
317
318 n_boot <- 100
319 importance_matrix1 <- matrix(NA, nrow = n_boot, ncol = length(rf_model$importance[, 1]))
320 importance_matrix2 <- matrix(NA, nrow = n_boot, ncol = length(rf_model$importance[, 1]))
321 importance_matrix3 <- matrix(NA, nrow = n_boot, ncol = length(rf_model$importance[, 1]))
322
323 colnames(importance_matrix1) <- sort(rownames(rf_model$importance))
324 colnames(importance_matrix2) <- sort(rownames(rf_model$importance))
325 colnames(importance_matrix3) <- sort(rownames(rf_model$importance))
326
327 for (i in 1:n_boot) {
328 boot_idx <- sample(1:nrow(roc.df), replace = TRUE)
329 boot_data <- roc.df[boot_idx, ]
330 rf_boot <- randomForest(
331 disease_status ~ BQCV_scaled + ABPV_scaled + DWVA_scaled + KBV_scaled +
332 LSV1_scaled + SBV_scaled + DWVB_scaled + Napis_scaled + Ncera_scaled +
333 Mplut_scaled + Plarv_scaled + breeder + builder,
334 data = boot_data,
335 importance = TRUE,
336 ntree = 1000
337 )
338 #------
339 importance.df.tmp <- randomForest::importance(rf_boot) %>% as.data.frame() %>% arrange(factor(rownames(.)))
340 importance_matrix1[i, ] <- importance.df.tmp[, "MeanDecreaseAccuracy"]
341 importance_matrix2[i, ] <- importance.df.tmp[, "MeanDecreaseGini"]
342 #------
343 }
344
345 # Compute 95% CI for each variable and combine into a data frame
346 importance_ci_df <- data.frame(
347 measure = c(rep("MeanDecreaseAccuracy", length(colnames(importance_matrix1))), rep("MeanDecreaseGini", length(colnames(importance_matrix2)))),
348 Variable = c(colnames(importance_matrix1), colnames(importance_matrix2)),
349 Mean = c(colMeans(importance_matrix1), colMeans(importance_matrix2)),
350 SD = c(apply(importance_matrix1, 2, sd), apply(importance_matrix2, 2, sd)),
351 CI_lower = c(t(apply(importance_matrix1, 2, function(x) quantile(x, probs = c(0.025, 0.975))))[, 1], t(apply(importance_matrix2, 2, function(x) quantile(x, probs = c(0.025, 0.975))))[, 1]),
352 CI_upper = c(t(apply(importance_matrix1, 2, function(x) quantile(x, probs = c(0.025, 0.975))))[, 2], t(apply(importance_matrix2, 2, function(x) quantile(x, probs = c(0.025, 0.975))))[, 2])) %>% arrange(Mean) %>%
353 mutate(Variable = factor(Variable, levels=unique(Variable))) %>% mutate(Variable = gsub("_scaled", "", Variable))
354 #----------------------------Plots
355 print(importance_ci_df %>% arrange(measure))## measure Variable Mean SD CI_lower CI_upper
## 1 MeanDecreaseAccuracy ABPV 1.14238296 1.57051676 -1.616809092 4.2602330
## 2 MeanDecreaseAccuracy KBV 3.47400937 1.75113153 0.000000000 6.4718860
## 3 MeanDecreaseAccuracy Ncera 3.49247723 2.41473296 -1.102045175 8.3498972
## 4 MeanDecreaseAccuracy DWVA 4.25388563 1.53611889 1.192487309 7.0151773
## 5 MeanDecreaseAccuracy LSV1 4.27690687 2.25631476 -0.007431712 8.3147952
## 6 MeanDecreaseAccuracy SBV 4.94779696 1.92705236 1.757144107 8.6059124
## 7 MeanDecreaseAccuracy Plarv 5.43314248 2.12927715 0.974078796 8.9930927
## 8 MeanDecreaseAccuracy DWVB 5.49469579 2.26273893 2.024002737 10.1008990
## 9 MeanDecreaseAccuracy Napis 6.72719089 2.12815061 2.279841130 10.2344228
## 10 MeanDecreaseAccuracy breeder 8.67713499 3.14122914 2.309612238 15.2462876
## 11 MeanDecreaseAccuracy builder 9.96082259 3.23738370 4.087658721 15.4986756
## 12 MeanDecreaseAccuracy BQCV 18.01569677 2.50974710 13.152280654 22.4367543
## 13 MeanDecreaseAccuracy Mplut 24.48538149 2.01218251 20.431877842 27.7826520
## 14 MeanDecreaseGini ABPV 0.05035977 0.05020640 0.002172080 0.1818321
## 15 MeanDecreaseGini KBV 0.12729563 0.08294934 0.018229399 0.3475086
## 16 MeanDecreaseGini Ncera 0.23153378 0.11192730 0.042543860 0.4478525
## 17 MeanDecreaseGini DWVA 0.28545591 0.11643120 0.122642633 0.5290488
## 18 MeanDecreaseGini LSV1 0.29240834 0.17767174 0.111147284 0.5956527
## 19 MeanDecreaseGini SBV 0.34851666 0.15936021 0.116299233 0.6524892
## 20 MeanDecreaseGini DWVB 0.43592082 0.24162010 0.177619071 1.2456566
## 21 MeanDecreaseGini Plarv 0.46869322 0.20254659 0.209562376 0.9978756
## 22 MeanDecreaseGini Napis 0.62406089 0.23539972 0.302746602 1.1607892
## 23 MeanDecreaseGini breeder 1.69510317 0.50805360 1.006328763 3.0258485
## 24 MeanDecreaseGini builder 1.89226762 0.59520243 0.945118538 3.1065781
## 25 MeanDecreaseGini BQCV 2.73735599 0.66659305 1.505622791 3.9913831
## 26 MeanDecreaseGini Mplut 4.48171606 0.48423923 3.475197822 5.4249952Figure 1E
356 fig1E <- ggplot(importance_ci_df %>% filter(measure=="MeanDecreaseAccuracy") %>%
357 arrange(measure) %>% mutate(Variable = factor(Variable, levels=unique(Variable))),
358 aes(x = Variable, y = Mean)) +
359 geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0) + geom_point(size = 2, shape=21, fill="white") +
360 coord_flip() + theme_minimal(base_size = 10) + labs(x = "Variable", y = "%IncMSE (Mean Decrease Accuracy)") +
361 theme(panel.border=element_rect(colour="black", fill=NA),
362 panel.grid = element_blank(),
363 axis.title.y=element_blank())
364 fig1E
Figure 1 (Full panel)
365 #---------------------------------------------
366 # 5) Combine all plots
367 #---------------------------------------------
368
369 plot.mid <- cowplot::plot_grid(fig1B, labels=c("B"), label_size=10)
370 plot.bottom <- cowplot::plot_grid(fig1C, fig1D, fig1E, ncol=3, labels=c("C", "D", "E"), label_size=10)
371
372 fig1.final <- cowplot::plot_grid(plot.mid, plot.bottom, ncol=1, rel_heights=c(0.555, 0.45))
373 fig1.final
5. Microbiota profiling
Microbiota profiling reveals disease-associated community
restructuring
Having identified the major pathogen-associated
signals, we next asked whether disease was accompanied by broader
disruption of the queen pupal microbiota. 16S rRNA gene profiling was
used to examine community composition, diversity, bacterial abundance,
M. plutonius lineage structure, and the taxa that distinguish
healthy from diseased individuals.
Here, we load the necessary libraries and import the ASV
table derived from 16S rRNA gene sequencing
374 library(dplyr)
375 library(ggplot2)
376 library(stringr)
377 library(readr)
378 library(eoffice)
379 library(zCompositions)
380 library(compositions)
381 library(ggbeeswarm)
382 library(reshape2)
383 library(heatmap3)
384 library(phyloseq)
385 library(vegan)
386 library(ALDEx2)
387 library(ggdendro)
388 library(cowplot)
389 library(dendsort)
390 library(dendextend)
391 library(cluster)
392 library(factoextra)
393 library(ggtree)
394 library(ggtreeExtra)
395 library(WdStar)
396
397 #---------------------------------------------
398 # Import data tables
399 #---------------------------------------------
400 d.b <- read.table("data/V3V4_ASV_count_table.txt",
401 sep="\t", quote="", check.names=F, header=T, row.names=1, comment.char="") %>%
402 mutate(tax.vector= str_split_fixed(rownames(.), ";", 2)[,2]) %>%
403 `rownames<-`(str_split_fixed(rownames(.), ";", 2)[,1]) %>%
404 `colnames<-`(str_split_fixed(colnames(.), "_S", 2)[,1]) %>%
405 filter(grepl("Bacteria;", .$tax.vector)) %>%
406 `rownames<-`(paste(rownames(.), .$tax.vector, sep=";")) %>%
407 dplyr::select(!tax.vector) %>%
408 {. ->> tmp} %>%
409 cbind(.,
410 str_split_fixed(rownames(.), ";", 8)) %>%
411 `colnames<-`(c(colnames(tmp), "SV", "domain", "phylum", "class", "order", "family", "genus", "species")) %>%
412 `rownames<-`(paste(.$SV, .$domain, .$phylum, .$class, .$order, .$family, .$genus, .$species, sep=";")) %>% dplyr::select(-SV,-domain,-phylum,-class,-order,-family,-genus,-species) %>%
413 {. ->> d.b.unfilt}
414
415 d.b.species <- d.b %>% cbind(., str_split_fixed(rownames(.), ";", 8)[,2:8] %>%
416 `colnames<-`(c("domain", "phylum", "class", "order", "family", "genus", "species"))) %>%
417 group_by(species) %>% mutate_at(vars(!matches("domain|phylum|class|order|family|genus|species")), ~sum(.)) %>%
418 dplyr::slice(1) %>% ungroup() %>% tibble::column_to_rownames("species") %>%
419 dplyr::select(!matches("domain|phylum|class|order|family|genus|species"))
420
421 d.b.genus <- d.b %>% cbind(., str_split_fixed(rownames(.), ";", 8)[,2:8] %>%
422 `colnames<-`(c("domain", "phylum", "class", "order", "family", "genus", "species"))) %>%
423 group_by(genus) %>% mutate_at(vars(!matches("domain|phylum|class|order|family|genus|species")), ~sum(.)) %>%
424 dplyr::slice(1) %>% ungroup() %>% tibble::column_to_rownames("genus") %>%
425 dplyr::select(!matches("domain|phylum|class|order|family|genus|species"))
426
427 #---------------------------------------------
428 # Import metadata
429 #---------------------------------------------
430 metadata <- read.table("data/metadata.txt", sep="\t",
431 quote="",
432 check.names=F,
433 header=T,
434 row.names=1,
435 comment.char="") %>%
436 tibble::rownames_to_column("sample_id")
437
438
439 ################################
440 ################################ RAW (for phyloseq)
441 ################################
442 #================================
443 #Null cutoff for raw data
444 filt.cutoff.raw= 0
445 prev.cutoff.raw=0
446 #================================
447
448 d.b.filt.raw <- d.b %>%
449 tibble::rownames_to_column() %>% filter(!grepl("Chloroplast|Mitochondria",rowname)) %>%
450 tibble::column_to_rownames()
451
452 d.b.clr.raw <- d.b.filt.raw %>% select_if(colSums(.) != 0) %>%
453 #FILTER BY ABUNDANCE CUTOFF
454 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > filt.cutoff.raw) %>%
455 mutate_all(~.+1) %>% #Simple zero inflation
456 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
457 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>% #CONVERT TO CLR
458 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
459 arrange(sample_id) %>% arrange(disease_status)
460
461 d.b.counts.raw <- d.b.filt.raw %>% select_if(colSums(.) != 0) %>%
462 #FILTER BY ABUNDANCE CUTOFF
463 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > filt.cutoff.raw) %>%
464 as.data.frame() %>%
465 sjmisc::rotate_df() %>%
466 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>% arrange(sample_id) %>% arrange(disease_status)
467
468 d.b.prev.raw <- d.b.filt.raw %>% select_if(colSums(.) != 0) %>%
469 #FILTER BY ABUNDANCE CUTOFF
470 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > filt.cutoff.raw) %>%
471 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
472 as.data.frame() %>%
473 sjmisc::rotate_df() %>%
474 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
475 arrange(sample_id) %>% arrange(disease_status) %>%
476 {d.b.freq.raw <<- .} %>%
477 mutate_at(vars(matches("s__|g__")), ~ifelse(. >= prev.cutoff.raw, 1, 0))
478
479
480 ################################
481 ################################ SPECIES-level
482 ################################
483
484 d.b.filt.species <- d.b.species %>%
485 tibble::rownames_to_column() %>% filter(!grepl("Chloroplast|Mitochondria",rowname)) %>%
486 tibble::column_to_rownames()
487
488 d.b.clr.species <- d.b.filt.species %>% select_if(colSums(.) != 0) %>%
489 #FILTER BY ABUNDANCE CUTOFF
490 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
491 mutate_all(~.+1) %>% #Simple zero inflation
492 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
493 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>% #CONVERT TO CLR
494 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
495 arrange(sample_id) %>% arrange(disease_status)
496
497 d.b.counts.species <- d.b.filt.species %>% select_if(colSums(.) != 0) %>%
498 #FILTER BY ABUNDANCE CUTOFF
499 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
500 as.data.frame() %>%
501 sjmisc::rotate_df() %>%
502 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>% arrange(sample_id) %>% arrange(disease_status)
503
504 d.b.prev.species <- d.b.filt.species %>% select_if(colSums(.) != 0) %>%
505 #FILTER BY ABUNDANCE CUTOFF
506 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
507 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
508 as.data.frame() %>%
509 sjmisc::rotate_df() %>%
510 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
511 arrange(sample_id) %>% arrange(disease_status) %>%
512 {d.b.freq.species <<- .} %>%
513 mutate_at(vars(matches("s__|g__")), ~ifelse(. >= 0.00001, 1, 0))
514
515 ################################
516 ################################ GENUS-level
517 ################################
518
519 d.b.filt.genus <- d.b.genus %>%
520 tibble::rownames_to_column() %>% filter(!grepl("Chloroplast|Mitochondria",rowname)) %>%
521 tibble::column_to_rownames()
522
523 d.b.clr.genus <- d.b.filt.genus %>% select_if(colSums(.) != 0) %>%
524 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
525 #FILTER BY ABUNDANCE CUTOFF
526 mutate_all(~.+1) %>% #Simple zero inflation
527 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
528 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>% #CONVERT TO CLR
529 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
530 arrange(sample_id) %>% arrange(disease_status)
531
532 d.b.counts.genus <- d.b.filt.genus %>% select_if(colSums(.) != 0) %>%
533 #FILTER BY ABUNDANCE CUTOFF
534 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
535 as.data.frame() %>%
536 sjmisc::rotate_df() %>%
537 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>% arrange(sample_id) %>% arrange(disease_status)
538
539 d.b.prev.genus <- d.b.filt.genus %>% select_if(colSums(.) != 0) %>%
540 #FILTER BY ABUNDANCE CUTOFF
541 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
542 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
543 as.data.frame() %>%
544 sjmisc::rotate_df() %>%
545 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
546 arrange(sample_id) %>% arrange(disease_status) %>%
547 {d.b.freq.genus <<- .} %>%
548 mutate_at(vars(matches("s__|g__")), ~ifelse(. >= 0.00001, 1, 0))
549
550
551 ################################
552 ################################ ASV-level
553 ################################
554
555 d.b.filt.asv <- d.b %>%
556 tibble::rownames_to_column() %>% filter(!grepl("Chloroplast|Mitochondria",rowname)) %>%
557 tibble::column_to_rownames()
558
559 d.b.clr.asv <- d.b.filt.asv %>% select_if(colSums(.) != 0) %>%
560 #FILTER BY ABUNDANCE CUTOFF
561 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
562 mutate_all(~.+1) %>% #Simple zero inflation
563 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
564 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>% #CONVERT TO CLR
565 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
566 arrange(sample_id) %>% arrange(disease_status)
567
568 d.b.counts.asv <- d.b.filt.asv %>% select_if(colSums(.) != 0) %>%
569 #FILTER BY ABUNDANCE CUTOFF
570 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
571 as.data.frame() %>%
572 sjmisc::rotate_df() %>%
573 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>% arrange(sample_id) %>% arrange(disease_status)
574
575 d.b.prev.asv <- d.b.filt.asv %>% select_if(colSums(.) != 0) %>%
576 #FILTER BY ABUNDANCE CUTOFF
577 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
578 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
579 as.data.frame() %>%
580 sjmisc::rotate_df() %>%
581 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
582 arrange(sample_id) %>% arrange(disease_status) %>%
583 {d.b.freq.asv <<- .} %>%
584 mutate_at(vars(matches("s__|g__")), ~ifelse(. >= 0.00001, 1, 0))Composition and clustering
Taxonomic profiles and hierarchical clustering provide an overview of microbiota composition across individual queen pupae. Samples organize into distinct community states that broadly reflect disease status while retaining variation associated with breeder and builder origin.
585 ################################
586 ################################ Dendrogram
587 ################################
588 dend.in <- d.b.clr.asv
589
590 #---Absolute
591 d.otu <- dend.in %>% tibble::column_to_rownames("sample_id") %>% dplyr::select(matches("s__|g__")) %>%
592 sjmisc::rotate_df() %>% arrange(desc(rowSums(.))) %>% filter(grepl("s__|g__", rownames(.))) %>% sjmisc::rotate_df()
593
594 d.meta <- dend.in %>% tibble::column_to_rownames("sample_id") %>% dplyr::select(!matches("s__|g__"))
595
596 #Cluster samples and taxa
597 cluster.tax <- hclust(dist(t(d.otu), method="euclidian"), method="ward.D2")
598 cluster.samples <- hclust(dist(d.otu, method="euclidian"), method="ward.D2")
599
600 #Dendrogram for sample and taxa clustering
601 g1.dendro.samples <- ggdendro::ggdendrogram(cluster.samples, labels = TRUE, rotate=FALSE)
602 g2.dendro.tax<- ggdendro::ggdendrogram(cluster.tax, labels = TRUE, rotate=TRUE) + scale_y_reverse()
603
604 #Dendrogram with labels coloured by cluster
605 kcluster.plot <- fviz_nbclust(d.otu, kmeans, method = "silhouette") +
606 ggtitle("Silhouette Method for Optimal k")
607 kcluster.plot
608 # Plot the dendrogram with color-coded labels
609 plot.dendro <- ggplot() +
610 geom_segment(data = dend_data$segments,
611 aes(x = x, y = y, xend = xend, yend = yend)) +
612 geom_text(data = labels_df, aes(x = x, y = y,
613 label = label,
614 color = cluster),
615 hjust = 1, vjust = 0.5,
616 size = 3, angle = 90) + # Rotate labels
617 scale_color_manual(values = c("#DD8B8B", "#00B050")) + # Adjust colors as needed
618 theme_minimal() +
619 theme(axis.text.x.bottom=element_blank(), axis.title.x=element_blank()) +
620 labs(title = "Dendrogram Clustering (Hierarchical-Based)") +
621 coord_cartesian(clip = "off") + # Allow text to be drawn outside the plot area
622 theme(plot.margin = ggplot2::margin(10, 50, 10, 10)) # Expand right margin for labels
623 plot.dendro
624 #::::::::::::::::::::::::::::
625 # ggtree plot version
626 #::::::::::::::::::::::::::::
627 tree.df <- fortify(dend_data) %>%
628 mutate(cluster=dplyr::recode(!!!setNames(labels_df$cluster, labels_df$label),label))
629
630 tree.g <- ggtree(tree.df) +
631 geom_rootedge(rootedge = 0.1) +
632 geom_tippoint(aes(colour = cluster)) +
633 scale_colour_manual(values=c("1" = scales::alpha("#DD8B8B", .9),
634 "2" = scales::alpha("#00B050", .8))) +
635 coord_flip() + scale_x_reverse() + scale_y_reverse() + theme_tree() +
636 ggplot2::theme(legend.position = "none",
637 axis.text.x = element_blank(),
638 axis.title.x = element_blank(),
639 axis.text.y = element_blank(),
640 axis.ticks.x = element_blank())
641
642 #Rotate branches for clearer visualization
643 tree.g.rotate <- tree.g %>% ggtree::rotate(50) %>% ggtree::rotate(56)%>% ggtree::rotate(51)
644 tree.g.rotate
Figure 2A
645 ################################
646 ################################ Bar plot
647 ################################
648 bar.in <- d.b.freq.genus
649
650 d.otu <- bar.in %>% tibble::column_to_rownames("sample_id") %>% dplyr::select(matches("s__|g__")) %>%
651 sjmisc::rotate_df() %>% arrange(desc(rowSums(.))) %>% filter(grepl("s__|g__", rownames(.)))
652
653 d.meta <- bar.in %>% tibble::column_to_rownames("sample_id") %>% dplyr::select(!matches("s__|g__"))
654
655 #Top 29
656 mer.xxx.sub.match <- rownames(d.otu %>% arrange(desc(rowSums(.))) %>% .[1:29,]) #top 30
657 taxa.match <- d.otu %>% filter(., rownames(d.otu) %in% mer.xxx.sub.match) %>% arrange(desc(rowSums(.)))
658 remainder <- 1-colSums(taxa.match, dims=1)
659 p<-rbind(taxa.match,remainder)
660 rownames(p)[nrow(p)]<-"remainder"
661 p.names <- colnames(p)
662 y.bac <- as.matrix(p)
663
664 #--------------------------
665 # ggplot2 Bar plot
666 #--------------------------
667 metadata <- metadata %>% mutate(cluster_group=dplyr::recode(!!!setNames(labels_df$cluster ,labels_df$label), sample_id)) %>%
668 mutate(cluster_order=dplyr::recode(!!!setNames(labels_df$x ,labels_df$label), sample_id)) %>% arrange(cluster_order) %>%
669 mutate(cluster_group=ifelse(cluster_group==1, "Diseased", "Healthy"))
670
671 y.bac.mer <- as.data.frame(y.bac) %>% sjmisc::rotate_df() %>% tibble::rownames_to_column() %>%
672 merge(metadata, ., by.x="sample_id", by.y="rowname", sort=FALSE, all=FALSE) %>%
673 arrange(sample_id) %>% arrange(cluster_group)
674
675 d.melt <- reshape2::melt(y.bac.mer, id.vars=c(colnames(metadata))) %>%
676 mutate(sample_id = factor(sample_id, levels=ggtree::get_taxa_name(tree.g.rotate))) %>%
677 mutate(value_abs = value * qpcr_16S_adj)
678
679 plot.bar <-ggplot(d.melt, aes(x=sample_id, y=value, fill=variable)) +
680 geom_col(position="stack", width=0.9, alpha=1, color=NA, linewidth=0) + scale_y_reverse() +
681 scale_fill_manual(values=c("#7F7F7F", "lightskyblue3", "gold", "darkseagreen1","steelblue",
682 "red2", "pink","darkolivegreen3","#663366","darkorange3",
683 "darkorchid3","#C0C0C0","yellow","#00B9FF","#E8FF10",
684 "#F740B7","#9999FF", "beige", "#00E8FF", "#BEBADA", "#FDB462",
685 "#66C2A5", "#FFFFB3", "#20FF00", "#E41A1C", "#00B9FF",
686 "#A65628","#B900FF", "#E5D8BD", "#2000FF", "#CBD5E8")) +
687 theme(panel.background=element_blank(), axis.text.x=element_text(angle=90, size=6, vjust=0.5),
688 axis.title.x=element_blank(),
689 legend.text=element_text(size=6), legend.key.size=unit(3, "mm")) +
690 ylab("16S rRNA proportion") +
691 guides(fill = guide_legend(title = "Species", ncol = 1, reverse=TRUE)) # Ensures the legend is in one column
692
693 #----------------------------------
694 #Strip legend #1 - Cluster
695 #----------------------------------
696 annotation_plot1 <- ggplot(data = d.melt %>% distinct(sample_id, .keep_all=TRUE),
697 aes(x=sample_id, y=1, fill = cluster_group)) +
698 geom_tile(height = 0.1, aes(), colour="white", linewidth=0.3) +
699 scale_y_continuous(position = "right") +
700 scale_fill_manual( values=c("Diseased" = scales::alpha("#DD8B8B", .9), "Healthy" = scales::alpha("#00B050", .8)),
701 labels=c("Diseased", "Healthy")) +
702 theme_void() + theme(legend.position = "right") + ylab("Group") +
703 guides(fill = guide_legend(title = "Group")) +
704 theme(legend.text=element_text(size=6), legend.key.size=unit(3, "mm"),
705 axis.title.y.right = element_text(size = 8, hjust=0, face="bold",
706 margin = ggplot2::margin(l = 10, r = 10)))
707
708 #----------------------------------
709 #Strip legend #2 - Breeder ID
710 #----------------------------------
711 breeder_colours <- setNames(colorRampPalette(RColorBrewer::brewer.pal(9,"Pastel1"))(length(unique(d.melt$breed))),
712 unique(d.melt$breed))
713
714 annotation_plot2 <- ggplot(data = d.melt %>% distinct(sample_id, .keep_all=TRUE),
715 aes(x=sample_id, y=1, fill = breed)) +
716 geom_tile(height = 0.1, aes(), colour="white", linewidth=0.3) +
717 scale_y_continuous(position = "right") +
718 scale_fill_manual(values=breeder_colours)+
719 theme_void() + theme(legend.position = "right") + ylab("Breeder") +
720 guides(fill = guide_legend(title = "Breeder")) +
721 theme(legend.text=element_text(size=6), legend.key.size=unit(3, "mm"),
722 axis.title.y.right = element_text(size = 8, hjust=0, face="bold",
723 margin = ggplot2::margin(l = 10, r = 10)))
724
725 #----------------------------------
726 #Strip legend #3 - Builder ID
727 #----------------------------------
728 builder_colours <- setNames(colorRampPalette(RColorBrewer::brewer.pal(8,"Set2"))(length(unique(d.melt$build))),
729 unique(d.melt$build))
730
731 annotation_plot3 <- ggplot(data = d.melt %>% distinct(sample_id, .keep_all=TRUE),
732 aes(x=sample_id, y=1, fill = build)) +
733 geom_tile(height = 0.1, aes(), colour="white", linewidth=0.3) +
734 scale_y_continuous(position = "right") +
735 scale_fill_manual(values=colorRampPalette(RColorBrewer::brewer.pal(8,"Set2"))(length(unique(d.melt$build)))) +
736 theme_void() + theme(legend.position = "right") + ylab("Builder") +
737 guides(fill = guide_legend(title = "Builder")) +
738 theme(legend.text=element_text(size=6), legend.key.size=unit(3, "mm"),
739 axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
740
741 #---------------------------------------------
742 # COMBINE bar plots + tiles with patchwork
743 #---------------------------------------------
744 library(patchwork)
745
746 fig2.a <- tree.g.rotate +
747 annotation_plot1 + guides(fill = guide_legend(ncol=2, reverse=TRUE)) +
748 annotation_plot2 + guides(fill = guide_legend(ncol=3, reverse=TRUE)) +
749 annotation_plot3 + guides(fill = guide_legend(ncol=3, reverse=TRUE)) +
750 # Sidebar 5% of plot width
751 plot.bar + plot_layout(heights = c(0.15, 0.03, 0.03, 0.03, 1), guides = "collect") +
752 guides(fill = guide_legend(ncol=2, reverse=TRUE))
753
754 fig2.a
Phyloseq object conversion
755 #---------------------------------------------
756 # Convert to phyloseq object
757 #---------------------------------------------
758 sub.asv <- d.b.counts.raw %>% ungroup() %>%
759 tibble::column_to_rownames("sample_id") %>%
760 dplyr::select(matches("g__|s__")) %>% dplyr::select_if(colSums(.) != 0)
761
762 sub.meta <- d.b.counts.raw %>% ungroup() %>%
763 mutate(cluster_group=dplyr::recode(!!!setNames(labels_df$cluster ,labels_df$label), sample_id)) %>%
764 mutate(cluster_order=dplyr::recode(!!!setNames(labels_df$x ,labels_df$label), sample_id)) %>%
765 arrange(cluster_order) %>%
766 mutate(cluster_group=ifelse(cluster_group==1, "Diseased", "Healthy")) %>%
767 tibble::column_to_rownames("sample_id") %>%
768 dplyr::select(!matches("g__|s__"))
769
770 OTU.phyloseq = otu_table(as.matrix(sub.asv), taxa_are_rows = FALSE)
771 tax.in <- str_split_fixed(colnames(sub.asv), ";", 8)
772 colnames(tax.in) <- c("SV",
773 "domain",
774 "phylum",
775 "class",
776 "order",
777 "family",
778 "genus",
779 "species")
780 rownames(tax.in) <- colnames(sub.asv)
781 TAX.phyloseq = tax_table(as.matrix(tax.in))
782 META.phyloseq = sample_data(data.frame(sub.meta))
783 sample_names(META.phyloseq) <- sample_names(OTU.phyloseq)
784 ps <-phyloseq(OTU.phyloseq, TAX.phyloseq, META.phyloseq)
785 ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 168 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 9 sample variables ]
## tax_table() Taxonomy Table: [ 168 taxa by 8 taxonomic ranks ]786 #ASV counts table
787 reactable::reactable(otu_table(ps),
788 width = 1000,
789 wrap = FALSE,
790 striped = TRUE,
791 highlight = TRUE,
792 defaultPageSize = 15,
793 compact = TRUE,
794 defaultColDef = reactable::colDef(
795 align = "left",
796 headerStyle = list(textAlign = "left")),
797 theme = reactable::reactableTheme(
798 style = list(fontSize = "11px"),
799 cellPadding = "3px 6px",
800 headerStyle = list(fontSize = "10px")))801 #Sample data table
802 reactable::reactable(sample_data(ps),
803 width = 1000,
804 wrap = FALSE,
805 striped = TRUE,
806 highlight = TRUE,
807 defaultPageSize = 15,
808 compact = TRUE,
809 defaultColDef = reactable::colDef(
810 align = "left",
811 headerStyle = list(textAlign = "left")),
812 theme = reactable::reactableTheme(
813 style = list(fontSize = "11px"),
814 cellPadding = "3px 6px",
815 headerStyle = list(fontSize = "10px")))816 #Taxonomy table
817 reactable::reactable(tax_table(ps),
818 width = 1000,
819 wrap = FALSE,
820 striped = TRUE,
821 highlight = TRUE,
822 defaultPageSize = 15,
823 compact = TRUE,
824 defaultColDef = reactable::colDef(
825 align = "left",
826 headerStyle = list(textAlign = "left")),
827 theme = reactable::reactableTheme(
828 style = list(fontSize = "11px"),
829 cellPadding = "3px 6px",
830 headerStyle = list(fontSize = "10px")))Diversity metrics
A variety of diversity metrics were assessed to compare between healthy and diseased communities to determine whether disease is associated with loss of within-sample microbial diversity.
831 ####################################
832 #################################### Alpha diversity
833 ####################################
834
835 # Set ggplot2 themes for diversity plots:
836 diversity.theme <- theme_classic() +
837 theme(legend.position="none",
838 plot.title = element_text(size=6,
839 color="black"),
840 axis.title.x=element_blank(),
841 axis.title.y=element_text(size=8),
842 axis.line = element_blank(),
843 panel.border = element_rect(color = "black",
844 fill = NA,
845 linewidth = 0.5))
846
847 #---------------------Calculate diversity metrics with phyloseq and microbiome R packages
848 diversity.stats <- cbind(sub.meta,
849 estimate_richness(ps, split = TRUE, measures = c("Observed",
850 "Chao1",
851 "ACE",
852 "Shannon",
853 "Simpson",
854 "InvSimpson",
855 "Fisher")),
856 microbiome::dominance(ps, index = "all", relative = TRUE, aggregate = FALSE))
857
858 #---------------------- Shannon
859 pwc1 <- diversity.stats %>% ungroup() %>%
860 rstatix::wilcox_test(Shannon ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
861 rstatix::add_xy_position()
862
863 alpha.1 <- ggplot(diversity.stats, aes(x=cluster_group, y=Shannon)) +
864 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
865 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
866 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
867 "Healthy" = scales::alpha("#00B050", 0.4))) +
868 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
869 "Healthy" = scales::alpha("#00B050", 0.8))) +
870 ylab("Shannon Diversity") +
871 ggpubr::stat_pvalue_manual(pwc1, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
872 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
873 guides(fill = guide_legend(title = "Group")) +
874 ggtitle("Shannon") + diversity.theme
875
876 #----------------------Observed
877 pwc2 <- diversity.stats %>% ungroup() %>%
878 rstatix::wilcox_test(Observed ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
879 rstatix::add_xy_position()
880
881 alpha.2 <- ggplot(diversity.stats, aes(x=cluster_group, y=Observed)) +
882 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
883 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
884 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
885 "Healthy" = scales::alpha("#00B050", 0.4))) +
886 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
887 "Healthy" = scales::alpha("#00B050", 0.8))) +
888 ylab("Observed") +
889 ggpubr::stat_pvalue_manual(pwc2, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
890 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
891 guides(fill = guide_legend(title = "Group")) +
892 ggtitle("Observed") + diversity.theme
893
894 #----------------------Dbp
895 pwc3 <- diversity.stats %>% ungroup() %>%
896 rstatix::wilcox_test(dbp ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
897 rstatix::add_xy_position()
898
899 alpha.3 <- ggplot(diversity.stats, aes(x=cluster_group, y=dbp)) +
900 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
901 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
902 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
903 "Healthy" = scales::alpha("#00B050", 0.4))) +
904 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
905 "Healthy" = scales::alpha("#00B050", 0.8))) +
906 ylab("dbp Diversity") +
907 ggpubr::stat_pvalue_manual(pwc3, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
908 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
909 guides(fill = guide_legend(title = "Group")) +
910 ggtitle("Dbp") +
911 theme_classic() + diversity.theme
912
913 #----------------------Gini
914 pwc4 <- diversity.stats %>% ungroup() %>%
915 rstatix::wilcox_test(gini ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
916 rstatix::add_xy_position()
917
918 alpha.4 <- ggplot(diversity.stats, aes(x=cluster_group, y=gini)) +
919 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
920 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
921 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
922 "Healthy" = scales::alpha("#00B050", 0.4))) +
923 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
924 "Healthy" = scales::alpha("#00B050", 0.8))) +
925 ylab("Gini index") +
926 ggpubr::stat_pvalue_manual(pwc4, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
927 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
928 guides(fill = guide_legend(title = "Group")) +
929 ggtitle("gini") + diversity.theme
930
931 #---------------------- Absolute
932 pwc5 <- diversity.stats %>% ungroup() %>%
933 rstatix::wilcox_test(absolute ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
934 rstatix::add_xy_position()
935
936 alpha.5 <- ggplot(diversity.stats, aes(x=cluster_group, y=absolute)) +
937 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
938 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
939 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
940 "Healthy" = scales::alpha("#00B050", 0.4))) +
941 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
942 "Healthy" = scales::alpha("#00B050", 0.8))) +
943 ylab("Absolute") +
944 ggpubr::stat_pvalue_manual(pwc5, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
945 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
946 guides(fill = guide_legend(title = "Group")) +
947 ggtitle("Absolute") + diversity.theme
948
949 #----------------------Core abund
950 pwc6 <- diversity.stats %>% ungroup() %>%
951 rstatix::wilcox_test(core_abundance ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
952 rstatix::add_xy_position()
953
954 alpha.6 <- ggplot(diversity.stats, aes(x=cluster_group, y=core_abundance)) +
955 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
956 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
957 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
958 "Healthy" = scales::alpha("#00B050", 0.4))) +
959 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
960 "Healthy" = scales::alpha("#00B050", 0.8))) +
961 ylab("Core abundance") +
962 ggpubr::stat_pvalue_manual(pwc6, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
963 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
964 guides(fill = guide_legend(title = "Group")) +
965 ggtitle("Core abundance") + diversity.theme
966
967 #----------------------Simpson
968 pwc7 <- diversity.stats %>% ungroup() %>%
969 rstatix::wilcox_test(simpson ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
970 rstatix::add_xy_position()
971
972 alpha.7 <- ggplot(diversity.stats, aes(x=cluster_group, y=simpson)) +
973 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
974 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
975 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
976 "Healthy" = scales::alpha("#00B050", 0.4))) +
977 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
978 "Healthy" = scales::alpha("#00B050", 0.8))) +
979 ylab("Simpson dominance") +
980 ggpubr::stat_pvalue_manual(pwc7, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
981 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
982 guides(fill = guide_legend(title = "Group")) +
983 ggtitle("Simpson dominance") + diversity.theme
984
985 #----------------------Dmn
986 pwc8 <- diversity.stats %>% ungroup() %>%
987 rstatix::wilcox_test(dmn ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
988 rstatix::add_xy_position()
989
990 alpha.8 <- ggplot(diversity.stats, aes(x=cluster_group, y=dmn)) +
991 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
992 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
993 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
994 "Healthy" = scales::alpha("#00B050", 0.4))) +
995 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
996 "Healthy" = scales::alpha("#00B050", 0.8))) +
997 ylab("Dmn") +
998 ggpubr::stat_pvalue_manual(pwc8, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
999 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
1000 guides(fill = guide_legend(title = "Group")) +
1001 ggtitle("Dmn") + diversity.theme
1002
1003 #Combine all diversity metrics
1004 cowplot::plot_grid(alpha.1, alpha.2, alpha.3, alpha.4, alpha.5, alpha.6, alpha.7, alpha.8, ncol=4)
Figure 2B
1005 #Shannon diversity for final plot
1006 fig2.b <- alpha.1
1007 fig2.b
Figure 2C
1008 #Simpson dominance for final plot
1009 fig2.c <- alpha.7
1010 fig2.c
PERMANOVA tests
Differences in overall microbiota composition were tested using complementary distance metrics, including Aitchison distance (Euclidean distance on CLR-transformed ASV profiles), Bray–Curtis dissimilarity, and Jaccard distance. PERMANOVA models were used to evaluate the independent contributions of disease-associated clustering, breeder, and builder effects using 999 permutations. WdS and associated post hoc tests were additionally applied as complementary multivariate tests of differences between healthy- and disease-associated community states.
1011 ######################
1012 ###################### # Aitchison
1013 ######################
1014
1015 # Subset data
1016 sub.otu <- d.b.clr.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1017 dplyr::select(matches("s__")) %>% dplyr::select_if(colSums(.) != 0)
1018 sub.meta <- d.b.clr.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1019 dplyr::select(!matches("s__")) %>% mutate(breed_build=paste0(breed, "___", build)) %>%
1020 mutate_at(vars(matches("^breed$|^build$|breed_build")), ~as.factor(.))
1021
1022 #-------------------
1023 # PERMANOVA
1024 #-------------------
1025 set.seed(500)
1026 adonis.aitchison <- vegan::adonis2(sub.otu ~ breed + build + disease_status,
1027 data = sub.meta, permutations=999, method = "euclidean", by = "margin")
1028 adonis.aitchison## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = sub.otu ~ breed + build + disease_status, data = sub.meta, permutations = 999, method = "euclidean", by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## breed 6 1660.4 0.16673 1.2427 0.205
## build 9 3483.0 0.34976 1.7379 0.026 *
## disease_status 1 776.0 0.07792 3.4846 0.001 ***
## Residual 9 2004.1 0.20125
## Total 29 9958.5 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11029 #-------------------
1030 # WdS
1031 #-------------------
1032 #WdS test
1033 set.seed(500)
1034 WdS.aitchison <- WdS.test(vegan::vegdist(sub.otu, method = "euclidean"),
1035 as.factor(sub.meta$disease_status), nrep = 999)
1036 WdS.aitchison##
## Distance-based Multivariate Welch ANOVA
##
## data: euclidean distance matrix with 30 observations and grouping factor with 2 levels
## WdS = 4.5589, dfb = 1, number of permutations = 999, p-value = 0.001
## sample estimates:
## effect size estimator of variance, omega-squared (ω²)
## 0.10605041037 ######################
1038 ###################### # Bray-Curtis
1039 ######################
1040 # Subset data
1041 sub.otu <- d.b.counts.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1042 dplyr::select(matches("s__")) %>% dplyr::select_if(colSums(.) != 0)
1043
1044 sub.meta <- d.b.counts.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1045 dplyr::select(!matches("s__")) %>% mutate(breed_build=paste0(breed, "___", build)) %>%
1046 mutate_at(vars(matches("^breed$|^build$|breed_build")), ~as.factor(.))
1047
1048 #-------------------
1049 # PERMANOVA
1050 #-------------------
1051 set.seed(500)
1052 adonis.bray <- vegan::adonis2(sub.otu ~ breed + build + disease_status,
1053 data = sub.meta, permutations=999, method = "bray", by = "margin")
1054 adonis.bray## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = sub.otu ~ breed + build + disease_status, data = sub.meta, permutations = 999, method = "bray", by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## breed 6 0.9340 0.14213 1.3649 0.208
## build 9 1.3778 0.20965 1.3422 0.223
## disease_status 1 0.9187 0.13980 8.0554 0.002 **
## Residual 9 1.0265 0.15620
## Total 29 6.5715 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11055 #-------------------
1056 # WdS
1057 #-------------------
1058 #WdS test
1059 set.seed(500)
1060 WdS.bray <- WdS.test(vegan::vegdist(sub.otu, method = "bray"),
1061 as.factor(sub.meta$disease_status), nrep = 999)
1062 WdS.bray##
## Distance-based Multivariate Welch ANOVA
##
## data: bray distance matrix with 30 observations and grouping factor with 2 levels
## WdS = 13.454, dfb = 1, number of permutations = 999, p-value = 0.001
## sample estimates:
## effect size estimator of variance, omega-squared (ω²)
## 0.29334931063 ######################
1064 ###################### # Jaccard
1065 ######################
1066
1067 # Subset data
1068 sub.otu <- d.b.counts.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1069 dplyr::select(matches("s__")) %>% dplyr::select_if(colSums(.) != 0)
1070 sub.meta <- d.b.counts.asv %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1071 dplyr::select(!matches("s__")) %>%
1072 mutate(breed_build=paste0(breed, "___", build)) %>%
1073 mutate_at(vars(matches("^breed$|^build$|breed_build")), ~as.factor(.))
1074
1075 #-------------------
1076 # PERMANOVA
1077 #-------------------
1078 set.seed(500)
1079 adonis.jaccard <- vegan::adonis2(sub.otu ~ breed + build + disease_status,
1080 data = sub.meta, permutations=999, method = "jaccard", by = "margin")
1081 adonis.jaccard## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = sub.otu ~ breed + build + disease_status, data = sub.meta, permutations = 999, method = "jaccard", by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## breed 6 1.2604 0.15571 1.2652 0.216
## build 9 1.9404 0.23973 1.2986 0.187
## disease_status 1 0.8770 0.10835 5.2823 0.003 **
## Residual 9 1.4943 0.18461
## Total 29 8.0943 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11082 #-------------------
1083 # WdS
1084 #-------------------
1085 #WdS test
1086 set.seed(500)
1087 WdS.jaccard <- WdS.test(vegan::vegdist(sub.otu, method = "jaccard"),
1088 as.factor(sub.meta$disease_status), nrep = 999)
1089 WdS.jaccard##
## Distance-based Multivariate Welch ANOVA
##
## data: jaccard distance matrix with 30 observations and grouping factor with 2 levels
## WdS = 9.9269, dfb = 1, number of permutations = 999, p-value = 0.001
## sample estimates:
## effect size estimator of variance, omega-squared (ω²)
## 0.2293253Total bacterial abundance
Absolute bacterial loads were compared to determine whether the disease-associated microbiota shift reflects an overall increase in bacterial biomass.
1090 #####################
1091 ##################### Boxplots of CC3/13 vs CC12
1092 #####################
1093
1094 #Merge counts and metadata together, then melt dataframe
1095 d.mer <- merge(metadata, d.b %>% sjmisc::rotate_df(), by.x="sample_id", by.y="row.names", sort=FALSE, all=TRUE)
1096 d.mer.melt <- reshape2::melt(d.mer, id.vars=c(colnames(metadata)))
1097
1098 #Convert counts using qPCR data to get estimated absolute abundance
1099 d.b.freq <- d.b %>% mutate_all(~.+1) %>%
1100 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>% #FILTER BY ABUNDANCE CUTOFF
1101 apply(., 2, function(x) x/sum(x)) %>% as.data.frame() %>% #CONVER TO FREQ
1102 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=TRUE) %>%
1103 mutate(combined_names = (paste(.$sample_id, .$group, .$breed, .$build, sep="_")))
1104
1105 d.b.freq.melt <- reshape2::melt(d.b.freq, id.vars=c(colnames(metadata), "combined_names")) %>%
1106 filter(grepl("zotu", variable)) %>%
1107 mutate_at(vars(value), as.numeric) %>%
1108 filter(grepl("s__Melissococcus_plutonius$", variable)) %>%
1109 mutate(value_abs = value*(10^qpcr_16S_adj)) %>%
1110 mutate(value_abs_other = (1-value)*(10^qpcr_16S_adj)) %>%
1111 mutate(value_abs_sum = value_abs + value_abs_other) %>%
1112 mutate(CC3_CC12_ratio= 10^qpcr_CC3_adj/(10^qpcr_CC3_adj+ 10^qpcr_CC12_adj)) %>%
1113 mutate(value_abs_CC3= log10(value_abs * CC3_CC12_ratio)) %>%
1114 mutate(value_abs_CC12= log10(value_abs * (1-CC3_CC12_ratio))) %>%
1115 mutate(value_rel_CC3= value * CC3_CC12_ratio) %>%
1116 mutate(value_rel_CC12= value * (1-CC3_CC12_ratio)) %>%
1117 mutate(value_abs = log10(value_abs)) %>%
1118 mutate(value_abs_other = log10(value_abs_other)) %>%
1119 mutate(value_abs_sum = log10(value_abs_sum))
1120
1121 d.b.freq.melt <- reshape2::melt(d.b.freq, id.vars=c(colnames(metadata), "combined_names")) %>%
1122 filter(grepl("zotu", variable)) %>%
1123 mutate_at(vars(value), as.numeric) %>%
1124 filter(grepl("s__Melissococcus_plutonius$", variable)) %>%
1125 mutate(value_abs = value*(qpcr_16S_adj)) %>%
1126 mutate(value_abs_other = (1-value)*(qpcr_16S_adj)) %>%
1127 mutate(value_abs_sum = value_abs + value_abs_other) %>%
1128 mutate(CC3_CC12_ratio= 10^qpcr_CC3_adj/(10^qpcr_CC3_adj+ 10^qpcr_CC12_adj)) %>%
1129 mutate(value_abs_CC3= value_abs * CC3_CC12_ratio) %>%
1130 mutate(value_abs_CC12= value_abs * (1-CC3_CC12_ratio)) %>%
1131 mutate(value_rel_CC3= value * CC3_CC12_ratio) %>%
1132 mutate(value_rel_CC12= value * (1-CC3_CC12_ratio))Figure 2D
Total bacteria inlcuding M. plutonius
Absolute
bacterial loads were compared to determine whether the
disease-associated microbiota shift reflects an overall increase in
bacterial biomass.
1133 #----------------BOX PLOT (total bacteria inlcuding M. plutonius)
1134 pwc7 <- d.b.freq.melt %>%
1135 rstatix::wilcox_test(value_abs_sum ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1136 rstatix::add_xy_position()
1137
1138 fig2.d <- ggplot(d.b.freq.melt, aes(x=cluster_group, y=value_abs_sum)) +
1139 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
1140 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
1141 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
1142 "Healthy" = scales::alpha("#00B050", 0.4))) +
1143 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1144 "Healthy" = scales::alpha("#00B050", 0.8))) +
1145 ylab("Log10 gene copies") +
1146 ylim(0, 7.5) +
1147 ggtitle("Total bacteria") +
1148 ggpubr::stat_pvalue_manual(pwc7, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
1149 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
1150 guides(fill = guide_legend(title = "Group")) + diversity.theme
1151
1152 fig2.d
Figure 2E
Total bacteria NOT inlcuding M. plutonius
Removing
the contribution of M. plutonius allows the remaining bacterial
community to be examined independently, helping distinguish pathogen
expansion from a generalized increase in bacterial abundance.
1153 pwc8 <- d.b.freq.melt %>%
1154 rstatix::wilcox_test(value_abs_other ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1155 rstatix::add_xy_position()
1156
1157 fig2.e <- ggplot(d.b.freq.melt, aes(x=cluster_group, y=value_abs_other)) +
1158 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
1159 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
1160 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
1161 "Healthy" = scales::alpha("#00B050", 0.4))) +
1162 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1163 "Healthy" = scales::alpha("#00B050", 0.8))) +
1164 ylab("Log10 gene copies") +
1165 ylim(0, 7.5) +
1166 ggtitle("Total bacteria NOT including M. plutonius") +
1167 ggpubr::stat_pvalue_manual(pwc8, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
1168 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
1169 guides(fill = guide_legend(title = "Group")) + diversity.theme
1170
1171 fig2.e
Figure 2F
Total loads of M. plutonius using CC3-specific primers
Clonal-complex-specific qPCR was used to quantify
CC3-associated M. plutonius and determine whether this lineage
contributes to the disease-associated increase in pathogen burden.
1172 pwc9 <- d.b.freq.melt %>%
1173 rstatix::wilcox_test(value_abs_CC3 ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1174 rstatix::add_xy_position()
1175
1176 boxplot3 <- ggplot(d.b.freq.melt, aes(x=cluster_group, y=value_abs_CC3)) +
1177 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
1178 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
1179 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
1180 "Healthy" = scales::alpha("#00B050", 0.4))) +
1181 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1182 "Healthy" = scales::alpha("#00B050", 0.8))) +
1183 ylab("Log10 gene copies") +
1184 ylim(0, 7.5) +
1185 ggtitle("Melissococcus plutonius - CC3") +
1186 ggpubr::stat_pvalue_manual(pwc9, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
1187 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
1188 guides(fill = guide_legend(title = "Group")) + diversity.theme
1189
1190 boxplot3
1191 fig2.f <- boxplot3Figure 2G
Total loads of M. plutonius using CC12-specific primers
CC12-specific quantification determines whether the lineage
recovered from affected queen pupae is preferentially associated with
disease and helps resolve the strain-level basis of the M.
plutonius signal.
1192 pwc10 <- d.b.freq.melt %>%
1193 rstatix::wilcox_test(value_abs_CC12 ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1194 rstatix::add_xy_position()
1195
1196 boxplot4 <- ggplot(d.b.freq.melt, aes(x=cluster_group, y=value_abs_CC12)) +
1197 geom_boxplot(aes(fill=cluster_group), outliers=FALSE) +
1198 geom_point(aes(fill=cluster_group), size=1, alpha=0.5, position = position_jitter(width=0.2, height =0.03)) +
1199 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.5),
1200 "Healthy" = scales::alpha("#00B050", 0.4))) +
1201 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1202 "Healthy" = scales::alpha("#00B050", 0.8))) +
1203 ylab("Log10 gene copies") +
1204 ylim(0, 7.5) +
1205 ggtitle("Melissococcus plutonius - CC12") +
1206 ggpubr::stat_pvalue_manual(pwc10, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
1207 vjust=0, bracket.nudge.y=0, bracket.shorten=0.1) +
1208 guides(fill = guide_legend(title = "Group")) + diversity.theme
1209
1210 boxplot4
1211 fig2.g <- boxplot4Figure 2H
Here, ordination of community composition tests whether healthy and diseased queen pupae harbour compositionally distinct microbiota and identifies the taxa contributing most strongly to this separation.
1212 ################################
1213 ################################ PCoA PLOT
1214 ################################
1215 pcoa.in <- d.b.clr.genus %>%
1216 mutate(cluster_group=dplyr::recode(!!!setNames(labels_df$cluster ,labels_df$label), sample_id)) %>%
1217 mutate(cluster_order=dplyr::recode(!!!setNames(labels_df$x ,labels_df$label), sample_id)) %>%
1218 arrange(cluster_order) %>%
1219 mutate(cluster_group=ifelse(cluster_group==1, "Diseased", "Healthy"))
1220
1221 sub.asv <- pcoa.in %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1222 dplyr::select(matches("s__|g__")) %>% dplyr::select_if(colSums(.) != 0)
1223 sub.asv <- -min(sub.asv) + 1 + sub.asv
1224 sub.meta <- pcoa.in %>% ungroup() %>% tibble::column_to_rownames("sample_id") %>%
1225 dplyr::select(!matches("s__|g__"))
1226
1227 #No formula for regular PCoA using Bray-Curtis distance
1228 cap <- capscale(sub.asv ~ 1, distance = "bray")
1229
1230
1231 # Extract sample scores and species (taxa) scores
1232 VALUES.g2 <- merge(as.data.frame(scores(cap, display = "sites")),
1233 sub.meta,
1234 by="row.names",
1235 sort=FALSE, all=TRUE) %>%
1236 dplyr::rename(PC1="MDS1", PC2="MDS2")
1237
1238 LOADINGS2 <- as.data.frame(scores(cap, display = "species")) %>%
1239 dplyr::rename(PC1="MDS1", PC2="MDS2") %>%
1240 mutate(Variable=rownames(.)) %>%
1241 mutate(magnitude=sqrt(PC1^2 + PC2^2)) %>%
1242 arrange(desc(magnitude)) %>%
1243 #Plot top 15 taxa with greatest over magnitude of influence on ordination
1244 mutate(Variable=ifelse(row_number()<=15, Variable, NA))
1245
1246 #Print table of loads and magnitude values
1247 #-------------------------------------------------
1248 reactable::reactable(
1249 LOADINGS2,
1250 width = 1000,
1251 striped = TRUE,
1252 highlight = TRUE,
1253 defaultPageSize = 15,
1254 compact = TRUE,
1255 defaultColDef = reactable::colDef(
1256 align = "left",
1257 headerStyle = list(textAlign = "left")),
1258 theme = reactable::reactableTheme(
1259 style = list(fontSize = "11px"),
1260 cellPadding = "3px 6px",
1261 headerStyle = list(fontSize = "10px")))1262 #-------------------------------------------------
1263
1264 base_plot <- ggplot(VALUES.g2, aes(x=PC1, y=PC2, color=cluster_group,
1265 label=cluster_group,fill=cluster_group)) +
1266 stat_ellipse(mapping = NULL, data = NULL, geom = "polygon", position = "identity",
1267 type = "norm", level = 0.68, segments = 51,
1268 na.rm = FALSE, show.legend = NA, inherit.aes = TRUE) +
1269 geom_point(aes(), size=2, alpha=rep(0.85)) +
1270 xlab(paste0("Axis 1 (", round((cap$CA$eig / sum(cap$CA$eig) * 100)[1], 2), "%)")) +
1271 ylab(paste0("Axis 2 (", round((cap$CA$eig / sum(cap$CA$eig) * 100)[2], 2), "%)")) +
1272 geom_segment(data = LOADINGS2, aes(x = 0, y = 0, xend = (PC1 * 1), yend = (PC2* 1)),
1273 arrow = arrow(length = unit(5/20, "picas")),
1274 alpha=0.5, color = "black", inherit.aes = FALSE, linewidth=0.15) +
1275 ggrepel::geom_text_repel(data = LOADINGS2,aes(x =PC1 * 1.2, y = PC2 * 1.2, label = Variable),
1276 size = 2.5, alpha = 0.4, box.padding = 0.1, # Adjust padding around text
1277 point.padding = 0.2, # Space from original point
1278 max.overlaps = 100, inherit.aes = FALSE) +
1279 #---Colours
1280 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.2),
1281 "Healthy" = scales::alpha("#00B050", 0.05))) +
1282 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1283 "Healthy" = scales::alpha("#00B050", 0.8))) +
1284 labs(color = "Group") +
1285 theme(panel.border = element_rect(color ="black", fill = NA, size = 0.5, linetype = 1))+
1286 theme(axis.text.x=element_text(angle=0, hjust=1, size=8),
1287 axis.text.y=element_text(size=8),
1288 axis.title.x = element_text(size=10),
1289 axis.title.y = element_text(size=10)) +
1290 theme(panel.background = element_blank(),
1291 panel.grid.major = element_blank(),
1292 panel.grid.minor = element_blank())
1293
1294 #Pull out legend for plotting
1295 legend.a <- base_plot + guides(fill ="none",
1296 colour =guide_legend(title = "Group",
1297 override.aes = list(size = 3)))
1298 legend.a.plot <- cowplot::plot_grid(ggplotGrob(
1299 legend.a
1300 )$grobs[[
1301 which(sapply(ggplotGrob(legend.a)$grobs,
1302 function(x) x$name) == "guide-box")
1303 ]])
1304 type = c("density", "histogram", "boxplot",
1305 "violin", "densigram")
1306
1307 pb <- ggplot_build(base_plot)
1308
1309 # Extract axis limits
1310 xlim <- pb$layout$panel_params[[1]]$x.range
1311 ylim <- pb$layout$panel_params[[1]]$y.range
1312
1313 #===================
1314 # Lower - Boxplot 1
1315 #===================
1316 VALUES.g2$cluster_group <- factor(VALUES.g2$cluster_group,
1317 levels=c("Diseased", "Healthy"))
1318
1319 pwc.box.1 <- VALUES.g2 %>%
1320 rstatix::wilcox_test(PC1 ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1321 rstatix::add_xy_position()
1322
1323 print(pwc.box.1)## # A tibble: 1 × 16
## estimate .y. group1 group2 n1 n2 statistic p conf.low conf.high
## <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -0.708 PC1 Diseas… Healt… 17 13 6 5.01e-7 -0.854 -0.471
## # ℹ 6 more variables: method <chr>, alternative <chr>, y.position <dbl>,
## # groups <named list>, xmin <dbl>, xmax <dbl>1324 # study boxplot axis 1
1325 box.1 <- ggplot(VALUES.g2, aes(y=PC1, x=cluster_group)) +
1326 xlab(paste0('group\n','')) +
1327 geom_boxplot(aes(colour=cluster_group, fill=cluster_group), outlier.size = 0.1) +
1328 geom_point(size=0.7, alpha=0.5, position=position_jitter(width=0.07, height=-0.07)) +
1329 ggpubr::stat_pvalue_manual(pwc.box.1, label="p", hide.ns=FALSE, size=3, bracket.size=0.1,
1330 vjust=1, bracket.nudge.y=-0.0) +
1331 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.2),
1332 "Healthy" = scales::alpha("#00B050", 0.05))) +
1333 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1334 "Healthy" = scales::alpha("#00B050", 0.8))) +
1335 ylim(xlim) +
1336 theme(axis.ticks.y = element_blank(),
1337 panel.background = element_rect(fill='white',
1338 color = 'black'),
1339 axis.text.y = element_blank(),
1340 axis.title.x = element_blank(),
1341 axis.title.y = element_blank(),
1342 panel.grid = element_blank()) +
1343 coord_flip()
1344
1345 #===================
1346 # Right - Boxplot 1
1347 #===================
1348 # study boxplot axis 1
1349 pwc.box.2 <- VALUES.g2 %>%
1350 rstatix::wilcox_test(PC2 ~ cluster_group, p.adjust.method = "BH", detailed=TRUE) %>%
1351 rstatix::add_xy_position()
1352
1353 print(pwc.box.2)## # A tibble: 1 × 16
## estimate .y. group1 group2 n1 n2 statistic p conf.low conf.high
## <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 0.104 PC2 Diseased Healthy 17 13 120 0.711 -0.301 0.429
## # ℹ 6 more variables: method <chr>, alternative <chr>, y.position <dbl>,
## # groups <named list>, xmin <dbl>, xmax <dbl>1354 box.2 <- ggplot(VALUES.g2, aes(y=PC2, x=cluster_group)) +
1355 xlab(paste0('group\n','')) +
1356 geom_boxplot(aes(colour=cluster_group, fill=cluster_group), outlier.size = 0.1) +
1357 geom_point(size=0.7, alpha=0.5, position=position_jitter(width=0.07, height=-0.07)) +
1358 ggpubr::stat_pvalue_manual(pwc.box.2, label="p", hide.ns=FALSE, size=3,
1359 bracket.size=0.1, vjust=0, bracket.nudge.y=0.1) +
1360 ylim(c(-1.2, 1.1)) +
1361 ylim(ylim) +
1362 scale_fill_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.2),
1363 "Healthy" = scales::alpha("#00B050", 0.05))) +
1364 scale_colour_manual(values=c("Diseased" = scales::alpha("#DD8B8B", 0.9),
1365 "Healthy" = scales::alpha("#00B050", 0.8))) +
1366 theme(axis.ticks.x = element_blank(),
1367 panel.background = element_rect(fill='white', color = 'black'),
1368 axis.text.x = element_blank(),
1369 axis.title.x = element_blank(),
1370 axis.title.y = element_blank(),
1371 panel.grid = element_blank())1372 #===================
1373 # COMBINED
1374 #===================
1375
1376 library(patchwork)
1377 bray.plot.left <- (box.1 + theme(legend.position="none")) +
1378 (base_plot + theme(legend.position="none")) +
1379 plot_layout(heights = c(0.15, 0.85), guides = "collect")
1380
1381 bray.plot.right <- legend.a.plot + theme(plot.margin = ggplot2::margin(0,0,0,0, "cm")) +
1382 (box.2 + theme(legend.position="none")) +
1383 plot_layout(heights = c(0.15, 0.85), guides = "collect")
1384
1385 bray.plot.final <- {bray.plot.left | bray.plot.right} +
1386 plot_layout(widths = c(0.85, 0.15) ,
1387 guides = "collect")
1388
1389 bray.plot.final
1390 fig2.h <- bray.plot.finalHeatmap
Here, we scale relative abundance across samples (columns) per genus (rows) to highlight relevant clustering patterns between disease-associated taxa (e.g., Melissococcus, Enterococcus) and health-associated taxa (e.g., Ralstonia, Caulobacter, Mesorhizobium, Bradyrhizobium)
1391 ################################
1392 ################################ Heatmap
1393 ################################
1394
1395 heat.in <- d.b.clr.genus
1396 d.b.prev.otu <- heat.in %>% tibble::column_to_rownames("sample_id") %>%
1397 dplyr::select(matches("g__|s__")) %>%
1398 apply(., 2, function(x) scales::rescale(x, to=c(-1,1)))
1399
1400 d.b.prev.meta <- heat.in %>% tibble::column_to_rownames("sample_id") %>%
1401 dplyr::select(!matches("g__|s__"))
1402
1403 #Cluster taxa and samples
1404 cluster.tax <- hclust(dist(t(d.b.prev.otu), method="euclidian"), method="ward.D2")
1405 cluster.tax <- dendsort(cluster.tax, type="average")
1406 cluster.samples <- hclust(dist(d.b.prev.otu, method="euclidian"), method="ward.D2")
1407 cluster.samples <- dendsort(cluster.samples, type="average")
1408
1409 #Dendrograms for taxa and sample clustering
1410 g1.dendro.samples <- ggdendro::ggdendrogram(cluster.samples, labels = TRUE, rotate=FALSE)
1411 g2.dendro.tax<- ggdendro::ggdendrogram(cluster.tax, labels = FALSE, rotate=TRUE) +
1412 scale_y_reverse() +theme_void()
1413
1414 #Melt dataframe
1415 count.table.melt <- melt(heat.in, id.vars=c("sample_id", colnames(d.b.prev.meta))) %>%
1416 group_by(variable) %>%
1417 mutate(value=scales::rescale(value, to=c(-1,1))) %>% ungroup() %>%
1418 mutate(variable= factor(variable, levels=(cluster.tax$labels[cluster.tax$order]))) %>%
1419 mutate(sample_id_factor = factor(sample_id, levels=ggtree::get_taxa_name(tree.g.rotate)))
1420
1421 k2.heatmap <- ggplot(data = count.table.melt, aes(x=sample_id_factor, y=variable, fill=value)) +
1422 geom_tile() +
1423 scale_color_gradientn(colours=colorRampPalette(c("#000004FF","#000004FF","#451077FF", "#CD4071FF", "#FCFDBFFF"))(100),
1424 aesthetics = "fill") + #, limits=c(-4,4)
1425 scale_y_discrete(position = "right") +
1426 theme(axis.text.x=element_text(size=8,
1427 angle=90,
1428 vjust=0.5),
1429 axis.title.x=element_blank(),
1430 axis.title.y=element_blank(),
1431 axis.text.y=element_text(size=6,
1432 hjust=1))
1433
1434 k2.heatmap.combined <- tree.g.rotate +
1435 annotation_plot1 + guides(fill = guide_legend(ncol=1, reverse=TRUE)) +
1436 annotation_plot2 + guides(fill = guide_legend(ncol=1, reverse=TRUE)) +
1437 annotation_plot3 + guides(fill = guide_legend(ncol=1, reverse=TRUE)) +
1438 k2.heatmap + plot_layout(heights = c(0.05, 0.03,0.03,0.03, 1), guides = "collect") +
1439 guides(fill = guide_legend(ncol=1, reverse=TRUE, title="Scaled rel. abund."))
1440
1441 side.tree <- (plot_spacer() / g2.dendro.tax + plot_layout(heights = c(0.10, 0.9)))
1442
1443 k2.heatmap.combined.full <- (side.tree | k2.heatmap.combined) + plot_layout(widths = c(0.08, 0.88))
1444 k2.heatmap.combined.full
ALDEx2 statistics
Compositional differential-abundance analysis identifies the bacterial taxa most strongly associated with healthy or diseased microbiota, moving from community-level separation to the specific organisms underlying that shift. The following are the relevant statistics used to generate Figure 2I.
1445 #--------------------------
1446 # ALDEx2 statistics
1447 #--------------------------
1448
1449 d.b.filt.aldex <- d.b.species %>%
1450 tibble::rownames_to_column() %>% filter(!grepl("Chloroplast|Mitochondria",rowname)) %>%
1451 tibble::column_to_rownames() %>%
1452 tibble::rownames_to_column("tax") %>%
1453 #Correct M. plutonius taxonomy
1454 mutate(tax= gsub("s__Melissococcus_uncl", "s__Melissococcus_plutonius", tax)) %>%
1455 group_by(tax) %>% mutate_at(vars(!matches("tax")), ~sum(.)) %>% dplyr::slice(1) %>% ungroup() %>%
1456 tibble::column_to_rownames("tax")
1457
1458 d.b.counts.aldex <- d.b.filt.aldex %>% select_if(colSums(.) != 0) %>%
1459 #FILTER BY ABUNDANCE CUTOFF
1460 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.0001) %>%
1461 as.data.frame() %>%
1462 sjmisc::rotate_df() %>%
1463 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1464 arrange(sample_id) %>% arrange(cluster_group)
1465
1466 #---------------------------------Organize input for ALDEx2
1467 aldex.in <- d.b.counts.aldex
1468 d.aldex.otu <- aldex.in %>% tibble::column_to_rownames("sample_id") %>%
1469 dplyr::select(!matches(colnames(metadata))) %>%
1470 sjmisc::rotate_df() %>% arrange(desc(rowSums(.)))
1471 d.aldex.meta <- aldex.in %>% tibble::column_to_rownames("sample_id")
1472
1473 #---------------------------------
1474 #For reproducibility
1475 set.seed(500)
1476 d.aldex.meta$cluster_order## [1] 5 4 6 11 3 13 12 9 10 2 8 7 1 14 16 15 17 27 20 25 29 22 23 18 30
## [26] 24 19 28 26 211477 aldex.clr <- aldex.clr(d.aldex.otu, d.aldex.meta$cluster_group, mc.samples = 128,
1478 verbose = TRUE, useMC=TRUE, denom = "all")
1479 aldex.clr.eff <- aldex.effect(aldex.clr, verbose = FALSE,
1480 include.sample.summary = FALSE, useMC=TRUE, CI=FALSE)
1481 aldex.clr.ttest <- aldex.ttest(aldex.clr, verbose = TRUE, paired.test=FALSE)## |------------(25%)----------(50%)----------(75%)----------|1482 aldex.clr.merge <- cbind(aldex.clr.ttest, aldex.clr.eff) %>%
1483 mutate(labels=rownames(.)) %>%
1484 cbind(str_split_fixed(rownames(.), ";", 8)[,1:8], .) %>%
1485 dplyr::rename(sv=`1`,domain = `2`,phylum=`3`,class=`4`,
1486 order=`5`,family=`6`,genus=`7`,species=`8`) %>%
1487 mutate(sv_species=paste0(sv, "___", species))Supp Data 1C
1488 write_tsv(aldex.clr.merge, "figures/Supp_Data_1C.tsv")Figure 2I
Effect size plots visualizing ALDEx2-derived statistics.
1489 #::::::::::::::::::::::::
1490 # Effect bars
1491 #::::::::::::::::::::::::
1492 aldex.clr.merge.g <- aldex.clr.merge %>% arrange(desc(abs(effect))) %>%
1493 .[1:22,] %>%
1494 arrange(effect) %>%
1495 #Exclude for graphing
1496 #filter(!grepl("Lactobacillus_apis|Lactobacillus_panisapium", sv)) %>%
1497 mutate_at(vars(species,sv, sv_species), ~factor(., levels=unique(.)))
1498
1499 fig2.i <-ggplot(aldex.clr.merge.g, aes(x=sv_species, y=effect)) + geom_bar(stat="identity") +
1500 coord_flip() + theme_classic() +
1501 theme(axis.title.y=element_blank(),
1502 axis.text.y = element_text(size=8)) + ylab("ALDEx2 median effect size") +
1503 ylim(-2.6,2.6) +
1504 geom_hline(yintercept = c(-1,0,1),
1505 linetype = "longdash") #+
1506 fig2.i
Figure 2 (Full panel)
Collectively, these analyses show that disease is associated not simply with the presence of M. plutonius, but with a broader restructuring of the queen pupal microbiota characterized by altered diversity, community composition, and specific disease- and health-associated taxa.
1507 ############################################################
1508 # Combine plots
1509 ############################################################
1510
1511 fig2.left<- cowplot::plot_grid(fig2.a, labels=c("A"), label_size=10)
1512 fig2.right <- cowplot::plot_grid(fig2.b, fig2.c, fig2.d, fig2.e, fig2.f, fig2.g, ncol=2,
1513 labels=c("B", "C", "D", "E", "F", "G"), label_size=10)
1514
1515 fig2.top <- cowplot::plot_grid(fig2.left, fig2.right, ncol=2, rel_widths=c(0.6,0.4))
1516 fig2.bottom <- cowplot::plot_grid(fig2.h,
1517 fig2.i + theme(axis.text.y = element_text(size=6),
1518 axis.title.x=element_text(size=7)),
1519 ncol=2, rel_widths=c(0.6,0.4), labels=c("H", "I"), label_size=10)
1520
1521 fig2.final <- cowplot::plot_grid(fig2.top, fig2.bottom, ncol=1, nrow=2, rel_heights=c(0.5, 0.5))
1522 fig2.final
6. Functional profiling
Functional profiling reveals disease-associated shifts in
microbiome potential
We next asked whether the taxonomic
restructuring of the microbiota was accompanied by corresponding changes
in functional potential. Shotgun metagenomic profiles were therefore
compared across multiple annotation frameworks to identify functions and
pathways enriched in healthy or diseased queen pupae.
1523 library(dplyr)
1524 library(Biostrings)
1525 library(stringr)
1526 library(reshape2)
1527 library(ggplot2)
1528 library(heatmap3)
1529 library(ggheatmap)
1530 library(readr)
1531 library(ggdendro)
1532 library(cowplot)
1533 library(dendsort)
1534 library(ggiraph)
1535 library(ggordiplots)
1536 library(Maaslin2)
1537 library(ggrepel)
1538 library(ALDEx2)
1539 library(ape)
1540 library(ggtree)
1541 library(ggtreeExtra)
1542
1543 #--------------------------------------
1544 # Import Metacerberus roll up tables
1545 #--------------------------------------
1546
1547 input.list <- unique(paste0(str_split_fixed(list.files("data/metacerberus_results/final/rollup"), "_", 5)[,3], "_",
1548 str_split_fixed(list.files("data/metacerberus_results/final/rollup"), "_", 5)[,4]))
1549 rollup.list <- c("AMRFinder",
1550 "CAZy",
1551 "COG",
1552 "GVDB",
1553 "KOFam_all_FOAM",
1554 "KOFam_all_KEGG",
1555 "KOFam_prokaryote_FOAM",
1556 "KOFam_prokaryote_KEGG",
1557 "PFAM",
1558 "PGAP",
1559 "PHROG",
1560 "PVOG",
1561 "TIGRFAM",
1562 "VOG")
1563
1564 metadata <- readr::read_tsv("data/metadata.txt")
1565
1566 #--------------------------------------
1567 #Source the import function script
1568 #--------------------------------------
1569 source("scripts/function_rollup_collector.R")
1570 collector.rollup <- function_rollup_import(input.list = input.list, rollup.list= rollup.list)
1571
1572 #-------------------------------------------------------
1573 # Loop for ALDEx2 differential abundance tests
1574 #-------------------------------------------------------
1575 collector.aldex.list <- list()
1576
1577 for(i in rollup.list){
1578 # Load data / filter / format
1579 d.b <- collector.rollup[[i]] %>% distinct(concat_column, .keep_all=TRUE) %>%
1580 tibble::column_to_rownames("concat_column")
1581 #-----------------------
1582 d.b.filt <- d.b %>%
1583 tibble::rownames_to_column() %>%
1584 tibble::column_to_rownames() %>%
1585 sjmisc::rotate_df() %>% merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1586 tibble::column_to_rownames("sample_id") %>%
1587 dplyr::select(!matches(colnames(metadata))) %>% sjmisc::rotate_df() #%>% filter(rowMeans(. > 1) > 0.3)
1588
1589 d.b.filt2 <- d.b.filt %>% select_if(colSums(.) != 0) %>%
1590 #FILTER BY ABUNDANCE CUTOFF
1591 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.001) %>%
1592 mutate_all(~.+1) %>% #Simple zero inflation
1593 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
1594 {tmp <<- .} %>%
1595 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>% #CONVERT TO CLR
1596 {tmp2 <<- .} %>%
1597 sjmisc::rotate_df() %>%
1598 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1599 arrange(sample_id) %>% arrange(disease_status) %>% {d.b.clr <<- .}
1600
1601 d.b.counts <- d.b.filt %>% select_if(colSums(.) != 0) %>%
1602 #FILTER BY ABUNDANCE CUTOFF
1603 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.001) %>%
1604 as.data.frame() %>%
1605 sjmisc::rotate_df() %>%
1606 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1607 arrange(sample_id) %>% arrange(disease_status)
1608
1609 d.b.prev <- d.b.filt %>% select_if(colSums(.) != 0) %>%
1610 #FILTER BY ABUNDANCE CUTOFF
1611 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.001) %>%
1612 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
1613 as.data.frame() %>%
1614 sjmisc::rotate_df() %>%
1615 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1616 arrange(sample_id) %>% arrange(disease_status) %>%
1617 {d.b.freq <<- .} %>%
1618 mutate_at(vars(matches("__")), ~ifelse(as.numeric(.) >= 0.0000000001, 1, 0))
1619
1620 # ALDEx2 differential relative abundance test
1621 set.seed(500)
1622 d.aldex.otu <- d.b.counts %>% tibble::column_to_rownames("sample_id") %>%
1623 dplyr::select(!matches(colnames(metadata))) %>% sjmisc::rotate_df()
1624 d.aldex.meta <- d.b.counts %>% tibble::column_to_rownames("sample_id") %>%
1625 dplyr::select(matches(colnames(metadata)))
1626
1627 aldex.in <- aldex.clr(d.aldex.otu, d.aldex.meta$disease_status, mc.samples = 128,
1628 verbose = TRUE, useMC=TRUE, denom = "all")
1629 aldex.clr.eff <- aldex.effect(aldex.in, verbose = FALSE,
1630 include.sample.summary = FALSE, useMC=FALSE, CI=FALSE)
1631 aldex.clr.ttest <- aldex.ttest(aldex.in, verbose = TRUE, paired.test=FALSE)
1632 aldex.clr.merge <- cbind(aldex.clr.ttest, aldex.clr.eff) %>% mutate(labels=rownames(.))
1633
1634 collector.aldex.list[[i]] <- aldex.clr.merge %>% mutate(category=paste(i))
1635 }
1636
1637 #----------------------------
1638 # Combine ALDEx2 results
1639 #----------------------------
1640 collector.aldex.list.combined <- dplyr::bind_rows(collector.aldex.list) %>%
1641 filter(!grepl("KOFam_all", category)) %>%
1642 mutate(category=gsub("KOFam_prokaryote_", "", category))Effect plot (Interactive)
Differential-abundance analyses across complementary functional databases provide a broad view of the metagenomic features associated with each disease state. This reveals that the transition from health to disease extends beyond taxonomic turnover to coordinated changes in microbial functional potential.
1643 #----------------------------
1644 # Effect size plot
1645 #----------------------------
1646
1647 source("scripts/function_interactive_strip_plot.R")
1648
1649 aldex.plot2 <- interactive_strip_plot(df= collector.aldex.list.combined,
1650 eff="effect", pvalue="wi.ep", group="category", labs="labels",
1651 effect_cut=1,
1652 pval_cut=0.05,
1653 jitter.height=0.25,
1654 point.size=0.1,
1655 text.size=0,
1656 min.x=NULL,
1657 max.x=NULL,
1658 plot.labs=TRUE) +
1659 xlab("Median effect size (log2)") +
1660 guides(colour = guide_legend(override.aes = list(size = 4), title="Key"))+
1661 theme(legend.position = "right",
1662 legend.key = element_blank(),
1663 axis.text.y=element_text(size=7),
1664 axis.line =element_blank(),
1665 axis.ticks.x=element_line(),
1666 axis.ticks.length = unit(0.1, "cm"),
1667 panel.background=element_rect(fill=NA,
1668 colour="black",
1669 linewidth=0.5),
1670 panel.border = element_blank())
1671
1672 ggiraph::girafe(ggobj = aldex.plot2,
1673 options = list(ggiraph::opts_hover(css = "fill: yellow; stroke: yellow; r: 3px;" )))ALDEx2 statistics
Output the entire ALDEx2 statistics results table to TSV.
1674 reactable::reactable(collector.aldex.list.combined,
1675 width = 1000,
1676 wrap = FALSE,
1677 striped = TRUE,
1678 highlight = TRUE,
1679 resizable = TRUE,
1680 defaultPageSize = 15,
1681 compact = TRUE,
1682 defaultColDef = reactable::colDef(
1683 align = "left",
1684 headerStyle = list(textAlign = "left")),
1685 theme = reactable::reactableTheme(
1686 style = list(fontSize = "11px"),
1687 cellPadding = "3px 6px",
1688 headerStyle = list(fontSize = "10px")))Supp Data 1D
1689 write_tsv(collector.aldex.list.combined, "figures/Supp_Data_1D.tsv")Marginal density plot
1690 #----------------------------
1691 # Marginal density plot
1692 #----------------------------
1693 top_density <- ggplot(collector.aldex.list.combined %>%
1694 mutate(group=ifelse(effect < 0, "BAD", "GOOD")),
1695 aes(x = effect, fill = group, colour=group)) +
1696 geom_density(aes(y = after_stat(count)), adjust=2, bw=0.05, position = "identity") +
1697 scale_fill_manual(values=c("BAD" = scales::alpha("#DD8B8B", 0.3),
1698 "GOOD" = scales::alpha("#00B050", 0.2)),
1699 labels=c("BAD" = "Diseased",
1700 "GOOD" = "Healthy")) +
1701 scale_colour_manual(values=c("BAD" = scales::alpha("#DD8B8B", 0.9),
1702 "GOOD" = scales::alpha("#00B050", 0.8)),
1703 labels=c("BAD" = "Diseased",
1704 "GOOD" = "Healthy")) +
1705 theme_void() +
1706 xlim(range(ggplot_build(aldex.plot2)$layout$panel_params[[1]]$x.range))+
1707 theme(legend.position = "right") + theme(plot.tag = element_blank(),
1708 axis.line = element_blank(),
1709 axis.title = element_blank(),
1710 axis.text = element_blank(),
1711 axis.ticks = element_blank(),
1712 axis.ticks.length = unit(0, "pt"),
1713 plot.margin = ggplot2::margin(0,0,0,0,"cm"),
1714 panel.spacing = unit(0, "cm"))
1715 top_density
Ortholog counts
1716 #--------------------------------------------------------
1717 # Count number of significant features in each group
1718 #--------------------------------------------------------
1719 #Calculate total number of orthologs enriched in Diseased group
1720 nrow(collector.aldex.list.combined %>% filter(wi.ep <0.05 & effect <0))## [1] 301721 #Calculate total number of orthologs enriched in Healthy group
1722 nrow(collector.aldex.list.combined %>% filter(wi.ep <0.05 & effect >0))## [1] 184Figure 3A
1723 collector.aldex.list.combined.bar <- collector.aldex.list.combined %>%
1724 mutate(group=ifelse(effect < 0, "BAD", "GOOD")) %>%
1725 mutate(effect2=ifelse(wi.ep <0.05, 1, 0)) %>% group_by(category, group) %>%
1726 mutate(counts=ifelse(sum(effect2)==0, 0.1, sum(effect2))) %>% dplyr::slice(1) %>% ungroup() %>%
1727 mutate(group = factor(group, levels=c("GOOD", "BAD")))
1728
1729 bar.counts <- ggplot(collector.aldex.list.combined.bar,
1730 aes(x = factor(category, levels=rev(unique(category))), y=counts)) +
1731 geom_bar(aes(fill=group), stat="identity", width=0.5, position=position_dodge()) +
1732 scale_fill_manual(values=c("BAD" = scales::alpha("#DD8B8B", 0.9),
1733 "GOOD" = scales::alpha("#00B050", 0.5)),
1734 labels=c("BAD" = "Diseased",
1735 "GOOD" = "Healthy")) +
1736 scale_y_continuous(breaks=c(0,30)) +
1737 coord_flip() + theme_void() +
1738 theme(legend.position="none",
1739 axis.text.x=element_text(size=8, angle=0),
1740 axis.ticks.x=element_line(),
1741 axis.ticks.length = unit(0.1, "cm"),
1742 axis.line.x=element_line(linewidth=0.25))
1743
1744 #-------------------------------
1745 # Combine Fig 3A
1746 #-------------------------------
1747
1748 library(patchwork)
1749 #For interactive
1750 fig3.a.combined.1 <- top_density + aldex.plot2 + plot_layout(heights = c(0.1, 1))
1751 fig3.a.combined.2 <- plot_spacer() + bar.counts + plot_layout(heights = c(0.1, 1))
1752 fig3.a <- {fig3.a.combined.1 | fig3.a.combined.2} + plot_layout(widths = c(1, 0.2), guides = "collect")
1753
1754 ggiraph::girafe(ggobj = fig3.a,
1755 options = list(ggiraph::opts_hover(css = "fill: yellow; stroke: yellow; r: 2px;" )))1756 #For saving as static image
1757 fig3.a.combined.1 <- (top_density + theme(legend.position="none")) +
1758 (aldex.plot2+ theme(legend.position="none")) + plot_layout(heights = c(0.1, 1))
1759 fig3.a.combined.2 <- plot_spacer() + bar.counts + plot_layout(heights = c(0.1, 1))
1760 fig3.a <- {fig3.a.combined.1 | fig3.a.combined.2} + plot_layout(widths = c(1, 0.2), guides = "collect")KEGG module profiles
KEGG module-level analysis resolves these broader functional differences into specific metabolic and cellular pathways and visualizes how disease-associated functions are distributed across individual samples.
1761 #Run script to parse EGGNOG KOs
1762 #--------------------------------------------------------------------
1763 source("scripts/function_eggnog_parse_KEGG_KOs.R")
1764 #Outputs 2 files:
1765 # "data/eggnog_results/eggnog_annotations_parsed.tsv"
1766 # "data/eggnog_results/eggnog_annotations_parsed_modules.tsv"
1767 #--------------------------------------------------------------------
1768
1769 kegg.names <- readr::read_tsv("/home/brendan/.conda/envs/metacerberus/lib/python3.12/site-packages/meta_cerberus/DB/KEGG.tsv",
1770 show_col_types=FALSE, col_names=TRUE) %>% distinct(ID, .keep_all=TRUE)
1771
1772 eggnog.functions <- read_tsv("data/eggnog_results/eggnog_annotations_parsed.tsv") %>%
1773 `colnames<-`(gsub("_trim_nohost", "", colnames(.))) %>% as.data.frame(.) %>%
1774 `rownames<-`(paste(.$KO,.$KO,.$KO, .$KO, .$KO, .$KO,.$KO, .$KO, sep=";")) %>% dplyr::select(-KO)
1775
1776 eggnog.modules <- read_tsv("data/eggnog_results/eggnog_annotations_parsed_modules.tsv") %>%
1777 `colnames<-`(gsub("_trim_nohost", "", colnames(.))) %>% as.data.frame(.) %>%
1778 filter(!is.na(module_accession)) %>%
1779 mutate(pathway_class2 = str_split_fixed(pathway_class, "; ", 3)[,2]) %>%
1780 dplyr::relocate("pathway_class2", .after="pathway_class") %>%
1781 `rownames<-`(paste(.$module_accession,.$module_accession,.$module_accession,
1782 .$module_accession, .$module_accession, .$module_accession,
1783 .$module_accession, .$module_accession, sep=";")) %>%
1784 mutate_at(vars(matches("BC")), ~as.numeric(.))
1785
1786 index.module.module <- setNames(eggnog.modules$pathway_name , eggnog.modules$module_accession)
1787 index.module.class <- setNames(eggnog.modules$pathway_class , eggnog.modules$module_accession)
1788 index.module.class2 <- setNames(eggnog.modules$pathway_class2 , eggnog.modules$module_accession)
1789
1790 d.b <- eggnog.modules %>% dplyr::select(matches("BCMy"))
1791
1792 #----------------
1793
1794 metadata <- readr::read_tsv("data/metadata.txt")
1795
1796 ################################
1797 ################################ Load data / filter / format
1798 ################################
1799 d.b.filt <- d.b %>%
1800 tibble::rownames_to_column() %>%
1801 tibble::column_to_rownames()
1802
1803 d.b.filt2 <- d.b.filt %>% dplyr::select_if(colSums(.) != 0) %>%
1804 #FILTER BY ABUNDANCE CUTOFF
1805 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.000001) %>%
1806 mutate_all(~.+1) %>% #Simple zero inflation
1807 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
1808 {tmp <<- .} %>%
1809 #CONVERT TO CLR
1810 apply(., 2, function(x) log2(x) - mean(log2(x))) %>% as.data.frame() %>%
1811 {tmp2 <<- .} %>%
1812 sjmisc::rotate_df() %>%
1813 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1814 arrange(sample_id) %>% {d.b.clr <<- .}
1815
1816 d.b.counts <- d.b.filt %>% select_if(colSums(.) != 0) %>%
1817 #FILTER BY ABUNDANCE CUTOFF
1818 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.000001) %>%
1819 as.data.frame() %>%
1820 sjmisc::rotate_df() %>%
1821 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>% arrange(sample_id)
1822
1823 d.b.prev <- d.b.filt %>% select_if(colSums(.) != 0) %>%
1824 #FILTER BY ABUNDANCE CUTOFF
1825 filter(apply(apply(., 2, function(x) x/sum(x)), 1, mean) > 0.000001) %>%
1826 apply(., 2, function(x) x/sum(x)) %>% #CONVER TO FREQ
1827 as.data.frame() %>%
1828 sjmisc::rotate_df() %>%
1829 merge(metadata, ., by.x="sample_id", by.y="row.names", all=FALSE) %>%
1830 arrange(sample_id) %>% #arrange(group) %>%
1831 {d.b.freq <<- .} %>%
1832 mutate_at(vars(!matches(colnames(metadata))), ~ifelse(. >= 0.000001, 1, 0))MaAsLin2 statistics
1833 #-----------------------------------------
1834 # MAASLIN to get significant functions
1835 #-----------------------------------------
1836 set.seed(500)
1837 d.mer.counts <- d.b.clr %>% `rownames<-`(.$sample_id) %>%
1838 dplyr::select(!matches(colnames(metadata))) %>%
1839 mutate_all(., ~scales::rescale(., to=c(-1,1)))
1840 d.mer.meta <- d.b.clr %>% `rownames<-`(.$sample_id) %>%
1841 dplyr::select(matches(colnames(metadata)))
1842
1843 maaslin.results <- Maaslin2(d.mer.counts, d.mer.meta,
1844 output='figures/FIGURES___05_maaslin2_output',
1845 transform = "none",
1846 min_prevalence = 0, min_abundance = -50000,
1847 plot_heatmap=FALSE, plot_scatter=FALSE, standardize = FALSE,
1848 analysis_method="LM", normalization="none",
1849 fixed_effects= c("disease_status"),
1850 #random_effects=c("breed", "build"),
1851 reference=c("disease_status,Diseased"))
1852
1853 collector.df.bind <- maaslin.results$results %>%
1854 mutate(col_group = case_when(qval < 0.05 & abs(coef) > 1 ~ "positive_pval_coef",
1855 qval < 0.05 ~ "positive_pval",
1856 TRUE ~ "group1")) %>%
1857 arrange(feature) %>% mutate(feature_rename=str_split_fixed(feature, "[.]", 8)[,1]) %>%
1858 mutate(module = dplyr::recode(!!!index.module.module, feature_rename)) %>%
1859 mutate(class = dplyr::recode(!!!index.module.class , feature_rename)) %>%
1860 mutate(module=paste0(feature_rename, ": ", module))
1861
1862 collector.df.bind.save <- collector.df.bind %>%
1863 dplyr::select(-name,-col_group,-feature_rename) %>%
1864 mutate(feature=str_split_fixed(feature,"\\.", 2)[,1])Supp Data 1E
1865 readr::write_tsv(collector.df.bind.save, file="figures/Supp_Data_1E.tsv")KEGG module volcano plot
1866 #Get most significant fetures for plotting
1867 #median(collector.df.bind$qval)
1868 sig.list <- (collector.df.bind %>% arrange(pval) %>% filter(pval < 0.05))$feature_rename
1869
1870 #--------------------------------------
1871 # Volcano plot of maaslin2 statistics
1872 #--------------------------------------
1873
1874 volcano.g <- ggplot(collector.df.bind,
1875 aes(x = coef, y = log10(pval), colour=col_group)) +
1876 geom_point_interactive(aes(color = col_group,
1877 tooltip = paste0(
1878 "KEGG ID: ", .data[["feature_rename"]],
1879 "\nModule name: ", str_split_fixed(.data[["module"]], " ", 2)[,2],
1880 "\nModule hierarchy: ", .data[["class"]],
1881 "\nCoefficient: ", round(.data[["coef"]], 3),
1882 "\np-value: ", sprintf("%.2e", .data[["pval"]]),
1883 "\np-value: ", sprintf("%.2e", .data[["qval"]])),
1884 data_id =feature_rename ),
1885 size=1.5,
1886 alpha=0.7,
1887 hover_nearest = TRUE,
1888 position= position_jitter(width=0.05,height=0.1)) +
1889 scale_color_manual(values=c(positive_pval_coef = "blue",
1890 positive_pval = "red3",
1891 group1="grey60"),
1892 labels=c(group1 ="ns",
1893 positive_pval = 'qval < 0.05',
1894 positive_pval_coef = 'qval < 0.05 + |coef| > 1')) +
1895 geom_vline(xintercept = c(0), linetype="dashed", color="black", alpha=0.9) +
1896 labs(x = "Log Fold Change (coef)", y = "-Log10 p-value") + scale_y_reverse() +
1897 geom_text_repel(aes(label = ifelse(qval < 0.01 & coef >0 | qval < 0.01 & coef < 0, module, "")),
1898 max.overlaps=50, size = 2) +
1899 theme_minimal() + # theme(legend.position = "none") +
1900 guides(colour=guide_legend(title="Key")) +
1901 ggplot2::theme(legend.text=element_text(size=8),
1902 legend.key.size=unit(4, "mm"),
1903 panel.border = element_rect(color = "black",
1904 fill = NA,
1905 size = 0.5),
1906 panel.grid=element_blank())
1907
1908 volcano.labels <- cowplot::ggdraw() + cowplot::draw_label("Enriched in\nDiseased group",
1909 size = 8, alpha = 1,
1910 colour=scales::alpha("#DD8B8B", 1)) |
1911 cowplot::ggdraw() + cowplot::draw_label("Enriched in\nHealthy group",
1912 size = 8, alpha = 1, colour=scales::alpha("#00B050", 0.8))
1913 volcano.g.final <- (volcano.labels / volcano.g) + plot_layout(heights =c(0.1, 1))
1914
1915 ggiraph::girafe(ggobj = volcano.g.final,
1916 options = list(ggiraph::opts_hover(css = "fill: yellow; stroke: yellow; r: 3px;" )))KEGG module heatmap
1917 #----------------------
1918 # KEGG module heatmap
1919 #----------------------
1920 d.b.clr.filt <- d.b.clr %>% `colnames<-`(gsub(";", "_", colnames(.)))
1921 d.mer.meta <- d.b.clr.filt %>% dplyr::select(matches(colnames(metadata)))
1922 d.mer.counts <- d.b.clr.filt %>% tibble::column_to_rownames("sample_id") %>%
1923 dplyr::select(!matches(colnames(metadata))) %>% sjmisc::rotate_df() %>%
1924 tibble::rownames_to_column() %>% rowwise() %>%
1925 mutate(mean_BAD = mean(c_across((d.mer.meta %>%
1926 filter(grepl("Diseased",disease_status)))$sample_id))) %>%
1927 mutate(mean_GOOD = mean(c_across((d.mer.meta %>%
1928 filter(grepl("Healthy",disease_status)))$sample_id))) %>%
1929 mutate(mean_diff = (mean_GOOD - mean_BAD)) %>% ungroup() %>%
1930 tibble::column_to_rownames() %>% arrange(desc(abs(mean_diff))) %>%
1931 filter(str_split_fixed(rownames(.), "_", 2)[,1] %in% sig.list[1:50])
1932
1933 df.heat <- d.b.clr %>% arrange(disease_status) %>% tibble::column_to_rownames("sample_id") %>%
1934 dplyr::select(gsub("_",";", rownames(d.mer.counts))) %>%
1935 `colnames<-`(str_split_fixed(colnames(.), ";", 2)[,1])
1936 #----------scale and then add back in rownames
1937 df.heat <- apply(df.heat , 2, function(x) scale(x)) %>% as.data.frame() %>% `rownames<-`(rownames(df.heat))
1938 #----------
1939
1940 df.heat.melt <- d.b.clr %>% arrange(disease_status) %>%
1941 tibble::column_to_rownames("sample_id") %>%
1942 dplyr::select(disease_status, breed, build) %>% cbind(., df.heat)
1943
1944 #-----------------------
1945 cluster.tax <- hclust(dist(t(df.heat), method="euclidian"), method="ward.D2")
1946 cluster.samples <- hclust(dist(df.heat, method="euclidian"), method="ward.D2")
1947 cluster.samples <- dendsort(cluster.samples, type="average")
1948 #---------------------------------------------------------------------------------------------------------
1949 cluster.samples.x <- as.phylo(cluster.samples)
1950 g1.dendro.samples.gg <- ggtree(cluster.samples.x, ladderize=TRUE) +
1951 coord_flip() + scale_x_reverse()
1952
1953 tip_labels_in_order <- g1.dendro.samples.gg$data$label[g1.dendro.samples.gg$data$isTip][order(g1.dendro.samples.gg$data$y[g1.dendro.samples.gg$data$isTip])]
1954
1955 #---------------------------------------------------------------------------------------------------------
1956
1957 #Dendrograms
1958 g1.dendro.samples <- ggdendro::ggdendrogram(cluster.samples, labels = FALSE, rotate=FALSE)+ theme_void()
1959 g2.dendro.tax <- ggdendro::ggdendrogram(cluster.tax, labels = FALSE, rotate=TRUE) + scale_y_reverse() +theme_void()
1960
1961 #Main heatmap
1962 count.table.melt <- melt(df.heat.melt %>% tibble::rownames_to_column("sample_id"),
1963 id.vars=c("sample_id", "disease_status", "breed", "build")) %>%
1964 merge(., collector.df.bind, by.x="variable", by.y="feature_rename", sort=FALSE, all=FALSE) %>%
1965 dplyr::rename(value=value.x) %>%
1966 mutate(variable=str_split_fixed(variable, ";", 2)[,1]) %>%
1967 mutate(variable_factor = factor(variable, levels=cluster.tax$labels[cluster.tax$order])) %>%
1968 arrange(variable_factor) %>%
1969 mutate(module = dplyr::recode(!!!index.module.module, variable)) %>%
1970 mutate(module=paste0(variable, ": ", module)) %>%
1971 mutate(module_factor = factor(module, levels=unique(module))) %>%
1972 mutate(class = dplyr::recode(!!!index.module.class, variable)) %>%
1973 mutate(class_factor = factor(class, levels=unique(class))) %>%
1974 mutate(class2 = dplyr::recode(!!!index.module.class2, variable)) %>%
1975 mutate(class2_factor = factor(class2, levels=unique(class2))) %>%
1976 mutate(sample_id_factor= factor(sample_id,
1977 levels=cluster.samples$labels[cluster.samples$order]))
1978
1979 col.index <- setNames(colorRampPalette(RColorBrewer::brewer.pal(8,"Dark2"))(length(unique(sort(count.table.melt$class2)))),
1980 unique(sort(count.table.melt$class2)))
1981
1982 count.table.melt <- count.table.melt %>% mutate(class_colours = dplyr::recode(!!!col.index, class2))
1983
1984 mas.heatmap <- ggplot(data = count.table.melt,
1985 aes(x=sample_id_factor,
1986 y=module_factor,
1987 fill=value)) +
1988 geom_tile_interactive(aes(tooltip = paste0("KEGG ID: ", .data[["variable"]],
1989 "\nModule group: ", .data[["class2"]],
1990 "\nCoefficient: ", round(.data[["coef"]], 3),
1991 "\np-value: ", sprintf("%.2e", .data[["pval"]]),
1992 "\np-value: ", sprintf("%.2e", .data[["qval"]])),
1993 data_id =paste0(sample_id_factor, "_", variable)),
1994 hover_nearest = TRUE) +
1995 scale_color_gradientn(colours=(viridis::magma(100, alpha=0.85, begin=0.25)),
1996 aesthetics = "fill") +
1997 scale_y_discrete(position = "right") +
1998 theme(legend.text=element_text(size=8),
1999 legend.key.size=unit(3, "mm"),
2000 panel.background=element_blank(),
2001 axis.title.y=element_blank(),
2002 axis.text.x=element_text(size=8,
2003 angle=90,
2004 vjust=0.5),
2005 axis.title.x=element_blank(),
2006 axis.text.y=element_text(size=6, hjust=1,
2007 colour=(as.character((count.table.melt %>%
2008 distinct(module_factor, .keep_all=TRUE))$class_colours))))
2009
2010 #mas.heatmap
2011 legend.df <- count.table.melt %>% distinct(class2, class_colours) %>%
2012 arrange(class2) %>% mutate(class2= factor(class2, levels=unique(class2)))
2013 legend.plot <- ggplot(legend.df, aes(x=1, y=class2, color=class2)) +
2014 geom_point() + scale_color_manual(values = setNames(legend.df$class_colours,
2015 legend.df$class2)) +
2016 guides(color = guide_legend(override.aes = list(size = 5))) +
2017 theme_void() + theme(legend.position = "right",legend.title = element_blank())
2018
2019 #Pull out legend for plotting
2020 fig3.b.legend <- legend.plot + guides(fill ="none",
2021 colour =guide_legend(title = "Group",
2022 override.aes = list(size = 3))) +
2023 theme(legend.text=element_text(size=8),
2024 legend.key.size=unit(3, "mm"))
2025 fig3.b.legend.plot <- cowplot::plot_grid(ggplotGrob(
2026 fig3.b.legend
2027 )$grobs[[
2028 which(sapply(ggplotGrob(fig3.b.legend)$grobs,
2029 function(x) x$name) == "guide-box")
2030 ]])
2031
2032 #mas.heatmap | fig3.b.legend.plot + plot_layout(ncol=1, widths = c(1, 0.05))
2033 #cowplot::plot_grid(mas.heatmap, legend.plot, ncol=2, rel_widths=c(1,0.05))
2034
2035 #----------------------------------
2036 #Strip legend #1 - Disease status
2037 #----------------------------------
2038 category_colors <- c("Diseased" = scales::alpha("#DD8B8B", 0.9),
2039 "Healthy" = scales::alpha("#00B050", 0.8))
2040
2041 annotation_plot1 <- ggplot(data = count.table.melt %>% distinct(sample_id, .keep_all=TRUE) %>%
2042 mutate(sample_id=factor(sample_id, levels=unique(.$sample_id))),
2043 aes(x=sample_id_factor, y=1, fill = disease_status)) +
2044 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2045 scale_fill_manual(values = category_colors, name = "Category") + # Group color scale
2046 scale_y_continuous(position = "right") +
2047 theme_void() + theme(legend.position = "right") + ylab("Group") +
2048 theme(legend.text=element_text(size=8),
2049 legend.key.size=unit(3, "mm"),
2050 axis.title.y.right = element_text(size = 8,
2051 hjust=0,
2052 face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2053
2054 #----------------------------------
2055 #Strip legend #2 - Breeder ID
2056 #----------------------------------
2057 breeder_colours <- setNames(colorRampPalette(RColorBrewer::brewer.pal(9,"Pastel1"))(length(unique(count.table.melt$breed))),
2058 unique(count.table.melt$breed))
2059
2060
2061 annotation_plot2 <- ggplot(data = count.table.melt %>% distinct(sample_id, .keep_all=TRUE) %>%
2062 mutate(sample_id=factor(sample_id, levels=unique(.$sample_id))),
2063 aes(x=sample_id_factor, y=1, fill = breed)) +
2064 geom_tile(height = 0.1, aes()) +
2065 scale_fill_manual(values=breeder_colours, name = "Breeder ID") + #Breeder colour scale
2066 scale_y_continuous(position = "right") +
2067 theme_void() + theme(legend.position = "right") + ylab("Breed") +
2068 guides(fill = guide_legend(ncol=2)) +
2069 theme(legend.text=element_text(size=8),
2070 legend.key.size=unit(3, "mm"),
2071 axis.title.y.right = element_text(size = 8,
2072 hjust=0, face="bold",
2073 margin =ggplot2::margin(l = 10, r = 10)))
2074
2075 #----------------------------------
2076 #Strip legend #3 - Builder ID
2077 #----------------------------------
2078 builder_colours <- setNames(colorRampPalette(RColorBrewer::brewer.pal(8,"Set2"))(length(unique(count.table.melt$build))),
2079 unique(count.table.melt$build))
2080
2081 annotation_plot3 <- ggplot(data = count.table.melt %>% distinct(sample_id, .keep_all=TRUE) %>%
2082 mutate(sample_id=factor(sample_id, levels=unique(.$sample_id))),
2083 aes(x=sample_id_factor, y=1, fill = build)) +
2084 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2085 scale_fill_manual(values = builder_colours, name = "Builder ID") + # Builder color scale
2086 scale_y_continuous(position = "right") +
2087 theme_void() + theme(legend.position = "right") + ylab("Build") +
2088 guides(fill = guide_legend(ncol=2)) +
2089 theme(legend.text=element_text(size=8),
2090 legend.key.size=unit(3, "mm"),
2091 axis.title.y.right = element_text(size = 8,
2092 hjust=0,
2093 face="bold",
2094 margin = ggplot2::margin(l = 10, r = 10))) Figure 3B
2095 #---------------------------------------------------------
2096 # Fig 3B - Combined heatmap + sidebars using Patchwork
2097 #---------------------------------------------------------
2098 fig3.b.side.dendro <- ((plot_spacer() +
2099 labs(tag = "B") +
2100 theme(plot.tag = element_text(face = "bold", size = 12),
2101 plot.tag.position = c(0, 1))) / plot_spacer() / plot_spacer() / plot_spacer() / g2.dendro.tax) +
2102 plot_layout(heights = c(0.01,0.01,0.01,0.01,1)) # guides = "collect"
2103
2104 fig3.b.middle <- (g1.dendro.samples /
2105 annotation_plot1 /
2106 annotation_plot2 /
2107 annotation_plot3 /
2108 mas.heatmap) + plot_layout(heights = c(0.05,0.03,0.03,0.03,1))
2109
2110 fig3.b.legend.plot.x <- (plot_spacer() / plot_spacer() / plot_spacer() / plot_spacer() / fig3.b.legend.plot) +
2111 plot_layout(heights = c(0.1,0.1,0.1,0.1,0.1)) # guides = "collect"
2112
2113 fig3.b.combined <- (fig3.b.side.dendro | fig3.b.middle | fig3.b.legend.plot.x) +
2114 plot_layout(widths = c(0.1, 1, 0.1), guides = "collect")
2115 fig3.b.combined
Figure 3 (Full panel)
Together, the functional analyses demonstrate that disease-associated microbiota are functionally distinct as well as taxonomically altered, supporting a broader transition in the ecological state of the queen-associated microbial community.
2116 #---------------------------------------------------------
2117 # Fig 3 (Full panel)
2118 #---------------------------------------------------------
2119 #Remake interactive plot for Fig 3A
2120 fig3.a.combined.1 <- top_density + aldex.plot2 + plot_layout(heights = c(0.1, 1)) # Sidebar 5% of plot width
2121 fig3.a.combined.2 <- plot_spacer() + bar.counts + plot_layout(heights = c(0.1, 1))
2122 fig3.a.combined <- {fig3.a.combined.1 | fig3.a.combined.2} + plot_layout(widths = c(1, 0.2), guides = "collect")
2123
2124 #Combine Fig 3A + 3B panels
2125 fig3 <- (fig3.a.combined.1 |
2126 fig3.a.combined.2 |
2127 fig3.b.side.dendro |
2128 fig3.b.middle |
2129 fig3.b.legend.plot.x) + plot_layout(widths = c(1, 0.2, 0.2, 1, 0.1) ,
2130 guides = "collect") + guides(fill = guide_legend(ncol=2))
2131
2132 fig3.final <- cowplot::ggdraw(fig3) +
2133 cowplot::draw_label("A", x = 0.005, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 12) +
2134 cowplot::draw_label("B", x = 0.25, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 12)
2135
2136 #Plot interactive full panel
2137 fig3.final.interactive <- ggiraph::girafe(ggobj = fig3.final,
2138 width_svg = 12, height_svg = 8,
2139 options = list(ggiraph::opts_hover(css = "fill: yellow; stroke: yellow; r: 2px;")))
2140
2141 fig3.final.interactive7. Infection assays
Experimental infection establishes strain-specific virulence in
queen larvae.
Associations observed during the outbreak
cannot by themselves establish causality. We therefore challenged larvae
experimentally with the reference M. plutonius strain and the
queen-derived qk206 isolate, alone or in combination with BQCV, to test
strain-specific virulence and the contribution of viral co-infection.
Import survival data for worker and queen larval infection
assays. Set ggplot theme for survival curve graphing.
2142 #---------------------------------
2143 # Load libraries
2144 #---------------------------------
2145
2146 library(survival)
2147 library(survminer)
2148 library(dplyr)
2149 library(readr)
2150 library(readxl)
2151 library(ggplot2)
2152
2153 #--------------------------------------
2154 # Read in worker survival data
2155 #--------------------------------------
2156 df.workers <- readr::read_tsv("data/infection_survival_data.tsv") %>% filter(grepl("worker", caste)) %>%
2157 dplyr::select(Day, `1_Control`, `3_ATCC`, `5_qk206`)
2158
2159 df.workers.melt <- reshape2::melt(df.workers, id.vars=c("Day")) %>%
2160 mutate(sample_id=paste0("sample_", row_number())) %>%
2161 filter(!is.na(value)) %>% dplyr::rename(group=variable) %>%
2162 arrange(group) %>%
2163 mutate(group=factor(group, levels=c("1_Control", "3_ATCC", "5_qk206")))
2164
2165 #--------------------------------------
2166 # Read in queen survival data
2167 #--------------------------------------
2168 df.queens <- readr::read_tsv("data/infection_survival_data.tsv") %>% filter(grepl("queen", caste)) %>% dplyr::select(-caste)
2169
2170 df.queens.melt <- reshape2::melt(df.queens, id.vars=c("Day")) %>% mutate(sample_id=paste0("sample_", row_number())) %>%
2171 filter(!is.na(value)) %>% dplyr::rename(group=variable) %>%
2172 mutate(group=factor(group, levels=c("1_Control",
2173 "2_BQCV",
2174 "3_ATCC",
2175 "4_ATCC+BQCV",
2176 "5_qk206",
2177 "6_qk206+BQCV")))
2178
2179 #---------------------------------
2180 # Set surival plotting themes
2181 #---------------------------------
2182 surv_x_scale <- scale_x_continuous(limits = c(0, 6),
2183 breaks = seq(0, 6, by = 1))
2184 surv_y_scale <- scale_y_continuous(limits = c(0, 1),
2185 breaks = seq(0, 1, by = 0.25),
2186 labels = c("0", "", "0.5", "", "1.00"))
2187
2188 surv.theme <- theme_classic() + theme(legend.position="none",
2189 plot.title = element_text(size=10, color="black"),
2190 axis.ticks=element_line(linewidth=0.3),
2191 axis.line = element_line(linewidth=0.3))Worker infection survival
Worker larval infections provide a comparative virulence model for the reference and queen-derived M. plutonius strains. Differences in survival establish that pathogenic potential varies substantially among bacterial isolates.
2192 #--------------------------------------
2193 # Worker larval infection statistics
2194 #--------------------------------------
2195 #Overall
2196 survdiff(Surv(Day, value) ~ group, data = df.workers.melt)## Call:
## survdiff(formula = Surv(Day, value) ~ group, data = df.workers.melt)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## group=1_Control 12 3 6.49 1.8759 2.808
## group=3_ATCC 34 16 16.87 0.0449 0.136
## group=5_qk206 11 8 3.64 5.2202 6.988
##
## Chisq= 8.3 on 2 degrees of freedom, p= 0.022197 #Equivalent to Mantel-Cox
2198 res1 <- survminer::pairwise_survdiff(Surv(Day, value) ~ group, data = df.workers.melt, p.adjust.method = "BH", rho = 0)
2199 #signif(res1$p.value, 6)
2200 res1##
## Pairwise comparisons using Log-Rank test
##
## data: df.workers.melt and group
##
## 1_Control 3_ATCC
## 3_ATCC 0.207 -
## 5_qk206 0.039 0.040
##
## P value adjustment method: BH2201 # Endpoint stats at D6
2202 fit1 <- survival::survfit(Surv(Day, value) ~ group, data = df.workers.melt )
2203 fit1.x <- summary(fit1, times = 6)
2204 fit1.df <- data.frame(group = fit1.x$strata,
2205 time = fit1.x$time,
2206 survival = fit1.x$surv,
2207 err=fit1.x$std.err,
2208 lower_95_CI = fit1.x$lower,
2209 upper_95_CI = fit1.x$upper)
2210 fit1.df## group time survival err lower_95_CI upper_95_CI
## 1 group=1_Control 6 0.7500000 0.12500000 0.5409964 1.0000000
## 2 group=3_ATCC 6 0.5294118 0.08560081 0.3856226 0.7268164
## 3 group=5_qk206 6 0.2727273 0.13428163 0.1039025 0.7158652Queen infection survival
Experimental infection of queen larvae directly tests whether the queen-derived qk206 isolate is sufficient to reproduce the lethal phenotype observed during the outbreak and whether BQCV modifies disease severity.
2211 #--------------------------------------
2212 # Queen larval infection statistics
2213 #--------------------------------------
2214 #Overall
2215 survdiff(Surv(Day, value) ~ group, data = df.queens.melt)## Call:
## survdiff(formula = Surv(Day, value) ~ group, data = df.queens.melt)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## group=1_Control 59 10 20.0 4.97376 6.77768
## group=2_BQCV 60 9 20.2 6.18173 8.44166
## group=3_ATCC 62 13 20.5 2.77314 3.80332
## group=4_ATCC+BQCV 61 20 19.7 0.00619 0.00842
## group=5_qk206 67 31 19.8 6.27466 8.60998
## group=6_qk206+BQCV 69 37 19.8 14.87173 20.46264
##
## Chisq= 40.3 on 5 degrees of freedom, p= 1e-072216 #Equivalent to Mantel-Cox
2217 res2 <- survminer::pairwise_survdiff(Surv(Day, value) ~ group, data = df.queens.melt, p.adjust.method = "BH", rho = 0)
2218 #signif(res1$p.value, 6)
2219 res2##
## Pairwise comparisons using Log-Rank test
##
## data: df.queens.melt and group
##
## 1_Control 2_BQCV 3_ATCC 4_ATCC+BQCV 5_qk206
## 2_BQCV 0.79633 - - - -
## 3_ATCC 0.60339 0.46962 - - -
## 4_ATCC+BQCV 0.08099 0.05135 0.21205 - -
## 5_qk206 0.00125 0.00072 0.00612 0.14306 -
## 6_qk206+BQCV 0.00013 0.00010 0.00064 0.02729 0.46962
##
## P value adjustment method: BH2220 # Endpoint stats at D6
2221 fit2 <- survival::survfit(Surv(Day, value) ~ group, data = df.queens.melt )
2222 fit2.x <- summary(fit2, times = 6)
2223 fit2.df <- data.frame(group = fit2.x$strata,
2224 time = fit2.x$time,
2225 survival = fit2.x$surv,
2226 err=fit2.x$std.err,
2227 lower_95_CI = fit2.x$lower,
2228 upper_95_CI = fit2.x$upper)
2229 fit2.df## group time survival err lower_95_CI upper_95_CI
## 1 group=1_Control 6 0.8305085 0.04884499 0.7400858 0.9319789
## 2 group=2_BQCV 6 0.8500000 0.04609772 0.7642862 0.9453265
## 3 group=3_ATCC 6 0.7903226 0.05169900 0.6952212 0.8984332
## 4 group=4_ATCC+BQCV 6 0.6721311 0.06010522 0.5640732 0.8008895
## 5 group=5_qk206 6 0.5373134 0.06091439 0.4302574 0.6710071
## 6 group=6_qk206+BQCV 6 0.4637681 0.06003468 0.3598430 0.5977076Figure 4A
Plot worker larvae survival curves
2230 #--------------------------------------
2231 # Plot survival curves - Workers
2232 #--------------------------------------
2233 pal1 <- c("grey70","#D3AE77","#FF9B9B")
2234
2235 fig4.a <- survminer::ggsurvplot(fit1,
2236 data = df.workers.melt,
2237 pval = TRUE,
2238 risk.table = TRUE,
2239 palette= pal1)
2240
2241 #----------------------------------------------------
2242 # Fig 4A - Interactive
2243 #----------------------------------------------------
2244 fig4.a.surv.df <- fig4.a$plot$data %>%
2245 mutate(surv_d6=ifelse(time==6, surv, 0)) %>%
2246 mutate(stderr_d6=ifelse(time==6, std.err, 0)) %>%
2247 group_by(strata) %>% mutate(surv_d6 = max(surv_d6)) %>% mutate(stderr_d6 = max(stderr_d6)) %>% ungroup()
2248
2249
2250 fig4.a.gg <- ggplot(fig4.a.surv.df, aes(x = time, y = surv,
2251 colour = strata, group = strata)) +
2252 ggiraph::geom_step_interactive(aes(data_id = strata,
2253 tooltip = paste0(
2254 "Treatment: ", group,
2255 "\nSurvival: ", round(surv_d6, 2), "±", round(stderr_d6, 2))),
2256 direction = "hv",
2257 hover_nearest = TRUE,
2258 linewidth = 0.8) +
2259 scale_colour_manual(values = pal1, name="Treatment") +
2260 ylab("Survival probability") + xlab("Time") +
2261 surv_x_scale + surv_y_scale + surv.theme + guides(colour = guide_legend(override.aes = list(linewidth = 0.7))) +
2262 theme(legend.background = element_rect(colour = "black", fill = NA, linewidth = 0.3),
2263 legend.title=element_text(size=8, face="plain"),
2264 legend.text=element_text(size=6), legend.key.height = unit(3, "mm"),
2265 legend.key.width = unit(10, "mm"), legend.position="inside",
2266 legend.justification = c(0, 0), legend.position.inside = c(0.05, 0.05))
2267
2268 #Format risk table
2269 #----------------------------------------------------
2270 fig4.a$table$layers[[1]]$aes_params$size <- 3
2271 fig4.a$table$theme$axis.text.y$size <- 8
2272 fig4.a$table$theme$axis.text.x$size <- 8
2273 fig4.a$table$theme$plot.title$size <- 10
2274 #----------------------------------------------------
2275
2276 fig4.a.gg.interactive <- cowplot::plot_grid(fig4.a.gg,
2277 fig4.a$table + theme(axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8)),
2278 nrow=2, rel_heights=c(0.7, 0.3))
2279
2280
2281 girafe(ggobj = fig4.a.gg.interactive, width_svg = 9, height_svg = 6, options = list(opts_hover(css = "stroke-width:3px;")))Figure 4B
Plot queen larvae survival curves
2282 #--------------------------------------
2283 # Plot survival curves - Queens
2284 #--------------------------------------
2285 pal2 <- c("grey70","#90D3C1","#D3AE77","#C5944E","#FF9B9B", "#FF0000")
2286 linetypes2 <- c("solid", "solid", "solid", "dotted", "solid", "dotted")
2287
2288 fig4.b <- survminer::ggsurvplot(fit2, data = df.queens.melt,
2289 pval = TRUE,
2290 risk.table = TRUE,
2291 palette= pal2,
2292 linetype=linetypes2)
2293
2294 #----------------------------------------------------
2295 # Fig 4B - Interacttive
2296 #----------------------------------------------------
2297 fig4.b.surv.df <- fig4.b$plot$data %>%
2298 mutate(surv_d6=ifelse(time==6, surv, 0)) %>%
2299 mutate(stderr_d6=ifelse(time==6, std.err, 0)) %>%
2300 group_by(strata) %>% mutate(surv_d6 = max(surv_d6)) %>% mutate(stderr_d6 = max(stderr_d6)) %>% ungroup()
2301
2302
2303 fig4.b.gg <- ggplot(fig4.b.surv.df, aes(x = time, y = surv,
2304 colour = strata, linetype = strata, group = strata)) +
2305 ggiraph::geom_step_interactive(aes(data_id = group,
2306 tooltip = paste0(
2307 "Treatment: ", group,
2308 "\nSurvival: ", round(surv_d6, 2), "±", round(stderr_d6, 2))),
2309 direction = "hv",
2310 hover_nearest = TRUE,
2311 linewidth = 0.8) +
2312 scale_colour_manual(values = pal2, name="Treatment") +
2313 scale_linetype_manual(values = linetypes2, name="Treatment") +
2314 ylab("Survival probability") + xlab("Time") +
2315 surv_x_scale + surv_y_scale + surv.theme + guides(colour = guide_legend(override.aes = list(linewidth = 0.7))) +
2316 theme(legend.background = element_rect(colour = "black", fill = NA, linewidth = 0.3),
2317 legend.title=element_text(size=8, face="plain"),
2318 legend.text=element_text(size=6), legend.key.height = unit(3, "mm"),
2319 legend.key.width = unit(10, "mm"), legend.position="inside",
2320 legend.justification = c(0, 0), legend.position.inside = c(0.05, 0.05))
2321
2322 #Format risk table
2323 #----------------------------------------------------
2324 fig4.b$table$layers[[1]]$aes_params$size <- 3
2325 fig4.b$table$theme$axis.text.y$size <- 8
2326 fig4.b$table$theme$axis.text.x$size <- 8
2327 fig4.b$table$theme$plot.title$size <- 10
2328 #----------------------------------------------------
2329
2330 fig4.b.gg.interactive <- cowplot::plot_grid(fig4.b.gg,
2331 fig4.b$table + theme(axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8)),
2332 nrow=2, rel_heights=c(0.7, 0.3))
2333
2334 girafe(ggobj = fig4.b.gg.interactive, width_svg = 9, height_svg = 6, options = list(opts_hover(css = "stroke-width:3px;")))Worker vs queen survival
2335 #---------------------------------------------
2336 # Fig 4A-4B - Combine with patchwork
2337 #---------------------------------------------
2338
2339 fig4.ab <- cowplot::ggdraw() +
2340 cowplot::draw_plot(fig4.a.gg | fig4.b.gg) +
2341 cowplot::draw_label("A", x = 0.02, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2342 cowplot::draw_label("Worker survival ", x = 0.28, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10) +
2343 cowplot::draw_label("B", x = 0.52, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2344 cowplot::draw_label("Queen survival ", x = 0.78, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10)
2345
2346 girafe(ggobj = fig4.ab, options = list(opts_hover(css = "stroke-width:3px;")))Pathogen loads
Pathogen quantification following experimental challenge confirms exposure and allows differences in M. plutonius burden among infection treatments to be related to the observed survival phenotypes.
2347 #---------------------------------------------
2348 # Import pathogen load data
2349 #---------------------------------------------
2350
2351 pathogen.df <- readr::read_tsv("data/infection_pathogen_loads.tsv") %>%
2352 ungroup() %>%
2353 mutate(bacterial = case_when(group %in% c("01_NTC", "02_BQCV") ~ "None",
2354 group %in% c("03_ATCC", "04_ATCC_BQCV") ~ "ATCC",
2355 group %in% c("05_qk206", "06_qk206_BQCV") ~ "qk206"),
2356 viral = case_when(group %in% c("02_BQCV", "04_ATCC_BQCV", "06_qk206_BQCV") ~ "BQCV",
2357 TRUE ~ "None"),
2358 bacterial = factor(bacterial, levels = c("None", "ATCC", "qk206")),
2359 viral = factor(viral, levels = c("None", "BQCV")))Figure 4C
2360 #-------------------------------
2361 # M. plutonius loads
2362 #-------------------------------
2363
2364 #Pairwise estimated marginal means
2365 emmeans.df <- pathogen.df %>% rstatix::emmeans_test(qPCR_DNA_Mplut ~ group, p.adjust.method = "BH", detailed = TRUE) %>%
2366 mutate(p.adj = round(p.adj, 4)) %>% mutate(p.adj = ifelse(p.adj==0, "<0.0001", p.adj)) %>%
2367 rstatix::add_xy_position() %>%
2368 mutate(y.position=y.position*0.5) %>%
2369 mutate(cld_group = paste0(group1, "-", group2))
2370
2371 cld.output.mplut <- setNames(emmeans.df$p.adj, emmeans.df$cld_group)
2372 cld.letters.mplut <- multcompView::multcompLetters(cld.output.mplut)$Letters
2373 pathogen.df <- pathogen.df %>% mutate(letters.mplut = dplyr::recode(group, !!!cld.letters.mplut))
2374
2375 fig4.c <- ggplot(pathogen.df, aes(x=group, y=qPCR_DNA_Mplut)) +
2376 stat_summary(fun = mean, geom = "col", width = 0.9, aes(fill = group)) +
2377 stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1), geom = "errorbar", width = 0.1, linewidth = 0.2) +
2378 ggiraph::geom_point_interactive(aes(data_id = sample_id,
2379 tooltip = paste0("Group: ", group,
2380 "\nIndividual: ", sample_id,
2381 "\nValue: ", round(qPCR_DNA_Mplut, 1))),
2382 position = ggbeeswarm::position_quasirandom(width = 0.2, varwidth = TRUE),
2383 alpha = 0.6, colour = "grey10", hover_nearest = TRUE, size = 1) +
2384 scale_fill_manual(values=c( "01_NTC"="grey70",
2385 "02_BQCV"="#90D3C1" ,
2386 "03_ATCC"="#FFF2CC",
2387 "04_ATCC_BQCV"= "#FFD966",
2388 "05_qk206"="#FF9B9B",
2389 "06_qk206_BQCV"= "#FF0000")) +
2390 ylab("M. plutonius (log10 copies)") +
2391 geom_text(data=pathogen.df, aes(y=max(qPCR_DNA_Mplut)*1.01,
2392 label = letters.mplut, group = group), vjust = 0, size=3, color="#2E2E2E",angle = 0) +
2393 guides(fill = guide_legend(title = "Group")) +
2394 ylim(0,8) +
2395 geom_hline(yintercept=2, linetype="dashed", lwd=0.7, color="black", alpha=0.5)+
2396 theme_classic() + theme(legend.position="none",
2397 plot.title = element_text(size=10, color="black"),
2398 axis.title.x=element_blank(),
2399 axis.ticks=element_line(linewidth=0.3),
2400 axis.line = element_line(linewidth=0.3),
2401 axis.text.x=element_text(angle=45, hjust=1))
2402
2403 girafe(ggobj = fig4.c, width=10, height=6, options = list(opts_hover(css = "stroke-width:3px;")))Figure 4D
BQCV abundance was quantified across treatment groups to determine how viral infection behaves during single and mixed infections and whether bacterial co-infection alters viral burden.
2404 #-------------------------------
2405 # BQCV loads
2406 #-------------------------------
2407 #Pairwise estimated marginal means
2408 emmeans.df <- pathogen.df %>% rstatix::emmeans_test(qPCR_RNA_BQCV ~ group, p.adjust.method = "BH", detailed = TRUE) %>%
2409 rstatix::add_xy_position() %>%
2410 mutate(y.position=y.position*0.5) %>%
2411 mutate(cld_group = paste0(group1, "-", group2))
2412
2413 cld.output.bqcv <- setNames(emmeans.df$p.adj, emmeans.df$cld_group)
2414 cld.letters.bqcv <- multcompView::multcompLetters(cld.output.bqcv)$Letters
2415 pathogen.df <- pathogen.df %>% mutate(letters.bqcv = dplyr::recode(group, !!!cld.letters.bqcv))
2416
2417 fig4.d <- ggplot(pathogen.df, aes(x=group, y=qPCR_RNA_BQCV)) +
2418 stat_summary(fun = mean, geom = "col", width = 0.9, aes(fill = group)) +
2419 stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1),
2420 geom = "errorbar", width = 0.1, linewidth = 0.2) +
2421 ggiraph::geom_point_interactive(aes(data_id = sample_id,
2422 tooltip = paste0("Group: ", group,
2423 "\nIndividual: ", sample_id,
2424 "\nValue: ", round(qPCR_DNA_Mplut, 1))),
2425 position = ggbeeswarm::position_quasirandom(width = 0.2, varwidth = TRUE),
2426 alpha = 0.6, colour = "grey10", hover_nearest = TRUE, size = 1) +
2427 scale_fill_manual(values=c( "01_NTC"="grey70",
2428 "02_BQCV"="#90D3C1" ,
2429 "03_ATCC"="#FFF2CC",
2430 "04_ATCC_BQCV"= "#FFD966",
2431 "05_qk206"="#FF9B9B",
2432 "06_qk206_BQCV"= "#FF0000")) +
2433 ylab("BQCV (log10 copies)") +
2434 geom_text(data=pathogen.df, aes(y=max(qPCR_RNA_BQCV)*1.01,
2435 label = letters.bqcv, group = group), vjust = 0, size=3, color="#2E2E2E",angle = 0) +
2436 guides(fill = guide_legend(title = "Group")) +
2437 ylim(0,8) +
2438 geom_hline(yintercept=2, linetype="dashed", lwd=0.7, color="black", alpha=0.5)+
2439 theme_classic() + theme(legend.position="none",
2440 plot.title = element_text(size=10, color="black"),
2441 axis.title.x=element_blank(),
2442 axis.ticks=element_line(linewidth=0.3),
2443 axis.line = element_line(linewidth=0.3),
2444 axis.text.x=element_text(angle=45, hjust=1))
2445
2446 girafe(ggobj = fig4.d, width=10, height=6, options = list(opts_hover(css = "stroke-width:3px;")))Weights of survivors
Larval weight provides a complementary measure of disease severity among individuals that survive infection, capturing developmental costs that are not apparent from survival alone.
2447 #---------------------------------------------
2448 # Import weight data
2449 #---------------------------------------------
2450
2451 weights.df <- readr::read_tsv("data/infection_weights.tsv") %>%
2452 ungroup() %>%
2453 mutate(bacterial = case_when(group %in% c("01_NTC", "02_BQCV") ~ "None",
2454 group %in% c("03_ATCC", "04_ATCC_BQCV") ~ "ATCC",
2455 group %in% c("05_qk206", "06_qk206_BQCV") ~ "qk206"),
2456 viral = case_when(group %in% c("02_BQCV", "04_ATCC_BQCV", "06_qk206_BQCV") ~ "BQCV",
2457 TRUE ~ "None"),
2458 bacterial = factor(bacterial, levels = c("None", "ATCC", "qk206")),
2459 viral = factor(viral, levels = c("None", "BQCV")))
2460
2461 #-------------------------------
2462 # Weight comparisons
2463 #-------------------------------
2464
2465 emmeans.df <- weights.df %>% rstatix::emmeans_test(weight ~ group, p.adjust.method = "BH", detailed = TRUE) %>%
2466 mutate(p.adj = round(p.adj, 4)) %>% mutate(p.adj = ifelse(p.adj==0, "<0.0001", p.adj)) %>%
2467 rstatix::add_xy_position() %>%
2468 mutate(y.position=y.position*0.5) %>%
2469 mutate(cld_group = paste0(group1, "-", group2))
2470
2471 cld.output.weight <- setNames(emmeans.df$p.adj, emmeans.df$cld_group)
2472 cld.letters.weight <- multcompView::multcompLetters(cld.output.weight)$Letters
2473 weights.df <- weights.df %>% mutate(letters.weight = dplyr::recode(group, !!!cld.letters.weight))
2474
2475 #Summary stats
2476 weight.stats <- weights.df %>% group_by(group) %>%
2477 dplyr::summarise(n = sum(!is.na(weight)),
2478 mean = mean(weight, na.rm = TRUE),
2479 SE = sd(weight, na.rm = TRUE) / sqrt(n))
2480 weight.stats## # A tibble: 6 × 4
## group n mean SE
## <chr> <int> <dbl> <dbl>
## 1 01_NTC 9 271. 16.2
## 2 02_BQCV 9 247. 8.99
## 3 03_ATCC 9 233. 13.1
## 4 04_ATCC_BQCV 9 173. 13.5
## 5 05_qk206 9 121. 15.4
## 6 06_qk206_BQCV 9 111. 14.5Figure 4E
2481 fig4.e <- ggplot(weights.df, aes(x=group, y=weight)) +
2482 stat_summary(fun = mean, geom = "col", width = 0.9, aes(fill = group)) +
2483 stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1),
2484 geom = "errorbar", width = 0.1, linewidth = 0.2) +
2485 ggiraph::geom_point_interactive(aes(data_id = sample_id,
2486 tooltip = paste0("Group: ", group,
2487 "\nIndividual: ", sample_id,
2488 "\nWeight: ", round(weight, 1), " mg")),
2489 position = ggbeeswarm::position_quasirandom(width = 0.2, varwidth = TRUE),
2490 alpha = 0.6, colour = "grey10", hover_nearest = TRUE, size = 1) +
2491 #ggbeeswarm::geom_quasirandom(varwidth = TRUE, width= 0.2, alpha = 0.6, color = "grey10", size = 1) +
2492 scale_fill_manual(values=c( "01_NTC"="grey70",
2493 "02_BQCV"="#90D3C1" ,
2494 "03_ATCC"="#FFF2CC",
2495 "04_ATCC_BQCV"= "#FFD966",
2496 "05_qk206"="#FF9B9B",
2497 "06_qk206_BQCV"= "#FF0000")) +
2498 ylab("Weight (mg)") +
2499 geom_text(data=weights.df, aes(y=max(weight)*1.01, label = letters.weight, group = group),
2500 vjust = 0, size=3, color="#2E2E2E",angle = 0) +
2501 guides(fill = guide_legend(title = "Group")) +
2502 theme_classic() + theme(legend.position="none",
2503 plot.title = element_text(size=10, color="black"),
2504 axis.title.x=element_blank(),
2505 axis.ticks=element_line(linewidth=0.3),
2506 axis.line = element_line(linewidth=0.3),
2507 axis.text.x=element_text(angle=45, hjust=1))
2508
2509 girafe(ggobj = fig4.e, width=10, height=6, options = list(opts_hover(css = "stroke-width:3px;")))Immune gene expression
Host immune-gene expression was examined to determine whether the different infection treatments generate distinct physiological responses in surviving larvae.
2510 #---------------------------------------------
2511 # Import immune gene expression data
2512 #---------------------------------------------
2513
2514 immune.df <- readr::read_tsv("data/infection_gene_expression.tsv") %>%
2515 ungroup() %>%
2516 mutate(bacterial = case_when(group %in% c("01_NTC", "02_BQCV") ~ "None",
2517 group %in% c("03_ATCC", "04_ATCC_BQCV") ~ "ATCC",
2518 group %in% c("05_qk206", "06_qk206_BQCV") ~ "qk206"),
2519 viral = case_when(group %in% c("02_BQCV", "04_ATCC_BQCV", "06_qk206_BQCV") ~ "BQCV",
2520 TRUE ~ "None"),
2521 bacterial = factor(bacterial, levels = c("None", "ATCC", "qk206")),
2522 viral = factor(viral, levels = c("None", "BQCV")))
2523
2524 #------------------------------------
2525 # Plot immune gene expression PCA
2526 #------------------------------------
2527 df.counts <- immune.df %>% tibble::column_to_rownames("sample_id") %>% dplyr::select(matches("RNA"))
2528 df.meta <- immune.df %>% tibble::column_to_rownames("sample_id")
2529
2530 d.pcx <- prcomp(df.counts, center=TRUE, scale=TRUE) #scale=TRUE to normalize expression for each immune gene
2531 d.mvar <- sum(d.pcx$sdev^2)
2532 # Calculate the PC1 and PC2 variance
2533 PC1 <- paste("PC1: ", round(sum(d.pcx$sdev[1]^2)/d.mvar, 3))
2534 PC2 <- paste("PC2: ", round(sum(d.pcx$sdev[2]^2)/d.mvar, 3))
2535 PC3 <- paste("PC3: ", round(sum(d.pcx$sdev[3]^2)/d.mvar, 3))
2536 #biplot(d.pcx, var.axes=T, scale=0, xlab=PC1, ylab=PC2, cex=c(0.5, 0.5))
2537
2538
2539 LOADINGS<- data.frame(Variables=rownames(d.pcx$rotation), d.pcx$rotation)
2540 VALUES<-merge(df.meta, d.pcx$x[,c(1,2,3)], by.x="row.names", by.y="row.names", all=F) %>% dplyr::rename(sample_id="Row.names")
2541
2542 VALUES.g <- VALUES
2543 hull_data <- VALUES.g %>%
2544 tidyr::drop_na() %>%
2545 group_by(group) %>%
2546 dplyr::slice(chull(PC1, PC2)) %>%
2547 mutate(linetype=ifelse(grepl("01_NTC|02_BQCV|03_ATCC$|05_qk206$", group), "solid", "dotted"))Figure 4F
2548 fig4.f <- ggplot(VALUES.g, aes(x=PC1, y=PC2, color=group, label=sample_id)) +
2549 scale_colour_manual(values=c("01_NTC"="grey70",
2550 "02_BQCV"="#90D3C1" ,
2551 "03_ATCC"="#D3AE77",
2552 "04_ATCC_BQCV"= "#FFD966",
2553 "05_qk206"="#FF9B9B",
2554 "06_qk206_BQCV"= "#FF0000")) +
2555 scale_fill_manual(values=c("01_NTC"="grey70",
2556 "02_BQCV"="#90D3C1" ,
2557 "03_ATCC"="#D3AE77",
2558 "04_ATCC_BQCV"= "#FFD966",
2559 "05_qk206"="#FF9B9B",
2560 "06_qk206_BQCV"= "#FF0000")) +
2561 scale_linetype_manual(values = c("01_NTC" = "solid",
2562 "02_BQCV" = "solid",
2563 "03_ATCC" = "solid",
2564 "04_ATCC_BQCV" = "dashed",
2565 "05_qk206" = "solid",
2566 "06_qk206_BQCV" = "dashed")) +
2567 ggiraph::geom_point_interactive(aes(data_id = sample_id,
2568 tooltip = paste0("Group: ", group,
2569 "\nIndividual: ", sample_id)),
2570 alpha = rep(0.85), hover_nearest = TRUE, size = 1.5) +
2571 geom_segment(data = LOADINGS, aes(x = 0, y = 0, xend = (PC1*4), yend = (PC2*4)),
2572 arrow = arrow(length = unit(5/20, "picas")),
2573 color = "grey40", inherit.aes = FALSE, size=0.1) +
2574 annotate("text", x = (LOADINGS$PC1*4.5), y = (LOADINGS$PC2*4.5),
2575 label = LOADINGS$Variables, size=3, alpha=0.4) +
2576 stat_ellipse(
2577 aes(group = group, colour = group, fill = group, linetype=group),
2578 geom = "polygon",
2579 alpha = 0.1,
2580 type = "norm",
2581 level = 0.68,
2582 segments = 51,
2583 linewidth = 0.6,
2584 show.legend = FALSE) +
2585 xlab(PC1) + ylab(PC2) + theme(legend.position="none") +
2586 theme(panel.border = element_blank(), # element_rect(color ="grey", fill = NA, size = 0.5, linetype = 1),
2587 axis.line = element_line(),
2588 axis.text.x=element_text(angle=0, hjust=1, size=6),
2589 axis.text.y=element_text(size=6),
2590 axis.title.x = element_text(size=8),
2591 axis.title.y = element_text(size=8),
2592 panel.background = element_blank(),
2593 panel.grid.major = element_blank(),
2594 panel.grid.minor = element_blank(),
2595 legend.position = "none")
2596
2597
2598 girafe(ggobj = fig4.f, width=10, height=6, options = list(opts_hover(css = "stroke-width:3px;")))Figure 4 (Full panel)
Taken together, the infection experiments move the study from association to causality. Virulence is strongly strain dependent, the queen-derived qk206 isolate can reproduce severe disease in queen larvae, and BQCV acts as a context-dependent modifier rather than a sufficient cause on its own.
2599 #---------------------------------------------------------
2600 # Fig 4 (Full panel)
2601 #---------------------------------------------------------
2602
2603 #Add labels for Fig 4C-D
2604 fig4.cd <- cowplot::ggdraw() + cowplot::draw_plot(fig4.c + fig4.d) +
2605 cowplot::draw_label("C", x = 0.02, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2606 cowplot::draw_label("M. plutonius loads ", x = 0.28, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10) +
2607 cowplot::draw_label("D", x = 0.52, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2608 cowplot::draw_label("BQCV loads ", x = 0.78, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10)
2609
2610 #Add labels for Fig 4E-F
2611 fig4.ef <- cowplot::ggdraw() + cowplot::draw_plot(fig4.e + fig4.f) +
2612 cowplot::draw_label("E", x = 0.02, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2613 cowplot::draw_label("Weights of survivors ", x = 0.28, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10) +
2614 cowplot::draw_label("F", x = 0.52, y = 0.995, hjust = 0, vjust = 1, fontface = "bold", size = 10) +
2615 cowplot::draw_label("Immune gene expression ", x = 0.78, y = 0.995, hjust = 0.5, vjust = 1, fontface = "bold", size = 10)
2616
2617 fig4.final <- fig4.ab / fig4.cd / fig4.ef + plot_layout(heights = c(1.1, 1, 1)) #, guides = "collect"
2618
2619 girafe(ggobj = fig4.final, width=9, height=12, options = list(opts_hover(css = "stroke: blue; stroke-width:3px;")))8. Genome comparisons
Genome comparisons place queen-derived M. plutonius within
the broader species population
The strong difference in
virulence between M. plutonius strains prompted us to investigate their
genomic relationships. Population-genetic and pangenomic analyses were
used to determine where the queen-derived isolates fall within known M.
plutonius diversity and whether they possess distinctive patterns of
gene content.
MLST goeBURST plot
The following code generates a minimum spanning tree using the PHYLOViZ (https://github.com/phyloviz) goeBURST algorithm - a java implementation of the eBURST algorithm proposed by Feil et al (https://doi.org/10.1128/jb.186.5.1518-1530.2004) that ensures an optimal solution for the placement of links between Sequence Types.
MLST profiles of M. plutonius strains were obtained from pubMLST database (https://pubmlst.org; accessed on 05-June-2026) and combined with the MLST profiles determined for the isolates in this study.
2620 library(dplyr)
2621 library(readr)
2622 library(stringr)
2623 library(igraph)
2624 library(ggraph)
2625 library(tidyr)
2626 library(ggraph)
2627 library(scatterpie)
2628 library(scales)
2629 library(V8)
2630 library(jsonlite)
2631
2632 #----------------------------------------
2633 # Load dataframe
2634 #----------------------------------------
2635 df.mlst <- readxl::read_xlsx("data/pubMLST_results_June5_formatted_for_phyloviz.xlsx", sheet="ALL_july2026") %>%
2636 filter(blanks<1) %>% arrange(ST) %>% mutate(caste=ifelse(is.na(caste), "unknown", caste))
2637
2638 df.mlst.g.CC3 <- df.mlst %>% group_by(ST) %>% dplyr::slice(1) %>%
2639 mutate(ST=paste0("ST", ST)) %>% filter(grepl("3|13", `clonal complex`))
2640 df.mlst.g.CC12 <- df.mlst %>% group_by(ST) %>% dplyr::slice(1) %>%
2641 mutate(ST=paste0("ST", ST)) %>% filter(grepl("12", `clonal complex`))
2642 df.mlst.g <- df.mlst %>% group_by(ST) %>% dplyr::slice(1) %>% mutate(ST=paste0("ST", ST))
2643
2644 #----------------------------------------------------------
2645 #Source the goeBURST algorithm and layout functions
2646 #----------------------------------------------------------
2647 source("scripts/function_goeBURST_algorithm.R")
2648 source("scripts/function_goeBURST_layout.R")
2649
2650 #----------------------------------------------------------
2651 #Download javascript code from PHYLOViZ github:
2652 #----------------------------------------------------------
2653 viva.file <- "scripts/vivagraph_julbra.js"
2654 if (!file.exists(viva.file)) {download.file(url = paste0("https://raw.githubusercontent.com/", "bfrgoncalves/Online-PhyloViZ/master/", "public/javascripts/dependencies/vivagraph_julbra.js"), destfile = viva.file, mode = "wb")}
2655
2656 #------------------------------------------------
2657 #Run goeBURST algorithm and apply layout
2658 #------------------------------------------------
2659 result <- goeBURST_algorithm(profiles=df.mlst.g, id_col = "ST",
2660 locus_cols=c("argE","galK","gbpB","purR"))
2661 set.seed(3)
2662 pv.layout <- goeBURST_layout(graph = result$graph,
2663 viva_file = viva.file,
2664 distance_attribute = "distance")
2665
2666 coords <- pv.layout$coordinates
2667 caste.columns <- unique(df.mlst$caste)
2668 caste.colours <- c("worker" = "#AEC7E8",
2669 "unknown" = "#DEEBF7",
2670 "drone" = "#FFBB78",
2671 "queen" = "#C59EE2")
2672
2673 index.sum <- df.mlst %>%
2674 mutate(ST_plot = paste0("ST", ST),
2675 caste = ifelse(is.na(caste) | caste == "", "unknown", caste)) %>%
2676 count(ST_plot, caste, name = "n")
2677
2678 node.pie.df <- index.sum %>% pivot_wider(names_from = caste, values_from = n, values_fill = 0)
2679
2680 #------------------------------------------------# Interactive graph
2681 xy <- as.matrix(coords[, c("x", "y")])
2682
2683 #tk.id <- igraph::tkplot(result$graph, layout = xy,
2684 # vertex.label = igraph::V(result$graph)$name,
2685 # vertex.size = 15, edge.color = "grey70")
2686 #Save coordinates
2687 #xy.new <- igraph::tk_coords(tk.id, norm = FALSE)
2688 #saveRDS(xy.new, "data/goeBURST_coords.RDS")
2689 #igraph::tk_close(tk.id)
2690 xy.new <- readRDS("data/goeBURST_coords.RDS")
2691 coords$x <- xy.new[, 1]
2692 coords$y <- xy.new[, 2]
2693
2694 #------------------------------------------------#Reproduce layout
2695 lay <- create_layout(result$graph, layout = "manual",
2696 x = coords$x,
2697 y = coords$y)
2698 layout.df <- as.data.frame(lay) %>%
2699 left_join(node.pie.df, by = c("name" = "ST_plot")) %>%
2700 mutate(across(all_of(caste.columns), ~ tidyr::replace_na(.x, 0)))
2701
2702 edge.ends <- igraph::ends(result$graph, igraph::E(result$graph),names = FALSE)
2703 xy <- as.matrix(coords[, c("x", "y")])
2704 edge.lengths <- sqrt(rowSums((xy[edge.ends[, 1], , drop = FALSE] - xy[edge.ends[, 2], , drop = FALSE])^2 ))
2705
2706 node.radius <- median(edge.lengths) * 0.08
2707 layout.df$r <- node.radius
2708
2709 node.radius <- median(edge.lengths) * 0.35
2710 layout.df$r <- node.radius
2711
2712 #------------------------------------------------
2713 fig5.b1 <- ggraph(lay) + geom_edge_link(colour = "grey70", linewidth = 0.6) +
2714 scatterpie::geom_scatterpie( data = layout.df, aes(x = x,y = y, r = r ),
2715 cols = caste.columns, colour = NA, linewidth = 0.3) +
2716 scale_fill_manual(values = caste.colours, breaks = caste.columns, drop = FALSE) +
2717 geom_text(data = layout.df, aes(x = x,y = y,label = name), size = 2) +
2718 coord_equal( clip = "off") + theme_void() +labs( fill = "Caste") +
2719 theme(legend.position="none") #+ coord_flip()
2720 fig5.b1
2721 #----------------------------------------
2722 # Donut plot of MLST-caste relationship
2723 #----------------------------------------
2724 donut.df <- df.mlst %>%
2725 mutate(caste = ifelse(is.na(caste) | trimws(caste) == "", "unknown", trimws(caste) )) %>%
2726 count(caste, name = "n") %>%
2727 mutate(proportion = n / sum(n), label = paste0(caste, "\n", scales::percent(proportion)))
2728
2729 total_n <- sum(donut.df$n)
2730
2731 fig5.b2 <- ggplot( donut.df, aes(x = 2, y = n, fill = caste)) +
2732 geom_col(width = 1, colour = "white") + coord_polar(theta = "y") + xlim(0.25, 2.5) +
2733 annotate("text", x = 0.25, y = sum(donut.df$n) / 2,
2734 label = paste0("Total\nN = ", total_n),
2735 size = 4, fontface = "bold", lineheight = 1.1) +
2736 scale_fill_manual( values = c("worker" = "#AEC7E8",
2737 "unknown" = "#DEEBF7",
2738 "drone" = "#FFBB78",
2739 "queen" = "#C59EE2")) +
2740 theme_void() +
2741 labs(fill = "Caste")
2742 fig5.b2
2743 #----------------------------
2744 # Fig 5B - Combine plots
2745 #----------------------------
2746 fig5.b.final <- fig5.b1 + patchwork::inset_element(fig5.b2, left = 0, bottom = 0.5, right = 0.5, top = 1)
2747 fig5.b.final
Pangenome analysis
Gene presence–absence profiles were compared across available M. plutonius genomes to identify conserved and variable regions of the pangenome and to determine whether the queen-derived isolates possess distinctive accessory-gene complements.
2748 library(dplyr)
2749 library(Biostrings)
2750 library(stringr)
2751 library(reshape2)
2752 library(ggplot2)
2753 library(heatmap3)
2754 library(ggheatmap)
2755 library(readr)
2756 library(ggdendro)
2757 library(cowplot)
2758 library(dendsort)
2759 library(ggiraph)
2760 library(ggordiplots)
2761 library(ggtree)
2762 library(ggtreeExtra)
2763 library(phyloseq)
2764 library(ape)
2765 library(phangorn)
2766
2767 metadata.genomes <- readr::read_csv("data/genome_information.csv")
2768 pano.mer.wide.meta2 <- readr::read_tsv("data/pangenome_matrix.tsv") %>% tibble::column_to_rownames("combined")
2769
2770 #-------------------------------------
2771 # Gene presence/absence comparisons
2772 #-------------------------------------
2773 df.heat.meta <- pano.mer.wide.meta2 %>% tibble::rownames_to_column("combined") %>%
2774 dplyr::rowwise() %>%
2775 mutate(mean_BCQMy = mean(c_across(matches("BCQMy204|BCQMy205|BCQMy206")))) %>%
2776 mutate(mean_BCQMy2 = mean(c_across(matches("BCQMy223|BCQMy224|BCQMy225")))) %>%
2777 mutate(mean_BCQMy_both = mean(c_across(matches("BCQMy2")))) %>%
2778 mutate(mean_atypical = mean(c_across(matches("_atypical")))) %>%
2779 mutate(mean_typical = mean(c_across(matches("_typical")))) %>%
2780 mutate(mean_others = mean(c_across(!matches("BCQMy|combined|mean_BCQMy")))) %>%
2781 mutate(mean_all = mean(c_across(!matches("mean_BCQMy|mean_others|combined")))) %>%
2782 mutate(var = var(c_across(!matches("mean_|combined")))) %>%
2783 mutate(var_all = var(c_across(!matches("mean_|var|combined")))) %>%
2784 ungroup() %>% tibble::column_to_rownames("combined") %>%
2785 mutate(mean_diff = (mean_BCQMy - mean_others)) %>%
2786 mutate(mean_diff_both = (mean_BCQMy_both - mean_others)) %>%
2787 mutate(mean_diff_typical = (mean_typical - mean_atypical))
2788
2789
2790 #Gather top orthologous genes and then dereplicate
2791 df.heat.meta1 <- df.heat.meta %>%
2792 arrange(desc(abs(mean_diff))) %>% filter(var!=0) %>% dplyr::slice(1:50)
2793 df.heat.meta2 <- df.heat.meta %>%
2794 arrange(desc(abs(mean_diff_both))) %>% filter(var!=0) %>% dplyr::slice(1:50)
2795 df.heat.meta3 <- df.heat.meta %>%
2796 arrange(desc(mean_diff_typical)) %>% filter(var!=0) %>% dplyr::slice(1:50)
2797 df.heat.meta4 <- df.heat.meta %>%
2798 filter(grepl("_vMP",rownames(.)))
2799 df.heat.meta5 <- df.heat.meta %>%
2800 filter(grepl("group_4440|group_4618|group_4600|group_5551",rownames(.)))
2801 df.heat.meta6 <- df.heat.meta %>% filter(grepl("group_4998|group_5369|group_4264|group_5508|group_4273|group_5017|group_4238|group_4575|group_3751|group_4083|group_506|group_4188|group_2237|group_3853|group_3654|group_5098|group_4934|group_4731|group_4156|group_4052|group_3741|group_5437|group_5425|group_5297|group_3773",rownames(.)))
2802
2803 #==================================================== For plotting
2804 genes.to.keep <- unique(c(rownames(df.heat.meta1),
2805 rownames(df.heat.meta2),
2806 rownames(df.heat.meta3),
2807 rownames(df.heat.meta4),
2808 rownames(df.heat.meta5),
2809 rownames(df.heat.meta6)))
2810 #====================================================
2811 gene.list <- str_split_fixed(rownames(df.heat.meta %>% filter(mean_diff==1)), "___", 2)[,1]
2812
2813 df.heat <- pano.mer.wide.meta2 %>% filter(rownames(.) %in% genes.to.keep) %>%
2814 sjmisc::rotate_df() %>% `colnames<-`(paste0(str_split_fixed(colnames(.), "___", 7)[,5], "___",
2815 str_split_fixed(colnames(.), "___", 7)[,1], "___",
2816 str_split_fixed(colnames(.), "___", 7)[,4], "___",
2817 str_split_fixed(colnames(.), "___", 7)[,3], "___",
2818 str_split_fixed(colnames(.), "___", 7)[,2]))
2819
2820 #Cluster samples/genes
2821 cluster.tax <- hclust(dist(t(df.heat), method="euclidian"), method="ward.D2")
2822 cluster.tax <- dendsort(cluster.tax, type="average") #min / average
2823 cluster.samples <- hclust(dist(df.heat, method="euclidian"), method="ward.D2")
2824
2825 #Reorder samples
2826 cluster.samples.phylo <- as.phylo(cluster.samples)
2827
2828 #Generate dendrograms for heatmap
2829 ggtree.samples <- ggtree(cluster.samples.phylo, ladderize=FALSE)
2830 side.tree <- ggtree.samples
2831 sample.order <- side.tree$data %>% dplyr::filter(isTip) %>% dplyr::arrange((y)) %>% dplyr::pull(label)
2832
2833 side.tree <- ggtree::rotate(ggtree.samples, 41)
2834 sample.order <- side.tree$data %>% dplyr::filter(isTip) %>% dplyr::arrange(y) %>% dplyr::pull(label)
2835
2836 #Reorder taxa
2837 cluster.tax.phylo <- as.phylo(cluster.tax)
2838 ggtree.taxa <- ggtree(cluster.tax.phylo, ladderize=FALSE)
2839 g2.dendro.tax <- ggtree.taxa
2840 tax.order <- g2.dendro.tax$data %>% dplyr::filter(isTip) %>% dplyr::arrange(y) %>% dplyr::pull(label)
2841
2842 #Base pangenome heatmap
2843 pangenome.melt <- melt(df.heat %>% tibble::rownames_to_column("combined"), id.vars=c("combined")) %>%
2844 mutate(variable_factor = factor(variable, levels=tax.order)) %>%
2845 mutate(combined_factor= factor(combined, levels=sample.order)) %>%
2846 arrange(combined_factor) %>%
2847 mutate(combined_factor=gsub("___", " (", combined_factor)) %>%
2848 mutate(combined_factor=paste0(combined_factor, ")")) %>%
2849 mutate(combined_factor= factor(combined_factor, levels=unique(combined_factor))) %>%
2850 arrange(variable_factor) %>%
2851 mutate(phage=str_split_fixed(variable_factor, "___", 2)[,1]) %>%
2852 mutate(phage=ifelse(phage=="NA", "not_phage", phage)) %>%
2853 mutate(genome=str_split_fixed(combined,"___", 3)[,1]) %>%
2854 mutate(strain=str_split_fixed(combined,"___", 3)[,2])
2855
2856 pangenome.melt.distinct <- pangenome.melt %>% distinct(variable_factor, .keep_all=TRUE)
2857 pangeome.heatmap <- ggplot(data = pangenome.melt, aes(x=variable_factor, y=combined_factor, fill=value)) +
2858 geom_tile(colour="white", linewidth=0.1) +
2859 geom_tile_interactive(aes(tooltip = paste0("Genome accession: ", .data[["genome"]],
2860 "\nStrain: ", .data[["strain"]],
2861 "\nOrtholog group: ", .data[["variable"]],
2862 "\nOrtholog presence : ", ifelse(value==1, "Present", "Absent"),
2863 "\nPhage-associated : ", .data[["phage"]]),
2864 data_id =paste0(combined, "_", variable)),
2865 colour="white", linewidth=0.1, hover_nearest = TRUE) +
2866 scale_color_gradientn(colours=colorRampPalette(c("grey90","#8FAADC"))(100),
2867 aesthetics = "fill") +
2868 scale_y_discrete(position = "right") +
2869 theme(axis.title.x=element_blank(),
2870 axis.ticks.x=element_blank(),
2871 axis.text.x=element_blank(), #axis.text.x=element_text(size=2, angle=90, hjust=1),
2872 axis.text.y=element_text(size=7, hjust=1))
2873
2874 #---------------------------------------------------------------------
2875 # Generate genome info indices
2876 #---------------------------------------------------------------------
2877 metadata.genomes <- metadata.genomes %>% #readr::read_csv("../metadata.genomes.txt") %>%
2878 mutate(plasmid=case_when(pmp1=="Y" & pmp19=="Y" ~ "Both",
2879 pmp1=="Y" & pmp19=="N" ~ "pMP1",
2880 pmp1=="N" & pmp19=="Y" ~ "pMP19",
2881 ))
2882
2883 metadata.genomes.index <- setNames(metadata.genomes$location, metadata.genomes$genome)
2884 metadata.genomes.index2 <- setNames(metadata.genomes$typical_or_atypical, metadata.genomes$genome)
2885 metadata.genomes.index3 <- setNames(metadata.genomes$pmp19, metadata.genomes$genome)
2886 metadata.genomes.index4 <- setNames(metadata.genomes$oxytet_MIC, metadata.genomes$genome)
2887
2888 #-----------------------------------------
2889 #HORIZONTAL strip legend #1
2890 #-----------------------------------------
2891 pangenome.melt.horizontal <- pangenome.melt %>% arrange(variable_factor) %>%
2892 mutate(phage=str_split_fixed(variable_factor, "___", 2)[,1]) %>%
2893 mutate(phage=ifelse(phage=="NA", "not_phage", phage)) %>%
2894 distinct(variable_factor, .keep_all=TRUE) %>%
2895 mutate(combined_factor=factor(combined_factor, levels=unique(.$combined_factor)))
2896
2897 annotation_plot.horizontal <- ggplot(data = pangenome.melt.horizontal, aes(x=variable_factor, y=1, fill = phage)) +
2898 geom_tile(height = 0.1, aes()) +
2899 scale_fill_manual(values=c(vMP1="#61D6FF",
2900 vMP2="#7FB1DF",
2901 vMP3="#F8F695",
2902 vMP4="#F9CBBD",
2903 vMP5="#FAA0FC",
2904 not_phage="grey60")) +
2905 scale_y_continuous(position = "right") +# coord_flip()+
2906 theme_void() + theme(legend.position = "right") + ylab("Group") +
2907 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2908
2909 #-----------------------------------------
2910 #VERTICAL strip legend #1
2911 #-----------------------------------------
2912 strip.colours.1 <- c( "YES" = "#FAA0FC", "NO" = "#FFE699")
2913
2914 annotation_plot <- ggplot(data = pangenome.melt %>% arrange(combined_factor) %>% distinct(combined, .keep_all=TRUE) %>%
2915 mutate(group=ifelse(grepl("BCQMy", combined), "YES", "NO")) %>%
2916 mutate(combined=factor(combined, levels=unique(.$combined))),
2917 aes(x=combined, y=1, fill = group)) +
2918 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2919 scale_fill_manual(values=strip.colours.1) +
2920 scale_y_continuous(position = "right") + coord_flip()+
2921 theme_void() + theme(legend.position = "right") + ylab("Group") +
2922 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2923
2924 #-----------------------------------------
2925 #VERTICAL strip legend #2
2926 #-----------------------------------------
2927 pangenome.melt.location <- pangenome.melt %>% arrange(combined_factor) %>% distinct(combined, .keep_all=TRUE) %>%
2928 mutate(combined2=str_split_fixed(combined,"___", 2)[,1]) %>%
2929 mutate(location=dplyr::recode(!!!metadata.genomes.index, combined2)) %>%
2930 mutate(combined=factor(combined, levels=unique(.$combined)))
2931
2932 strip.colours.2 <- colorRampPalette(c(RColorBrewer::brewer.pal(12,"Set3")))(length(unique(pangenome.melt.location$location)))
2933
2934 annotation_plot2 <- ggplot(pangenome.melt.location,
2935 aes(x=combined, y=1, fill = location)) +
2936 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2937 scale_fill_manual(values=strip.colours.2) +
2938 scale_y_continuous(position = "right") + coord_flip()+
2939 theme_void() + theme(legend.position = "right") + ylab("Group") +
2940 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2941
2942 #-----------------------------------------
2943 #VERTICAL strip legend #3
2944 #-----------------------------------------
2945 strip.colours.3 <- c( "Atypical" = "#FFC000", "Typical" = "#595959")
2946 pangenome.melt.typical <- pangenome.melt %>% arrange(combined_factor) %>% distinct(combined, .keep_all=TRUE) %>%
2947 mutate(combined2=str_split_fixed(combined,"___", 2)[,1]) %>%
2948 mutate(type=dplyr::recode(!!!metadata.genomes.index2, combined2))%>%
2949 mutate(combined=factor(combined, levels=unique(.$combined)))
2950
2951 annotation_plot3 <- ggplot(data = pangenome.melt.typical ,
2952 aes(x=combined, y=1, fill = type)) +
2953 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2954 scale_fill_manual(values=strip.colours.3) +
2955 scale_y_continuous(position = "right") + coord_flip()+
2956 theme_void() + theme(legend.position = "right") + ylab("Group") +
2957 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2958
2959 #-----------------------------------------
2960 #VERTICAL strip legend #4
2961 #-----------------------------------------
2962 strip.colours.4 <- c( "N" = "#E5E5E5", "Y" = "#C00000")
2963
2964 annotation_plot4 <- ggplot(data = pangenome.melt %>% arrange(combined_factor) %>% distinct(combined, .keep_all=TRUE) %>%
2965 mutate(combined2=str_split_fixed(combined,"___", 2)[,1]) %>%
2966 mutate(plasmid=dplyr::recode(!!!metadata.genomes.index3, combined2)) %>%
2967 mutate(combined=factor(combined, levels=unique(.$combined))),
2968 aes(x=combined, y=1, fill = plasmid)) +
2969 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2970 scale_fill_manual(values=strip.colours.4) +
2971 #scale_fill_manual(values = category_colors, name = "Category") + # Group color scale
2972 scale_y_continuous(position = "right") + coord_flip()+
2973 theme_void() + theme(legend.position = "right") + ylab("Group") +
2974 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10)))
2975
2976 #-----------------------------------------
2977 #VERTICAL strip legend #4
2978 #-----------------------------------------
2979 pangenome.melt.abx <- pangenome.melt %>% arrange(combined_factor) %>% distinct(combined, .keep_all=TRUE) %>%
2980 mutate(combined2=str_split_fixed(combined,"___", 2)[,1]) %>%
2981 mutate(abx=dplyr::recode(!!!metadata.genomes.index4, combined2)) %>%
2982 mutate(abx= ifelse(is.na(abx), 0, abx)) %>%
2983 mutate(abx_group= case_when(as.numeric(abx) >0 & as.numeric(abx) <16 ~ "<16ug/mL",
2984 as.numeric(abx) >=16 ~ ">=16ug/mL",
2985 TRUE ~ "unknown")) %>%
2986 mutate(combined=factor(combined, levels=unique(.$combined)))
2987
2988 strip.colours.4 <- c( "unknown" = "#7D7F81",
2989 "<16ug/mL" = "#1F4468",
2990 ">=16ug/mL" = "#65AEE8")
2991
2992 annotation_plot5 <- ggplot(data = pangenome.melt.abx,
2993 aes(x=combined, y=1, fill = abx_group)) +
2994 geom_tile(height = 0.1, aes()) + # Adjust width to keep sidebar slim
2995 scale_fill_manual(values=strip.colours.4) +
2996 scale_y_continuous(position = "right") + coord_flip()+
2997 theme_void() + theme(legend.position = "right") + ylab("Group") +
2998 theme(axis.title.y.right = element_text(size = 8, hjust=0, face="bold", margin = ggplot2::margin(l = 10, r = 10))) Figure 5C
2999 #-------------------------------
3000 # Combine Fig 5C
3001 #-------------------------------
3002 fig5.c.combined1 <- side.tree + annotation_plot + annotation_plot2 + annotation_plot3 +
3003 annotation_plot4 +annotation_plot5 + pangeome.heatmap + plot_layout(widths = c(0.04,0.03,0.03, 0.03,0.03,0.03, 1)) # Sidebar 5% of plot width
3004 top.legend <- cowplot::ggdraw(g2.dendro.tax + coord_flip() + scale_x_reverse()) /(annotation_plot.horizontal + theme(legend.position="none")) + plot_layout(heights = c(0.7, 0.3))
3005
3006 fig5.c.combined2 <- (plot_spacer() + cowplot::ggdraw(top.legend) +
3007 plot_layout(widths = c(0.15, 0.85)) )/ fig5.c.combined1 + plot_layout(heights = c(0.1, 1))
3008 fig5.c.final <- fig5.c.combined2 + plot_layout(guides = "collect") & theme(legend.title = element_text(size = 8, face = "bold"),
3009 legend.text = element_text(size = 6),
3010 legend.key.size = unit(0.3, "cm"))
3011
3012 girafe(ggobj = fig5.c.final, width=12, height=7, options = list(opts_hover(css = "stroke: blue; stroke-width:3px;")))