HostDepletionKitComparison.Rmd

---
title: "Host Depletion Kit Comparisons"
author: "ErikaGanda"
output:
  html_document: default
  pdf_document: default
  word_document: default
editor_options: 
  chunk_output_type: console
---
```{r}
#Data File: OpenRefineHostDepletionData1.23.19
library(ggplot2)
library(dplyr)
library(emmeans)
library(multcompView)
library(tidyverse)
library(multcomp)
library(car)

#Read Data
AllData <- read.table("Inputs/OpenRefineHostDepletionData1.23.19.txt", sep="\t",  fill = TRUE, header=TRUE)

AllData$Category[AllData$Category  == "NTC"] <- "UTC"

SampleData <- AllData %>% filter(Category!="Standard") %>% filter(Category!="NEBdiluent")
dim (SampleData)

#Getting All variables and all levels within each

#sapply(SampleData, levels)
#SampleData %>% 
#     sapply(levels)

#Summary Statistics
SampleData.summary <- SampleData %>%
  group_by(Target.Name,Category,BioRep) %>%
  summarize(mean_CopyN_permLMilk=mean(LogCopiespermLofMilk,na.rm=T),
            st_dev=sd(LogCopiespermLofMilk,na.rm=T),
            n_missing=sum(is.na(LogCopiespermLofMilk)),
            n_total=n())
# write.table (SampleData.summary, "SampleData.summary.txt", sep="\t" ) 


#Summary Statistics by Kit
SampleData.summary <- SampleData %>%
  group_by(Target.Name,Category) %>%
  summarize(mean_CopyN_permLMilk=mean(LogCopiespermLofMilk,na.rm=T),
            st_dev=sd(LogCopiespermLofMilk,na.rm=T),
            n_missing=sum(is.na(LogCopiespermLofMilk)),
            n_total=n())
# write.table (SampleData.summary, "SampleData.summary.txt", sep="\t" ) 
```


### Total Bovine Copy Numbers
```{r}

#Filter Subset from Sample Data
Bovine <- SampleData %>% filter(Target.Name=="BovineUXT")
dim(Bovine)

#Summary Statistics
Bovine.summary <- Bovine %>%
  group_by(Category) %>%
  summarize(mean_CopyN_permLMilk=mean(LogCopiespermLofMilk,na.rm=T),
            st_dev=sd(LogCopiespermLofMilk,na.rm=T),
            n_missing=sum(is.na(LogCopiespermLofMilk)),
            n_total=n())
# write.table (Bovine.summary, "Bovine.summary.txt", sep="\t" )
```


```{r}
#Plot Raw Means and Standard Deviations
Colors <- c("20uM" = "#d7191c", "Molz"= "#a1d99b", "NEB"= "#feb24c","UTC"= "#595959", "Qia"="#2c7bb6")
ggplot(data=Bovine,
       mapping=aes(x=Category,y=LogCopiespermLofMilk, color=Category, ylab="Copy Numbers")) +
  ylab ("Log10 Copies / mL of Milk")+
     geom_boxplot(lwd=1)+
  theme_bw()+
 ggtitle("Bovine DNA Copy Numbers")+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
  scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))

ggplot(data=Bovine,
       mapping=aes(x=Category,y=LogCopiespermLofMilk, color=Category, shape=BioRep)) +
  ylab ("Log10 Copies / mL of Milk")+
    geom_point(position=position_jitterdodge(jitter.width=0, dodge.width = .6),size = 2,stroke = .7)+
  ggtitle("Bovine DNA Copy Numbers")+
  scale_shape_discrete(solid=F) +
  theme_bw()+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
  scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))

```

```{r}
# ERIKA: I guess this is needed
Bovine$Category <- ifelse(Bovine$Category == "20uM", "TwentyuM", Bovine$Category)
Bovine <- mutate(Bovine, Category = factor(Category))

# Statistical Analysis
m_Bovine.LogCopiespermLofMilk <- lm( LogCopiespermLofMilk ~ Category*BioRep -1, data=Bovine ) ## ERIKA: I am guessing you want to
                                                                                              ## remove the intercept (-1)
summary(m_Bovine.LogCopiespermLofMilk)

AIC(m_Bovine.LogCopiespermLofMilk) 

# Model Fit Plots

plot(x=predict(m_Bovine.LogCopiespermLofMilk),y=resid(m_Bovine.LogCopiespermLofMilk))
  
  # using ggplot2   
Colors <- c("TwentyuM" = "#d7191c", "Molz"= "#a1d99b", "NEB"= "#feb24c","UTC"= "#595959", "Qia"="#2c7bb6")
ggplot(m_Bovine.LogCopiespermLofMilk, aes(x=predict(m_Bovine.LogCopiespermLofMilk), y=resid(m_Bovine.LogCopiespermLofMilk), color=Category)) +
  geom_point()+
  theme_bw()+
  ggtitle("Bovine - Model Fit - Residuals vs Predicted")+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  geom_hline(yintercept = 1) +
  geom_hline(yintercept = -1)

  # qqplots
qqnorm(resid(m_Bovine.LogCopiespermLofMilk))
qqline(resid(m_Bovine.LogCopiespermLofMilk))

summary(m_Bovine.LogCopiespermLofMilk)

#All residuals are < 1. 
```

