figure5_sup14.tree_with_patients_data.Rmd

---
title: "Genome comparison of isolates from hospital and patient"
output:
  html_document:
    df_print: paged
---

### Genome comparison of isolates from hospital and patient

Load generic libraries
```{r message=FALSE, warning=FALSE}
source('configuration.r')
```

Load plot specific libraries
```{r message=FALSE}
library(foreach)
library(ggtree)
library(phytools)
library(magrittr)
library(readr)
```

Draw a species tree using 16S 
```{r}
tr <- midpoint.root(read.tree('../tables/trees/16s.tree'))
tr$tip.label <-  c("E. anophelis", "S. aureus", "S. epidermidis", "A. baumannii")

ggtree(tr, lwd=2, layout='circular') + 
  geom_tippoint(size=6, shape=18)  + 
  geom_treescale(linesize=2, offset=0.1, width=0.02, y=4, fontsize = 0)  #+
##geom_tiplab()

ggsave('../plots/fig5_16S.svg', width = 5, height = 5)
```

Function to plot tree
```{r}
strains.all <- read.table('../output_tables/strain_cluster_summary.tsv', head=TRUE, stringsAsFactors = FALSE)
colors <- pal_npg('nrc')(9)
get_tree_data <- function(species, diff_threshold=0.0001){
    strain.dat <- filter(strains.all, Species==species)
    dist.dat <- read.table(paste0('../tables/',species,'_opera_ms_mummer_dist.dat')) %>% 
      select(V1, V2, coverage=V3, ngsnp=V9, identity=V10)
    ## add strain information to dist.dat
    tmp <- 
      select(strain.dat, Nanopore_ID, strain) %>% 
      merge(dist.dat, ., by.x=1, by.y=1, all.x=TRUE) %>% ## annotate left
      merge(strain.dat, by.x=2, by.y=1, all.x=TRUE) %>%  ## merge right
      select(V1,V2,strain.x,strain.y,coverage,ngsnp,identity) %>% 
      filter(coverage>80, identity>(1-diff_threshold) * 100) 
    ## aggregated overlap data
    overlaps <- bind_rows(
      filter(tmp, !is.na(strain.x) & is.na(strain.y)) %>% 
        select(strain=strain.x, V1, V2, coverage,ngsnp,identity),
      filter(tmp, !is.na(strain.y) & is.na(strain.x)) %>% 
        select(strain=strain.y, V1, V2, coverage,ngsnp,identity)
    )
    ## get antibiotics profile
    distinct(strain.dat, strain, Antibiotics) %>% 
      mutate(detect=Antibiotics) %>% 
      pivot_wider(id_cols = strain, names_from = Antibiotics, values_from = detect) %>% 
      select(strain, AMP, CHLOR, KAN, STREP, TET) %>% 
      column_to_rownames("strain")  -> profile.antibiotics
    ## get tree
    tree.strains <- read.tree(paste0('../tables/trees/', species, '.parsnp.tree'))
    tree.strains$tip.label %<>% str_replace_all(c(".trimmed.fasta"="", "nanopore.cons.cluster_"="s", ".ref"=""))
    ## get size of nodes
    anno <- count(strain.dat, strain)
    list(tree=tree.strains, overlaps=overlaps, antibiotics=profile.antibiotics, anno=anno)
}
```

#### *Staphylococcus aureus*

Overlapping analysis
```{r}
dat <- get_tree_data('Staphylococcus_aureus')

## patient overlapping profile
patient <- distinct(dat$overlaps, strain, V2) %>% 
  filter(str_detect(V2, 'ADD')) %>% 
  count(strain, name = '99.99%') %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain, `99.99%`) 

## persistance after 1.5 year
persistance <- distinct(dat$overlaps, strain, V2) %>% 
  filter(str_detect(V2, 'barcode')) %>%  ## 3rd time point identified by keyword 'barcode' from qcat demuxer
  count(strain, name = '99.99%') %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain, `99.99%`) 
```

Full overlapping table
```{r}
arrange(dat$overlaps %>% filter(str_detect(V2, 'ADD')), desc(identity))
```

