joelewis101/blantyreESBL
Data and Code Repository For a Study Describing Antimicrobial Resistance in Blantyre, Malawi

Package index
Search the joelewis101/blantyreESBL package
Vignettes
Functions
46
Source code
18
Man pages
46
Home
/
GitHub
/
joelewis101/blantyreESBL
/

In joelewis101/blantyreESBL: Data and Code Repository For a Study Describing Antimicrobial Resistance in Blantyre, Malawi 

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

Introduction

This document contains reproducing analysis code which generates the tables and figures for the manuscript:


Genomic analysis of extended-spectrum beta-lactamase (ESBL) producing Escherichia coli colonising adults in Blantyre, Malawi reveals previously undescribed diversity


Joseph M Lewis^1,2,3,4^, , Madalitso Mphasa^1^, Rachel Banda^1^, Matthew Beale^4^, Jane Mallewa^5^, Catherine Anscombe^1,2^, Allan Zuza^1^, Adam P Roberts^2^, Eva Heinz*^2^, Nicholas Thomson*^4,6^, Nicholas A Feasey*^1,2^

  1. Malawi Liverpool Wellcome Research Programme, Blantyre, Malawi
  2. Liverpool School of Tropical Medicine, Liverpool, UK
  3. University of Liverpool, Liverpool, UK
  4. Wellcome Sanger Institute, Hinxton, UK
  5. Kamuzu University of Health Sciences, Malawi

  6. London School of Tropical Medicine and Hygiene, London, UK

  7. = contributed equally

library(tidytree)
library(ggtree)
library(ape)
library(phytools)
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(forcats)
library(stringr)
library(blantyreESBL)
library(lubridate)
library(ggnewscale)
library(kableExtra)
library(patchwork)
library(here)
library(scales)
library(ggstance)
library(blantyreESBL)
library(countrycode)
library(viridis)

write_figs <- FALSE

if (write_figs) {
  if (!dir.exists(here("figures"))) {dir.create(here("figures"))}
  if (!dir.exists(here("tables"))) {dir.create(here("tables"))}

    if (!dir.exists(here("figures/ecoli-genomics"))) {
    dir.create(here("figures/ecoli-genomics"))
  }

  if (!dir.exists(here("tables/ecoli-genomics"))) {
    dir.create(here("tables/ecoli-genomics"))
    }

}

Assembly stats and accession numbers

btESBL_sequence_sample_metadata %>%
  mutate(date_of_collection = data_date) %>% 
  select(accession,
         lane,
         supplier_name,
         pid,
         date_of_collection,
         N50) %>%
  kbl(caption =
        "Accession numbers and assembly statistics for included samples"
      ) %>%
  kable_classic(full_width = FALSE)

if (write_figs) {
  write_csv(
    btESBL_sequence_sample_metadata %>%
      mutate(date_of_collection = data_date) %>%
      filter(species == "E. coli") %>%
      select(
        accession,
        lane,
        supplier_name,
        pid,
        date_of_collection,
        N50,
        Cluster,
        ST,
        ecoli_phylogroup,
        ecoli_pathotype
      ) %>%
      left_join(
        btESBL_qrdr_mutations %>%
          filter(genus == "E. coli") %>%
          select(-genus) %>%
          semi_join(
            btESBL_CARD_qrdr_mutations,
            by = c(
              "variant" = "variant",
              "gene" = "gene"
            )
          ) %>%
          pivot_wider(,
            id_cols = lane,
            names_from = gene,
            values_from = variant,
            values_fn = function(x) paste0(x, collapse = ",")
          )
      ) %>%
      mutate(ecoli_pathotype = if_else(is.na(ecoli_pathotype),
        "none_identfied",
        ecoli_pathotype
      )) %>%
      mutate(across(
        matches("gyr|par"),
        ~ if_else(is.na(.x), "no", .x)
      )) %>%
      left_join(
        btESBL_amrgenes %>%
          select(-genus) %>%
          arrange(ref_seq) %>%
          pivot_wider(
            id_cols = lane, names_from = ref_seq,
            values_from = ref_seq,
            values_fn = length,
            values_fill = 0
          ) %>%
          mutate(across(where(is.numeric), ~ if_else(.x == 1, "yes", "no")))
      ),
    here("tables/ecoli-genomics/sample_metadata.csv")
  )
}

Phylogroup, MLST plasmid replicons

btESBL_sequence_sample_metadata %>% 
  filter(species == "E. coli") %>% 
  mutate(Phylogroup = ecoli_phylogroup,
         Pathotype = ecoli_pathotype) ->
  dassimEcoli_BTEcoli.accession


pgroup_cols <- c(hue_pal()(9)[c(1:5,7:8)], "gray")

names(pgroup_cols) <-
  sort(unique(dassimEcoli_BTEcoli.accession$Phylogroup))

dassimEcoli_BTEcoli.accession %>% 
  group_by(ST) %>% 
  mutate(n = n()) %>% 
  filter(n > 2) %>% 
  ggplot(aes(fct_rev(fct_infreq(ST)), fill = Phylogroup)) +
  geom_bar() +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "ST", y = "Number of Samples") +
  scale_fill_manual(values = pgroup_cols) -> 
  p1 
#  scale_fill_viridis_d(option = "cividis") -> p1

dassimEcoli_BTEcoli.accession %>% 
  ggplot(aes(fct_rev(fct_infreq(Phylogroup)), fill = Phylogroup)) +
  geom_bar() +
  theme_bw() +
  theme(legend.position = "none", 
        axis.text.x = element_text(angle = 45,
                                   hjust = 1)) +
#  scale_fill_viridis_d(option = "cividis") +
  scale_fill_manual(values = pgroup_cols) +
  coord_flip() +
  labs(x = "Phylogroup", y = "Number of samples") -> p2

btESBL_plasmidreplicons %>% 
  filter(species == "E. coli") %>% 
   mutate(ref_seq = gsub("_.+$","", ref_seq),
         ref_seq = gsub("\\.[0-9]","", ref_seq),
         ref_seq = case_when(
           grepl("Col", ref_seq) ~ "Col",
           grepl("rep", ref_seq) ~ "rep",
           grepl("^FIA", ref_seq) ~ "IncFIA",
           !grepl("Inc", ref_seq) ~ "Other",
           TRUE ~ ref_seq
          )) %>% 
  filter(grepl("Inc", ref_seq)) %>% 
  ggplot(aes(fct_rev(fct_infreq(ref_seq)))) +
  geom_bar() +
  coord_flip() +
  theme_bw() +
  labs(y = "Number of samples",
       x = "Plasmid replicon") -> plasm_rep_prev




p2 + p1 + plasm_rep_prev + plot_spacer() + 
  plot_annotation(tag_levels = "A") +
  plot_layout(ncol = 2) -> mlst_plot

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/pgroup_mlst_plasm_plot.pdf"),
    mlst_plot,
    width = 7,
    height = 8)
  ggsave(
    here("figures/ecoli-genomics/pgroup_mlst_plasm_plot.svg"),
    mlst_plot,
    width = 7,
    height = 8)
}
mlst_plot

Compare within- to between- patient diversity

btESBL_sequence_sample_metadata %>%
  filter(species == "E. coli") %>%
  group_by(pid) %>%
  summarise(
    n_sts_per_participant = length(unique(ST)),
    n_samples_per_participant = n()
  ) %>%
  mutate(
    n_samples_per_participant = factor(n_samples_per_participant,
      levels = 1:5
    ),
    n_sts_per_participant = factor(n_sts_per_participant, levels = 1:5)
  ) %>%
  count(n_sts_per_participant, n_samples_per_participant, .drop = FALSE) %>%
  filter(as.numeric(as.character(n_samples_per_participant)) >=
    as.numeric(as.character(n_sts_per_participant))) %>%
  ggplot(aes(n_samples_per_participant, n_sts_per_participant,
    label = n, fill = n, color = -n
  )) +
  geom_tile(color = NA) +
  geom_text() +
  scale_fill_viridis_c() +
  scale_color_gradient(low = "black", high = "white") +
  guides(color = "none") +
  labs(
    x = "No. of samples per participant",
    y = "No. of STs per participant",
    fill = "No. of\nParticipants"
  ) +
  theme_bw() -> samples_vs_st_plot


btESBL_snpdists_esco %>%
  rename(sample1 = sample) %>%
  pivot_longer(-sample1,
    names_to = "sample2",
    values_to = "snpdist"
  ) %>%
  left_join(
    btESBL_sequence_sample_metadata %>%
      filter(species == "E. coli") %>%
      transmute(
        pid1 = pid,
        ST1 = ST,
        sample1 = lane
      )
  ) %>%
  left_join(
    btESBL_sequence_sample_metadata %>%
      filter(species == "E. coli") %>%
      transmute(
        pid2 = pid,
        ST2 = ST,
        sample2 = lane
      )
  ) %>%
  rowwise() %>%
  mutate(filter_var = paste(sort(c(sample1, sample2)), collapse = ",")) -> df

df[!duplicated(df$filter_var), ] -> df

df %>%
  filter(sample1 != sample2) %>%
  filter(ST1 == ST2) %>%
  filter(!is.na(pid1) & !is.na(pid2)) %>%
  mutate(btwn_var = if_else(pid1 == pid2, "Within-participant, same ST",
    "Between-participant, same ST"
  )) %>%
  ggplot(aes(snpdist)) +
  # geom_density() +
  geom_histogram(bins = 50, color = "black", size = 0.4) +
  xlim(0, 200) +
  facet_wrap(~btwn_var, scales = "free") +
  theme_bw() +
  guides(fill = "none") +
  labs(
    x = "Core genome SNP distance",
    y = "Number of sample-pairs"
  ) +
  scale_fill_viridis_d(option = "cividis") -> snpdist_plot


((samples_vs_st_plot + plot_spacer() + plot_layout(widths = c(1.5,1))) /
  snpdist_plot) +
  plot_annotation(tag_levels = "A") -> within_participant_diversity_plot

within_participant_diversity_plot


if (write_figs) {

  ggsave(
    here("figures/ecoli-genomics/within_participant_diversity_plot.pdf"),
  within_participant_diversity_plot,
    width = 6,
    height = 6)

  ggsave(
    here("figures/ecoli-genomics/within_participant_diversity_plot.svg"),
    within_participant_diversity_plot,
    width = 6,
    height = 6)

}

Compare phylogroup to horesh collection

As per reviewer's comments - stratify Horesh collection by invasive vs carriage and ESBL vs not ESBL for comparison

# add esbl status to btESBL_ecoli_horesh_metadata

btESBL_ecoli_horesh_metadata %>%
  left_join(
    btESBL_ecoli_horesh_amr %>%
      left_join(
        select(
          btESBL_NCBI_phenotypic_bl,
          allele_name,
          class
        ),
        by = c("ref_seq" = "allele_name")
      ) %>%
      group_by(ID) %>%
      summarise(res = if_else(
        any(grepl("ESBL|AmpC|Carbapenemase", class)),
        "ESBL/AmpC/Carbapenemase",
        "No ESBL/AmpC/Carbapenemase"
      )),
    by = c("ID" = "ID")
  ) -> btESBL_ecoli_horesh_metadata