```{r}
## JDS Help
## Create vector of pairwise comparisons to test 
comparisons <- unique(Bovine$Category)
comparisons <- comparisons[which(comparisons != "UTC")]
comparisons <- t(combn(comparisons, 2))
comparisons <- paste(comparisons[,1], comparisons[,2], sep=" - ")
comparisons <- paste0(comparisons, " == 0")

## The multcomp package needs some updating
## This is all to deal with a poorly designed interface and internal logic in multcomp
tmp_vcov <- vcov(m_Bovine.LogCopiespermLofMilk)
tmp_vcov[is.na(tmp_vcov)] <- 1e-8 # this value doesn't matter
tmp_coef <- coef(m_Bovine.LogCopiespermLofMilk)
tmp_coef[is.na(tmp_coef)] <- 1e-8 # this value doesn't matter
comps <- glht(m_Bovine.LogCopiespermLofMilk, mcp(Category = comparisons),
              vcov.=tmp_vcov, coef.=tmp_coef)

## Now summarise results - using manual p-value correction
comps <- summary(comps, test=univariate())
comps$test$pvalues <- p.adjust(comps$test$pvalues, method="bonferroni")

## Now view results
print(comps)
```


## Total Bacterial DNA Copy Numbers


```{r}

#Filter Subset from Sample Data
TotalBacterialDNA <- SampleData %>% filter(Target.Name=="Total16S")
dim(TotalBacterialDNA)

#Summary Statistics
TotalBacterialDNA.summary <- TotalBacterialDNA %>%
  group_by(Category) %>%
  summarize(mean_CopyN_permLMilk=mean(LogCopiespermLofMilk,na.rm=T),
            st_dev=sd(LogCopiespermLofMilk,na.rm=T),
            n_missing=sum(is.na(LogCopiespermLofMilk)),
            n_total=n())
# write.table (TotalBacterialDNA.summary, "TotalBacterialDNA.summary.txt", sep="\t" )
```


```{r}
#Plot Raw Means and Standard Deviations
Colors <- c("20uM" = "#d7191c", "Molz"= "#a1d99b", "NEB"= "#feb24c","UTC"= "#595959", "Qia"="#2c7bb6")
ggplot(data=TotalBacterialDNA,
       mapping=aes(x=Category,y=LogCopiespermLofMilk, color=Category, ylab="Copy Numbers")) +
  ylab ("Log10 Copies / mL of Milk")+
     geom_boxplot(lwd=1)+
  theme_bw()+
 ggtitle("TotalBacterialDNA DNA Copy Numbers")+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
  scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))

ggplot(data=TotalBacterialDNA,
       mapping=aes(x=Category,y=LogCopiespermLofMilk, color=Category, shape=BioRep)) +
  ylab ("Log10 Copies / mL of Milk")+
    geom_jitter(width=0.25)+
  scale_shape_discrete(solid=F) +
    ggtitle("TotalBacterialDNA DNA Copy Numbers")+
  theme_bw()+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
  scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))

```

```{r}
# Statistical Analysis
m_TotalBacterialDNA.LogCopiespermLofMilk <- lm( LogCopiespermLofMilk ~ Category*BioRep, data=TotalBacterialDNA )

summary(m_TotalBacterialDNA.LogCopiespermLofMilk)

AIC (m_TotalBacterialDNA.LogCopiespermLofMilk)

# Model Fit Plots

plot(x=predict(m_TotalBacterialDNA.LogCopiespermLofMilk),y=resid(m_TotalBacterialDNA.LogCopiespermLofMilk))
  
  # using ggplot2
Colors <- c("20uM" = "#d7191c", "Molz"= "#a1d99b", "NEB"= "#feb24c","UTC"= "#595959", "Qia"="#2c7bb6")
ggplot(m_TotalBacterialDNA.LogCopiespermLofMilk, aes(x=predict(m_TotalBacterialDNA.LogCopiespermLofMilk), y=resid(m_TotalBacterialDNA.LogCopiespermLofMilk), color=Category)) +
  geom_point()+
  theme_bw()+
  ggtitle("TotalBacterialDNA - Model Fit - Residuals vs Predicted")+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  geom_hline(yintercept = 1) +
  geom_hline(yintercept = -1)

  # qqplots
qqnorm(resid(m_TotalBacterialDNA.LogCopiespermLofMilk))
qqline(resid(m_TotalBacterialDNA.LogCopiespermLofMilk))

summary(m_TotalBacterialDNA.LogCopiespermLofMilk)

#All residuals are < 1. 
```


