Section 3 AB-effects (resequenced)

3.1 Metadata

Sample	Genome	Group	Species	Treatment
DR4	KB1	1687	Enterococcus_faecalis	Ciprofloxacin
DR1	KB1	1681	Enterococcus_faecalis	Water
DR8	KB1	1691	Enterococcus_faecalis	Tetracyclin
DR11	KB1	1696	Enterococcus_faecalis	Vancomycin
DR3	I46	1681	Clostridium_innocuum	Water
DR6	I46	1687	Clostridium_innocuum	Ciprofloxacin
DR10	I46	1691	Clostridium_innocuum	Tetracyclin
DR12	I46	1696	Clostridium_innocuum	Vancomycin
DR7	YL58	1687	Blautia_coccoides	Ciprofloxacin
DR13	YL58	1696	Blautia_coccoides	Vancomycin

3.2 Process Lofeq vcf files

source("utils.R")
orf_shapes = c(coding = 15, `non-coding` = 3)

vcfToDataframe <- function(vcf.files, contig_mapping = contig_mapping, gff.df = gff.df) {
    require(vcfR)
    res <- list()
    for (file in vcf.files) {
        # message(file)
        library(data.table)
        vcf.content <- vcfR::read.vcfR(file, verbose = FALSE)
        vcf.fix <- as.data.frame(vcf.content@fix)  # contains chr, position and substitution informations
        vcf.info <- vcfR::INFO2df(vcf.content)  # get INFO field, contains DP, AF informations
        # no annotation since they are not yet processed with snpEff vcf.ann <-
        # data.frame(do.call('rbind', strsplit(as.character(vcf.info$ANN),'|',fixed = TRUE))) #
        # split ANN field, contains information if files are processed by snpEff there are variants
        if (nrow(vcf.fix) > 0) {
            dat <- as.data.frame(cbind(vcf.fix[, c(1, 2, 4, 5, 6)], vcf.info[, c(1, 2)]))
            dat$majorAF <- sapply(dat$AF, minorAfToMajorAf)  # transforms e.g. AF of 0.1 to 0.9, 0.9 stays 0.9 and 0.5 stays 0.5
            dat$genome <- contig_mapping[match(dat$CHROM, contig_mapping$contig), ]$genome  # map chr information to genome name e.g. NHMU01000001.1 -> i48
            dat$genome_hr <- translateGenomeIdToFullName(tolower(dat$genome))
            dat$mouse.id <- translateSampletoMouse(tools::file_path_sans_ext(basename(file)))
            dat$mouse.group <- translateMouseIdToTreatmentGroup(translateSampletoMouse(tools::file_path_sans_ext(basename(file))))
            dat$dp <- as.numeric(as.matrix(vcf.info$DP))
            
            dt.gff <- data.table(start = gff.df$start, end = gff.df$end, chr = as.character(as.matrix(gff.df$chr)), 
                feature = gff.df$product)
            colnames(dat)[1:2] <- c("chr", "start")
            dat$start <- as.integer(as.matrix(dat$start))
            dat$chr <- as.character(as.matrix(dat$chr))
            dat$end <- dat$start
            dat2 <- as.data.table(dat)
            setkey(dt.gff, chr, start, end)
            annotated <- foverlaps(dat2, dt.gff, type = "within", mult = "first")
            res[[tools::file_path_sans_ext(basename(file))]] <- annotated  # add vcf df to list
        } else {
            message("Skipping")
        }
    }
    df <- as.data.frame(do.call(rbind, res))  # merge list to df
    return(df)
}

gff.files <- Sys.glob("data/references/gff/*.gff")
gff.df <- NULL
for (gff.file in gff.files) {
    message(gff.file)
    gff <- rtracklayer::readGFF(gff.file)
    # subset since different columns are present on gff files
    relevant <- data.frame(start = gff$start, end = gff$end, type = as.character(as.matrix(gff$type)), 
        gene = as.character(as.matrix(gff$gene)), product = as.character(as.matrix(gff$product)), 
        chr = as.character(as.matrix(gff$seqid)))
    relevant$genome <- substr(basename(gff.file), 1, nchar(basename(gff.file)) - 4)
    gff.df <- rbind(gff.df, relevant)
}

## data/references/gff/I46.gff

## data/references/gff/KB1.gff

## data/references/gff/YL58.gff