bind_rows(
btESBL_ecoli_horesh_metadata %>% 
  transmute(Study = "Global",
            Phylogroup = Phylogroup),
btESBL_sequence_sample_metadata %>% 
  filter(species == "E. coli") %>% 
  transmute(Study = "This Study",
            Phylogroup = ecoli_phylogroup)
) %>% 
  mutate(Phylogroup =
           if_else(Phylogroup == "Not Determined",
                   "Unknown",
                   Phylogroup),
         Phylogroup = fct_infreq(Phylogroup)
  ) %>% 
  count(Study, Phylogroup, .drop = FALSE) %>% 
  group_by(Study) %>% 
  mutate(prop = n/sum(n)) %>% 
  ggplot(aes(fct_reorder(Phylogroup, prop),
             prop,
             group = Study,
             fill = Study)) +
  geom_col(position = "dodge") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Phylogroup",
       y = "Proportion") ->
  malawi_vs_global_pgroup

malawi_vs_global_pgroup

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/global_vs_malawi_pgroup.svg"),
    malawi_vs_global_pgroup,
    width = 5,
    height = 4
  )
  ggsave(
    here("figures/ecoli-genomics/global_vs_malawi_pgroup.pdf"),
    malawi_vs_global_pgroup,
    width = 5,
    height = 4
  )
}

bind_rows(
  btESBL_ecoli_horesh_metadata %>%
    transmute(
      Study = "Global",
      Phylogroup = Phylogroup,
      res = res,
      source = case_when(
        Isolation == "Unknown" ~ "Unknown",
        Isolation %in% c("Feces", "Rectal", "Rectal
                                                Swab") ~ "Stool",
        TRUE ~ "Invasive"
      )
    ),
  btESBL_sequence_sample_metadata %>%
    filter(species == "E. coli") %>%
    transmute(
      Study = "This Study",
      Phylogroup = ecoli_phylogroup,
      res = "ESBL/AmpC/Carbapenemase",
      source = "Stool"
    )
) %>%
  mutate(
    Phylogroup =
      case_when(
        Phylogroup == "Not Determined" ~ "Unknown",
        grepl("clade|Shigella", Phylogroup) ~ "Other",
        TRUE ~ Phylogroup
      )
  ) %>%
  filter(!is.na(res)) %>%
  #  mutate(across(everything(), as.factor)) %>% */
  count(Study, Phylogroup, res, source, .drop = FALSE) %>%
  group_by(Study, res, source) %>%
  mutate(prop = n / sum(n),
         source = factor(source,
                      levels = c("Stool","Invasive","Unknown")
         )) %>%
  ggplot(aes(Phylogroup,
    prop,
    group = Study,
    fill = Study
  )) +
  geom_col(position = "dodge") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "Phylogroup",
    y = "Proportion"
  ) +
  facet_grid(source ~ res) -> malawi_vs_global_pgroup_stratified

malawi_vs_global_pgroup_stratified

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/global_vs_malawi_pgroup_stratified.svg"),
    malawi_vs_global_pgroup_stratified,
    width = 6,
    height = 6
  )

  ggsave(
    here("figures/ecoli-genomics/global_vs_malawi_pgroup_stratified.pdf"),
    malawi_vs_global_pgroup_stratified,
    width = 6,
    height = 6
  )

}

Virulence

Only 2 aEPEC/EPEC, 10 EAEC.

Malawi phylogeny

dassimEcoli_BTEcoli.accession %>% 
  as.data.frame() -> dassimEcoli_BTEcoli.accession 

rownames(dassimEcoli_BTEcoli.accession ) <-
  dassimEcoli_BTEcoli.accession$lane

 dassimEcoli_BTEcoli.accession %>% 
  select(lane, ST) %>% 
  group_by(ST) %>% 
  mutate(n = n(),
         ST = if_else(n == 1, "Other", ST)) %>% 
  mutate(ST = if_else(!ST %in% c("Novel", "Other"), 
                      paste0("ST", ST), ST)) %>% 
  arrange(fct_infreq(ST)) %>% 
  pivot_wider(id_cols = lane, 
              names_from = ST,
              values_from = ST,
              values_fn = length,
              values_fill = 0) %>% 
  mutate(across(everything(), as.character)
         ) %>% 
  relocate(Other, .after = everything()) %>% 
  as.data.frame() ->
  mlst_onehot

rownames(mlst_onehot) <- mlst_onehot$lane

(
(
  ggtree(btESBL_coregene_tree_esco) %>%
    gheatmap(
      select(dassimEcoli_BTEcoli.accession, Phylogroup),
      width = 0.05,
      color = NA,
      font.size = 3,
      colnames_angle = 90,
      colnames_position = "top",
      colnames_offset_y = 3,
      hjust = 0
    ) +
    scale_fill_manual(
      values = hue_pal()(8),
      name = "Phylogroup") + 
   new_scale_fill()
) %>%
  gheatmap(select(mlst_onehot, -lane),
           font.size = 2.5,
           color = NA,
           colnames_angle = 90,
           offset = .0022,
           colnames_offset_y = 3,
           colnames_position = "top",
           hjust = 0)  +
    scale_fill_manual(values = c("lightgrey", "black"), guide = "none") +
    new_scale_fill()
  ) %>% 
  gheatmap(select(dassimEcoli_BTEcoli.accession, Pathotype),
           width = 0.05,
           font.size = 3,
           color = NA,
           colnames_angle = 90,
           offset = .001,
           colnames_offset_y = 3,
           colnames_position = "top",
           hjust = 0)  +
  scale_fill_manual(values = viridis_pal(option = "plasma")(6)[c(3,5)],
                    name = "Pathotype",
                    na.translate = FALSE) +
  ylim(NA, 540) + 
  annotate("text", x = 0.03, y = 520, label = "Sequence Type",
           size = 3) -> malawi_treeplot


malawi_treeplot

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/malawi_treeplot.svg"),
    malawi_treeplot,
    width = 9,
    height = 9
  )
  ggsave(
    here("figures/ecoli-genomics/malawi_treeplot.pdf"),
    malawi_treeplot,
    width = 9,
    height = 9
  )
}

AMR determinants

# add class ro which resistnce is conferred --------------------------------

quinolone <- "Par|Gyr|Par|Qnr|Qep|Nor|GyrA|GyrB|ParC|ParE"
tetracycline <- "Tet"
sulphonamide <- "Sul"
aminoglycoside <- "Str|Aad|Aac|Aph|Rmt|APH"
streptothricin <- "Sat"
macrolide <- "Mph|Mdf|Erm|Ere"
fosfomycin <- "Fos"
chloramphenicol <- "Cat|FloR|Cml"
trimethoprim <- "Dfr"
rifampicin <- "Arr"
ESBL <- "SHV_12"
penicillinase <- "OKP|SCO|LEN|LAP"
ampc <- "CMY"

btESBL_amrgenes %>% 
  semi_join(dassimEcoli_BTEcoli.accession, 
            by = c( "lane")) %>% 
  select(-genus) %>% 
  rename(gene = ref_seq) %>% 
  # add in QRDR mutations
  bind_rows(
    btESBL_qrdr_mutations %>% 
      filter(genus == "E. coli") %>% 
      select(-genus) %>% 
      semi_join(
        btESBL_CARD_qrdr_mutations, 
        by = c("variant" = "variant",
              "gene" =  "gene")
      ) %>% 
      select(gene, lane) %>% 
      unique() %>% 
      mutate(gene =
               gsub("(^.{1})", '\\U\\1',
                    gene,
                    perl = TRUE))
      ) %>% 
  filter(!grepl("AmpH|AMPH|AmpC|MrdA|MefB", gene)) %>% 
  # add in beta-lactamases
  left_join(
    select(btESBL_NCBI_phenotypic_bl, allele_name, class),
    by = c("gene" = "allele_name")) %>% 
  mutate(class = case_when(
    str_detect(gene, quinolone) ~ "Quinolone",
    str_detect(gene, tetracycline) ~ "Tetracycline",
    str_detect(gene, sulphonamide) ~ "Sulphonamide",
    str_detect(gene, aminoglycoside) ~ "Aminoglycoside",
    str_detect(gene, streptothricin) ~ "Streptothricin",
    str_detect(gene, macrolide) ~ "Macrolide",
    str_detect(gene, fosfomycin) ~ "Fosfomycin",
    str_detect(gene, chloramphenicol) ~ "Chloramphenicol",
    str_detect(gene, rifampicin) ~ "Rifampicin",
    str_detect(gene,trimethoprim) ~ "Trimethoprim",
    str_detect(gene,ESBL) ~ "ESBL",
    str_detect(gene,penicillinase) ~ "Penicillinase",
    str_detect(gene,ampc) ~ "AmpC",
    TRUE ~ class
  )) %>% 
  mutate(gene = if_else(gene == "TEM_95",
                        "TEM_1",
                        gene)) ->
  amr

### plot prevalence ----------------------------------------


amr %>% 
  group_by(class) %>% 
  mutate(n_class = length(class),
         n_genes_in_class = n_distinct(gene)) %>% 
  select(class, n_class,n_genes_in_class) %>% 
  unique() %>% 
  arrange(n_class) %>% 
  ungroup() %>% 
  mutate(end = cumsum(n_genes_in_class),
         start = lag(end, default = 0),
         textpos = start + 0.6 + 0.5 * (end - start)) -> annotate.df

amr %>% 
  group_by(class) %>% 
  mutate(n_class = length(class),
         n_genes_in_class = n_distinct(gene),
         gene = gsub("_", "-", gene)) %>% 
ggplot( aes(fct_reorder(fct_rev(fct_infreq(gene)), n_class), 
            fill = class)) + 
  geom_bar() +
  theme_bw() +
  coord_flip(ylim = c(-5,450),clip = "off", expand = FALSE, xlim = c(0, 76)) +
  annotate(geom = "segment",
             x = annotate.df$start + 0.5 + 0.2, 
             xend = annotate.df$end + 0.5 - 0.2, 
             y = -115, yend = -115 ) +
  annotate(geom = "text", y = -125,
           x = annotate.df$textpos,
           label = annotate.df$class, 
           size = 3, 
           hjust = 1) +
  labs(y = "Number of samples") + 
  theme(plot.margin = unit(c(0.2,0.5,0.2,4), "cm"),
        axis.title.y  = element_blank(),
        legend.position = "none") -> amrplot

amrplot


if (write_figs) {
  ggsave(
    filename = here("figures/ecoli-genomics/amr_plot.pdf"),
    plot = amrplot,
    width = 6,
    height = 10
  )
  ggsave(
    filename = here("figures/ecoli-genomics/amr_plot.svg"),
    plot = amrplot,
    width = 6,
    height = 10
  )
}

AMR mapped to tree

amr %>% 
  # filter(!gene %in% c("OqxA", "OqxB", "FosA")) %>% 
  # remove FosA as not in KpI which we are plottingh
  mutate(
         gene = gsub("_","-", gene)) %>% 
  group_by(class) %>% 
  mutate(n_class = length(class)) %>% 
  ungroup() %>% 
  group_by(gene) %>% 
  mutate(n_gene = length(gene)) %>% 
  arrange(desc(n_class), desc(n_gene)) %>% 
  select(-n_class, -n_gene) %>% 
  pivot_wider(id_cols = lane,
              values_from = class,
              names_from = gene) %>% 

  as.data.frame() ->
  amr.ariba.maptotree


amr %>% 
    # filter(!gene %in% c("OqxA", "OqxB", "FosA")) %>% 
    mutate(
           gene = gsub("_","-", gene)) %>% 
    group_by(class) %>% 
    mutate(n_class = length(class)) %>% 
    ungroup() %>% 
    group_by(gene) %>% 
    mutate(n_gene = length(gene)) %>% 
    arrange(desc(n_class), desc(n_gene)) %>% pull(class) %>% 
  unique() -> class_order