```{r}
# Check Tukey-adjusted pairwise comparison of kit estimates
# JDS help!

```

### Manuscript Figures
```{r}
ggplot() +
  geom_point(data=TotalBacterialDNA,aes(x=Category,y=LogCopiespermLofMilk,color=Category, shape=BioRep), position=position_jitterdodge(jitter.width=0, dodge.width = .5),show.legend = F) +
  scale_color_manual(values=Colors)+
#  geom_point(data=df2_TotalBacterialDNA.model,aes(x=Category,y=emmean,fill=Category),show.legend = F)+
#  geom_errorbar(data=df2_TotalBacterialDNA.model,aes(x=Category,ymin=lower.CL,ymax=upper.CL),width=0.5)+  
#  geom_text(data =data.frame(m_TotalBacterialDNA.LogCopiespermLofMilk_cld),aes(x=Category,label=`.group`, y=emmean), nudge_y = .7, nudge_x = -0.0, fontface = "bold") +
ylim(1.7, 7)+
xlab("Kit")+
ylab ("Log10 Copies / mL of Milk")+
scale_shape(solid = F)+
  theme_bw()+
  ggtitle("Total Bacterial DNA Copy Numbers")+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))

#ggsave("PanelA-BovBactKitComparison.TIFF", width = 3, height = 2.5 , units = "in", dpi = 600)

ggplot() +
  geom_point(data=Bovine,aes(x=Category,y=LogCopiespermLofMilk,color=Category, shape=BioRep), position=position_jitterdodge(jitter.width=0, dodge.width = .6), show.legend = F) +
  scale_color_manual(values=Colors)+
  #geom_point(data=df2_Bovine.model,aes(x=Category,y=emmean,fill=Category), show.legend = F)+
  #geom_errorbar(data=df2_Bovine.model,aes(x=Category,ymin=lower.CL,ymax=upper.CL),width=0.5)+  
  #geom_text(data =data.frame(m_Bovine.LogCopiespermLofMilk_cld),aes(x=Category,label=`.group`, y=emmean),nudge_y = .7,nudge_x = -0.0, fontface = "bold") +
ylim(1.7, 7)+
xlab("Kit")+
ylab ("Log10 Copies / mL of Milk")+
  theme_bw()+
  ggtitle("Bovine Copy Numbers")+
  scale_shape(solid = F)+
  scale_color_manual(values=Colors)+
  theme(panel.grid.major = element_blank())+
  theme(panel.grid.minor = element_blank())+
  theme(panel.background = element_blank())+
  theme(panel.border = element_blank())+
  theme(axis.line = element_line())+
  theme(axis.text.x=element_text(angle=25,vjust=0.5))+
scale_x_discrete(limits=c("UTC", "NEB", "20uM", "Qia", "Molz"))
#ggsave("PanelB-BovBactKitComparison.TIFF", width = 3, height = 2.5 , units = "in", dpi = 600)


ggplot()+
geom_boxplot(data=SampleData,aes(x=Cat.Target, y=LogCopiespermLofMilk, fill=Target.Name), show.legend = F)+
 scale_fill_manual(values = c("lightgray", "white"))+
  geom_point(data=SampleData,aes(x=Cat.Target,y=LogCopiespermLofMilk,color=Category, shape=BioRep), size = 2.5, stroke = .5, position=position_jitterdodge(jitter.width=0.05, dodge.width = 1), show.legend = F)+
  scale_shape(solid = F)+
  theme_bw()+
  aes(colour=factor(Category))+
  scale_colour_manual(values = Colors)+
  theme(legend.title=element_blank())+
  xlab(NULL)+
  ylab(NULL)+
  theme(panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  panel.border = element_blank(),
  axis.line = element_line(),
  axis.text.x=element_text(angle=90,vjust=0.5))+
  #ylim(1.5, 7.5)+
  scale_x_discrete(limits=c("UTC.Bacteria","UTC.Bovine"," "," ","NEB.Bacteria","NEB.Bovine"," "," ","PMA.Bacteria","PMA.Bovine"," "," ","Qia.Bacteria","Qia.Bovine"," "," ","Molz.Bacteria","Molz.Bovine"," "," "))

#ggsave("PanelC-BovBactKitComparison.TIFF", width = 7.5, height = 3.5 , units = "in", dpi = 600)

```


```{r}
sessionInfo()
```