contig_mapping <- read.csv2("data/contig_mapping_new_ref.csv", sep = ";", header = T, stringsAsFactors = F)  # this file contains contig names of the 12 OligoMM genomes
vcf.files <- Sys.glob("out_reseq/all_vcf/*.vcf")
vcf.samples <- suppressWarnings(vcfToDataframe(vcf.files, contig_mapping, gff.df = gff.df))

vcf.samples$feature <- as.character(as.matrix(vcf.samples$feature))
vcf.samples[which(is.na(vcf.samples$feature)), ]$feature <- "outside ORFs"
vcf.samples[which(vcf.samples$feature == "hypothetical protein"), ]$feature <- "hypothetical proteins"

vcf.samples$start <- NULL
vcf.samples$end <- NULL
vcf.samples$i.end <- NULL
colnames(vcf.samples)[3] <- "POS"
vcf.samples$ref_size <- nchar(as.character(as.matrix(vcf.samples$REF)))
vcf.samples$alt_size <- nchar(as.character(as.matrix(vcf.samples$ALT)))
vcf.samples$alteration <- paste(as.character(vcf.samples$REF), "->", as.character(vcf.samples$ALT))
vcf.samples$alteration_type <- "SNP"
vcf.samples[which(vcf.samples$ref_size < vcf.samples$alt_size), ]$alteration_type <- "insertion"
vcf.samples[which(vcf.samples$ref_size > vcf.samples$alt_size), ]$alteration_type <- "deletion"
saveRDS(vcf.samples, file = "data/rds/reseq.rds")

Save as table

dat <- readRDS("data/rds/reseq.rds")
write.table(vcf_samples, file = "results/tables/omm_antibiotic_resequencing_variants_long.csv", 
    sep = ";", col.names = T, row.names = F, quote = F)

3.3 Plot AF frequency

dat <- readRDS("data/rds/reseq.rds")
p1 <- ggplot(dat, aes(AF, fill = genome)) + geom_histogram()
p1 <- p1 + facet_grid(mouse.id + mouse.group ~ genome + genome_hr)
p1 <- p1 + theme_minimal() + xlab("AF") + ylab("occurence")
plotly::ggplotly(p1)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 3.1: AF of resequenced strains

p <- ggplot(dat, aes(majorAF, fill = genome)) + geom_histogram()
p <- p + facet_grid(mouse.id + mouse.group ~ genome + genome_hr)
p <- p + theme_minimal() + xlab("AF") + ylab("occurence")
plotly::ggplotly(p)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 3.2: major AF of resequenced strains

