timo_0.01_SNVDiversity.Rmd

---
title: "R Notebook"
output: html_notebook
---

```{r}
library("tidyr")
library('ggplot2')
library('dplyr')
library("glue")

wkdir = "~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01"
setwd(wkdir)
savedir = "~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01/Output_Figures"

source("~/Desktop/GitHub/Obesity/NewExtractions/H9N2/FD_functions.R")
```

```{r}
diet = c("Obese","Lean","Control")
dietColors = c("#FF9933","#66CCFF","#606060")
names(dietColors) = diet
DietcolScale_fill <- scale_fill_manual(name = "grp",values = dietColors)
DietcolScale <- scale_colour_manual(name = "grp",values = dietColors)
```

# Specifying thresholds and plotting variables
```{r}
cov_cut = 200
freq_cut = 0.01
pvalcut  = 0.05

ntlist = c("A","C","G","T")
SEGMENTS = c('H9N2_PB2','H9N2_PB1','H9N2_PA','H9N2_HA','H9N2_NP','H9N2_NA','H9N2_MP','H9N2_NS')
```

#Loading metadata
This includes titer and Ct values when applicable. ND indicates qPCR was run with a negative result; 0 indicates plaque assay or HAI was run with a negative result. NA for any values indicate that data was missing. Sacrificed indicates there was no data at that time point because the ferret had already been sacrficied for pathology. 
```{r}
metafile = "~/Desktop/GitHub/Obesity/NewExtractions/H9N2/H9_Metadata.csv"

meta = read.csv(file=metafile,header=T,sep=",",na.strings = c(''))
meta = filter(meta, resequenced == "yes")

meta$Ct_Mgene = as.numeric(meta$Ct_Mgene)
meta$titer = as.numeric(meta$titer)
meta$log10_titer = as.numeric(meta$log10_titer)

meta$inf_route = factor(meta$inf_route, levels = c("Index","Contact","Aerosol","Control"))
```

# Loading in coverage file & segment size information
```{r}
cov = read.csv("./avg_coverage/H9N2.coverage.csv", header = TRUE, sep = ",")

seg_sizes = "../SegmentSize.csv"
#manually changed NS size from 838 to 864 to match the longest length found in the cov data rather than the Segemnt size value from the .gtf file - think about this
sizes = read.csv(file=seg_sizes,header=T,sep=",",na.strings = c(''))
GenomeSize = (sizes %>% filter(segment == 'H9N2_GENOME'))$SegmentSize

cov$segment = factor(cov$segment, levels = SEGMENTS)
```

# Checking if data passes thresholds 
```{r}
cov_check = CoverageAcross(cov,cov_cut,40,sizes, wkdir)
```

Merging coverage check info with the rest of the metadata
```{r}
meta = merge(meta, cov_check, by.x = c("sample"), by.y = c("name"), all.y = TRUE)

nrow(meta)
count(meta,quality)
```

# Loading in variant files
```{r}
varfile = "./varfiles/H9N2.VariantsOnly.0.01.200.csv"

# read and rearrange the data
vars = read.csv(file=varfile,header=T,sep=",",na.strings = c(''))
vars$name = vars$sample
```

# Rearranging variant dataframe
```{r}
vdf = ArrangeVarWRep(vars)
# already have replicate data in the varfiles from running CompareReps.v2.py script
vdf = vdf[!duplicated(vdf), ] %>% droplevels()
nrow(vdf)
```

# Adding metadata
```{r}
vdf = merge(vdf,meta, by = c("sample","segment"))
vdf = vdf[!duplicated(vdf), ] %>% droplevels()

vdf$segment = factor(vdf$segment, levels = SEGMENTS)
```

# How many unique ferrets? How many unique good quality ferrets?
```{r}
select(vdf, ferretID) %>% unique() %>% count()
#total number of unique ferrets that were sequenced

select(vdf, ferretID,quality) %>% unique() %>% filter(quality == "good") %>% nrow()
select(vdf, ferretID,quality) %>% unique() %>% filter(quality == "good")
#number of ferrets with at least one sample passing cov_check QC

filter(vdf, quality == "bad") %>% select(sample) %>% unique()
#samples that failed cov_check QC
# mostly late time points
# d02 time points: 1417 is an aerosol, most other aerosol samples weren't resequenced 
# 1913 is contact lean,1981 and 1986 and 2235 are contact obese 
```

```{r}
nrow(vdf)
vdf = filter(vdf, quality == "good") %>% unique()
nrow(vdf)

vdf = filter(vdf, sample != "2244_d12") %>% unique()
#technically passes due to the low coverage threhold needed to include W17 but has incomplete segments that aren't the Pol genes

good_names = c(levels(factor(vdf$sample)))
```

# Filtering variant df with frequency cutoffs
```{r}
nrow(vdf)
vdf = filter(vdf, minorfreq1 >= freq_cut & 
               minorfreq2 >= freq_cut & 
               minor %in% ntlist &
               major %in% ntlist) %>% 
            droplevels()
# based on MAF study, reps and 0.01% cutoff was best combo
#filter each replicate separately rather than using the average
nrow(vdf)
# does not eliminate any variants here

vdf = filter(vdf, inf_route == "Index" | inf_route == "Contact" | inf_route == "Control") %>% unique()
# ignoring aerosol for now

vdf = filter(vdf, !(ferretID == "2232" & inf_route == "Index")) %>% filter(!(ferretID == "2239" & inf_route == "Index"))
# these two ferrets are listed as both index and contacts but I'm only using the contact info as the second infection wasn't successful
```