rownames(amr.ariba.maptotree) <- amr.ariba.maptotree$lane

colz = hue_pal()(13)
names(colz) <- sort(unique(amr %>% 
                             filter(class != "Fosfomycin") %>% 
                             pull(class)))

ggtree(btESBL_coregene_tree_esco ) %>% 
#ggtree(tree) %>% 
    gheatmap(
        select(amr.ariba.maptotree,-lane),
        width = 5, 
        color = NA, 
        font.size = 3, 
        colnames_angle = 90, 
        colnames_position = "top", 
        colnames_offset_y = 0,
        hjust = 0,
        offset = 0
    ) +
  ylim(NA, 500) +
  scale_fill_manual(name = "Class", values = colz, 
                    breaks = class_order,
                    na.translate = FALSE) ->
  malawi_tree_with_amr

malawi_tree_with_amr


if (write_figs) {

  ggsave(
    here("figures/ecoli-genomics/malawi_treeplot_with_amr.pdf"),
    malawi_tree_with_amr,
    width = 12,
    height = 7)
  ggsave(
    here("figures/ecoli-genomics/malawi_treeplot_with_amr.svg"),
    malawi_tree_with_amr,
    width = 12,
    height = 7)

}

Malawi isolates in a global context

Phylogeny

 #prepare metadata ----------
# lookup list of horesh merged clusters 



purrr::map_df(
  unique(grep("_", btESBL_ecoli_global_popPUNK_clusters$Cluster, value = TRUE)),
              function(x)
                data.frame(mergeclust = x,
                           splts = strsplit(x, "_")[[1]])) %>%
  mutate(splts = as.numeric(splts)) -> mergclust_lookup

bind_rows(
  dassimEcoli_BTEcoli.accession %>%
    transmute(
      lane = lane,
      ST = ST,
      Phylogroup = Phylogroup,
      study = "dassim",
      Country = "Malawi",
      Continent = "Africa",
      Year = year(data_date),
      res = "ESBL/AmpC/Carbapenemase",
      source = "Colonising"
    ) %>%
    left_join(
      btESBL_ecoli_global_popPUNK_clusters,
      by = "lane"
    ),
  btESBL_ecoli_musicha_metadata %>%
    transmute(
      lane = lane,
      Cluster = Cluster,
      ST = if_else(is.na(ST),
        "Novel",
        as.character(ST)
      ),
      Phylogroup = phylogroup,
      study = "musicha",
      Country = "Malawi",
      Continent = "Africa",
      Year = Year
    ),
  btESBL_ecoli_horesh_metadata %>%
    mutate(
      assembly_name_recode =
        gsub("\\..*$", "", Assembly_name)
    ) %>%
    left_join(mergclust_lookup, by = c("PopPUNK" = "splts")) %>%
    transmute(
      lane = gsub(
        "#",
        "_",
        if_else(
          !is.na(name_in_presence_absence),
          name_in_presence_absence,
          assembly_name_recode
        )
      ),
      Cluster = if_else(!is.na(mergeclust),
        mergeclust,
        as.character(PopPUNK)
      ),
      ST = gsub("~", "", ST),
      Phylogroup = Phylogroup,
      study = "horesh",
      Country = Country,
      Continent = Continent,
      Year = case_when(
        Year == "Unknown" ~ NA_real_,
        grepl("s", Year) ~ as.numeric(gsub("s", "", Year)),
        grepl("-", Year) ~ as.numeric(gsub("-.*$", "", Year)),
        TRUE ~ as.numeric(str_trim(Year))
      ),
      Pathotype = Pathotype,
      res = res,
      source = case_when(
        Isolation == "Unknown" ~ "Unknown",
        Isolation %in% c("Feces", "Rectal", "Rectal
                                                Swab") ~ "Colonising",
        TRUE ~ "Infecting"
      )
    )
) %>%
  mutate(
    Continent =
      case_when(
        grepl("America", Continent) ~ "Americas",
        is.na(Continent) ~ "Unknown",
        TRUE ~ Continent
      )
  ) %>%
  mutate(
    cluster_recode =
      if_else(
        !grepl(
          paste(c(
            "\\<1", 2:1153, "1184\\>"
          ), collapse = "\\>|\\<"),
          Cluster
        ) &
          !grepl("_", Cluster),
        NA_character_,
        Cluster
      ),
    cluster_recode = factor(
      cluster_recode,
      levels =
        unique(cluster_recode)[
          order(
            as.numeric(
              gsub("_.*$", "", unique(cluster_recode))
            )
          )
        ]
    ),
    Phylogroup =
      case_when(
        Phylogroup %in% c(
          "Unknown", "Not Determined",
          "E or cladeI"
        ) ~
        "Unknown",
        is.na(Phylogroup) ~ "Unknown",
        TRUE ~ Phylogroup
      )
  ) -> df_all
df_all <- as.data.frame(df_all)

write_csv(df_all, here("tables/ecoli-genomics/horesh_globaltree_metadata.csv"))

rownames(df_all) <- df_all$lane

#define color palettes

pgroup_cols <- c(brewer_pal(type = "qual", palette = 8)(10),
                            "white")
  names(pgroup_cols) <- sort(unique(df_all$Phylogroup))

continent_cols <- c(viridis_pal(option = "plasma")(5),
                    "white")
names(continent_cols) <- sort(unique(df_all$Continent))

country_cols <-
  c("Malawi" =  "grey15",
    "Not Malawi" = "lightgrey")

btESBL_ecoli_globaltree -> horesh_context_tree

popPUNK_clusts_inplots <- 
  df_all %>% 
   filter(lane %in% 
           c(tree_subset(
             btESBL_ecoli_globaltree,
             3555, levels_back = 0)$tip.label,
             tree_subset(
               btESBL_ecoli_globaltree,
               5404, levels_back = 0)$tip.label,
             tree_subset(
               btESBL_ecoli_globaltree,
               5652, levels_back = 0)$tip.label)
  ) %>% 
  select(cluster_recode) %>% 
  unique() %>% 
  arrange(cluster_recode) %>% 
  pull(cluster_recode)


popPUNK_cols <- c(
  "dodgerblue2", "#E31A1C", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "black", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown"
)


#   c(
#   "palegreen2",
#   "#E31A1C",
#   "green4",
#   "#6A3D9A",
#   "#FF7F00",
#   "steelblue4",
#   "gold1",
#   "skyblue2",
#   "dodgerblue2",
#   "#FDBF6F",
#   "gray70",
#   "maroon",
#   "orchid1",
#   "darkturquoise",
#   "darkorange4",
#   "brown",
#   "white"
# )



names(popPUNK_cols) <- df_all %>% 
  filter(lane %in% 
           c(tree_subset(
             btESBL_ecoli_globaltree,
             3565, levels_back = 0)$tip.label,
             tree_subset(
               btESBL_ecoli_globaltree,
               5404, levels_back = 0)$tip.label,
             tree_subset(
               btESBL_ecoli_globaltree,
               5652, levels_back = 0)$tip.label)
  ) %>% 
  select(cluster_recode) %>% 
  unique() %>% 
  arrange(cluster_recode) %>% 
  pull(cluster_recode)

popPUNK_cols <- popPUNK_cols[!is.na(names(popPUNK_cols))]



offset_add <- 0

(
    (
      (
        (
          ggtree(btESBL_ecoli_globaltree, size = 0.2) %<+%
            (select(df_all, lane, ST) %>%
               mutate(
                 ST = if_else(ST == 167,
                              as.character("hilight"), NA_character_)
               ))  %>%
            gheatmap(
              select(df_all, `Year`),
              font.size = 4,
              width = 0.03,
              colnames_position = "top",
              color = NA,
              colnames_offset_y = 5,
              colnames_angle = 90,
              offset = 0.0006 + offset_add,
              hjust = 0
            )  +
            scale_fill_viridis_c(name = "Year", na.value = "white",
                                 guide = guide_colorbar(order = 2)) +
            new_scale_fill()
        ) %>%
          gheatmap(
            select(df_all, Continent),
            font.size = 4,
            width = 0.03,
            colnames_position = "top",
            color = NA,
            colnames_offset_y = 5,
            colnames_angle = 90,
            offset = 0.0012 + offset_add,
            hjust = 0
          ) +
          scale_fill_manual(values = 
                              continent_cols,
                            name = "Continent",
                            guide = guide_legend(order = 3)) +
          new_scale_fill()
      ) %>%
        gheatmap(
          select(df_all, Country) %>%
            mutate(Country = if_else(
              Country == "Malawi",
              "Malawi", "Not Malawi"
            )),
          font.size = 4,
          width = 0.03,
          colnames_position = "top",
          color = NA,
          colnames_offset_y = 5,
          colnames_angle = 90,
          offset = 0.0018 + offset_add,
          hjust = 0
        ) +
        scale_fill_manual(values = country_cols,
                          name = "Country",
                          guide = guide_legend(order = 4)) +
        new_scale_fill()
    ) %>%
      gheatmap(
        select(df_all, Phylogroup) %>% 
          mutate(Phylogroup = if_else(
            is.na(Phylogroup) ,"Unknown",
            Phylogroup)),
        width = 0.03,
        color = NA,
        font.size = 4,
        colnames_angle = 90,
        colnames_position = "top",
        colnames_offset_y = 3,
        hjust = 0,
        offset = 0 + offset_add
      ) +
      scale_fill_manual(values = pgroup_cols,
      name = "Phylogroup",
      na.translate = FALSE,
      guide = guide_legend(order = 1)) +
      new_scale_fill()
  ) +
    ylim(NA, 3550) +# + geom_tippoint(aes(color = ST)) + 
    geom_treescale(y = 2366, x = 0.004,offset = 30)  -> p

  p + geom_hilight( node = 3565, alpha = 0.3,
                    fill = "steelblue",
                    color = NA) +  # phylogroup A
    geom_cladelabel(node = 3565, label = "C", barsize = 0,offset = 0) +
    geom_highlight(node = 5404, alpha = 0.3, fill = "steelblue",
                   color = NA, extend = 0.0002) + # st 410 and context
    geom_cladelabel(node = 5404, label = "B", barsize = 0,offset = 0) +
    geom_highlight(node = 5652, alpha = 0.3, fill = "steelblue",
                   color = NA) + # ST131
    geom_cladelabel(node = 5652, label = "D", barsize = 0,offset = 0) -> p

  # phylo a subreee -----------


offset_add <- 0.0000 #0.0015
textsize <-  3
cladelab_offset <- 0.0001
tippoint_alpha = 0.7

(
  ggtree(
    tree_subset(
      btESBL_ecoli_globaltree, 3565, levels_back = 0)
    ) %<+% select(df_all, lane, cluster_recode) %>% 
    gheatmap(
    select(df_all, Country) %>%
      mutate(Country = if_else(
        Country == "Malawi",
        "Malawi", "Not Malawi")),
    font.size = 3,
    width = 0.04,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = offset_add,
    hjust = 0
    ) +
    scale_fill_manual(
      values = country_cols, 
      name = "Country",
      guide = "none") + 
      guides(fill = "none") + 
      new_scale_fill()
) + 
   ylim(NA, 100) +
   geom_tippoint(aes(color = cluster_recode), 
                na.rm = TRUE,
                size = 1,
                alpha = tippoint_alpha) +
  scale_color_manual(name = "PopPUNK\nCluster",
                       values = popPUNK_cols,
                       na.translate = FALSE) +
  geom_treescale(y = 66,
                 x = 0.00005,
                 offset  = 2,
                 width = 0.00001) -> p_phylo_a_subtree

p_phylo_a_subtree + 
  geom_cladelabel(node = 119, label = "ST167", 
                  fontsize = textsize,
                  offset = 0.00002,
                  offset.text = 0.000005) +
  geom_cladelabel(node = 99, label = "ST167", 
                  fontsize = textsize,
                  offset = 0.0001665,
                  offset.text = 0.000005) +
  theme(legend.position = "bottom")  +
  guides(
    color = guide_legend(
      nrow = 2,
      byrow = TRUE,
      override.aes = list(alpha = 1))
    ) ->
  plot_phyloa_subtree

# ST131 zoom

offset_add <- 0.0004
(
  ggtree(
    tree_subset(
      btESBL_ecoli_globaltree, 5652, levels_back = 0)
    ) %<+% select(df_all, lane, cluster_recode) %>%
    gheatmap(
      select(df_all, Country) %>%
      mutate(Country = if_else(
        Country == "Malawi",
        "Malawi", "Not Malawi")),
      font.size = 3,
      width = 0.07,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 0,
      colnames_angle = 90,
      offset =  offset_add,
      hjust = 0
    ) +
    scale_fill_manual(values = country_cols, 
                        name = "Country",
                        guide = "none") +
    new_scale_fill()
) +
  ylim(NA, 125) +
  geom_tippoint(aes(color = cluster_recode), 
                na.rm = TRUE, 
                size = 1,
                alpha = tippoint_alpha) +
  scale_color_manual(name = "PopPUNK\nCluster",
                     values = popPUNK_cols,
                     na.translate = FALSE) +
  geom_treescale(y = 83,
                 x = 0.0001,
                 offset  = 3.6,
                 width = 0.0001) -> p_st131_subtree

p_st131_subtree + 
  geom_cladelabel(node = 104, label = "ST131", 
                  fontsize = textsize,
                  offset = 0.0001,
                  offset.text = 0.00002) +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(
    nrow = 2,
    byrow = TRUE,
    override.aes = list(alpha = 1)) 
    ) -> plot_st131tree