Test for significance of MDR
```{r}
test.dat <- rownames_to_column(dat$antibiotics, 'strain') %>% 
  mutate(MDR=rowSums(!is.na(dat$antibiotics)) > 2) %>% 
  select(strain, MDR) %>% 
  merge(patient,all=TRUE, by=1) %>% 
  merge(persistance, all=TRUE, by=1)

p1 <- binom.test(filter(test.dat, MDR) %>% pull(`99.99%.x`) %>% sum(na.rm = T), 
           sum(test.dat$`99.99%.x`, na.rm = T), 
           sum(test.dat$MDR)/nrow(test.dat)
           )$p.value
p2 <- binom.test(filter(test.dat, MDR, !is.na(`99.99%.x`)) %>% pull(`99.99%.y`) %>% sum(na.rm = T), 
           sum(test.dat$`99.99%.y`, na.rm = T), 
           filter(test.dat, MDR, !is.na(`99.99%.x`)) %>% nrow / nrow(test.dat) ## MDR + overlap with patient
           )$p.value

print(paste0('Test for 99.99% overlapping with patient: ', p1))
print(paste0('Test for 99.99% overlapping with time point 3: ', p2))
```


Plot the tree
```{r fig.height=6, fig.width=6, message=FALSE}
d1 <- rownames_to_column(dat$antibiotics) %>% melt(id.vars='rowname') %>% mutate(`Antibiotic resistance`=value)
d2 <- melt(patient) %>% mutate(N=value) 
d3 <- melt(persistance) %>% mutate(N=value)

p <- ggtree(dat$tree, lwd=1, color=colors[3]) %<+% select(d2, strain, N) %<+% dat$anno + 
  geom_tippoint(aes(size=n,color=is.na(N))) +
  geom_tiplab(hjust=1, fontface='bold', offset=0.004, align = TRUE, color=colors[3]) +
  geom_treescale(x=0.002, y=15, color = colors[3], linesize=1, width=0.002) + 
  scale_color_manual(values=c('black', 'grey'),name="Patient overlap", labels=c("Yes", "No"))

p1 <- p + 
  geom_facet(panel='Antibiotic\nresistance profile', data=d1, geom=geom_tile, 
               mapping=aes(x = as.numeric(variable),  fill = `Antibiotic resistance`), color='white', lwd=1) + 
  #geom_facet(panel='Htun et al\nDerivative', data=d2, geom=geom_point, 
  #           mapping=aes(x = as.numeric(variable), size=N), color=colors[3]) +
  # geom_facet(panel='Timepoint 3\nDerivative', data=d3, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[3]) +
  scale_fill_manual(values=colors[c(1,5,2,6,4)], na.translate=FALSE) +
  scale_size_continuous(limits = c(1,150), trans='sqrt', breaks =c(1,10,20,50,100,150) ) +
  theme(legend.position = 'right', strip.background =element_rect(fill=NA, color=NA), 
        strip.text = element_text(size=10, face='bold'))

p1
legend <- cowplot::get_legend(p1)
p1 <- p1 + theme(legend.position = 'none')
        
tree1 <- facet_widths(p1, c(2,1))
tree1
```

#### *Elizabethkingia anophelis*

Overlapping analysis
```{r}
dat <- get_tree_data('Elizabethkingia_anophelis')

## patient overlapping profile
patient <- 
  filter(dat$overlaps, str_detect(V2, 'ADD')) %>% 
  mutate(anno=ifelse(!str_detect(V2, 'WEB'), 'Teo et al', 'This study')) %>%
  group_by(strain, V2, anno) %>% summarise(identity=max(identity)) %>% ungroup() %>% 
  group_by(strain, anno) %>% 
  summarise('99.99%'=n(), '99.999%'=sum(identity>99.999)) %>% ungroup() %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain:`99.999%`)

## persistance after 1.5 year
persistance <- filter(dat$overlaps, str_detect(V2, 'barcode')) %>% 
  group_by(strain, V2) %>% summarise(identity=max(identity)) %>% ungroup() %>% 
  group_by(strain) %>% 
  summarise('99.99%'=n(), '99.999%'=sum(identity>99.999)) %>% ungroup() %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain:`99.999%`)
```

Full overlapping table
```{r}
arrange(dat$overlaps %>% filter(str_detect(V2, 'ADD')), desc(identity)) 
```


