fig_5D.Rmd

---
title: "Fig_5D"
author: "Max Schubert"
date: "6/18/2019"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

library(data.table)
#library(stringr)
library(plyr)
library(gridExtra)
require(ggplot2)
require(data.table)
require(scales)
library(ggpubr)
library(plotly)
#dan's suggestion for plotting:
library(ggrepel)

parent_folder = here::here('fig_5D')
source(here::here('constants.R'))
```


###### trimmed version
```{r get_trimmed_data}
#the new stuff, annotated
trimmed_table = "t1000_t1500_trimmedvariants.csv"
trimmed_path = paste(parent_folder, trimmed_table, sep = "/")
trimmed_dt = data.table(fread(trimmed_path))

#ditch columns with na's
#cutting down table to only a few columns
trimmed_dt = Filter(function(x)!all(is.na(x)), trimmed_dt) #kills columns with all NA
trimmed_dt = trimmed_dt[,c( "UID", "POSITION",
                                        "REF", "ALT", "AF", "IS_HET",
                                        "EXPERIMENT_SAMPLE_LABEL", "DP",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE")] #keeps only specified columns
#throw out alleles with low allelic fraction
trimmed_dt = subset(trimmed_dt, AF > .90)

#compute summary stats
trimmed_dt[,mean_dp := mean(DP) , by=c("EXPERIMENT_SAMPLE_LABEL")]
trimmed_dt[,sd_dp := sd(DP) , by=c("EXPERIMENT_SAMPLE_LABEL")]
#compute z scores
trimmed_dt[, z_score := (DP- mean_dp)/sd_dp]
```
```{r fix variant effect annotations}
# frameshifts, stops, and alterations of start codons are nonsense
trimmed_dt[(INFO_EFF_EFFECT == "frameshift_variant") |
             (INFO_EFF_EFFECT == "stop_gained"),
           INFO_EFF_EFFECT := "Nonsense variant"]
# upstream, downstream, intergenic will all be intergenic
trimmed_dt[(INFO_EFF_EFFECT == "upstream_gene_variant") |
             (INFO_EFF_EFFECT == "downstream_gene_variant") |
             (INFO_EFF_EFFECT == "intergenic_region"),
           INFO_EFF_EFFECT := "Intergenic Variant"]
#discard annotations that only occur rarely and/or in rare SNPs
trimmed_dt= trimmed_dt[INFO_EFF_EFFECT != "splice_acceptor_variant+intron_variant"]
trimmed_dt= trimmed_dt[INFO_EFF_EFFECT != "intragenic_variant"]
#rename start_lost
trimmed_dt[(INFO_EFF_EFFECT == "start_lost"),
           INFO_EFF_EFFECT := "Start Codon Variant"]
#rename missense and synonymous to be more aesthetically pleasing
trimmed_dt[(INFO_EFF_EFFECT == "missense_variant"),
           INFO_EFF_EFFECT := "Missense Variant"]
trimmed_dt[(INFO_EFF_EFFECT == "synonymous_variant"),
           INFO_EFF_EFFECT := "Synonymous Variant"]
```


```{r plot_trimmed}
###### new version with zscores
fig4d_z <- ggplot(trimmed_dt[EXPERIMENT_SAMPLE_LABEL == "t1000.cutadapt.trim.fastq.bz2"],
       aes(x=POSITION, y=z_score, label = INFO_EFF_GENE, color = INFO_EFF_EFFECT)) + 
  geom_point(size = 3, alpha = 0.8, stroke = 0) +
  theme_linedraw() +
  geom_text(vjust = -1.1, data = subset(trimmed_dt[EXPERIMENT_SAMPLE_LABEL == "t1000.cutadapt.trim.fastq.bz2"],
                                      z_score > 2), show.legend = FALSE) +
  scale_colour_brewer(palette = "Set1") +
  theme(
    axis.text.x = element_text(size = fig_font_size),
    axis.text.y = element_text(size = fig_font_size),
    axis.title.y = element_text(size = fig_title_font_size),
    axis.title.x = element_text(size = fig_title_font_size),
    strip.background=element_blank(),
    panel.grid.minor=element_blank(),
    panel.grid.major.x=element_blank(),
    panel.border=element_blank(),
    axis.line=element_line(),
    legend.title = element_blank()) +
    #legend.position = c(.8,.8)) +
  scale_y_continuous(breaks = c(0,1,2,3,4,5,6),
                     limits = c(-0.7, 6)) +
  labs(
    x = "Genome position",
    y = "Enrichment of SNP (Z-score)")

fig4d_z
```

```{r save_fig}
ggsave(here::here('fig5D_wide.pdf'), plot=fig4d_z, width=8, 
       height=2.5, units='in', dpi=300, device = cairo_pdf)
```










########## cruft below
#todo delete
```{r get data}
#can update when i get the new stuff
test_table = "variants_fromT1000T1500.csv"

test_path = paste(parent_folder, test_table, sep = "/")
test_file_dt = data.table(fread(test_path))

#ditch columns with na's
#cutting down table to only a few columns
test_file_dt = Filter(function(x)!all(is.na(x)), test_file_dt) #kills columns with all NA
test_file_trimmed_dt = test_file_dt[,c( "UID", "POSITION",
                                        "REF", "ALT", "AF", "IS_HET",
                                        "EXPERIMENT_SAMPLE_LABEL", "DP",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE")] #keeps only specified columns
#throw out alleles with low allelic fraction
test_file_trimmed_dt = subset(test_file_trimmed_dt, AF > .90)


#compute summary stats
summary <- test_file_trimmed_dt[,
                        list(
                          N= .N,
                          mean_dp = mean(DP),
                          sd= sd(DP)),
                          #se= sd(DP) / sqrt(.N)),
                        by=c("EXPERIMENT_SAMPLE_LABEL")]

#compute SD bins above mean
T1000 <- test_file_trimmed_dt[EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt"]

meandp1000 = summary[EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt"]$mean_dp
sddp1000 = summary[EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt"]$sd
T1000[, z_score := (DP- meandp1000)/sddp1000]
```

```{r plot_1000}

###### new version with zscores
ggplot(T1000, aes(x=POSITION, y=z_score, label = INFO_EFF_GENE)) + 
  geom_point(size = 1, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = -1, data = subset(T1000, z_score > 1.5)) +
  theme(text = element_text(size=18),
        panel.grid.minor.y=element_blank(),
        panel.grid.minor.x=element_blank(),) +
  scale_y_continuous(breaks = c(0,1,2,3,4,5,6),
                     limits = c(-0.7, 7))
```

######## new version with manual followup on SNP calling
```{r get_more_data}
#the new stuff, annotated
annotated_table = "T1000_var_10_27.csv"
annotated_path = paste(parent_folder, annotated_table, sep = "/")
new_dt = data.table(fread(annotated_path))

#ditch columns with na's
#cutting down table to only a few columns
new_dt = Filter(function(x)!all(is.na(x)), new_dt) #kills columns with all NA
new_trimmed_dt = new_dt[,c( "UID", "POSITION", "manual_examination",
                                        "REF", "ALT", "AF", "IS_HET",
                                        "EXPERIMENT_SAMPLE_LABEL", "DP",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE")] #keeps only specified columns
#throw out alleles with low allelic fraction
new_trimmed_dt = subset(new_trimmed_dt, AF > .90)
#throw out alleles that failed manual inspection (only looked at top hits)
new_inspected_dt = new_trimmed_dt[! manual_examination %in% FALSE]

#compute summary stats
new_summary <- new_inspected_dt[,
                        list(
                          N= .N,
                          mean_dp = mean(DP),
                          sd= sd(DP)),
                          #se= sd(DP) / sqrt(.N)),
                        by=c("EXPERIMENT_SAMPLE_LABEL")]

#compute SD bins above mean
#T1000 <- new_inspected_dt[EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt"]

meandp1000_new = new_summary$mean_dp
sddp1000_new = new_summary$sd
new_inspected_dt[, z_score := (DP- meandp1000_new)/sddp1000_new]
```
```{r plot_1000_new}

###### new version with zscores
ggplot(new_inspected_dt, aes(x=POSITION, y=z_score, label = INFO_EFF_GENE)) + 
  geom_point(size = 1, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = -1, data = subset(new_inspected_dt, z_score > 2)) +
  theme(text = element_text(size=18),
        panel.grid.minor.y=element_blank(),
        panel.grid.minor.x=element_blank(),) +
  scale_y_continuous(breaks = c(0,1,2,3,4,5,6),
                     limits = c(-0.7, 7))
```


#8/29/2019 trying new stuff
```{r plot_new}

#convert counts to relative measure of frequency
test_file_trimmed_dt[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
test_file_trimmed_dt[, snp_freq := DP/snp_obs]

ggplot(test_file_trimmed_dt, aes(x=POSITION, y=snp_freq, label = INFO_EFF_GENE)) +
  geom_point(size = 2, alpha = 0.5) +
  #geom_col(position = position_dodge(), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(test_file_trimmed_dt, snp_freq > .02 )) +
  theme(text = element_text(size=18)) + #+
  #coord_polar()
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=2, ncol=1, scales = 'free_y') #+
  #coord_polar() # +
  #geom_text_repel() +
  #labs(
  #  x = "Mutations, by genomic coordinate",
  #  y = "Depth of coverage",
  #  title = "Loci conferring TMP resistance in sheared gDNA RLR",
  #  color = "Sample")
```


```{r alternative_plot}
ggplot(test_file_trimmed_dt[EXPERIMENT_SAMPLE_LABEL == "T1500_cutadapt"], aes(x=POSITION, y=DP, color=INFO_EFF_EFFECT, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(test_file_trimmed_dt, DP > 1000 & EXPERIMENT_SAMPLE_LABEL == "T1500_cutadapt")) +
  theme(text = element_text(size=18)) #+
  #coord_polar()
  #facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=2, ncol=1, scales = 'free_y') +
  #coord_polar() # +
  #geom_text_repel() +
  #labs(
  #  x = "Mutations, by genomic coordinate",
  #  y = "Depth of coverage",
  #  title = "Loci conferring TMP resistance in sheared gDNA RLR",
  #  color = "Sample")
```


###########
get data from plasmid libraries to use as a "before" sample for deseq2 analysis


#############
### use libraries before treatment with TMP as a "before" or "untreated" sample to determine FDR of SNPs for TMP1000 experiment

```{r get_before_data_38AB_40AB}
before38_table = "38A_38B_before.csv"
before40_table = "40A_40B_before.csv"

before38_path = paste(parent_folder, before38_table, sep = "/")
before40_path = paste(parent_folder, before40_table, sep = "/")
before38_dt = data.table(fread(before38_path))
before40_dt = data.table(fread(before40_path))

#ditch columns with all NA
before38_dt = Filter(function(x)!all(is.na(x)), before38_dt)
before40_dt = Filter(function(x)!all(is.na(x)), before40_dt)

#append both into one table
before_all_dt = rbind(before38_dt, before40_dt)

```

```{r plot_before_data_38}

before_sample = c("38A_induced", "38B_induced")
before_into_plot = before_all_dt[EXPERIMENT_SAMPLE_LABEL %in% before_sample]

#convert counts to relative measure of frequency
before_into_plot[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
before_into_plot[, snp_freq := DP/snp_obs]

ggplot(before_into_plot, aes(x=POSITION, y=snp_freq, label = INFO_EFF_GENE)) +
  geom_point(size = 2, alpha = 0.5) +
  #geom_col(position = position_dodge(), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(before_into_plot, snp_freq > .02 )) +
  theme(text = element_text(size=18)) +
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL)

```

```{r plot_before_data_40}

before_sample = c("40A_induced", "40B_induced")
before_into_plot = before_all_dt[EXPERIMENT_SAMPLE_LABEL %in% before_sample]

#convert counts to relative measure of frequency
before_into_plot[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
before_into_plot[, snp_freq := DP/snp_obs]

ggplot(before_into_plot, aes(x=POSITION, y=snp_freq, label = INFO_EFF_GENE)) +
  geom_point(size = 2, alpha = 0.5) +
  #geom_col(position = position_dodge(), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(before_into_plot, snp_freq > .02 )) +
  theme(text = element_text(size=18)) +
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL)

```


###########
## input genomic sequencing data

```{r get_data_genomes}
table_488489 = "488_489_variants.csv"
path_488489 = paste(parent_folder, table_488489, sep = "/")
dt_488489 = data.table(fread(path_488489))
dt_488489 = Filter(function(x)!all(is.na(x)), dt_488489) #ditch columns with all NA
```
```{r plot_genome_SNPs}
#before_sample = c("40A_induced", "40B_induced")
#before_into_plot = before_all_dt[EXPERIMENT_SAMPLE_LABEL %in% before_sample]
#convert counts to relative measure of frequency
#before_into_plot[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
#before_into_plot[, snp_freq := DP/snp_obs]

ggplot(dt_488489[EXPERIMENT_SAMPLE_LABEL == 488], aes(x=POSITION, y=DP, label = INFO_EFF_GENE)) +
  geom_point(size = 2, alpha = 0.5) +
  #geom_col(position = position_dodge(), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(dt_488489, DP > 100 )) +
  theme(text = element_text(size=18)) 

```





############
## try liquid selection data

```{r load_liquid_data}
liquid_table = "e1e2e3e4_liquid_TMP800.csv"
liquid_path = paste(parent_folder, liquid_table, sep = "/")
liquid_dt = data.table(fread(liquid_path))

#ditch columns with all NA
liquid_dt = Filter(function(x)!all(is.na(x)), liquid_dt)
```

```{r plot_liquid_data}

#convert counts to relative measure of frequency
liquid_dt[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
liquid_dt[, snp_freq := DP/snp_obs]

ggplot(liquid_dt, aes(x=POSITION, y=snp_freq, label = INFO_EFF_GENE)) +
  geom_point(size = 2, alpha = 0.5) +
  #geom_col(position = position_dodge(), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(liquid_dt, snp_freq > .02 )) +
  theme(text = element_text(size=18)) +
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL)

```



#############

# G10 data


```{r get_G10_data}
#can update when i get the new stuff
G10_table = "G10_new_variants.csv"
G10_path = paste(parent_folder, G10_table, sep = "/")
G10_dt = data.table(fread(G10_path))

#ditch columns with na's
#cutting down table to only a few columns
G10_dt = Filter(function(x)!all(is.na(x)), G10_dt) #kills columns with all NA
G10_dt_trimmed = G10_dt[,c( "UID", "POSITION",
                                        "REF", "ALT",
                                        "EXPERIMENT_SAMPLE_LABEL", "DP",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE")] #keeps only specified columns

```

```{r get_muts}
#load 38a and 38B induced
induced_all = G10_dt[EXPERIMENT_SAMPLE_LABEL %in% c("38A_induced", "38B_induced", "40A_induced", "40B_induced")]

induced_all[, total_induced_depth := sum(DP), by = c("UID")]
induced_all = induced_all[,c( "UID", "POSITION",
                                        "REF", "ALT",
                                        "EXPERIMENT_SAMPLE_LABEL", "DP",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE", "total_induced_depth")]

induced_snps = induced_all[EXPERIMENT_SAMPLE_LABEL %in% c("38A_induced")]
induced_snps= induced_snps[,c( "UID", "POSITION",
                                        "REF", "ALT",
                                        "INFO_EFF_CONTEXT","INFO_EFF_EFFECT",
                                        "INFO_EFF_AA","INFO_EFF_GT",
                                        "INFO_EFF_CLASS", "INFO_EFF_GENE", "total_induced_depth")]
```


```{r plot_40}
dt40plus = G10_dt_trimmed[EXPERIMENT_SAMPLE_LABEL %in% c("40A", "40B")]

#convert counts to relative measure of frequency
dt40plus[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
dt40plus[, snp_freq := DP/snp_obs]



plot_40 <- ggplot(dt40plus, aes(x=POSITION, y=snp_freq, color=EXPERIMENT_SAMPLE_LABEL, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +
  #scale_y_continuous(trans = log10_trans() ) + #limit=c(0.0005, 8), 
                     #breaks = c(.001, .01, .1, .3, 1, 3, 10, 100, 1000)) +
  geom_text(vjust = 0, data = subset(dt40plus_filter, snp_freq >= .01)) +
  theme(text = element_text(size=18),
        legend.position = 'none') +
  #coord_polar()
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=1, ncol=2)
  #coord_polar() # +
  #geom_text_repel() +
  #labs(
  #  x = "Mutations, by genomic coordinate",
  #  y = "Depth of coverage",
  #  title = "Loci conferring TMP resistance in sheared gDNA RLR",
  #  color = "Sample")

plot_40
```

#38 are the samples with beta recombinase
```{r plot_38}
dt38plus = G10_dt_trimmed[EXPERIMENT_SAMPLE_LABEL %in% c("38Aplus", "38Bplus")]

#convert counts to relative measure of frequency
dt38plus[, snp_obs := sum(DP), by= 'EXPERIMENT_SAMPLE_LABEL']
dt38plus[, snp_freq := DP/snp_obs]

#remove folA snp and see what's left:
#dt40plus_filter = dt40plus[UID != "aa76968f"]
#dt40plus_filter = dt40plus_filter[snp_freq > .001]

#dt40plus_faves = dt40plus_filter[INFO_EFF_GENE %in% c('cstA', 'igaA', 'uvrD')]

plot_38 <- ggplot(dt38plus, aes(x=POSITION, y=snp_freq, color=EXPERIMENT_SAMPLE_LABEL, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +
  #scale_y_continuous(trans = log10_trans() ) + #limit=c(0.0005, 8), 
                     #breaks = c(.001, .01, .1, .3, 1, 3, 10, 100, 1000)) +
  geom_text(vjust = 0, data = subset(dt38plus, snp_freq >= .01)) +
  theme(text = element_text(size=18),
        legend.position = 'none') +
  #coord_polar()
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=1, ncol=2)
  #coord_polar() # +
  #geom_text_repel() +
  #labs(
  #  x = "Mutations, by genomic coordinate",
  #  y = "Depth of coverage",
  #  title = "Loci conferring TMP resistance in sheared gDNA RLR",
  #  color = "Sample")
plot_38
```

```{r}
#plot both in grid
grid.arrange(
  grobs = list(plot_38, plot_40),
  #widths = c(2, 1, 1),
  layout_matrix = rbind(c(1, 1),
                        c(2, 2))
)
```

### what's in there besides fola?

```{r no_folA}
#remove folA snp and see what's left:
dt40plus_filter = dt40plus[UID != "aa76968f"]
#dt40plus_filter = dt40plus_filter[snp_freq > .001]

plot_40_filter <- ggplot(dt40plus_filter, aes(x=POSITION, y=snp_freq, color=EXPERIMENT_SAMPLE_LABEL, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +

  geom_text(vjust = 0, data = subset(dt40plus_filter, snp_freq >= .01)) +
  theme(text = element_text(size=18),
        legend.position = 'none') +
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=1, ncol=2)


dt38plus_filter = dt38plus[UID != "aa76968f"]
plot_38_filter <- ggplot(dt38plus_filter, aes(x=POSITION, y=snp_freq, color=EXPERIMENT_SAMPLE_LABEL, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +
  scale_y_continuous(limit=c(0, .25) )+
                     #breaks = c(.001, .01, .1, .3, 1, 3, 10, 100, 1000)) +
  geom_text(vjust = 0, data = subset(dt38plus, snp_freq >= .01)) +
  theme(text = element_text(size=18),
        legend.position = 'none') +
  #coord_polar()
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=1, ncol=2)


#plot both in grid
grid.arrange(
  grobs = list(plot_38_filter, plot_40_filter),
  #widths = c(2, 1, 1),
  layout_matrix = rbind(c(1, 1),
                        c(2, 2))
)


#closer look at 40B
ggplot(dt40plus_filter[EXPERIMENT_SAMPLE_LABEL == "40B"], aes(x=POSITION, y=snp_freq, color=INFO_EFF_EFFECT, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +

  geom_text(vjust = 0, data = subset(dt40plus_filter, snp_freq >= .002)) +
  theme(text = element_text(size=18))
        
        
        legend.position = 'none') #+
  facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=1, ncol=2)


```





###### old version
ggplot(test_file_trimmed_dt[EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt"], aes(x=POSITION, y=DP, color=EXPERIMENT_SAMPLE_LABEL, label = INFO_EFF_GENE)) + 
  geom_col(position = position_dodge(width = 0.9), size = 2, alpha = 0.5) +
  theme_bw() +
  geom_text(vjust = 0, data = subset(test_file_trimmed_dt, DP > 1000 & EXPERIMENT_SAMPLE_LABEL == "T1000_cutadapt")) +


  #line at mean
  
  #line at 
  
  #+
  #coord_polar()
  #facet_wrap(~EXPERIMENT_SAMPLE_LABEL, nrow=2, ncol=1, scales = 'free_y') +
  #coord_polar() # +
  #geom_text_repel() +
  #labs(
  #  x = "Mutations, by genomic coordinate",
  #  y = "Depth of coverage",
  #  title = "Loci conferring TMP resistance in sheared gDNA RLR",
  #  color = "Sample")