# st410 subtree



offset_add <- 0.00007
# offset_multiple <- 0.00008

(
  ggtree(
    tree_subset(
      btESBL_ecoli_globaltree, 5404, levels_back = 0)
    ) %<+% select(df_all, lane, cluster_recode) %>%
    gheatmap(
    select(df_all, Country) %>%
      mutate(Country = if_else(
        Country == "Malawi",
        "Malawi", "Not Malawi")),
    font.size = 3,
    width = 0.055,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = offset_add,
    hjust = 0,
    ) +
    scale_fill_manual(values = country_cols, 
                      name = "Country",
                      guide = "none") +
    new_scale_fill()
) +
  ylim(NA, 60) +
  geom_tippoint(aes(color = cluster_recode),
                na.rm = TRUE,
                size = 1,
                alpha = tippoint_alpha) +
  scale_color_manual(name = "PopPUNK\nCluster",
                     values = popPUNK_cols,
                     na.translate = FALSE) +
  geom_treescale(y = 40,
                 x = offset_add/2,
                 offset  = 1,
                 width = 0.00001) -> p_st410_subtree

p_st410_subtree +
  geom_cladelabel(48, label = "ST410", 
                  fontsize = textsize,
                  offset = offset_add/5,
                  offset.text = 0.000002) +
  theme(legend.position = "bottom")  +
  guides(
    color = guide_legend(
      nrow = 2, byrow = TRUE,
         override.aes = list(alpha = 1)
      )
    ) -> plot_st410tree

(plot_st410tree + plot_phyloa_subtree + plot_st131tree) + 
  plot_annotation(tag_levels = "A") + plot_layout(guides = "collect") & 
  theme(legend.position = "bottom")   -> subtreesplot

p / (subtreesplot) + plot_annotation(tag_levels = "A") + 
  plot_layout(guides = "collect", heights = c(2.5,1)) -> global_tree_plot

# height 12, width 9 loooks gooooooood
global_tree_plot

if (write_figs) {

  ggsave(
    here("figures/ecoli-genomics/global_treeplot.pdf"),
    global_tree_plot,
    width = 9,
    height = 13)
  ggsave(
    here("figures/ecoli-genomics/global_treeplot.svg"),
    global_tree_plot,
    width = 9,
    height = 12)

}

PopPUNK clusters

df_all %>%
  filter(study %in% c("dassim", "horesh")) %>%
  mutate(res = if_else(grepl("No", res), "No ESBL", "ESBL")) %>%
  group_by(study) %>%
  mutate(n_tot_study = n()) %>%
  group_by(Cluster, study) %>%
  summarise(
    n_cluster = n(),
    prop_cluster = n_cluster / unique(n_tot_study)
  ) %>%
  group_by(Cluster) %>%
  mutate(n_studies_cluster_is_in = length(Cluster)) %>%
  filter(
    grepl(
      paste(c("\\<1", 2:49, "50\\>"), collapse = "\\>|\\<"),
      Cluster
    ) |
      n_studies_cluster_is_in > 1 |
      study == "dassim"
  ) %>%
  filter(!(study == "dassim" & n_cluster <= 3 & n_studies_cluster_is_in == 1)) %>%
  ungroup() %>%
  mutate(
    prop_cluster =
      if_else(
        study == "dassim",
        -1 * prop_cluster,
        prop_cluster
      ),
    Cluster = factor(
      Cluster,
      levels =
        unique(Cluster)[
          order(
            as.numeric(
              gsub("_.*$", "", unique(Cluster))
            )
          )
        ]
    )
  ) %>%
  mutate(
    study = case_when(
      study == "dassim" ~ "This Study",
      TRUE ~ "Horesh et al."
    ),
    study = factor(study, levels = c("This Study", "Horesh et al."))
  ) %>%
  ggplot(aes(x = prop_cluster, y = fct_rev(Cluster), fill = study)) +
  geom_col(color = "black", size = 0.2) +
  theme_bw() +
  labs(
    y = "PopPUNK cluster", fill = "",
    x = "Proportion of samples"
  ) +
  scale_x_continuous(labels = c(0.1, 0.0, 0.1, 0.2, 0.3)) +
  scale_fill_manual(values = c(
    viridis_pal(option = "cividis")(8)[c(2)],
    "grey60"
  )) +
  theme(
    axis.text.y = element_text(size = 5),
    axis.text.x = element_text(size = 7),
    legend.position = "bottom"
  ) -> p1

df_all %>%
  filter(study %in% c("dassim", "horesh")) %>%
  mutate(
    study = case_when(
      study == "dassim" ~ "This Study",
      TRUE ~ "Horesh Collection"
    ),
    study = factor(study, levels = c("This Study", "Horesh Collection"))
  ) %>%
  mutate(Cluster = factor(
    Cluster,
    levels =
      unique(Cluster)[
        order(
          as.numeric(
            gsub("_.*$", "", unique(Cluster))
          )
        )
      ]
  )) %>%
  ggplot(aes(Cluster, color = study, group = study)) +
  stat_ecdf(pad = TRUE, geom = "step") +
  scale_x_discrete(
    breaks = as.character(c(1, seq(from = 100, to = 1400, by = 100))),
    expand = c(0.01, 0.01)
  ) +
  theme_bw() +
  labs(
    y = "Cumulative proportion",
    x = "PopPUNK cluster", color = ""
  ) +
  geom_vline(
    xintercept = "1184",
    alpha = 0.8,
    linetype = "dotted"
  ) +
  scale_color_manual(
    values = c(
      viridis_pal(option = "cividis")(8)[c(2)],
      "grey60"
    ),
    guide = "none"
  ) +
  theme(
    axis.text.y = element_text(size = 7),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 5),
    legend.position = "none"
  ) -> p2

  (
    p1 +
      (p2 + plot_spacer() + plot_layout(ncol = 1))
  ) + plot_layout(ncol = 2, guides = "collect") +
    plot_annotation(tag_levels = "A") &
    theme(legend.position = "bottom") -> pp_cluster_plots

  pp_cluster_plots

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/global_poppunk_clusterplot.pdf"),
    pp_cluster_plots,
    width = 6,
    height = 6
  )
  ggsave(
    here("figures/ecoli-genomics/global_poppunk_clusterplot.svg"),
    pp_cluster_plots,
    width = 6,
    height = 6
  )
}

# get clusters to include


df_all %>%
  mutate(res = if_else(grepl("No", res), "No ESBL", "ESBL")) %>%
  filter(study %in% c("dassim", "horesh")) %>%
  group_by(study) %>%
  mutate(n_tot_study = n()) %>%
  group_by(Cluster, study) %>%
  summarise(
    n_cluster = n(),
    prop_cluster = n_cluster / unique(n_tot_study)
  ) %>%
  group_by(Cluster) %>%
  mutate(n_studies_cluster_is_in = length(Cluster)) %>%
  filter(
    grepl(
      paste(c("\\<1", 2:49, "50\\>"), collapse = "\\>|\\<"),
      Cluster
    ) |
      n_studies_cluster_is_in > 1 |
      study == "dassim"
  ) %>%
  filter(!(study == "dassim" & 
           n_cluster <= 3 & n_studies_cluster_is_in == 1)) %>%
  ungroup() %>%
  pull(Cluster) %>%
  unique() -> clusters_to_include

df_all %>%
  mutate(res = if_else(grepl("No", res), "No ESBL", "ESBL")) %>%
  filter(study %in% c("dassim", "horesh")) %>%
  group_by(study) %>%
  filter(source != "Unknown") %>%
  mutate(n_tot_study = n()) %>%
  group_by(Cluster, study, res, source) %>%
  summarise(
    n_cluster = n(),
    prop_cluster = n_cluster / unique(n_tot_study)
  ) %>%
  filter(Cluster %in% clusters_to_include) %>%
  ungroup() %>%
  mutate(
    prop_cluster =
      if_else(
        study == "dassim",
        -1 * prop_cluster,
        prop_cluster
      ),
    Cluster = factor(
      Cluster,
      levels =
        unique(Cluster)[
          order(
            as.numeric(
              gsub("_.*$", "", unique(Cluster))
            )
          )
        ]
    )
  ) %>%
  mutate(
    var = case_when(
      study == "dassim" ~ "This Study\nColonising",
      study == "horesh" & grepl("Colonising", source) ~ "Horesh et al.\nStool",
      TRUE ~ "Horesh et al.\nInvasive"
    ),
    var = factor(var, levels = sort(unique(var))[c(3, 2, 1)])
  ) %>%
  ggplot(aes(x = prop_cluster, y = fct_rev(Cluster), fill = var)) +
  geom_col(color = "black", size = 0.2) +
  theme_bw() +
  labs(
    y = "PopPUNK cluster", fill = "",
    x = "Proportion of samples"
  ) +
  # scale_x_continuous(labels = c(0.1,0.0, 0.1, 0.2, 0.3)) +
  scale_fill_manual(
    values = c(
      viridis_pal(option = "cividis")(8)[c(2, 8, 6)]
    ),
    breaks = c(
      "This Study\nColonising",
      "Horesh et al.\nInvasive",
      "Horesh et al.\nStool"
    )
  ) +
  scale_color_manual(
    values = c(
      viridis_pal(option = "cividis")(8)[c(2, 8, 6)]
    ),
    breaks = c(
      "This Study\nColonising",
      "Horesh et al.\nInvasive",
      "Horesh et al.\nStool"
    )
  ) +
  theme(
    axis.text.y = element_text(size = 5),
    axis.text.x = element_text(size = 7),
    legend.position = "bottom"
  ) +
  facet_wrap(. ~ res) +
  scale_x_continuous(
    breaks = c(-0.1, 0, 0.1),
    labels = c(0.1, 0.0, 0.1)
  ) -> a

