---
title: "PhyloSeq/DeSeq2 Ovex and Sham analysis"
author: "Christine Olmstead"
date: "2/8/2020"
output: html_document
---

```{r setup}
    library(dplyr)
    library(ggplot2)
    library(vegan)
    library(stringr)
    library(svglite)
    library(extrafont)
    font_import()
    loadfonts(device = "win")
    library(DESeq2); packageVersion("DESeq2")
    library(phyloseq); packageVersion("phyloseq")
    library(biomformat); packageVersion("biomformat")
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown



Beta Diversity analysis (unweighetd and weighted UniFrac)
```{r}

otu <- read.csv("sv.sham.ovex.csv")
tax <- read.csv("tax.sham.ovex.csv")
combined_data <-left_join(otu, tax, by = NULL, copy = FALSE)

write.csv(combined_data,"merged.sv.sham.ovex.csv", row.names = FALSE)

taxa <- import_biom("table.sham.ovex.biom", treefilename = "tree.nwk")
map.file <-import_qiime_sample_data("ovex_mapping_sham_ovex_sample_removed.tsv")
OVX <- merge_phyloseq(taxa, map.file)
colnames(tax_table(OVX)) <- c("kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species" )

#subset data into treatment group
sham <-  subset_samples(OVX, Treatment == "Sham")  
ovex <-  subset_samples(OVX, Treatment == "Ovex")

#Plot sham and ovex abundance  by Phylum
plot_bar(sham, fill = "Phylum", title = "Sham Abundance")
plot_bar(ovex, fill = "Phylum", title = "Ovex Abundance")


#subset sample to include only T1-T5 (post-ovx)
OVX.post <-  subset_samples(OVX, time != "0")

#-----------------------unweighted uniFrac Beta diversity------------------------------

# PCoA plot using the unweighted UniFrac as distance
uw.unifrac.dist = phyloseq::distance(OVX.post, method="unifrac", weighted=F)
uw.ordination = ordinate(OVX.post, method="PCoA", distance=uw.unifrac.dist)

#PCoA plot to visualize unweighted uniFrac distance matrix
uwUniFrac <- plot_ordination(OVX.post, uw.ordination, color="Treatment", shape ="Treatment") + theme_bw() + theme(text = element_text(size=12, family = "Arial"),panel.grid.major = element_blank(),panel.grid.minor = element_blank(), legend.justification=c(0,0), legend.position=c(0,0), legend.background = element_blank(), legend.text = element_text(size = 12), legend.key.width = unit(0.55, "lines"), legend.key.height = unit(0.55, "lines"), legend.title=element_text(size=0),
        legend.box.background = element_rect(colour = "black"), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.margin=unit(c(1,5,1,5),"cm"),  plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) + scale_color_manual(values = c("Ovex" = "#ff6600", "Sham" = "#006600")) + geom_point(size = 2) + geom_vline(xintercept=c(0,0), linetype="dotted") + geom_hline(yintercept=c(0,0), linetype="dotted")  + labs(title = "PCoA - PC1 vs PC2 (unweighted)", subtitle = "Post-Treatment (Week 1-5)", x = "PC1- Percent variation explained 17%", y = "PC2- Percent variation explained 9%") +   annotate("label", x=-.195, y=.15, label = str_wrap("ADONIS p=0.036", 5), label.r = unit(0, 'lines'))

#PERMANOVA 
adonis(uw.unifrac.dist ~ sample_data(OVX.post)$Treatment)
 
ggsave(file="PCoA_unweighted.svg", plot=uwUniFrac, width=8, height=6)


#-----------------------weighted uniFrac Beta diversity------------------------------

# PCoA plot using the weighted UniFrac as distance
w.unifrac.dist = phyloseq::distance(OVX.post, method="unifrac", weighted=T)
w.ordination = ordinate(OVX.post, method="PCoA", distance=w.unifrac.dist)


#PCoA plot to visualize unweighted uniFrac distance matrix
uwUniFrac <- plot_ordination(OVX.post, w.ordination, color="Treatment", shape ="Treatment") + theme_bw() + theme(text = element_text(size=12, family = "Arial"),panel.grid.major = element_blank(),panel.grid.minor = element_blank(), legend.justification=c(0,0), legend.position=c(0,0), legend.background = element_blank(), legend.text = element_text(size = 12), legend.key.width = unit(0.55, "lines"), legend.key.height = unit(0.55, "lines"), legend.title=element_text(size=0),
        legend.box.background = element_rect(colour = "black"), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.margin=unit(c(1,5,1,5),"cm"),  plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) + scale_color_manual(values = c("Ovex" = "#ff6600", "Sham" = "#006600")) + geom_point(size = 2) + geom_vline(xintercept=c(0,0), linetype="dotted") + geom_hline(yintercept=c(0,0), linetype="dotted")  + labs(title = "PCoA - PC1 vs PC2 (weighted)", subtitle = "Post-Treatment (Week 1-5)", x = "PC1- Percent variation explained 89%", y = "PC2- Percent variation explained 4%") +   annotate("label", x=-.176, y=.05, label = str_wrap("ADONIS p=0.179", 5), label.r = unit(0, 'lines'))

#PERMANOVA
adonis(w.unifrac.dist ~ sample_data(OVX.post)$Treatment)

ggsave(file="PCoA_weighted.svg", plot=uwUniFrac, width=8, height=6)


```
Beta Diversity at Time 0

```{r}

#subset data to include only time 0
OVX.T0 <-  subset_samples(OVX, time == "0")

#----------------------------unweighed-------------------------------

# Calculate the unweighted uniFrac distance matrix
uw.unifrac.dist.T0 = phyloseq::distance(OVX.T0, method="unifrac", weighted=F)
uw.ordination.T0 = ordinate(OVX.T0, method="PCoA", distance=uw.unifrac.dist.T0)

# PCoA plot using the unweighted uniFrac as distance matrix
plot_ordination(OVX.T0, uw.ordination.T0 , color="Treatment", title = "Unweighted UniFrac Time 0") + theme(aspect.ratio=1)

#PERMANOVA
adonis(uw.unifrac.dist.T0 ~ sample_data(OVX.T0)$Treatment)

#------------------------------weighted----------------------

#Calculate the weighted uniFrac distance matrix
w.unifrac.dist.T0 = phyloseq::distance(OVX.T0, method="unifrac", weighted=T)
w.ordination.T0 = ordinate(OVX.T0, method="PCoA", distance=w.unifrac.dist.T0)

# PCoA plot using the weighted uniFrac as distance matrix
plot_ordination(OVX.T0, w.ordination.T0, color="Treatment", title = "Weighted UniFrac Time 0") + theme(aspect.ratio=1)

#PERMANOVA
adonis(w.unifrac.dist.T0 ~ sample_data(OVX.T0)$Treatment)

```

l