pairedAFplot <- function(dat = dat, case = "Water", control = "Ciprofloxacin") {
    require(reshape2)
    require(ggplot2)
    
    dat.pair <- dat[which(dat$mouse.group == case | dat$mouse.group == control), ]
    
    # dat.pair <- dat.pair[which(dat.pair$genome == genome),]
    
    # dat.pair <- dat.pair[which(!is.na(dat.pair$feature)),]
    data.wide <- dcast(dat.pair, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")
    data.wide$vartype <- ifelse(as.character(data.wide$feature) != "outside ORFs", "coding", 
        "non-coding")
    data.wide <- data.wide[which(data.wide$vartype != 0), ]
    
    colnames(data.wide)[6] <- "sample_1"
    colnames(data.wide)[7] <- "sample_2"
    data.wide[is.na(data.wide)] <- 0
    
    p <- ggplot(data.wide, aes(x = sample_2, y = sample_1, label = feature, color = genome, 
        shape = vartype))
    p <- p + geom_point(fill = NA) + theme_classic()  # + scale_color_manual(values = omm_colors)
    p <- p + theme(aspect.ratio = 1)
    p <- p + geom_hline(yintercept = 0.5, linetype = 2, color = "grey50")
    p <- p + geom_vline(xintercept = 0.5, linetype = 2, color = "grey50")
    p <- p + scale_shape_manual(values = orf_shapes)
    p <- p + theme(aspect.ratio = 1)
    p <- p + theme(strip.background = element_blank())
    p <- p + geom_abline(intercept = 0, slope = 1)
    p <- p + xlab(case) + ylab(control)
    p <- p + theme(panel.border = element_rect(colour = "black", fill = NA, size = 1))
    return(p)
}

suppressWarnings(ggplotly(pairedAFplot(dat, case = "Water", control = "Ciprofloxacin")))

Figure 3.3: Scatter plot of AF of two samples. Only variants with functional annotation are shown

suppressWarnings(ggplotly(pairedAFplot(dat, case = "Water", control = "Tetracyclin")))

Figure 3.4: Scatter plot of AF of two samples. Only variants with functional annotation are shown

suppressWarnings(ggplotly(pairedAFplot(dat, case = "Water", control = "Vancomycin")))

Figure 3.5: Scatter plot of AF of two samples. Only variants with functional annotation are shown

dat <- readRDS("data/rds/reseq.rds")
dat$ref_size <- nchar(as.character(as.matrix(dat$REF)))
dat$alt_size <- nchar(as.character(as.matrix(dat$ALT)))
dat$alteration <- paste(as.character(dat$REF), "->", as.character(dat$ALT))

# no grid version for html
geneCentricPlot2 <- function(dat = dat, genome = "KB1", type = "SNP") {
    require(reshape2)
    require(ggplot2)
    dat.subset <- dat[which(dat$genome == genome), ]
    if (type == "SNP") {
        dat.subset <- dat.subset[which(dat.subset$ref_size == 1 & dat.subset$alt_size == 1), 
            ]
    } else if (type == "insertion") {
        dat.subset <- dat.subset[which(dat.subset$ref_size == 1 & dat.subset$alt_size > 1), 
            ]
    } else {
        # deletion
        dat.subset <- dat.subset[which(dat.subset$ref_size > 1 & dat.subset$alt_size == 1), 
            ]
    }
    dat.subset$ID <- paste0(dat.subset$POS, "|", dat.subset$REF, "->", dat.subset$ALT)
    p <- ggplot(dat.subset, aes(x = feature, y = AF, color = mouse.group, label = alteration))
    p <- p + geom_point(size = 2, fill = NA)
    p <- p + coord_flip() + theme_classic()
    p <- p + ylim(c(0, 1))
    p <- p + theme(strip.background = element_blank()) + theme(strip.text.y = element_text(angle = 0), 
        strip.placement = "outside")
    p <- p + xlab("") + ylab("allele frequency")
    return(p)
}

3.4 KB1

3.4.1 SNP

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "KB1", type = "SNP")))

Figure 3.6: AF by function of one genome

3.4.2 insertion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "KB1", type = "insertion")))

Figure 3.7: AF by function of one genome

3.4.3 deletion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "KB1", type = "deletion")))

Figure 3.8: AF by function of one genome

3.5 I46

3.5.1 SNP

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "I46", type = "SNP")))

Figure 3.9: AF by function of one genome

3.5.2 insertion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "I46", type = "insertion")))

Figure 3.10: AF by function of one genome

3.5.3 deletion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "I46", type = "deletion")))

Figure 3.11: AF by function of one genome

3.6 YL58

3.6.1 SNP

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "YL58", type = "SNP")))

Figure 3.12: AF by function of one genome

3.6.2 insertion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "YL58", type = "insertion")))

Figure 3.13: AF by function of one genome

3.6.3 deletion

suppressWarnings(ggplotly(geneCentricPlot2(dat, genome = "YL58", type = "deletion")))

Figure 3.14: AF by function of one genome

3.7 mutation bias?

3.7.1 SNP

dat <- readRDS("data/rds/reseq.rds")
dat$ref_size <- nchar(as.character(as.matrix(dat$REF)))
dat$alt_size <- nchar(as.character(as.matrix(dat$ALT)))
dat <- dat[which(dat$ref_size == 1 & dat$alt_size == 1), ]
dat$alteration <- paste(as.character(dat$REF), "->", as.character(dat$ALT))
dat$dummy <- 1
dat.agg <- aggregate(dummy ~ genome + mouse.group + alteration, dat, sum)

p <- ggplot(dat.agg, aes(x = mouse.group, y = dummy, label = alteration, color = genome))
p <- p + geom_text() + facet_grid(~genome) + ylab("occurrence")
p <- p + stat_summary(fun.y = mean, geom = "line", aes(group = alteration), alpha = 0.1)

## Warning: `fun.y` is deprecated. Use `fun` instead.

p <- p + theme_classic()  # + scale_color_manual(values = omm_colors)
ggplotly(p)

3.7.2 insertions

dat <- readRDS("data/rds/reseq.rds")
dat$ref_size <- nchar(as.character(as.matrix(dat$REF)))
dat$alt_size <- nchar(as.character(as.matrix(dat$ALT)))
dat <- dat[which(dat$alt_size > dat$ref_size), ]
dat$alteration <- paste(as.character(dat$REF), "->", as.character(dat$ALT))
dat$dummy <- 1
dat.agg <- aggregate(dummy ~ genome + mouse.group + alteration, dat, sum)

