---
title: 'Case study: AMR analysis in the context of a diphtheria outbreak in Europe'
author: "Luca Freschi"
date: "2025-06-24"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Phylogenies and heatmaps, side by side

We first load the phylogenetic tree using the `treeio` package:

```{r}
library(treeio)
tree_377 = read.tree("t377b_rooted.nwk")
tree_377
```

We plot the tree with `ggtree`:

```{r fig.width=25, fig.height=20}
library(ggtree)
p_tree_377 <- ggtree(tree_377) + 
    geom_tiplab(size=4) +
    geom_treescale(y=20, offset=2)
p_tree_377
```

Then we load the AMR data and we reshape the data frame:
```{r}
# I load the AMR data
m_amr_pred_genes <- read.csv("predicted_amr_data_genes_377_cgMLST.tsv", sep="\t", row.names=1)

# Filtering the data frame
m_amr_pred_genes_filtered <- m_amr_pred_genes[colSums(m_amr_pred_genes > 0) > 0]

# I change the order in which the columns are displayed as well as the labels
right_order <- c('erm.X.', 'tet.W.', 'ant.2....Ia', 'blaOXA.2', 'aadA1',
                   'dfrA1', 'aadA5', 'cmx', 'aph.6..Id', 'aph.3....Ib', 'aph.3...Ia', 'tet.33.','sul1')
m_amr_pred_genes_filtered <- m_amr_pred_genes_filtered[, right_order]

colnames(m_amr_pred_genes_filtered) <- c("erm(X)","tet(W)", "ant(2\")-Ia","blaOXA-2", "aadA1",
      "dfrA1","aadA5","cmx","aph(6)-Id","aph(3\")-Ib", "aph(3\")-Ia","tet(33)","sul1"      
           )

m_amr_pred_genes_filtered[] <- lapply(m_amr_pred_genes_filtered, as.character)

head(m_amr_pred_genes_filtered,2)
```

We plot the tree and the heatmap side by side:
```{r, fig.width=30, fig.height=35}
p_tree_377_heatmap <- gheatmap(p_tree_377, 
         m_amr_pred_genes_filtered, 
         offset=0.2, 
         width=0.2, 
         font.size=4, 
         colnames_angle=-45, 
         hjust=0,
         color="black"
        ) +
    scale_fill_manual(
        breaks=c("0","1"), 
        values=c("thistle", "purple"), 
        labels=c("Absent", "Present"),
        name="Resistance genes"
    ) +
    vexpand(0.2,-1)
p_tree_377_heatmap
```
Finally, we save the output in pdf:
```{r}
ggsave("res_genes_cgMLST_377.pdf", p_tree_377_heatmap,height=30,width=27)
```