```{r}
vdf_cons = filter(vdf, major != stocknt) %>% 
  select(ferretID, DPI, segment, ntpos, aapos, stocknt, stockaa, major, majoraa, majorfreq, nonsyn) %>% unique()
vdf_cons$ntvar = paste0(vdf_cons$stocknt,vdf_cons$ntpos,vdf_cons$major)
vdf_cons$aavar = paste0(vdf_cons$stockaa,vdf_cons$aapos,vdf_cons$majoraa)
vdf_cons$var = vdf_cons$major
vdf_cons$var_aa = vdf_cons$majoraa
vdf_cons$freq = vdf_cons$majorfreq
vdf_cons = select(vdf_cons, ferretID, DPI, segment, ntpos, aapos, stocknt, stockaa, var, var_aa, freq, nonsyn,ntvar,aavar)

vdf_mins = filter(vdf, major == stocknt) %>% 
  select(ferretID, DPI, segment, ntpos, aapos, stocknt, stockaa, minor, minoraa, minorfreq, nonsyn) %>% unique()
vdf_mins$ntvar = paste0(vdf_mins$stocknt,vdf_mins$ntpos,vdf_mins$minor)
vdf_mins$aavar = paste0(vdf_mins$stockaa,vdf_mins$aapos,vdf_mins$minoraa)
vdf_mins$var = vdf_mins$minor
vdf_mins$var_aa = vdf_mins$minoraa
vdf_mins$freq = vdf_mins$minorfreq
vdf_mins = select(vdf_mins, ferretID, DPI, segment, ntpos, aapos, stocknt, stockaa, var, var_aa, freq, nonsyn,ntvar,aavar)

vdf_out = rbind(vdf_cons,vdf_mins) %>% unique()
write.table(vdf_out, file = "SNVs_Over_Thresholds.csv", sep = ",", row.names = FALSE)
```

# Filtering variant df by timo binocheck
```{r}
#vdf$binocheck = factor(vdf$binocheck, levels = c("False","R1","R2","True"))
#vdf = filter(vdf, binocheck != "False") %>% unique()
#nrow(vdf)

# binocheck is highly dependent on the allele frequency threshold used and also relatively conservative
# as a result, ignore this in favor of found in both replicates across ferrets and cohorts - this is more indicative of a real variant than binocheck
```

# How many variants per sample?
```{r}
nums = vdf %>% group_by(sample, diet, inf_route) %>% tally()
range(nums$n)
mean(nums$n)
median(nums$n)
nums = separate(nums, sample, into = c("ferretID","DPI"))
nums$samp = paste0(nums$DPI,"_",nums$ferretID)

SNVS_per_sample = ggplot(nums, aes(x = samp, y = n)) +
  geom_col(aes(fill = diet)) +
  PlotTheme1 +
  facet_grid(diet~DPI+inf_route) +
  DietcolScale_fill
print(SNVS_per_sample)
ggsave(SNVS_per_sample, filename = "SNVS_per_sample.pdf", path = savedir, width = 10, height = 5)

# this would only be for the minor variants
#ggplot(vdf, aes(x = sample, y = minorfreq)) +
#  geom_boxplot(outlier.shape = NA) +
#  PlotTheme1
```

# Loading transmission information
```{r}
transmission_info = "/Users/marissaknoll/Desktop/GitHub/Obesity/NewExtractions/H9N2/TransmissionPairs.csv"
pairs = read.csv(transmission_info, header = T)
```

# Plotting number of successful ferrets
```{r}
fercount = select(vdf,sample,ferretID,DPI,diet,inf_route)
fercount = fercount[!duplicated(fercount), ]  %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,DPI) %>% 
  tally()

fercount = separate(fercount,sample,into = c("ferretID","DPI"))
fercount2 = merge(fercount, pairs, by = c("ferretID"))
  
p1 = fercount2 %>% unique() %>% 
    ggplot(., aes(x= DPI, y = pair_numbers, fill = diet)) + 
    geom_tile(color = 'black') + 
    PlotTheme3 +
    DietcolScale_fill + 
    facet_grid(pair_diets~inf_route, scales = 'free', space = 'free')
print(p1)
ggsave("ferrets_tileplot.pdf", p1, path = savedir,)
```

# Pulling out consensus changes
```{r}
con_change = filter(vdf, stocknt != major) %>%
  filter(major %in% ntlist)
con_change = con_change[!duplicated(con_change), ]
con_change$var = paste0(con_change$ferretID,"_",con_change$segment,"_",
                        con_change$major,"_",con_change$ntpos,"_",con_change$minor)

consensus = unique(con_change$var)
length(consensus)

```

# Pulling out minor variants
```{r}
vdf$var = paste0(vdf$ferretID,"_",vdf$segment,"_",vdf$major,"_",vdf$ntpos,"_",vdf$minor)

minorvdf = filter(vdf, !(var %in% consensus))
minorvdf = minorvdf[!duplicated(minorvdf), ]
nrow(vdf) - nrow(minorvdf)
```

# Distribution of variants by AF