p <- ggplot(dat.agg, aes(x = mouse.group, y = dummy, label = alteration, color = genome))
p <- p + geom_text() + facet_grid(~genome) + ylab("occurrence")
p <- p + stat_summary(fun.y = mean, geom = "line", aes(group = alteration), alpha = 0.1)

## Warning: `fun.y` is deprecated. Use `fun` instead.

p <- p + theme_classic()  # + scale_color_manual(values = omm_colors)
ggplotly(p)

3.7.3 deletions

dat <- readRDS("data/rds/reseq.rds")
dat$ref_size <- nchar(as.character(as.matrix(dat$REF)))
dat$alt_size <- nchar(as.character(as.matrix(dat$ALT)))
dat <- dat[which(dat$alt_size < dat$ref_size), ]
dat$alteration <- paste(as.character(dat$REF), "->", as.character(dat$ALT))
dat$dummy <- 1
dat.agg <- aggregate(dummy ~ genome + mouse.group + alteration, dat, sum)

p <- ggplot(dat.agg, aes(x = mouse.group, y = dummy, label = alteration, color = genome))
p <- p + geom_text() + facet_grid(~genome) + ylab("occurrence")
p <- p + stat_summary(fun.y = mean, geom = "line", aes(group = alteration), alpha = 0.1)

## Warning: `fun.y` is deprecated. Use `fun` instead.

p <- p + theme_classic()  # + scale_color_manual(values = omm_colors)
ggplotly(p)

3.8 heatmap

3.8.1 KB1

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "KB1"), ]
rownames(data.wide) <- paste0(data.wide$POS, "-", data.wide$REF, "-", data.wide$ALT)
heat <- data.matrix(data.wide[, 6:9])
library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, column_names_gp = gpar(fontsize = 9), row_names_gp = gpar(fontsize = 3), 
    show_row_dend = F, show_row_names = F, show_column_dend = F)

3.8.2 I46

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "I46"), ]
rownames(data.wide) <- paste0(data.wide$POS, "-", data.wide$REF, "-", data.wide$ALT)
heat <- data.matrix(data.wide[, 6:9])
library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, column_names_gp = gpar(fontsize = 9), row_names_gp = gpar(fontsize = 9), 
    show_row_names = F, show_row_dend = F, show_column_dend = F)

3.8.3 YL58

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "YL58"), ]
rownames(data.wide) <- paste0(data.wide$POS, "-", data.wide$REF, "-", data.wide$ALT)
heat <- data.matrix(data.wide[, 6:9])
library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, column_names_gp = gpar(fontsize = 9), row_names_gp = gpar(fontsize = 9), 
    show_row_dend = F, show_row_names = F, show_column_dend = F)

3.9 Filter: present in all

3.9.1 KB1 (Enterococcusfaecalis)

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "KB1"), ]
rownames(data.wide) <- paste0(data.wide$POS, " (", data.wide$REF, " -> ", data.wide$ALT, ") ", 
    data.wide$feature)

heat <- data.matrix(data.wide[, 6:9])
heat <- heat[which(rowSums(!is.na(heat)) > 2), ]  # filter uncommon

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, column_names_gp = gpar(fontsize = 9), row_names_gp = gpar(fontsize = 9), 
    show_row_dend = F, row_names_max_width = unit(10, "cm"), show_row_names = T, show_column_dend = F)

3.9.2 YL58 (Blautia_coccoides)

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "YL58"), ]
rownames(data.wide) <- paste0(data.wide$POS, " (", data.wide$REF, " -> ", data.wide$ALT, ") ", 
    data.wide$feature)

heat <- data.matrix(data.wide[, 6:9])
heat <- heat[which(rowSums(!is.na(heat)) > 2), ]  # filter uncommon

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, row_names_max_width = unit(10, "cm"), column_names_gp = gpar(fontsize = 9), 
    row_names_gp = gpar(fontsize = 9), show_row_dend = F, show_row_names = T, show_column_dend = F)

3.9.3 I46 (Clostridium innocuum)