df_all %>%
  mutate(res = as.factor(if_else(grepl("No", res), "No ESBL", "ESBL"))) %>%
  filter(study %in% c("dassim", "horesh")) %>%
  #   filter(source != "Unknown")  %>%
  mutate(Cluster = factor(
    Cluster,
    levels =
      unique(Cluster)[
        order(
          as.numeric(
            gsub("_.*$", "", unique(Cluster))
          )
        )
      ]
  )) %>%
  mutate(
    var = case_when(
      study == "dassim" ~ "This Study\nColonising",
      study == "horesh" & grepl("Colonising", source) ~ "Horesh et al.\nStool",
      TRUE ~ "Horesh et al.\nInvasive"
    ),
    var = factor(var, levels = sort(unique(var))[c(3, 2, 1)])
  ) %>%
  count(Cluster, var, res, .drop = FALSE) %>%
  filter(!(grepl("This Study", var) & res == "No ESBL")) %>%
  group_by(var, res) %>%
  mutate(
    cumsum = cumsum(n),
    sum = sum(n),
    prop = cumsum / sum
  ) %>%
  ggplot(aes(Cluster, prop, color = var, fill = var, group = var)) +
  geom_step() +
  theme_bw() +
  labs(
    y = "Cumulative proportion",
    x = "PopPUNK cluster", color = ""
  ) +
  geom_vline(
    xintercept = "1184",
    alpha = 0.8,
    linetype = "dotted"
  ) +
  scale_color_manual(
    values = c(
      viridis_pal(option = "cividis")(8)[c(2, 8, 6)]
    ),
    breaks = c(
      "This Study\nColonising",
      "Horesh et al.\nInvasive",
      "Horesh et al.\nStool"
    ),
    guide = "none"
  ) +
  theme(
    axis.text.y = element_text(size = 7),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 5),
    legend.position = "none"
  ) +
  facet_wrap(. ~ res, ncol = 1) +
  scale_x_discrete(
    breaks = as.character(c(1, seq(from = 100, to = 1400, by = 100))),
    expand = c(0.01, 0.01)
  ) -> b


(p1 + (plot_spacer() + p2 + plot_layout(ncol = 1)) + plot_layout(
  guides =
    "auto", widths = c(1.7, 1)
)) +
  a + (b + plot_layout(ncol = 1)) +
  plot_annotation(tag_levels = "A") -> pp_cluster_plots_stratified

pp_cluster_plots_stratified

if (write_figs) {
  ggsave(
    here("figures/ecoli-genomics/global_poppunk_clusterplot_strat.pdf"),
    pp_cluster_plots_stratified,
    width = 8,
    height = 12
  )
  ggsave(
    here("figures/ecoli-genomics/global_poppunk_clusterplot_strat.svg"),
    pp_cluster_plots_stratified,
    width = 8,
    height = 12
  )
}

df_all %>%
  filter(study %in% c("dassim", "horesh")) %>%
  group_by(study, Cluster, ST) %>%
  summarise(n_ST = n()) %>%
  group_by(study, Cluster) %>%
  mutate(n_cluster = sum(n_ST))  %>%
  summarise(
    cluster_size = n_cluster,
            STprop = n_ST/n_cluster,
            STpropstr =
              paste0("ST", ST," (",
              formatC(n_ST / n_cluster, 2, format = "f"), ")"),
    )  ->
  df_clusttab

df_clusttab$STpropstr <- 
  factor(df_clusttab$STpropstr, levels = unique(
    df_clusttab$STpropstr[order(df_clusttab$STprop,
                                decreasing = TRUE)])
    )

df_clusttab %>% 
  filter(STprop >= 0.01) %>% 
  pivot_wider(
    id_cols = c(Cluster),
   names_from = study,
   values_from = c(cluster_size,STpropstr),
   values_fn = list( STpropstr = function(x) paste(sort(x),collapse = ";"),
                     cluster_size = unique),
   values_fill = list(STpropstr = "-",
                      cluster_size = 0)
     ) %>% 
  arrange(desc(cluster_size_dassim), desc(cluster_size_horesh)) %>% 
  filter(cluster_size_dassim > 0) %>% 
  left_join(
    df_all %>%
    filter(study %in% c("dassim", "horesh")) %>%
    group_by(Cluster) %>% 
    summarise(pgroup = names(sort(table(Phylogroup), decreasing=TRUE)[1]))
  ) ->
  df_clusttab


df_all %>%
  filter(study %in% c("dassim", "horesh")) %>%
  group_by(study, Cluster, Continent) %>%
  summarise(n_continent = n()) %>%
  group_by(study, Cluster) %>%
  mutate(n_cluster = sum(n_continent))  %>%
  summarise(
    cluster_size = n_cluster,
            contprop = n_continent/n_cluster,
            contpropstr =
              paste0( Continent," (",
              formatC(n_continent / n_cluster, 2, format = "f"), ")"),
    )  ->
  df_clust.cont

df_clust.cont$contpropstr <- 
  factor(df_clust.cont$contpropstr, levels = unique(
    df_clust.cont$contpropstr[order(df_clust.cont$contprop,
                                decreasing = TRUE)])
    )

left_join(df_clusttab,
df_clust.cont %>% 
  pivot_wider(
    id_cols = c(Cluster),
   names_from = study,
   values_from = contpropstr,
   values_fn = function(x) paste(sort(x),collapse = ";"),
   values_fill = "-",
  )
) %>% 
  arrange(desc(cluster_size_dassim), desc(cluster_size_horesh)) ->
  df_clusttab

df_all %>%
  filter(study %in% c("dassim", "horesh")) %>%
  mutate(Pathotype = gsub(" \\(predicted\\)","", Pathotype)) %>% 
  group_by(study, Cluster, Pathotype) %>%
  summarise(n_path = n()) %>%
  group_by(study, Cluster) %>%
  mutate(n_cluster = sum(n_path))  %>%
  summarise(
    cluster_size = n_cluster,
            pathprop = n_path/n_cluster,
            pathpropstr =
              paste0( Pathotype," (",
              formatC(n_path / n_cluster, 2, format = "f"), ")"),
    )  ->
  df_clust.path

df_clust.path$pathpropstr <- 
  factor(df_clust.path$pathpropstr, levels = unique(
    df_clust.path$pathpropstr[order(df_clust.path$pathprop,
                                decreasing = TRUE)])
    )


left_join(df_clusttab,
df_clust.path %>% 
  pivot_wider(
    id_cols = c(Cluster),
   names_from = study,
   values_from = pathpropstr,
   values_fn = function(x) paste(sort(x),collapse = ";"),
   values_fill = "-",
  ) %>% 
  transmute(
    Cluster = Cluster,
    horesh_pathotype = horesh,
)) ->
  df_clusttab



df_clusttab %>%
  select(Cluster,
         pgroup,
         cluster_size_dassim,
         STpropstr_dassim,
         cluster_size_horesh,
         STpropstr_horesh,
         horesh,
         horesh_pathotype
         ) %>%
  kbl(caption =
        "Cluster assigment of samples from this study compared to Horesh collection",
      col.names = c("Cluster","Phylogroup",
                    rep(c(
                    "n isolates",
                    "STs"),2),"Location", "Pathotype")
      ) %>%
  kable_classic(full_width = FALSE) %>% 
  add_header_above(c(" " = 2, "This Study" = 2, "Horesh Collection" = 4))

if (write_figs) {
 write_csv(df_clusttab %>% 
             transmute(
               Cluster = Cluster,
               Phylogroup = pgroup,
               n_this_study = cluster_size_dassim,
               STs_this_study = STpropstr_dassim,
               n_horesh = cluster_size_horesh,
               STs_horesh = STpropstr_horesh,
               location_horesh = horesh,
               pathotype_horesh = horesh_pathotype),
           here("tables/ecoli-genomics/cluster_stats.csv")
 )

}

# bits and bobs to calculate for manuscript

# how many clusters?

dassimEcoli_BTEcoli.accession %>% 
  group_by(Cluster) %>% 
  tally() %>% 
  nrow()

dassimEcoli_BTEcoli.accession %>% 
  group_by(Cluster) %>% 
  tally() %>% 
  summarise(
    median_size = median(n),
    lq = quantile(n, 0.25),
    uq = quantile(n, 0.75)
  )


dassimEcoli_BTEcoli.accession %>% 
  group_by(Cluster) %>% 
  tally() %>% 
  arrange(-n)

dassimEcoli_BTEcoli.accession %>% 
  group_by(Cluster, ST) %>% 
  tally() %>% 
  arrange(-n)

df_all %>% 
  filter(study == "dassim") %>% 
  mutate(numeric_cluster = as.numeric(
    gsub("_.*$","",Cluster)),
    clust_50orbelow =
           case_when(
             numeric_cluster <= 50 ~ "top 50",
             numeric_cluster <=1184 ~ "not top 50 but in horsh cluster",
             TRUE ~ "new cluster")
  ) %>% 
  group_by(clust_50orbelow) %>% 
  tally() %>% 
  mutate(
    tot = sum(n),
    prop = n/sum(n))

btESBL_ecoli_horesh_metadata %>% 
  mutate(clust_50orbelow =
           case_when(
             PopPUNK <= 50 ~ "top 50",
             PopPUNK <=1184 ~ "not top 50 but in horsh cluster",
             TRUE ~ "new cluster")
         ) %>% 
  group_by(clust_50orbelow) %>% 
  tally() %>% 
  mutate(
    tot = sum(n),
    prop = n/sum(n)) 

amr %>% 
  select(lane, class) %>% 
  unique() %>% 
  ungroup() %>% 
  mutate(n_tot = length(unique(lane))) %>% 
  group_by(class) %>% 
  summarise(n = n(),
            n_tot = max(n_tot),
            prop = n/n_tot
  )


  df_all %>% 
    filter(study == "dassim", ST != "Novel") %>% 
    group_by(ST) %>% 
    summarise(n = n()) %>% 
    ungroup() %>% 
    summarise(
      m = median(n),
      min = min(n),
      lq = quantile(n, 0.25),
      uq = quantile(n, 0.75),
      max = max(n)
    )

  #QRDR stuff    

  btESBL_amrgenes %>% 
  semi_join(dassimEcoli_BTEcoli.accession, 
            by = c( "lane")) %>% 
  select(-genus) %>% 
  rename(gene = ref_seq) %>% 
  # add in QRDR mutations
  bind_rows(
    btESBL_qrdr_mutations %>% 
      filter(genus == "E. coli") %>% 
      select(-genus) %>% 
      semi_join(
        btESBL_CARD_qrdr_mutations, 
        by = c("variant" = "variant",
              "gene" =  "gene")
      ) %>% 
      unique() %>% 
      mutate(gene =
               gsub("(^.{1})", '\\U\\1',
                    gene,
                    perl = TRUE))
      ) %>% 
  filter(!grepl("AmpH|AMPH|AmpC|MrdA|MefB", gene)) %>% 
  # add in beta-lactamases
  left_join(
    select(btESBL_NCBI_phenotypic_bl, allele_name, class),
    by = c("gene" = "allele_name")) %>% 
  mutate(class = case_when(
    str_detect(gene, quinolone) ~ "Quinolone",
    str_detect(gene, tetracycline) ~ "Tetracycline",
    str_detect(gene, sulphonamide) ~ "Sulphonamide",
    str_detect(gene, aminoglycoside) ~ "Aminoglycoside",
    str_detect(gene, streptothricin) ~ "Streptothricin",
    str_detect(gene, macrolide) ~ "Macrolide",
    str_detect(gene, fosfomycin) ~ "Fosfomycin",
    str_detect(gene, chloramphenicol) ~ "Chloramphenicol",
    str_detect(gene, rifampicin) ~ "Rifampicin",
    str_detect(gene,trimethoprim) ~ "Trimethoprim",
    str_detect(gene,ESBL) ~ "ESBL",
    str_detect(gene,penicillinase) ~ "Penicillinase",
    str_detect(gene,ampc) ~ "AmpC",
    TRUE ~ class
  )) %>% 
    filter(class == "Quinolone") -> quinolone.amr

  quinolone.amr %>% 
    mutate(gene = paste0(gene, 
                         if_else(
                           !is.na(variant), 
                           variant, ""))) %>% 
    select(lane,gene) %>% 
    unique() %>% 
    pivot_wider(id_cols = lane,
                names_from = gene,
                values_from = gene,
                values_fn = length,
                values_fill = 0) %>% 
    summarise(across(!contains("lane"), sum))