```{r}
min_index = filter(minorvdf, inf_route == "Index")

minorvdf_0.03 = filter(min_index, minorfreq < 0.03) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>% 
  mutate(cat = "1-2.99%")

minorvdf_0.05 = filter(min_index, minorfreq < 0.05 & minorfreq >= 0.03) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>% 
  mutate(cat = "3-4.99%")

minorvdf_0.1 = filter(min_index, minorfreq < 0.1 & minorfreq >= 0.05) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>% 
  mutate(cat = "5-9.99%")

minorvdf_0.3 = filter(min_index, minorfreq < 0.3 & minorfreq >= 0.1) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>% 
  mutate(cat = "10-29.99%")

minorvdf_0.5 = filter(min_index, minorfreq < 0.5 & minorfreq >= 0.3) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>%
  mutate(cat = "30-49.99%")

minors = rbind(minorvdf_0.03,minorvdf_0.05, minorvdf_0.1, minorvdf_0.3, minorvdf_0.5)
minors$cat = factor(minors$cat, levels = c("1-2.99%","3-4.99%","5-9.99%","10-29.99%","30-49.99%"))
mins = separate(minors, sample, into = c("ferretID","DPI"))
```

```{r}
con_change_1.0 = filter(con_change, inf_route == "Index" & ferretID != "2232") %>% 
  unique() %>% 
  group_by(sample,diet,inf_route) %>% 
  tally() %>%
  mutate(cat = "50-99.99%")

majors = con_change_1.0
majors$cat = factor(majors$cat, levels = c("50-99.99%"))
majs = separate(majors, sample, into = c("ferretID","DPI"))
```

```{r}
days = c("d02","d04","d06")
mins_in = filter(mins, inf_route == "Index" & DPI %in% days)
majs_in = filter(majs, inf_route == "Index" & DPI %in% days)

freq_dist = ggplot() +
  geom_boxplot(aes(x = cat, y = n, color = diet), mins_in, outlier.shape = NA) +
  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), mins_in, width = 0.2, alpha = 0.5) +
  geom_boxplot(aes(x = cat, y = n, color = diet), majs_in, outlier.shape = NA) +
  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), majs_in, width = 0.2, alpha = 0.5) +
  PlotTheme1 +
  DietcolScale
print(freq_dist)
ggsave(freq_dist, file = "freq_dist.pdf", path = savedir, width = 10, height = 7)

# This woudl be to look at all ferrets, not just index
#ggplot() +
#  geom_boxplot(aes(x = cat, y = n, color = diet), minors, outlier.shape = NA) +
#  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), minors, width = 0.2) +
#  geom_boxplot(aes(x = cat, y = n, color = diet), majors, outlier.shape = NA) +
#  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), majors, width = 0.2) +
#  PlotTheme1 +
#  DietcolScale
```

# Remake this plot separated out by nonsyn and syn
```{r}
minorvdf_0.03 = filter(minorvdf, minorfreq < 0.03) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>% 
  mutate(cat = "1-2.99%")

minorvdf_0.05 = filter(minorvdf, minorfreq < 0.05 & minorfreq >= 0.03) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>% 
  mutate(cat = "3-4.99%")

minorvdf_0.1 = filter(minorvdf, minorfreq < 0.1 & minorfreq >= 0.05) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>% 
  mutate(cat = "5-9.99%")

minorvdf_0.3 = filter(minorvdf, minorfreq < 0.3 & minorfreq >= 0.1) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>% 
  mutate(cat = "10-29.99%")

minorvdf_0.5 = filter(minorvdf, minorfreq < 0.5 & minorfreq >= 0.3) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>%
  mutate(cat = "30-49.99%")

minors = rbind(minorvdf_0.03,minorvdf_0.05, minorvdf_0.1, minorvdf_0.3, minorvdf_0.5)
minors$cat = factor(minors$cat, levels = c("1-2.99%","3-4.99%","5-9.99%","10-29.99%","30-49.99%"))
mins = separate(minors, sample, into = c("ferretID","DPI"))
```

```{r}
con_change_1.0 = filter(con_change, majorfreq < 1.0 & majorfreq >= 0.5) %>% 
  unique() %>% 
  group_by(sample,diet,inf_route,nonsyn) %>% 
  tally() %>%
  mutate(cat = "50-99.99%")

majors = con_change_1.0
majors$cat = factor(majors$cat, levels = c("50-99.99%"))
majs = separate(majors, sample, into = c("ferretID","DPI"))
```

```{r}
days = c("d02","d04","d06")
mins_in = filter(mins, inf_route == "Index" & DPI %in% days) %>% filter(!(is.na(nonsyn)))
mins_in$type = paste0(mins_in$nonsyn,"_",mins_in$cat)
majs_in = filter(majs, inf_route == "Index" & DPI %in% days)%>% filter(!(is.na(nonsyn)))
majs_in$type = paste0(majs_in$nonsyn,"_",majs_in$cat)

freq_dist = ggplot() +
  geom_boxplot(aes(x = cat, y = n, color = diet), mins_in, outlier.shape = NA) +
  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), mins_in, width = 0.2, alpha = 0.5) +
  geom_boxplot(aes(x = cat, y = n, color = diet), majs_in, outlier.shape = NA) +
  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), majs_in, width = 0.2, alpha = 0.5) +
  #facet_grid(~DPI) +
  PlotTheme1 +
  DietcolScale
print(freq_dist)
ggsave(freq_dist, file = "freq_dist.pdf", path = savedir, width = 10, height = 7)

ggplot() +
  geom_boxplot(aes(x = nonsyn, y = n, color = diet), mins_in, outlier.shape = NA) +
  geom_jitter(aes(x = nonsyn, y = n, color = diet, group = diet), mins_in, width = 0.2, alpha = 0.5) +
  #geom_boxplot(aes(x = cat, y = n, color = diet), majs_in, outlier.shape = NA) +
  #geom_jitter(aes(x = cat, y = n, color = diet, group = diet), majs_in, width = 0.2, alpha = 0.5) +
  facet_grid(~cat) +
  PlotTheme1 +
  DietcolScale
print(freq_dist)
ggsave(freq_dist, file = "freq_dist.pdf", path = savedir, width = 10, height = 7)

# This woudl be to look at all ferrets, not just index
#ggplot() +
#  geom_boxplot(aes(x = cat, y = n, color = diet), minors, outlier.shape = NA) +
#  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), minors, width = 0.2) +
#  geom_boxplot(aes(x = cat, y = n, color = diet), majors, outlier.shape = NA) +
#  geom_jitter(aes(x = cat, y = n, color = diet, group = diet), majors, width = 0.2) +
#  PlotTheme1 +
#  DietcolScale
```