dat <- readRDS("data/rds/reseq.rds")
# dat <- dat[which(dat$mouse.group == 'Water' | dat$mouse.group == 'Tetracyclin' ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "I46"), ]
rownames(data.wide) <- paste0(data.wide$POS, " (", data.wide$REF, " -> ", data.wide$ALT, ") ", 
    data.wide$feature)

heat <- data.matrix(data.wide[, 6:9])
heat <- heat[which(rowSums(!is.na(heat)) > 2), ]  # filter uncommon

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))
Heatmap(heat, name = "AF", col = col_fun, row_names_max_width = unit(10, "cm"), column_names_gp = gpar(fontsize = 9), 
    row_names_gp = gpar(fontsize = 9), show_row_dend = F, show_row_names = T, show_column_dend = F)

3.9.4 Filter: not in water but in AB

3.9.5 KB1 (Enterococcusfaecalis)

dat <- readRDS("data/rds/reseq.rds")
#dat <- dat[which(dat$mouse.group == "Water" | dat$mouse.group == "Tetracyclin" ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature  ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "KB1"),]
rownames(data.wide) <- paste0(data.wide$POS, " (",data.wide$REF, " -> ",data.wide$ALT, ") ", data.wide$feature)

heat <- data.wide[,6:9]
heat[is.na(heat)] <- 0
heat <- heat[which(heat$Water < 0.5 & (heat$Ciprofloxacin > 0.5 | heat$Tetracyclin > 0.5 | heat$Vancomycin > 0.5) ),] # filter water
heat <- data.matrix(heat)    

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))

Heatmap(heat, name = "AF", col = col_fun,
 column_names_gp = gpar(fontsize = 9),
  row_names_max_width  = unit(15, "cm"),

 row_names_gp = gpar(fontsize = 9),
 show_row_dend = F,
 show_row_names = T,
 show_column_dend = F
)

3.9.6 YL58 (Blautia coccoides)

dat <- readRDS("data/rds/reseq.rds")
#dat <- dat[which(dat$mouse.group == "Water" | dat$mouse.group == "Tetracyclin" ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature  ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "YL58"),]
rownames(data.wide) <- paste0(data.wide$POS, " (",data.wide$REF, " -> ",data.wide$ALT, ") ", data.wide$feature)

heat <- data.wide[,6:9]
heat[is.na(heat)] <- 0
heat <- heat[which(heat$Water < 0.5 & (heat$Ciprofloxacin > 0.5 | heat$Tetracyclin > 0.5 | heat$Vancomycin > 0.5) ),] # filter water
heat <- data.matrix(heat)    

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))

Heatmap(heat, name = "AF", col = col_fun,
 column_names_gp = gpar(fontsize = 9),
row_names_max_width  = unit(15, "cm"),

 row_names_gp = gpar(fontsize = 9),
 show_row_dend = F,
 show_row_names = T,
 show_column_dend = F
)

3.9.7 I46 (Clostridium innocuum)

dat <- readRDS("data/rds/reseq.rds")
#dat <- dat[which(dat$mouse.group == "Water" | dat$mouse.group == "Tetracyclin" ),]
data.wide <- dcast(dat, genome + POS + REF + ALT + feature  ~ mouse.group, value.var = "AF")

## Warning in dcast(dat, genome + POS + REF + ALT + feature ~ mouse.group, : The dcast
## generic in data.table has been passed a data.frame and will attempt to redirect to the
## reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. Please do this redirection yourself like reshape2::dcast(dat). In the
## next version, this warning will become an error.

data.wide <- data.wide[which(data.wide$genome == "I46"),]
rownames(data.wide) <- paste0(data.wide$POS, " (",data.wide$REF, " -> ",data.wide$ALT, ") ", data.wide$feature)

heat <- data.wide[,6:9]
heat[is.na(heat)] <- 0
heat <- heat[which(heat$Water < 0.5 & (heat$Ciprofloxacin > 0.5 | heat$Tetracyclin > 0.5 | heat$Vancomycin > 0.5) ),] # filter water
heat <- data.matrix(heat)    

library(circlize)
library(ComplexHeatmap)
heat[is.na(heat)] <- 0
col_fun = colorRamp2(c(0, 0.5, 1), c("white", "yellow", "red"))

Heatmap(heat, name = "AF", col = col_fun,
 column_names_gp = gpar(fontsize = 9),
  row_names_max_width  = unit(10, "cm"),

 row_names_gp = gpar(fontsize = 9),
 show_row_dend = F,
 show_row_names = T,
 show_column_dend = F
)