quinolone.amr %>%
  group_by(lane) %>%
  summarise(qrdr = any(grepl("Gyr|Par", gene)),
            qnr = any(grepl("Qnr", gene)),
            qep = any(grepl("Qep", gene))) %>%
  ungroup() %>%
  summarise(qrdr_alone = sum(qrdr & !qnr & !qep),
            qnr_alone = sum(qnr & !qrdr & !qep),
            qrdr_qnr = sum(qrdr & qnr),
            qrdr_qep = sum(qrdr & qep),
            qnr_qep = sum(qnr & qep)
  )

btESBL_plasmidreplicons %>% 
  filter(species == "E. coli") %>% 
   mutate(ref_seq = gsub("_.+$","", ref_seq),
         ref_seq = gsub("\\.[0-9]","", ref_seq),
         ref_seq = case_when(
           grepl("Col", ref_seq) ~ "Col",
           grepl("rep", ref_seq) ~ "rep",
           grepl("^FIA", ref_seq) ~ "IncFIA",
           !grepl("Inc", ref_seq) ~ "Other",
           TRUE ~ ref_seq
          )) %>% 
  filter(grepl("Inc", ref_seq)) %>%  
  arrange(fct_infreq(ref_seq)) %>% 
  pivot_wider(id_cols = lane,
              values_from = ref_seq,
              names_from = ref_seq,
              values_fn = length,
              values_fill = 0) %>% 
  mutate(across(where(is.numeric), ~ if_else(.x >= 1, "Present", "Absent"))) %>% 
  as.data.frame() ->
  btESBL_plasmidreplicons_onehot

rownames(btESBL_plasmidreplicons_onehot) <- btESBL_plasmidreplicons_onehot$lane



ggtree(btESBL_coregene_tree_esco) %>% 
  gheatmap(select(btESBL_plasmidreplicons_onehot, -lane),
           font.size = 2.5,
           color = NA,
           colnames_angle = 90,
           offset = 0,
           colnames_offset_y = 3,
           colnames_position = "top",
           hjust = 0) +
  ylim(NA,550) +
  labs(fill = "Plasmid\nReplicon") +
  scale_fill_manual(values = c("white", "darkgrey"), 
                     na.translate = FALSE) -> plas_replicons_map_to_tree
  plas_replicons_map_to_tree

if (write_figs) {

  ggsave(
    here("figures/ecoli-genomics/ecoli_plasmid_replicon_maptotree.pdf"),
    plas_replicons_map_to_tree,
    width = 6,
    height = 6)
  ggsave(
    here("figures/ecoli-genomics/ecoli_plasmid_replicon_maptotree.svg"),
    plas_replicons_map_to_tree,
    width = 6,
    height = 6)

}

ST410 trees

btESBL_ecoli_st410_metadata %>% 
  filter(accession %in% btESBL_ecoli_globalst410_tree$tip.label) %>% 
  bind_rows(
    btESBL_sequence_sample_metadata %>% 
      filter(lane %in% btESBL_ecoli_globalst410_tree$tip.label) %>% 
      transmute(accession = lane,
                `Collection Year` = year(data_date),
                host = "Homo sapiens",
                source = "Stool",
                Country = "Malawi")
  ) %>% 
  rename(
    Year = `Collection Year`) %>% 
  as.data.frame() -> st410_metadata

if (write_figs) {
st410_metadata %>% 
  filter(Country != "Malawi") %>% 
  transmute(accession = accession,
            country = Country,
            year = Year) %>% 
  write_csv(here("tables/ecoli-genomics/SUPP_T_st410_accession_numbers_and_metadata.csv"))
}

st410_metadata$malawi <- st410_metadata$Country == "Malawi"


st410_metadata$Continent <- 
  countrycode(
    sourcevar = st410_metadata$Country,
    origin = "country.name",
    destination = "un.region.name"
  )

st410_metadata %>%
  mutate(
    Continent = case_when(
      Country %in% c("Africa", "Asia", "Europe") ~
      Country,
      Country == "Blantyre" ~ "Africa",
      TRUE ~ Continent
    )
  ) ->
st410_metadata

rownames(st410_metadata) <- st410_metadata$accession


btESBL_ecoli_st410_amr %>%
  pivot_wider(
    names_from = gene,
    values_from = gene,
    values_fn = length,
    values_fill = 0
  ) %>%
  mutate(across(everything(), as.character)) %>%
  as.data.frame() ->
amr.ariba410


rownames(amr.ariba410) <- amr.ariba410$name




btESBL_ecoli_st410_amr %>%
  bind_rows(
    btESBL_amrgenes %>%
      filter(genus == "E. coli") %>%
      select(-genus) %>%
      transmute(
        name = lane,
        gene = ref_seq
      )
  ) %>%
  filter(name %in% btESBL_ecoli_globalst410_tree$tip.label) %>%
  left_join(
    select(btESBL_NCBI_phenotypic_bl, allele_name, class),
    by = c("gene" = "allele_name")
  ) %>%
  mutate(class = case_when(
    grepl("CMY", name) ~ "ESBL/AmpC",
    class == "ESBL" ~ "ESBL/AmpC",
    TRUE ~ class
  )) %>%
  filter(class %in% c("ESBL/AmpC", "Carbapenemase")) %>%
  mutate(gene = gsub("_", "-", gene)) %>%
  arrange(class, fct_infreq(gene)) %>%
  pivot_wider(
    id_cols = "name",
    names_from = gene,
    values_from = class
  ) %>%
  as.data.frame() ->
esbl_cpe.ariba410
rownames(esbl_cpe.ariba410) <- esbl_cpe.ariba410$name

btESBL_ecoli_st410_plasmids %>%
  bind_rows(
    btESBL_plasmidreplicons %>%
      mutate(
        ref_seq = gsub("\\..*$", "", ref_seq),
        ref_seq = gsub("_.*$", "", ref_seq),
        ref_seq = gsub("^FIA", "IncFIA", ref_seq)
      ) %>%
      filter(species == "E. coli") %>%
      select(-species) %>%
      rename(name = lane)
  ) %>%
  filter(name %in% btESBL_ecoli_globalst410_tree$tip.label) %>%
  filter(grepl("Inc", ref_seq)) %>%
  unique() %>%
  arrange(fct_infreq(ref_seq)) %>%
  pivot_wider(
    id_cols = name,
    names_from = ref_seq,
    values_from = ref_seq,
    values_fn = function(x) {
      return("Present")
    },
    values_fill = NA_character_
  ) %>%
  as.data.frame() -> st410_plasm_onehot
rownames(st410_plasm_onehot) <- st410_plasm_onehot$name


(
(  
  (ggtree(
    btESBL_ecoli_globalst410_tree) %<+% 
     (select(st410_metadata, accession, Country) %>% 
        mutate(Country = if_else(Country == "Malawi", "Malawi",NA_character_)))  %>%
     gheatmap(
       select(st410_metadata, `Year`),
       font.size = 4,
       width = 0.03,
       colnames_position = "top",
       color = NA,
       colnames_offset_y = 5,
       colnames_angle = 90,
       hjust = 0
     )  + 
     scale_fill_viridis(name = "Year", na.value = "white") +
     new_scale_fill()) %>% 
  gheatmap(
    select(st410_metadata, Continent),
    font.size = 4,
    width = 0.03,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 5,
    colnames_angle = 90,
    offset = 0.000008,
    hjust = 0
  ) + 
  scale_fill_viridis_d(option = "plasma", name = "Continent", na.translate =
                       FALSE) + 
  geom_tippoint(aes(color = Country), na.rm = TRUE, size = 0.5) + 
  scale_color_manual(values  = "red", na.translate = FALSE) +
  new_scale_fill() 
) %>% 
  gheatmap(
    select(esbl_cpe.ariba410, -name),
    font.size = 3,
    width = 0.6,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 5,
    colnames_angle = 90,
    offset = 0.000016,
    hjust = 0
  ) + scale_fill_manual(values = viridis_pal(option = "magma")(5)[c(2,3)],
                        na.translate = FALSE,
                        name = "Beta-\nlactamase") +
  new_scale_fill()
) + # %>% 
  # gheatmap(
  #   select(st410_plasm_onehot, -name),
  #   font.size = 3,
  #   width = 0.5,
  #   colnames_position = "top",
  #   color = NA,
  #   colnames_offset_y = 5,
  #   colnames_angle = 90,
  #   offset = 0.0026,
  #   hjust = 0
  # ) + scale_fill_manual(values = viridis_pal(option = "cividis")(7)[2],
  #                       na.translate = FALSE,
  #                       name = "Plasmid\nreplicon") +
  ylim(NA, 420) +
  geom_treescale(x = 0.00005, y = 250, offset = 5) +
   geom_highlight(node = 391 ,alpha = 0.3) +
  geom_cladelabel(
    node = 391,
    label = "B4/H24RxC",
    angle = 90,
    offset = 0,
    offset.text = -0.000032,
    barsize = 0,
    hjust = 0.5,
    fontsize = 3,
    color = "steelblue"
  )  -> st410_globaltree_plot

if (write_figs) {

ggsave(here("figures/ecoli-genomics/st410_globaltreeplot.svg"),
       st410_globaltree_plot, width = 9, height = 10)

ggsave(here("figures/ecoli-genomics/st410_globaltreeplot.pdf"),
       st410_globaltree_plot, width = 9, height = 10)
}
st410_globaltree_plot

tree_subset(
  btESBL_ecoli_globalst410_tree, 390, levels_back = 0)  -> st410_subtree

esbl_cpe.ariba410 %>%
  subset(name %in% st410_subtree$tip.label) %>%
  as.data.frame() ->
esbl_cpe.ariba410_subset

  esbl_cpe.ariba410_subset[
    c(TRUE, apply(esbl_cpe.ariba410_subset[-1], 2, function(x) !all(is.na(x))))
  ] ->
  esbl_cpe.ariba410_subset