# Testing for significance
```{r}
o = filter(minors, inf_route == "Index" & cat == "1-2.99%" & nonsyn == "syn" & diet == "Obese")
l = filter(minors, inf_route == "Index" & cat == "1-2.99%" & nonsyn == "syn" & diet == "Lean")
t.test(o$n,l$n)

# Significance: none
```

# Tallying SNVs
```{r}
# can make these groupings whatever you want

# count the number of SNVs per sample
group_list_seg = c('ferretID','segment',"DPI","diet","inf_route","cohort") # counts across each segment 
group_list_gen = c('ferretID',"DPI","diet","inf_route","cohort") # Counts across entire genome

seg_count = TallyIt(vdf, group_list_seg, "snv_count")
gen_count = TallyIt(vdf, group_list_gen, "snv_count") 

nonsyn_count_gen = c('ferretID',"DPI","diet","inf_route","cohort", "nonsyn")
ngen_count = TallyIt(vdf, nonsyn_count_gen, "snv_count") %>% filter(!is.na(nonsyn))
```

```{r}
# INCLUDING SEGMENTS WITH NO SNVS - but only using those that passed seq cutoff
reseq_seg = select(meta,ferretID,segment,DPI,diet,inf_route,cohort, quality) %>% 
  filter(quality == "good") %>% 
  unique()
seg_count = merge(seg_count,reseq_seg, all= TRUE)
seg_count = seg_count[!duplicated(seg_count), ]

seg_count$snv_count[is.na(seg_count$snv_count)] = 0
seg_count = filter(seg_count, !is.na(ferretID))
```

```{r}
reseq_gen = select(meta,ferretID,DPI,diet,inf_route,cohort,quality) %>% 
  filter(quality == "good") %>% 
  unique()
gen_count = merge(gen_count,reseq_gen, all = TRUE)
gen_count = gen_count[!duplicated(gen_count), ]

gen_count$snv_count[is.na(gen_count$snv_count)] = 0
gen_count = filter(gen_count, !is.na(ferretID))
```

# Average Number of Variants per Sample
```{r}
gen_count_avg = group_by(gen_count, DPI, diet, inf_route) %>%
  mutate(avgSNV = mean(snv_count), sdSNV = sd(snv_count))

seg_count_avg = group_by(seg_count, DPI, diet) %>%
  mutate(avgSNV = mean(snv_count), sdSNV = sd(snv_count))
```

```{r}
gen_count_segment_plot = filter(gen_count, inf_route == "Index" | inf_route == "Control") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock") %>%
ggplot(. , aes(x = diet, y = snv_count)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, aes(color = diet)) +
  facet_grid(~DPI) +
  ylab("SNV richness per sample") +
  PlotTheme1 +
  DietcolScale
print(gen_count_segment_plot)
#ggsave("gen_count_segment_plot.pdf",gen_count_segment_plot,path=savedir, width = 10, height = 5)

gen_count2 = gen_count %>% group_by(DPI,diet,inf_route) %>% mutate(avg_snv = mean(snv_count)) %>% ungroup() %>%
  filter(inf_route == "Index" | inf_route == "Control") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock") %>%
ggplot(. , aes(x = factor(DPI, levels = c("Stock","d02","d04","d06")), color = diet)) +
  #geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(y = snv_count), size = 1, width = 0.2, alpha = 0.5) +
  geom_point(aes(y = avg_snv), size = 3) +
  geom_line(aes(y = avg_snv, group = diet), size = 1.5) +
  ylab("SNV richness") +
  xlab("DPI") +
  PlotTheme1 +
  DietcolScale
print(gen_count2)
ggsave("gen_count2.pdf",gen_count2,path=savedir, width = 10, height = 7)

ngen_count$cat = paste0(ngen_count$diet,"_",ngen_count$nonsyn)

nonsyn_count = ngen_count %>% group_by(DPI,diet,inf_route,nonsyn) %>% mutate(avg_snv = mean(snv_count)) %>% ungroup() %>%
  filter(inf_route == "Index" | inf_route == "Control") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock") %>%
ggplot(. , aes(x = factor(DPI, levels = c("Stock","d02","d04","d06")), color = diet)) +
  geom_jitter(aes(y = snv_count, shape = nonsyn), size = 1, width = 0.2, alpha = 0.5) +
  geom_point(aes(y = avg_snv), size = 3) +
  geom_line(aes(y = avg_snv, group = cat, linetype = nonsyn), size = 1.5) +
  #geom_boxplot(aes(y = snv_count), outlier.shape = NA) +
  #geom_jitter(aes(y = snv_count, color = diet), width = 0.2) +
  ylab("SNV richness") +
  xlab("DPI") +
  #facet_grid(~diet) +
  PlotTheme1 +
  DietcolScale
print(nonsyn_count)
ggsave(nonsyn_count, file = "nonsyn_count.pdf", path = savedir, width = 12, height = 9)

counts = filter(ngen_count,inf_route == "Index", nonsyn == "nonsyn") %>% filter(DPI == "d02" | DPI == "d04" | DPI == "d06")
aov.counts = aov(snv_count ~ diet, counts)
summary(aov.counts)

o = filter(ngen_count, inf_route == "Index" & nonsyn == "syn" & diet == "Obese" & DPI == "d06")
l = filter(ngen_count, inf_route == "Index" & nonsyn == "syn" & diet == "Lean" & DPI == "d06")
t.test(o$snv_count, l$snv_count)
```

