This script plots Figure 1d: Host-microbiome associations of aspirin.

Scatterplot shows confidently deconfounded associations of aspirin usage in serum metabolome, host phenotype and microbiome features (Cliff’s delta), plotted against the effect size for the comparison between patients and healthy controls separated by clinical group. Thus, quadrants of the graph correspond to whether treatment signature direction is concordant with disease signature (indicating a potential severity marker) or discordant (indicating a possible direct effect of the medication). A subset of features is highlighted for interpretation and plotted as a bar graph.

Required file in the input_data folder:

Supplementary Table 6: Features of microbiome, host and metabolome impacted by different drug groups and drug compounds. Results of drug group (or drug compound according to the ATC classification) assessment for its impact on host and microbiome features for each patient group. Compound comparison with Maier et al., Nature 2018, tab shows microbiome features negatively impacted by the drug treatment (for the ATC-level compounds) in at least one patient group, and bacterial species whose growth was inhibited by the same drug in the in vitro experiment.

Load necessary libraries.

library (ggplot2)
library (cowplot)
library (ggrepel)
library (readxl)

Load input data from Supplementary Table 6.

fileFolder = 'C:/Users/mazimmer/Documents/Projects/metaCardis/SupplementaryTables/SupllTablesUpdate20201102/Submission/'
data <- read_excel(paste0(fileFolder,
                             "Supplementary_Table_6_2019-09-13434.xlsx"), 
                      sheet = "Drug group effect")
#Select only aspirine associations from the full table. 
data <- data[data$Effector=='Aspirine intake',]

Select a subset of data for plotting.

data$Number <- rownames (data)

samedata <- subset (data, Congruence == "Same")
oppositedata <- subset (data, Congruence == "Opposite")

# do not plot unannotated metabolome features
subdata <- subset (data, Congruence == "Opposite"  & ! (`Feature space` %in% c ("Serum, unannotated", "Urine, unannotated")))

# filter features by short name (to remove e.g. features measured by multiple methods)
subdata$Shortened.Name <- sub("\\(.*", "", subdata$`Feature display name`)
subdata <- subdata [! duplicated (subdata$Shortened.Name, subdata$Effector),]

Assemble graph for plotting.

p <- ggplot (data, aes (x = `Effect size`, y = `Comparison D`)) + 
  geom_vline (xintercept = 0) + 
  geom_hline (yintercept = 0) + 
  ylab ("Disease effect size (Cliff's delta)") + 
  xlab ("Drug effect size (Cliff's delta)") + 
  geom_point (aes (color = as.factor (sign (data$`Effect size`))), size = 5) + 
  theme (legend.position = "none", 
         axis.text = element_text (size = 20), 
         axis.title = element_text (size = 20)) + 
  scale_color_manual (values = c ("#d73027", "#4575b4")) + 
  geom_text_repel (data = subdata, aes (label = Number), color = "black", 
                   size = 3, min.segment.length = 0, force = 3, box.padding = 2,
                   point.padding = 1, segment.alpha = 0.5,
                   max.overlaps = 100) + 
  scale_x_continuous (breaks = c (-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0)) + 
  scale_y_continuous (breaks = c (-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0)) + 
  coord_cartesian (xlim = c (-1.0, 1.0), ylim = c (-1.0, 1.0)) + 
  ggtitle('Aspirin')

Select features with associations opposite from disease for plotting.

tsdata <- subset (subdata, Congruence == "Opposite")
tsdata$Number = factor (tsdata$Number, levels = unique (tsdata$Number [order (tsdata$`Effect size`, tsdata$Number)]), ordered = TRUE)
tsdataPos <- subset (tsdata, `Effect size` > 0)
tsdataNeg <- subset (tsdata, `Effect size` < 0)

Assemble the bar plot.

pDrug <- ggplot (tsdata, aes (x = `Effect size`, y = Number)) + 
  geom_segment (aes (xend = `Effect size`, x = 0, y = Number, yend = Number, 
                     color = as.factor (sign (tsdata$`Effect size`))), size = 4) +
  scale_x_continuous (breaks = c (-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0)) +
  coord_cartesian (xlim = c (-1.0, 1.0)) +
  theme (axis.title.x = element_blank (), axis.text.x = element_blank (),
         axis.ticks.x = element_blank (), axis.title.y = element_blank (),
         axis.text.y = element_blank (), axis.ticks.y = element_blank (),
         legend.position = "none") +
  geom_text (data = tsdataPos, aes (label = paste0 (tsdataPos$Number, ": ",
                                                    tsdataPos$`Feature display name`),
                                    y = Number, x = -0.01), hjust = 1, angle = 0, size = 3) +
  geom_text (data = tsdataNeg, aes (label = paste0 (tsdataNeg$Number, ": ",
                                                    tsdataNeg$`Feature display name`),
                                    y = Number, x = 0.01), hjust = 0, angle = 0, size = 3) +
  scale_color_manual (values = c ("#d73027", "#4575b4")) 

Combine two plots in a grid.

pGrid <- plot_grid (ncol = 1, p, pDrug, align = "v", axis = "bt", rel_heights = c (1, 2))

Display plot.

print (pGrid)

Uncomment to print plot to file.

#pdf(file = "fig1d_aspirine_associations.pdf", width=6, height=18)
#print (pGrid)
#dev.off ()
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] readxl_1.3.1  ggrepel_0.9.1 cowplot_1.1.1 ggplot2_3.3.3
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        highr_0.9         cellranger_1.1.0  bslib_0.2.4      
##  [5] compiler_3.6.1    pillar_1.6.0      jquerylib_0.1.4   tools_3.6.1      
##  [9] digest_0.6.27     jsonlite_1.7.2    evaluate_0.14     lifecycle_1.0.0  
## [13] tibble_3.1.0      gtable_0.3.0      pkgconfig_2.0.3   rlang_0.4.10     
## [17] DBI_1.1.1         yaml_2.2.0        xfun_0.22         withr_2.4.2      
## [21] stringr_1.4.0     dplyr_1.0.5       knitr_1.33        generics_0.1.0   
## [25] sass_0.3.1        vctrs_0.3.6       grid_3.6.1        tidyselect_1.1.0 
## [29] glue_1.4.2        R6_2.5.0          fansi_0.4.2       rmarkdown_2.7    
## [33] farver_2.1.0      purrr_0.3.4       magrittr_2.0.1    scales_1.1.1     
## [37] ellipsis_0.3.1    htmltools_0.5.1.1 assertthat_0.2.1  colorspace_2.0-0 
## [41] utf8_1.2.1        stringi_1.5.3     munsell_0.5.0     crayon_1.4.1