btESBL_ecoli_st410_plasmids %>%
  bind_rows(
    btESBL_plasmidreplicons %>%
      mutate(
        ref_seq = gsub("\\..*$", "", ref_seq),
        ref_seq = gsub("_.*$", "", ref_seq),
        ref_seq = gsub("^FIA", "IncFIA", ref_seq)
      ) %>%
      filter(species == "E. coli") %>%
      select(-species) %>%
      rename(name = lane)
  ) %>%
  filter(grepl("Inc", ref_seq)) %>%
  filter(name %in% st410_subtree$tip.label) %>%
  unique() %>%
  group_by(ref_seq) %>%
  mutate(n = n()) %>%
  filter(n > 2) %>%
  ungroup() %>%
  select(-n) %>%
  arrange(fct_infreq(ref_seq)) %>%
  pivot_wider(
    id_cols = name,
    names_from = ref_seq,
    values_from = ref_seq,
    values_fn = function(x) {
      return("Present")
    },
    values_fill = NA_character_
  ) %>%
  as.data.frame() -> st410_plasm_onehot_subset
rownames(st410_plasm_onehot_subset) <- st410_plasm_onehot_subset$name

(
  (
    (ggtree(st410_subtree, size = 0.3) %<+% (select(st410_metadata, accession, Country) %>% mutate(Country = if_else(Country == "Malawi", "Malawi", NA_character_))) %>%
      gheatmap(
        select(st410_metadata, `Year`),
        font.size = 4,
        width = 0.06,
        colnames_position = "top",
        color = NA,
        colnames_offset_y = 0,
        colnames_angle = 90,
        hjust = 0
      ) +
      scale_fill_viridis(
        name = "Year", limits = c(1978, 2020),
        na.value = "white"
      ) +
      new_scale_fill()) %>%
      gheatmap(
        select(st410_metadata, Continent),
        font.size = 4,
        width = 0.06,
        colnames_position = "top",
        color = NA,
        colnames_offset_y = 0,
        colnames_angle = 90,
        offset = 0.0000015,
        hjust = 0
      ) +
      scale_fill_viridis_d(
        option = "plasma", name = "Continent",
        na.translate = FALSE
      ) +
      geom_point2(aes(subset = as.numeric(label) < 95 & !isTip), size = 0.5) +
      geom_tippoint(aes(color = Country), na.rm = TRUE, size = 1) +
      scale_color_manual(values = "red", na.translate = FALSE) +
      new_scale_fill()
  ) %>%
    gheatmap(
      select(esbl_cpe.ariba410_subset, -name),
      font.size = 2.5,
      width = 0.4,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 0,
      colnames_angle = 90,
      offset = 0.0000031,
      hjust = 0
    ) + scale_fill_manual(
      values = viridis_pal(option = "magma")(5)[c(2, 3)],
      na.translate = FALSE,
      name = "Beta-\nlactamase"
    ) +
    new_scale_fill()
) %>%
  gheatmap(
    select(st410_plasm_onehot_subset, -name),
    font.size = 2.5,
    width = 0.45,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = 0.0000125,
    hjust = 0
  ) + scale_fill_manual(
    values = viridis_pal(option = "cividis")(7)[2],
    na.translate = FALSE,
    name = "Plasmid\nreplicon"
  ) +
  ylim(NA, 100) +
  geom_treescale(x = 0.000004, y = 66, offset = 2, width = 1e-6) ->
st410_subtree_plot

btESBL_ecoli_st167_plasmids %>%
  bind_rows(
    btESBL_plasmidreplicons %>%
      mutate(
        ref_seq = gsub("\\..*$", "", ref_seq),
        ref_seq = gsub("_.*$", "", ref_seq),
        ref_seq = gsub("^FIA", "IncFIA", ref_seq)
      ) %>%
      filter(species == "E. coli") %>%
      select(-species) %>%
      rename(name = lane)
  ) %>%
  filter(grepl("Inc", ref_seq)) %>%
  filter(name %in% btESBL_ecoli_globalst167_tree$tip.label) %>%
  unique() %>%
  pivot_wider(
    id_cols = name,
    names_from = ref_seq,
    values_from = ref_seq,
    values_fn = function(x) {
      return("Present")
    },
    values_fill = NA_character_
  ) %>%
  as.data.frame() -> st167_plasm_onehot

rownames(st167_plasm_onehot) <- st167_plasm_onehot$name





btESBL_ecoli_st167_metadata %>%
  filter(accession %in% btESBL_ecoli_globalst167_tree$tip.label) %>%
  bind_rows(
    btESBL_sequence_sample_metadata %>%
      filter(lane %in% btESBL_ecoli_globalst167_tree$tip.label) %>%
      transmute(
        accession = lane,
        `Collection Year` = year(data_date),
        `Source Niche` = "Human",
        Country = "Malawi"
      )
  ) %>%
  rename(Year = `Collection Year`) %>%
  as.data.frame() -> st167_metadata

if (write_figs) {
  st167_metadata %>%
    filter(Country != "Malawi") %>%
    transmute(
      accession = accession,
      source = `Source Niche`,
      country = Country,
      year = Year
    ) %>%
    write_csv(
      here("tables/ecoli-genomics/st167_accession_numbers_and_metadata.csv")
    )
}

st167_metadata$malawi <- st167_metadata$Country == "Malawi"

st167_metadata$Continent <-
  countrycode(
    sourcevar = st167_metadata$Country,
    origin = "country.name",
    destination = "un.region.name"
  )
rownames(st167_metadata) <- st167_metadata$accession




btESBL_ecoli_st167_amr %>% 
  bind_rows(
    btESBL_amrgenes %>%
      filter(genus == "E. coli") %>%
      select(-genus) %>%
      transmute(
        name = lane,
        gene = ref_seq
      )
  ) %>%
  filter(name %in% btESBL_ecoli_globalst167_tree$tip.label) %>%
  pivot_wider(names_from = gene, 
              values_from = gene,
              values_fn = length,
              values_fill = 0) %>%
  mutate(across(everything(), as.character)) %>% 
  as.data.frame() ->
  amr.ariba

rownames(amr.ariba) <- amr.ariba$name



amr.ariba %>% 
  rename(sample = name) %>% 
  pivot_longer(-sample) %>% 
  filter(value == 1) %>% 
  left_join(
    select(btESBL_NCBI_phenotypic_bl,
           allele_name, 
           class),
    by = c("name" = "allele_name")
  ) %>%  
  mutate(class = case_when(
    grepl("CMY", name) ~ "ESBL/AmpC",
    class == "ESBL" ~ "ESBL/AmpC",
    TRUE ~ class
        )) %>% 
  filter(class %in% c("ESBL/AmpC", "Carbapenemase") & value == 1) %>%
  mutate(name = gsub("_", "-", name)) %>% 
  arrange(class, fct_infreq(name)) %>% 
  pivot_wider(id_cols = "sample",
              names_from = name,
              values_from = class) %>% 
  as.data.frame() ->
  esbl_cpe.ariba

rownames(esbl_cpe.ariba) <- esbl_cpe.ariba$sample

(
(
(ggtree(btESBL_ecoli_globalst167_tree) %<+% 
   (select(st167_metadata, accession, Country) %>% 
      mutate(Country = if_else(Country == "Malawi", "Malawi",NA_character_)))  %>%
    gheatmap(
      select(st167_metadata, `Year`),
      font.size = 4,
      width = 0.03,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 0,
      colnames_angle = 90,
      hjust = 0
    )  + 
    scale_fill_viridis(name = "Year", na.value = "white") +
    new_scale_fill()) %>% 
  gheatmap(
    select(st167_metadata, Continent),
    font.size = 4,
    width = 0.03,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = 0.000003,
    hjust = 0
  ) + 
  scale_fill_viridis_d(option = "plasma", name = "Continent", na.translate =
                       FALSE) + 
  geom_tippoint(aes(color = Country), na.rm = TRUE, size = 1) + 
  scale_color_manual(values  = "red", na.translate = FALSE) +
    new_scale_fill() 
) %>% 
  gheatmap(
  select(esbl_cpe.ariba, -sample),
  font.size = 3,
  width = (ncol(esbl_cpe.ariba) -1)*0.03,
  colnames_position = "top",
  color = NA,
  colnames_offset_y = 0,
  colnames_angle = 90,
  offset = 0.000006,
  hjust = 0
) + scale_fill_manual(values = viridis_pal(option = "magma")(5)[c(2,3)],
                      na.translate = FALSE,
                      name = "Beta-\nlactamase") +
  new_scale_fill()
)  +  
  ylim(NA, 150) +
    geom_highlight(node = 142, alpha = 0.3, fill = "steelblue",
                   color = NA) + 
  geom_treescale(x = 0.000015, y = 160*0.666, offset = 2, width = 0.00001) ->
  st167_globaltree_plot

st167_globaltree_plot

if (write_figs) {

ggsave(here("figures/ecoli-genomics/st167_globaltreeplot.svg"),
       st167_globaltree_plot, width = 8, height = 10)

ggsave(here("figures/ecoli-genomics/st167_globaltreeplot.pdf"),
       st167_globaltree_plot, width = 8, height = 10)
}

tree_subset(
  btESBL_ecoli_globalst167_tree, 142, levels_back = 0) -> st167_subtree

esbl_cpe.ariba %>%
  subset(sample %in% st167_subtree$tip.label) %>%
  as.data.frame() ->
esbl_cpe.ariba167_subset

esbl_cpe.ariba167_subset[
  c(TRUE, apply(esbl_cpe.ariba167_subset[-1], 2, function(x) !all(is.na(x))))
] ->
esbl_cpe.ariba167_subset


btESBL_ecoli_st167_plasmids %>%
  bind_rows(
    btESBL_plasmidreplicons %>%
      mutate(
        ref_seq = gsub("\\..*$", "", ref_seq),
        ref_seq = gsub("_.*$", "", ref_seq),
        ref_seq = gsub("^FIA", "IncFIA", ref_seq)
      ) %>%
      filter(species == "E. coli") %>%
      select(-species) %>%
      rename(name = lane)
  ) %>%
  filter(grepl("Inc", ref_seq)) %>%
  filter(name %in% st167_subtree$tip.label) %>%
  unique() %>%
  group_by(ref_seq) %>%
  mutate(n = n()) %>%
  filter(n > 2) %>%
  ungroup() %>%
  select(-n) %>%
  arrange(fct_infreq(ref_seq)) %>%
  pivot_wider(
    id_cols = name,
    names_from = ref_seq,
    values_from = ref_seq,
    values_fn = function(x) {
      return("Present")
    },
    values_fill = NA_character_
  ) %>%
  as.data.frame() -> st167_plasm_onehot_subset

rownames(st167_plasm_onehot_subset) <- st167_plasm_onehot_subset$name