By segment
```{r}
seg_count$diet_segment = paste0(seg_count$diet,"_",seg_count$segment)
seg_count$diet_segment = factor(seg_count$diet_segment, 
                                levels =c ('Control_H9N2_PB2','Lean_H9N2_PB2','Obese_H9N2_PB2',
                                           'Control_H9N2_PB1','Lean_H9N2_PB1','Obese_H9N2_PB1',
                                           'Control_H9N2_PA','Lean_H9N2_PA','Obese_H9N2_PA',
                                           'Control_H9N2_HA','Lean_H9N2_HA','Obese_H9N2_HA',
                                           'Control_H9N2_NP','Lean_H9N2_NP','Obese_H9N2_NP',
                                           'Control_H9N2_NA','Lean_H9N2_NA','Obese_H9N2_NA',
                                           'Control_H9N2_MP','Lean_H9N2_MP','Obese_H9N2_MP',
                                           'Control_H9N2_NS','Lean_H9N2_NS','Obese_H9N2_NS'))

snv_count_segment_plot = filter(seg_count, inf_route == "Index" | inf_route == "Control") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock") %>%
ggplot(. , aes(x = diet_segment, y = snv_count)) +
  geom_boxplot(outlier.shape = NA, position = "dodge", aes(color = diet)) +
  #geom_jitter(width = 0.2, aes(color = diet)) +
  facet_grid(~DPI) +
  ylab("SNV richness per sample per segment") +
  PlotTheme1 +
  DietcolScale
print(snv_count_segment_plot)
ggsave("snv_count_segment_plot.pdf",snv_count_segment_plot,path=savedir, width = 10, height = 5)
# 03/07 now includes consensus changes (using vdf instead of minorvdf)
```

# T test for segment comparisons
```{r}
o = filter(seg_count, inf_route == "Index" & diet == "Obese" & DPI == "d06" & segment == "H9N2_PB2")
l = filter(seg_count, inf_route == "Index" & diet == "Lean" & DPI == "d06" & segment == "H9N2_PB2")
t.test(o$snv_count, l$snv_count)

# d06: NS are sig. diff
# nothing else is sig. diff

o_gen = filter(gen_count, inf_route == "Index" & diet == "Obese" & DPI == "d06")
l_gen = filter(gen_count, inf_route == "Index" & diet == "Lean" & DPI == "d06")
t.test(o_gen$snv_count, l_gen$snv_count)
```

# Calculating Shannon Entropy
```{r}
vdf = ShannonPos(vdf)
vdf$SegmentSize = as.numeric(vdf$SegmentSize)
vdf$shannon_perkb = (vdf$segment_shan/(vdf$SegmentSize/1000))
vdf$normalized_shannon = (vdf$shannon/GenomeSize)
```

```{r}
# shannon_ntpos = shannon entropy at that nt pos - should always be between 0 and 1 for each sample
# segment_shan = sum of all nt_pos per segment for each sample
# shannon = sum of all segment_shan across genome for each sample
# shannon_perkb = segment shannon per kb (segment specific) for each sample
# normalized_shannon = shannon divided by genome size (can make per kb by dividing by 1000) for each sample
```

# Average Shannon Entropy per Site per Sample (using normalized_shannon)
```{r}
shan_g = ungroup(vdf) %>%
  select(ferretID, DPI, diet, inf_route, cohort,normalized_shannon) %>%
  unique()

shan_g = merge(shan_g,reseq_gen %>% unique(), by = c("ferretID", "DPI", "diet", "inf_route","cohort"), all= TRUE) %>%
  filter(inf_route == "Index" | inf_route == "Control") %>% 
  unique()
shan_g$normalized_shannon[is.na(shan_g$normalized_shannon)] = 0

dim(shan_g)
shan_g <- shan_g[complete.cases(shan_g), ] 
dim(shan_g)

shan_gen_plot = ggplot(filter(shan_g, DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock"), 
       aes(x = diet, y = normalized_shannon/1000)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, aes(color = diet)) +
  facet_grid(~DPI) +
  xlab("Normalized Shannon entropy per kB") +
  PlotTheme1 +
  DietcolScale
print(shan_gen_plot)
ggsave("shan_gen_plot.pdf", shan_gen_plot, path = savedir, width = 10, height = 5)

# 03/07 now includes consensus changes (using vdf instead of minorvdf)

shan_gen_plot2 = shan_g %>% group_by(DPI,diet,inf_route) %>% mutate(avg_shan = mean(normalized_shannon)) %>% ungroup() %>%
  filter(inf_route == "Index" | inf_route == "Control") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06" | DPI == "Stock") %>%
ggplot(. , aes(x = factor(DPI, levels = c("Stock","d02","d04","d06")), color = diet)) +
  #geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(y = normalized_shannon), size = 1, width = 0.2, alpha = 0.5) +
  geom_point(size = 3, aes(y = avg_shan)) +
  geom_line(aes(y = avg_shan, group = diet), linewidth = 1.5) +
  #facet_grid(~DPI) +
  ylab("Shannon entropy") +
  xlab("DPI") +
  PlotTheme1 +
  DietcolScale
print(shan_gen_plot2)
ggsave("shan_gen_plot2.pdf", shan_gen_plot2, path = savedir, width = 10, height = 7)

counts = filter(shan_g,inf_route == "Index") %>% filter(DPI == "d02" | DPI == "d04" | DPI == "d06")
aov.counts = aov(normalized_shannon ~ diet, counts)
summary(aov.counts)
```