Plot the tree
```{r fig.height=2.5, fig.width=6, message=FALSE}
d1 <- rownames_to_column(dat$antibiotics) %>% melt(id.vars='rowname') %>% mutate(`Antibiotic resistance`=value)
d2 <- filter(patient, anno=='Teo et al') %>% melt() %>% 
  filter(value!=0) %>% 
  unite('variable', c(anno, variable)) %>% 
  mutate(N=value, variable=as.factor(variable))
d3 <- filter(patient, anno!='Teo et al') %>% melt() %>% 
  filter(value!=0) %>% 
  unite('variable', c(anno, variable)) %>% 
  mutate(N=value, variable=as.factor(variable))
d4 <- melt(persistance) %>% mutate(N=value)


p <- ggtree(midpoint.root(dat$tree), color=colors[1], lwd=1) %<+% select(d3, strain, N) %<+% dat$anno +
  geom_tippoint(aes(size=n, color=is.na(N))) +
  geom_tiplab(color=colors[1],hjust=1, fontface='bold', offset=0.001, align = TRUE) +
  geom_treescale(x=0, y=3.5, color = colors[1], linesize=1, width=0.001)  + 
  scale_color_manual(values=c('black', 'grey'),name="Patient overlap", labels=c("Yes", "No"))

p1 <- p + 
  #geom_facet(panel='Teo et al\nDerivative', data=d2, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[1]) +
  # geom_facet(panel='This study\nDerivative', data=subset(d3, variable=='This study_99.99%'), geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[1]) +
  # geom_facet(panel='This study\nTransfer', data=subset(d3, variable=='This study_99.999%'), geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[1])  + 
  # geom_facet(panel='Timepoint 3\nDerivative', data=subset(d4, variable=='99.99%'), geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[1]) +
  # geom_facet(panel='Timepoint 3\nTransfer', data=subset(d4, variable=='99.999%'), geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[1]) +
  geom_facet(panel='Antibiotic\nresistance profile', data=d1, geom=geom_tile, 
               mapping=aes(x = as.numeric(variable),  fill = `Antibiotic resistance`), color='white', lwd=1) + 
  scale_fill_manual(values=colors[c(1,5,2,6,4)],na.value='white',guide=F) +
  scale_size_continuous(limits = c(1,150), trans='sqrt', breaks =c(1,10,20,50,100,150) ) +
  theme(legend.position = 'none', strip.background =element_rect(fill=NA, color=NA), 
        strip.text = element_text(size=10, face='bold')) +
  ylim(0,5) 

tree2 <- facet_widths(p1, c(2,1))
tree2
```



#### *Acinetobacter baumannii*

Overlapping analysis
```{r}
dat <- get_tree_data('Acinetobacter_baumannii')

## patient overlapping profile
patient <- 
  filter(dat$overlaps, str_detect(V2, 'ADD')) %>% 
  mutate(anno=ifelse(!str_detect(V2, 'MBB'), 'Nasser et al', 'This study')) %>% 
  group_by(strain, V2, anno) %>% summarise(identity=max(identity)) %>% ungroup() %>% 
  group_by(strain, anno) %>% 
  summarise('99.99%'=n(), '99.999%'=sum(identity>99.999)) %>% ungroup() %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain:`99.999%`)
```

Full overlapping table
```{r}
arrange(dat$overlaps %>% filter(str_detect(V2, 'ADD')), desc(identity)) 
```

Test for significance of MDR
```{r}
test.dat <- 
  rownames_to_column(dat$antibiotics, 'strain') %>% 
  mutate(MDR=rowSums(!is.na(dat$antibiotics)) > 2) %>% 
  select(strain, MDR) %>% 
  merge(filter(patient, anno=='This study') %>% select(strain, `99.99%`), all=TRUE, by=1) %>% 
  merge(filter(patient, anno!='This study') %>% select(strain, `99.99%`), all=TRUE, by=1)

p1 <- binom.test(filter(test.dat, MDR) %>% select(`99.99%.x`, `99.99%.y`) %>% colSums(na.rm=T) %>% sum(),
           sum(c(test.dat$`99.99%.x`, test.dat$`99.99%.y`), na.rm = T), 
           sum(test.dat$MDR)/nrow(test.dat)
           )$p.value

print(paste0('Test for 99.99% overlapping with patients: ', p1))
```

Plot the tree

```{r fig.height=9, fig.width=8, message=FALSE}
d1 <- rownames_to_column(dat$antibiotics) %>% melt(id.vars='rowname') %>% mutate(`Antibiotic resistance`=value)
d2 <- filter(patient, anno=='Nasser et al') %>% melt() %>% 
  filter(value!=0) %>% 
  unite('variable', c(anno, variable)) %>% 
  mutate(N=value, variable=as.factor(variable))
d3 <- filter(patient, anno!='Nasser et al') %>% melt() %>% 
  filter(value!=0) %>% 
  unite('variable', c(anno, variable)) %>% 
  mutate(N=value, variable=as.factor(variable))

p <- ggtree(dat$tree, color=colors[4], lwd=1) %<+% select(d3, strain, N) %<+% dat$anno +
  geom_tippoint(aes(size=n, color=is.na(N))) +
  geom_tiplab(color=colors[4],hjust=1, fontface='bold', offset=0.003, align = TRUE) +
  geom_treescale(x=0, y=18, color = colors[4], linesize=1)  + 
  scale_color_manual(values=c('black', 'grey'),name="Patient overlap", labels=c("Yes", "No"))

p1 <- p + 
  # geom_facet(panel='Nasser et al\nDerivative', data=d2, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[4]) +  
  # geom_facet(panel='This study\nDerivative', data=d3, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[4]) +
  geom_facet(panel='Antibiotic\nresistance profile', data=d1, geom=geom_tile, 
               mapping=aes(x = as.numeric(variable),  fill = `Antibiotic resistance`), color='white', lwd=1) + 
  scale_fill_manual(values=colors[c(1,5,2,6,4)],na.value='white',guide=F) +
  scale_size_continuous(limits = c(1,150), trans='sqrt', breaks =c(1,10,20,50,100,150) ) +
  theme(legend.position = 'none', strip.background =element_rect(fill=NA, color=NA), 
        strip.text = element_text(size=10, face='bold'))
        
tree3 <- facet_widths(p1, c(2,1))
tree3
```