cont_cols <- viridis_pal(option = "plasma")(4)
names(cont_cols) <- c("Africa","Americas", "Asia","Europe")

  (
  (  
  (ggtree(st167_subtree, size = 0.3) %<+% (select(st167_metadata, accession, Country) %>% mutate(Country = if_else(Country == "Malawi", "Malawi",NA_character_)))  %>%
     gheatmap(
       select(st167_metadata, `Year`),
       font.size = 4,
       width = 0.06,
       colnames_position = "top",
       color = NA,
       colnames_offset_y = 0,
       colnames_angle = 90,
       hjust = 0
     )  + 
     scale_fill_viridis(name = "Year", limits = c(1978,2020),
                        na.value = "white") +
     new_scale_fill()) %>% 
    gheatmap(
      select(st167_metadata, Continent),
      font.size = 4,
      width = 0.06,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 0,
      colnames_angle = 90,
      offset = 0.000005,
      hjust = 0
    ) + 
    scale_fill_manual(values = cont_cols,
                      name = "Continent",
                      na.translate =
                      FALSE) + 
      geom_point2(aes(subset = as.numeric(label) < 95 & !isTip), size = 1) +
    geom_tippoint(aes(color = Country), na.rm = TRUE, size = 1) + 
    scale_color_manual(values  = "red", na.translate = FALSE) +
    new_scale_fill() 
) %>% 
  gheatmap(
    select(esbl_cpe.ariba167_subset, -sample),
    font.size = 2.5,
    width = 0.06*(ncol(esbl_cpe.ariba167_subset) - 1),
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = 0.00001,
    hjust = 0
  ) + scale_fill_manual(values = viridis_pal(option = "magma")(5)[c(2,3)],
                        na.translate = FALSE,
                        name = "Beta-\nlactamase") +
      new_scale_fill() 
    ) %>%   gheatmap(
    select(st167_plasm_onehot_subset, -name),
    font.size = 2.5,
    width = 0.045*(ncol(st167_plasm_onehot_subset) - 1) + 0.1,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 0,
    colnames_angle = 90,
    offset = 0.000029,
    hjust = 0
  ) + scale_fill_manual(values = viridis_pal(option = "cividis")(7)[2],
                        na.translate = FALSE,
                        name = "Plasmid\nreplicon") +
  ylim(NA, 60) +
  geom_treescale(x = 0.00001, y = 60*0.666, offset = 1, width = 1e-5) ->
    st167_subtree_plot

if (write_figs) {
  btESBL_ecoli_st131_metadata %>%
    filter(Country != "Malawi") %>%
    write_csv(here("tables/ecoli-genomics/SUPP_T_st131_accession_numbers_and_metadata.csv"))
}
btESBL_ecoli_st131_metadata$malawi <- btESBL_ecoli_st131_metadata$Country == "Malawi"


btESBL_ecoli_st131_metadata$Continent <-
  countrycode(
    sourcevar = btESBL_ecoli_st131_metadata$Country,
    origin = "country.name",
    destination = "un.region.name"
  )

btESBL_ecoli_st131_metadata$Continent[
  btESBL_ecoli_st131_metadata$Country == "Czech"
] <- "Europe"

btESBL_ecoli_st131_metadata$Continent[
  btESBL_ecoli_st131_metadata$Country == "Taiwan"
] <- "Asia"


rownames(btESBL_ecoli_st131_metadata) <- btESBL_ecoli_st131_metadata$Accession_number


btESBL_ecoli_st131_amr %>%
  pivot_wider(
    names_from = gene,
    values_from = gene,
    values_fn = length,
    values_fill = 0
  ) %>%
  mutate(across(everything(), as.character)) %>%
  as.data.frame() ->
amr.ariba131



rownames(amr.ariba131) <- amr.ariba131$name

btESBL_ecoli_st131_plasmids %>%
  filter(grepl("Inc", ref_seq)) %>%
  unique() %>%
  group_by(ref_seq) %>%
  mutate(n = n()) %>%
  filter(n >= 10) %>%
  ungroup() %>%
  select(-n) %>%
  arrange(fct_infreq(ref_seq)) %>%
  pivot_wider(
    id_cols = name,
    names_from = ref_seq,
    values_from = ref_seq,
    values_fn = function(x) {
      return("Present")
    },
    values_fill = NA_character_
  ) %>%
  as.data.frame() -> st131_plasm_onehot

rownames(st131_plasm_onehot) <- st131_plasm_onehot$name

esbl_cols <- viridis_pal(option = "magma")(5)[c(2, 3)]

names(esbl_cols) <- c("Carbapenemase", "ESBL/AmpC")

btESBL_ecoli_st131_amr %>%
  left_join(
    select(btESBL_NCBI_phenotypic_bl, allele_name, class),
    by = c("gene" = "allele_name")
  ) %>%
  mutate(class = case_when(
    grepl("CMY", name) ~ "ESBL/AmpC",
    class == "ESBL" ~ "ESBL/AmpC",
    TRUE ~ class
  )) %>%
  filter(class %in% c("ESBL/AmpC", "Carbapenemase")) %>%
  mutate(gene = gsub("_", "-", gene)) %>%
  arrange(class, fct_infreq(gene)) %>%
  pivot_wider(
    id_cols = "name",
    names_from = gene,
    values_from = class
  ) %>%
  as.data.frame() ->
#  select(-`TEM-15`) ->    # not in df?
esbl_cpe.ariba131
rownames(esbl_cpe.ariba131) <- esbl_cpe.ariba131$name

(
(
(ggtree(btESBL_ecoli_globalst131_tree, size = 0.2) %<+% 
   (select(btESBL_ecoli_st131_metadata, Accession_number, Country) %>% 
      mutate(Country = if_else(Country == "Malawi", "Malawi",NA_character_)))  %>%
    gheatmap(
      select(btESBL_ecoli_st131_metadata, `Year`),
      font.size = 4,
      width = 0.035,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 5,
      colnames_angle = 90,
      hjust = 0
    )  + 
     scale_fill_viridis_c(name = "Year", limits = c(1965,2020),
                          na.value = "white") +
    new_scale_fill()) %>% 
  gheatmap(
    select(btESBL_ecoli_st131_metadata, Continent),
    font.size = 4,
    width = 0.035,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 5,
    colnames_angle = 90,
    offset = 0.000005,
    hjust = 0
  ) + 
  scale_fill_viridis_d(option = "plasma", name = "Continent", 
                       na.translate = FALSE) + 
  geom_tippoint(aes(color = Country), na.rm = TRUE, size = 1) + 
  scale_color_manual(values  = "red", na.translate = FALSE) +
    new_scale_fill() 
) %>% 
  gheatmap(
  select(esbl_cpe.ariba131, -name),
  font.size = 2.5,
  width = (ncol(esbl_cpe.ariba131) -1)*0.03,
  colnames_position = "top",
  color = NA,
  colnames_offset_y = 5,
  colnames_angle = 90,
  offset = 0.0000103,
  hjust = 0
) + scale_fill_manual(values = esbl_cols,
                      na.translate = FALSE,
                      name = "Beta-\nlactamase") +
  new_scale_fill()
)   %>%   gheatmap(
    select(st131_plasm_onehot, -name),
    font.size = 2.5,
    width = 0.026*(ncol(st131_plasm_onehot) - 1) + 0.1,
    colnames_position = "top",
    color = NA,
    colnames_offset_y = 5,
    colnames_angle = 90,
    offset = 0.000059,
    hjust = 0
  ) + scale_fill_manual(values = viridis_pal(option = "cividis")(7)[2],
                        na.translate = FALSE,
                        name = "Plasmid\nreplicon") + 
  ylim(NA, 1000) +
  geom_treescale(x = 0.00002, y = 1000*0.666, offset = 5, width = 0.00001)   +
  theme(legend.position = "none") + 
      geom_point2(aes(subset = as.numeric(label) < 95 & !isTip), size = 0.5)  ->
  st131_globaltree_plot

st131_globaltree_plot

st410_subtree_plot + st167_subtree_plot  + st131_globaltree_plot + 
  plot_layout(widths = c(1,1.5),
              design = "AAAAACCCCCCCC
                        BBBBBCCCCCCCC") +
plot_annotation(tag_levels = "A") +
plot_layout(guides = "collect") -> st410_167_131_subtree_plot

st410_167_131_subtree_plot

if (write_figs) {

ggsave(here("figures/ecoli-genomics/st410_167_131_subtreeplot.svg"),
       st410_167_131_subtree_plot, width = 12, height = 8.2)

ggsave(here("figures/ecoli-genomics/st410_167_131_subtreeplot.pdf"),
       st410_167_131_subtree_plot, width = 12, height = 8.2)
}

# ----------------------------------------------------------------------
# plot subtree of ST131

tree_subset(btESBL_ecoli_globalst131_tree, 1090) -> st131_subtree


esbl_cpe.ariba131 %>%
  subset(name %in% st131_subtree$tip.label) %>%
  as.data.frame() ->
esbl_cpe.ariba131_subset

esbl_cpe.ariba131_subset[
  c(TRUE, apply(esbl_cpe.ariba131_subset[-1], 2, function(x) !all(is.na(x))))
] ->
esbl_cpe.ariba131_subset


st131_plasm_onehot %>%
  subset(name %in% st131_subtree$tip.label) %>%
  as.data.frame() -> st131_plasm_onehot_subset


st131_plasm_onehot_subset[
  c(TRUE, apply(st131_plasm_onehot_subset[-1], 2, function(x) !all(is.na(x))))
] ->

st131_plasm_onehot_subset

(
  (
    (ggtree(tree_subset(btESBL_ecoli_globalst131_tree, 1090)) %<+%
      (select(btESBL_ecoli_st131_metadata, Accession_number, Country) %>%
        mutate(Country = if_else(Country == "Malawi",
          "Malawi", NA_character_
        ))) %>%
      gheatmap(
        select(btESBL_ecoli_st131_metadata, `Year`),
        font.size = 4,
        width = 0.035,
        colnames_position = "top",
        color = NA,
        colnames_offset_y = 0,
        colnames_angle = 90,
        hjust = 0
      ) +
      scale_fill_viridis_c(
        name = "Year", limits = c(1965, 2020),
        na.value = "white"
      ) +
      new_scale_fill()) %>%
      gheatmap(
        select(btESBL_ecoli_st131_metadata, Continent),
        font.size = 4,
        width = 0.03,
        colnames_position = "top",
        color = NA,
        colnames_offset_y = 0,
        colnames_angle = 90,
        offset = 0.0000008,
        hjust = 0
      ) +
      scale_fill_viridis_d(
        option = "plasma", name = "Continent",
        na.translate = FALSE
      ) +
      geom_tippoint(aes(color = Country), na.rm = TRUE, size = 1) +
      scale_color_manual(values = "red", na.translate = FALSE) +
      new_scale_fill()
  ) %>%
    gheatmap(
      select(esbl_cpe.ariba131_subset, -name),
      font.size = 2.5,
      width = (ncol(esbl_cpe.ariba131_subset) - 1) * 0.03,
      colnames_position = "top",
      color = NA,
      colnames_offset_y = 0,
      colnames_angle = 90,
      offset = 0.0000015,
      hjust = 0
    ) + scale_fill_manual(
      values = esbl_cols,
      na.translate = FALSE,
      name = "Beta-\nlactamase"
    ) +
    new_scale_fill()
) %>% gheatmap(
  select(st131_plasm_onehot_subset, -name),
  font.size = 2.5,
  width = 0.026 * (ncol(st131_plasm_onehot_subset) - 1) + 0.1,
  colnames_position = "top",
  color = NA,
  colnames_offset_y = 0,
  colnames_angle = 90,
  offset = 0.0000032,
  hjust = 0
) + scale_fill_manual(
  values = viridis_pal(option = "cividis")(7)[2],
  na.translate = FALSE,
  name = "Plasmid\nreplicon"
) +
  ylim(NA, 80) +
  geom_treescale(x = 0.000001, y = 55, offset = 1, width = 0.000001) +
  geom_point2(aes(subset = as.numeric(label) < 95 & !isTip), size = 0.5) ->
st131_subtree_plot

st131_subtree_plot

if (write_figs) {

ggsave(here("figures/ecoli-genomics/st131_subtreeplot.svg"),
       st131_subtree_plot, width = 7, height = 6.5)

ggsave(here("figures/ecoli-genomics/st131_subtreeplot.pdf"),
       st131_subtree_plot, width = 7, height = 6.5)
}

Reproducibility - packages, R version and OS used

sessionInfo()


joelewis101/blantyreESBL documentation built on Nov. 1, 2023, 1:43 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close