# Test for significance
```{r}
o = filter(shan_g, DPI == "Stock" & diet == "Control")
l = filter(shan_g, DPI == "d06" & diet == "Lean")
t.test(o$normalized_shannon,l$normalized_shannon)
```

# Does higher titer correlate with more diversity?
```{r}
titers = select(meta, sample, ferretID, DPI, titer, log10_titer)
gen_count_titers = merge(gen_count, titers, by = c("ferretID","DPI")) %>% unique()

Titer_SNV_Corr = ggplot(filter(gen_count_titers, log10_titer > 1), aes(x = log10_titer, y = snv_count)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(~inf_route) +
  PlotTheme1
print(Titer_SNV_Corr)
ggsave(Titer_SNV_Corr, file = "Titer_SNV_Corr.pdf", path = savedir)

test_index = filter(gen_count_titers, inf_route == "Index")
test_contact = filter(gen_count_titers, inf_route == "Contact")

m = lm(test_contact$snv_count ~ test_contact$log10_titer)
summary(m)
```

# Do higher Ct values correlate with more diversity?
```{r}
ct_values = select(meta, sample, ferretID, DPI, Ct_Mgene, log10_titer)
gen_count_ct = merge(gen_count, ct_values, by = c("ferretID","DPI")) %>%
  filter(!(is.na(Ct_Mgene))) %>%
  filter(inf_route == "Index" | inf_route == "Contact") %>% 
  unique()

Ct_SNV_Corr = ggplot(gen_count_ct, aes(x = Ct_Mgene, y = snv_count)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(~inf_route) +
  PlotTheme1
print(Ct_SNV_Corr)
ggsave(Ct_SNV_Corr, file = "Ct_SNV_Corr.pdf", path = savedir, width = 10, height = 5)

test_index = filter(gen_count_ct, inf_route == "Index")
test_contact = filter(gen_count_ct, inf_route == "Contact")

m = lm(test_index$snv_count ~ test_index$Ct_Mgene)
summary(m)
```

# dNdS analysis
```{r}
# by ferret
#dNdS_ferret = minorvdf %>% 
#  ungroup() %>% 
#  group_by(ferretID,DPI,diet,inf_route) %>% 
#  count(nonsyn)

#dNdS_ferret = pivot_wider(dNdS_ferret,names_from = nonsyn, values_from = n)
#dNdS_ferret = select(dNdS_ferret, ferretID,DPI,nonsyn,syn)
#dNdS_ferret$dNdS = paste0(dNdS_ferret$nonsyn / dNdS_ferret$syn)
#dNdS_ferret$dNdS = as.numeric(dNdS_ferret$dNdS)

#dNdS_ferret = filter(dNdS_ferret, inf_route == "Index" | inf_route == "Contact")

#dNdS_ferret_plot = ggplot(dNdS_ferret, aes(x = DPI, y = dNdS, color = ferretID)) +
#  geom_point() +
#  geom_line(aes(group = ferretID)) +
#  facet_grid(~diet+inf_route) +
#  PlotTheme1
#print(dNdS_ferret_plot)
#ggsave("dNdS_ferret.pdf", dNdS_ferret_plot, path = savedir)
#ggsave("dNdS_ferret.png", dNdS_ferret_plot, path = savedir, width = 10, height = 5)
```

# Setting up AF dataframes for bottleneck calculation: Stock to index ferrets
```{r}
af_df = select(vdf, sample, segment, ntpos, stocknt, major, minor,minorfreq, cohort, inf_route, ferretID, diet, DPI) %>%
  ungroup() %>%
  unique()
af_df$var = paste0(af_df$segment,"_",af_df$major,"_",af_df$ntpos,"_",af_df$minor)

stock = filter(af_df, inf_route == "Control") %>% unique()
index = filter(af_df, inf_route == "Index") %>% unique()

samples = unique(index$sample)
  
stock_files_dir = "~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01/freqfiles/stock"
data = data.frame()
  
for(i in samples){
  print(i)
  
  n = filter(af_df, sample == i)
  
  c = unique(n$cohort)
  s = filter(stock, cohort %in% c)
  
  comp = merge(s, n, by = c("segment", "ntpos","stocknt","major","minor","cohort","var"), all.x = TRUE)
  comp$minorfreq.y[is.na(comp$minorfreq.y)] = 0
  
  df = select(comp,minorfreq.x,minorfreq.y)
  write.table(df, file = glue("{stock_files_dir}/Stock_{i}_freqs.csv"), row.names = FALSE, col.names = FALSE)
  
  var_af = select(comp, var, minorfreq.x,minorfreq.y) %>% 
    mutate(sample = i)
  data = rbind(data,var_af)
}
```

