Chapter 16 Antimicrobial resistance

genome_annotations_amr <- read_tsv("data/gene_annotations_amr.tsv.xz") %>%
    filter(resistance_type=="AMR")

genome_taxonomy_amr <- read_tsv("data/genome_taxonomy_amr.tsv")

16.1 Add genome column to the annotations dataframe

genome_annotations_amr <- genome_annotations_amr %>%
  mutate(genome = sub("@.*", "", gene))

16.2 Link metadata information to the genomes

gene_amr_metadata<-genome_annotations_amr %>%
    left_join(merged_data[c(1, 15, 16,23)], by="genome") %>%
    filter(time_point=="Wild") %>%
    left_join(genome_metadata[c(1, 3,6,7,8)], by = "genome")

16.3 AMR diversity

Amount of AMR genes per species

gene_amr_metadata %>%
  count(species.y, individual) %>%   # count AMR per individual
  group_by(species.y) %>%
  summarise(
    mean_AMR = mean(n),
    sd_AMR = sd(n),
    n_individuals = n()
  )%>%
tt()

species.y	mean_AMR	sd_AMR	n_individuals
Podarcis_liolepis	18627.56	6778.374	9
Podarcis_muralis	17735.22	5227.598	18

amr_indiv <- gene_amr_metadata %>%
  count(species.y, individual)

kruskal.test(n ~ species.y, data = amr_indiv)

    Kruskal-Wallis rank sum test

data:  n by species.y
Kruskal-Wallis chi-squared = 1.2804, df = 1, p-value = 0.2578

gene_amr_metadata %>%
  count(genome, phylum, species, resistance_target, species.y) %>%
  group_by(species.y) %>%
  ggplot(aes(y = n, x = species.y, color=species.y, fill=species.y)) +
      geom_col() +
   scale_color_manual(name="Species",
          breaks=c("Podarcis_muralis","Podarcis_liolepis"),
          labels=c("Podarcis muralis","Podarcis liolepis"),
          values=c('#008080', "#d57d2c"))+
  scale_fill_manual(name="Species",
          breaks=c("Podarcis_muralis","Podarcis_liolepis"),
          labels=c("Podarcis muralis","Podarcis liolepis"),
          values=c('#00808050', "#d57d2c50"))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_blank(),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Nº of AMR genes") 

Type of AMR genes per species

gene_amr_metadata %>%
  count(genome, phylum, species, resistance_target, species.y) %>%
  group_by(species.y) %>%
  ggplot(aes(y = n, x = species.y, color=resistance_target, fill=resistance_target)) +
      geom_col(position = "stack") +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_blank(),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8),
    legend.position = "bottom"
      ) +
  guides(fill = guide_legend(nrow = 3, byrow = TRUE)) +
  ylab("Nº of AMR genes") 

Nº of AMR genes per phylum

gene_amr_metadata %>%
  count(genome, phylum, species, resistance_target, species.y) %>%
  group_by(phylum) %>%
  ggplot(aes(y = n, x = phylum, color=species.y, fill=species.y)) +
      geom_col() +
   scale_color_manual(name="Species",
          breaks=c("Podarcis_muralis","Podarcis_liolepis"),
          labels=c("Podarcis muralis","Podarcis liolepis"),
          values=c('#008080', "#d57d2c"))+
  scale_fill_manual(name="Species",
          breaks=c("Podarcis_muralis","Podarcis_liolepis"),
          labels=c("Podarcis muralis","Podarcis liolepis"),
          values=c('#00808050', "#d57d2c50"))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_text(size=10),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle=45, hjust=1, vjust=1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Nº of AMR genes") +
  xlab("Phylum")

Nº of AMR genes per genome

gene_amr_metadata %>%
  count(genome, phylum, species, resistance_target, species.y) %>%
  group_by(phylum) %>%
  ggplot(aes(y = n, x = genome, color=phylum, fill=phylum)) +
      geom_col()+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_text(size=10),
    axis.title.y = element_text(size=10),
    axis.text.x = element_blank(),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Nº of AMR genes") +
  xlab("Genome")