### *Staphylococcus epidermidis*

Overlapping analysis
```{r}
dat <- get_tree_data('Staphylococcus_epidermidis')

se.meta <- read_tsv('../metadata/mrse_metadata.txt')

## patient overlapping profile
patient <- 
  filter(dat$overlaps, str_detect(V2, 'ADD')) %>% 
  filter(!str_detect(V2, 'SAMEA3712712|SAMEA3712767')) %>% ## removed because they were mutant
  filter(V2 %in% se.meta$pipeline) %>% 
  group_by(strain, V2) %>% summarise(identity=max(identity)) %>% ungroup() %>% 
  group_by(strain) %>% 
  summarise('99.99%'=n(), '99.999%'=sum(identity>99.999)) %>% ungroup() %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain:`99.999%`)

## persistance after 1.5 year
persistance <- filter(dat$overlaps, str_detect(V2, 'barcode')) %>% 
  group_by(strain, V2) %>% summarise(identity=max(identity)) %>% ungroup() %>% 
  group_by(strain) %>% 
  summarise('99.99%'=n(), '99.999%'=sum(identity>99.999)) %>% ungroup() %>% 
  merge(dat$antibiotics, all.y = TRUE, by.x=1, by.y=0) %>% 
  select(strain:`99.999%`)
```

Full overlapping table
```{r}
arrange(dat$overlaps %>% filter(str_detect(V2, 'ADD')), desc(identity))
```

Test for significance of MDR

```{r}
test.dat <- rownames_to_column(dat$antibiotics, 'strain') %>% 
  mutate(MDR=rowSums(!is.na(dat$antibiotics)) > 2) %>% 
  select(strain, MDR) %>% 
  merge(patient,all=TRUE, by=1) %>% 
  merge(persistance, all=TRUE, by=1)

p1 <- binom.test(filter(test.dat, MDR) %>% pull(`99.99%.x`) %>% sum(na.rm = T), 
           sum(test.dat$`99.99%.x`, na.rm = T), 
           sum(test.dat$MDR)/nrow(test.dat)
           )$p.value
p2 <- binom.test(filter(test.dat, MDR, !is.na(`99.99%.x`)) %>% pull(`99.99%.y`) %>% sum(na.rm = T), 
           sum(test.dat$`99.99%.y`, na.rm = T), 
           filter(test.dat, MDR, !is.na(`99.99%.x`)) %>% nrow / nrow(test.dat) ## MDR + overlap with patient
           )$p.value

print(paste0('Test for 99.99% overlapping with patient: ', p1))
print(paste0('Test for 99.99% overlapping with time point 3: ', p2))
```

```{r fig.height=10, fig.width=7, message=FALSE}
d1 <- rownames_to_column(dat$antibiotics) %>% melt(id.vars='rowname') %>% mutate(`Antibiotic resistance`=value)
d2 <- melt(patient) %>% 
  filter(value!=0) %>% 
  mutate(N=value, variable=as.factor(variable))
d3 <- melt(persistance) %>% filter(value!=0) %>% mutate(N=value)

p <- ggtree(dat$tree, color=colors[2], lwd=1) %<+% select(d2, strain, N) %<+% dat$anno + 
  geom_tippoint(aes(size=n, color=is.na(N))) +
  geom_tiplab(color=colors[2],hjust=1, fontface='bold', offset=0.002, align = TRUE) +
  geom_treescale(x=0, y=39, color = colors[2], linesize=1) + 
  scale_color_manual(values=c('black', 'grey'),name="Patient overlap", labels=c("Yes", "No"))
  
p1 <- p + 
  # geom_facet(panel='Lee et al\nDerivative', data=d2, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[2]) +
  # geom_facet(panel='Timepoint 3\nDerivative', data=d3, geom=geom_point, 
  #            mapping=aes(x = as.numeric(variable), size=N), color=colors[2]) +
  geom_facet(panel='Antibiotic\nresistance profile', data=d1, geom=geom_tile, 
               mapping=aes(x = as.numeric(variable),  fill = `Antibiotic resistance`), color='white', lwd=1) + 
  scale_fill_manual(values=colors[c(1,5,2,6,4)],na.value='white',guide=F) +
  scale_size_continuous(limits = c(1,150), trans='sqrt', breaks =c(1,10,20,50,100,150) ) +
  theme(legend.position = 'none', strip.background =element_rect(fill=NA, color=NA), 
        strip.text = element_text(size=10, face='bold'))
        
tree4 <- facet_widths(p1, c(2,1))
tree4
```