# Setting up AF dataframes for bottleneck calculation: Index to contact ferrets
```{r}
af_df = select(vdf, sample, segment, ntpos, stocknt, major, minor,minorfreq, cohort, inf_route, ferretID, diet, DPI) %>%
  ungroup() %>%
  unique()
af_df$var = paste0(af_df$segment,"_",af_df$major,"_",af_df$ntpos,"_",af_df$minor)
df = merge(af_df, pairs, by = c("ferretID"))

index = filter(df, inf_route == "Index") %>% unique()
first_time = c("1409_d04","1410_d02","1415_d02","1794_d02","1797_d02",
               "1913_d06","1914_d06","1980_d02","1981_d10","1986_d10","2231_d04","2232_d02","2239_d04","2243_d04")
contact = filter(df, sample %in% first_time) %>% unique()

samples = unique(index$sample)
  
ferret_dir = "~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01/freqfiles/transmission"
tran_df = data.frame()
  
for(i in samples){
  print(i)
  
  n = filter(index, sample == i)
  partner = unique(n$pair_numbers)
  c = filter(contact, pair_numbers %in% partner)
  
  if(nrow(c) > 0){
    
    s = unique(c$sample)
  
    comp = merge(n, c, by = c("segment", "ntpos","stocknt","major","minor","cohort","var"), all.x = TRUE)
    comp$minorfreq.y[is.na(comp$minorfreq.y)] = 0
  
    df = select(comp,minorfreq.x,minorfreq.y)
    write.table(df, file = glue("{ferret_dir}/{i}_{s}_freqs.csv"), row.names = FALSE, col.names = FALSE)
    
    transmission_af = select(comp, var, minorfreq.x,minorfreq.y, pair_numbers.x)
    tran_df = rbind(tran_df,transmission_af)
   
    }else(print("No transmission"))
}

```

# What happens to the stock variants in the index(af)?
```{r}
lowest_timepoint = c("1408_d02","1414_d04","1416_d02","1789_d02","1800_d02","1801_d02","1910_d02","1912_d02",
                     "1968_d02","1970_d02","1972_d02","1973_d02","1974_d02","1975_d04","1977_d02","1984_d02",
                     "2233_d04","2234_d04","2240_d04","2251_d04","2253_d04","2254_d02")
data = data %>% 
  filter(sample %in% lowest_timepoint) %>% 
  group_by(sample)

data_lost = filter(data, minorfreq.y == 0) %>% unique() %>% mutate(dynamics = "lost")
data_kept = filter(data, minorfreq.y > 0) %>% unique() %>% mutate(dynamics = "kept")
data = rbind(data_lost, data_kept) %>% unique() %>% arrange(sample)

df = group_by(data, sample) %>% 
  count(dynamics) %>% 
  pivot_wider(names_from = dynamics, values_from = n) %>%
  mutate(totalvar = kept +lost) %>%
  mutate(percent_lost = lost/totalvar)

meta_noseg = select(meta, sample, ferretID, DPI, cohort, inf_route, diet, STRAIN, quality) %>% unique()
df1 = merge(df, meta_noseg, by = c("sample"))

stock_lost = ggplot(df1, aes(x = diet, y = percent_lost)) +
  geom_boxplot() +
  geom_jitter(aes(color = diet), width = 0.2) +
  PlotTheme1 +
  DietcolScale +
  ylim(0,1) +
  ylab("Percentage of stock variants lost") +
  xlab("Diet")
print(stock_lost)
ggsave(stock_lost, file = "stock_lost.pdf", path = savedir)

o = filter(df1, diet == "Obese")
l = filter(df1, diet == "Lean")
t.test(o$percent_lost, l$percent_lost)
mean(o$percent_lost)
mean(l$percent_lost)
mean(df1$percent_lost)
```

# What happens to the variants in the stock in the index ferrets? (af)
```{r}
kept = filter(data, dynamics == "kept")
keep_vars = unique(kept$var)

vdf_stocks = filter(vdf, inf_route == "Index")
vdf_stocks$ntvar = paste0(vdf_stocks$segment,"_",vdf_stocks$major,"_",vdf_stocks$ntpos,"_",vdf_stocks$minor)
vdf_stocks = filter(vdf_stocks, ntvar %in% keep_vars)

vdf_stocks_in = filter(vdf_stocks, inf_route == "Index") %>% 
  filter(DPI == "d02" | DPI == "d04" | DPI == "d06")

stockvar_index = ggplot(vdf_stocks_in, aes(x = DPI, y = minorfreq, color = diet)) +
  geom_point() +
  geom_line(aes(group = ferretID)) +
  facet_grid(~ntvar) +
  PlotTheme1 +
  DietcolScale
print(stockvar_index)
ggsave(stockvar_index, file = "stockvar_index.pdf",path = savedir, width = 15, height = 5)
```

# Loading in data after running bottleneck code
```{r}
stock_index = read.csv("~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01/freqfiles/stock/all_stock_bottlenecks.csv")
stock_index = merge(stock_index, meta_noseg, by = c("sample","ferretID","DPI"))

transmission  = read.csv("~/Desktop/GitHub/Obesity/NewExtractions/H9N2/timo_0.01/freqfiles/transmission/transmission_bottlenecks.csv")
transmission = merge(transmission, meta_noseg, 
                     by.x = c("index_sample","index_ferretID","index_DPI"), 
                     by.y = c("sample","ferretID","DPI"))
transmission = merge(transmission, meta_noseg, 
                     by.x = c("contact_sample","contact_ferretID","contact_DPI"),
                     by.y = c("sample","ferretID","DPI"))
transmission$pair_diet = paste0(transmission$diet.x," > ",transmission$diet.y)
```

```{r}
ggplot(filter(stock_index, smallest_timepoint == "yes"), aes(x = diet, y = bottleneck_size)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, aes(color = diet)) +
  PlotTheme1 +
  DietcolScale +
  ggtitle("Bottleneck between stock and index ferrets")

o_bottle = filter(stock_index, diet == "Obese" & smallest_timepoint == "yes")
l_bottle = filter(stock_index, diet == "Lean" & smallest_timepoint == "yes")
t.test(o_bottle$bottleneck_size, l_bottle$bottleneck_size)

ggplot(filter(transmission, smallest_timepoint == "yes"), aes(x = pair_diet, y = bottleneck_size)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, aes(color = pair_diet)) +
  PlotTheme1 +
  DietcolScale +
  ggtitle("Bottleneck between index and conatct pairs")

o_bottle = filter(stock_index, diet == "Obese" & smallest_timepoint == "yes")
l_bottle = filter(stock_index, diet == "Lean" & smallest_timepoint == "yes")
t.test(o_bottle$bottleneck_size, l_bottle$bottleneck_size)
```

# How many SNVs are transmitted?
```{r}
index_contact_snv = filter(tran_df, minorfreq.y > 0) %>% 
  mutate(transmitted = "yes")
index_contact_nosnv = filter(tran_df, minorfreq.y == 0) %>% 
  mutate(transmitted = "no")

tran_df_count = rbind(index_contact_snv,index_contact_nosnv) %>%
  group_by(pair_numbers.x,transmitted) %>%
  tally() %>%
  pivot_wider(names_from = transmitted, values_from = n) %>%
  mutate(perc_transmitted = yes / no + yes) %>% 
  unique()

meta_small = select(meta, ferretID,cohort,inf_route,diet,STRAIN)
meta_small_pairs = merge(meta_small, pairs, by = c("ferretID"))

tran_df_count = merge(tran_df_count, meta_small_pairs, by.x = c("pair_numbers.x"), by.y = c("pair_numbers")) %>% 
  unique() %>%
  filter(inf_route == "Contact")
tran_df_count$perc_transmitted[is.na(tran_df_count$perc_transmitted)] = 0
```

```{r}
transmitted_snvs_diets = ggplot(tran_df_count, aes(x = pair_diets, y = perc_transmitted)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, aes(color = diet)) +
  PlotTheme1 +
  DietcolScale
print(transmitted_snvs_diets)
ggsave(transmitted_snvs_diets, file = "transmitted_snvs_diets.pdf", path = savedir)

transmitted_snvs_all = ggplot(tran_df_count, aes(x = STRAIN, y = perc_transmitted)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2) +
  PlotTheme1 +
  DietcolScale
print(transmitted_snvs_all)
ggsave(transmitted_snvs_all, file = "transmitted_snvs_all.pdf", path = savedir, height = 5, width = 5)
```

# Turnover rate
```{r}
multiday = ungroup(minorvdf) %>%
  filter(inf_route == "Index") %>%
  select(ferretID, DPI) %>% 
  unique() %>%
  group_by(ferretID) %>%
  tally() %>%
  filter(n > 1)

select_samps = unique(multiday$ferretID)

turnover = data.frame()
turnover2 = data.frame()
  
for(i in select_samps){
  
  print(i)
  df = filter(minorvdf, ferretID == i)
  days = count(df,DPI)
  
  df1 = df %>% group_by(ferretID,var) %>% tally()
  df1$days = df1$n
  
  df2 = df1 %>% group_by(days) %>% tally()
  df2$perc = (df2$n / nrow(df1))
  df2$ferretID = i
  
  turnover = rbind(turnover, df2)
  
  df_1day = filter(df2, days == "1")
  df_over1day = filter(df2, days != "1") %>%
    mutate(total_n = sum(n)) %>%
    select(ferretID,total_n) %>% unique() %>%
    mutate(days = "> 1") %>% 
    mutate(perc = total_n / nrow(df1))
  df_over1day$n = df_over1day$total_n 
  df_over1day = select(df_over1day,!(total_n))
  
  df3 = rbind(df_1day,df_over1day)
  turnover2 = rbind(turnover2,df3)
}
```

```{r}
turns = merge(turnover, meta_small, by = c("ferretID")) %>% unique
turns$days = as.character(turns$days)

ggplot(turns, aes(x = ferretID, y = perc, fill = factor(days, levels = c("4","3","2","1")))) +
  geom_col() +
  facet_grid(~diet) +
  PlotTheme1

turns2 = merge(turnover2, meta_small, by = c("ferretID")) %>% unique
turns2$days = as.character(turns2$days)

proportion_maintained = ggplot(turns2, aes(x = ferretID, y = perc, fill = days)) +
  geom_col() +
  facet_grid(~diet) +
  PlotTheme1
print(proportion_maintained)
ggsave(proportion_maintained, file = "proportion_maintained_index.pdf", path = savedir)


o = filter(turns2, diet == "Lean", days == "> 1")
l = filter(turns2, diet == "Obese", days == "> 1")
t.test(o$n, l$n)
# not significant counting each day separately or just by 1 vs >1
```