### Put four trees at the same scale
```{r fig.height=12, fig.width=13}
p1 <- cowplot::plot_grid(tree2, tree1, tree3, ncol=1, rel_heights = c(0.35,0.9,1))
cowplot::plot_grid(p1, tree4, legend, nrow=1, rel_widths = c(1,1,0.2))
ggsave('../plots/fig5_phylogenetic_trees.svg', height=12, width=13)
```


<!-- ### Overview tree -->
<!-- ```{r fig.height=15, fig.width=10} -->
<!-- full_mat <- read.table('../tables/genome_mash_dist.dat')/4 -->
<!-- idx <-  -->
<!--   str_detect(row.names(full_mat), 'Staphylococcus_aureus') | -->
<!--   str_detect(row.names(full_mat), 'Acinetobacter_baumannii') |  -->
<!--   str_detect(row.names(full_mat), 'Staphylococcus_epidermidis') | -->
<!--   str_detect(row.names(full_mat), 'Elizabethkingia_anophelis')  -->
<!-- full_mat <- full_mat[idx, idx] -->

<!-- meta <- read.table('../output_tables/merged_assembly_metadata.tsv', head=TRUE, sep='\t') %>%  -->
<!--   select(Genome_ID, Species, Nanopore_ID) %>% unite(id, Species, Nanopore_ID, remove = FALSE) -->
<!-- ## remove duplicates -->
<!-- idx <- rownames(full_mat) %in% meta$id -->
<!-- full_mat <- full_mat[idx, idx] -->
<!-- rownames(full_mat) <- colnames(full_mat) <- meta$Genome_ID[match(rownames(full_mat), meta$id)] -->

<!-- aux <- function(species){ ## function with side effect to fill in the mummer distance -->
<!--   dist.ab <- get_tree_data(species, get.dist = TRUE) -->
<!--   sel <- str_subset(rownames(dist.ab) ,'s_') -->
<!--   full_mat[sel, sel] <- dist.ab[sel, sel] -->
<!--   return(sel) -->
<!-- } -->
<!-- ## fill mash matrix with mummer distance -->
<!-- se <- aux('Staphylococcus_epidermidis') -->
<!-- sa <- aux('Staphylococcus_aureus') -->
<!-- ab <- aux('Acinetobacter_baumannii') -->
<!-- ea <- aux('Elizabethkingia_anophelis') -->

<!-- clustering <- hclust(as.dist(full_mat), method='single') -->
<!-- #anno <- meta[match(rownames(full_mat), meta$genomeID), ] -->
<!-- tree <- groupOTU(as.phylo(clustering), list(ea=ea, se=se, sa=sa, ab=ab)) -->

<!-- ## helper to identify internal clade for shading -->
<!-- # ggtree(tree, aes(color=group) , lwd=1.5) +  -->
<!-- #   geom_text2(aes(subset=!isTip, label=node))  -->

<!-- p <-  -->
<!--   ggtree(tree, aes(color=group) , lwd=1.5) +  -->
<!--   scale_color_manual(values=c("black", colors[c(4,1,3,2)])) + -->
<!--   geom_hilight(node=538, alpha=0.2, fill = colors[2]) + -->
<!--   geom_hilight(node=535, alpha=0.2, fill = colors[4]) + -->
<!--   geom_hilight(node=533, alpha=0.2, fill = colors[1]) + -->
<!--   geom_hilight(node=537, alpha=0.2, fill = colors[3]) + -->
<!--   theme(legend.position = 'left') -->
<!--   ##flip(p, 533, 534) %>% flip(535,536) %>% flip(539, 540) -->
<!-- p -->
<!-- ggsave("../plots/fig5_cladrogram.svg", height=15, width=10) -->
<!-- ``` -->


### Session informaton
```{r}
sessionInfo()
```