---
title: "deseq_sham_ovex"
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)
```

## Including Plots
Trimmed data for log2fold analysis time 4

```{r}
#import mapping file, phylogenetic tree, taxa 
biom <- import_biom("sham.ovex.trimmed.biom", treefilename = "tree.trimmed.nwk")
map <-import_qiime_sample_data("ovex_mapping_sham_ovex_sample_removed.tsv")

#format data for phyloseq
phylo <- merge_phyloseq(biom, map)
colnames(tax_table(phylo)) <- c("kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species" )

#convert time and treatement to factor type
sample_data(phylo)$Treatment <- as.factor(sample_data(phylo)$Treatment)
sample_data(phylo)$time <- as.factor(sample_data(phylo)$time)

#only include samples from time 5
phylo.week4 <- subset_samples(phylo, time =="4")

#convert phyloseq to deseq2
dSeq2 <- phyloseq_to_deseq2(phylo.week4, ~Treatment)

#geometric mean 
gm.mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(dSeq2), 1, gm.mean)
dSeq2= estimateSizeFactors(dSeq2, geoMeans = geoMeans)
dSeq2= DESeq(dSeq2, fitType="local")

alpha <- 0.05
log2.results <- results(dSeq2, contrast=c("Treatment", "Ovex", "Sham"), alpha=alpha)
log2.results <- log2.results[order(log2.results$padj, na.last = NA),]

#FDR p adjusted less than p = 0.05
log2.results.sig <- log2.results[(log2.results$padj < alpha), ]

log2.results.sig <- cbind(as(log2.results.sig, "data.frame"), as(tax_table(phylo)[rownames(log2.results.sig), ], "matrix"))

#only include taxa where log2 fold change is greater than +-2
log2.week4 <-  subset(log2.results.sig, log2FoldChange >=2 | log2FoldChange <=-2)

#create new df
feature.names <- c("Lachnospiraceae", "Turicibacter","S24-7", "Bacteroides", "Lachnospiraceae",  "RF32")
features <- data.frame(log2.week4[,2], feature.names)
names(features) <- c("log2FoldChange", "Feature")

#plot log2 fold change for ovex vs sham at week 4
week4 <- ggplot(features, aes(x = factor(Feature, levels = c( "Lachnospiraceae",  "S24-7", "Bacteroides","RF32","Turicibacter") ) ,y=log2FoldChange)) +
  geom_bar(position = position_dodge2(preserve = "single") ,stat="identity", width=2.2, colour = "black", fill = "black") +
  coord_flip() + theme_bw() +
theme(panel.border = element_blank(), axis.line.x = element_line(color="black", size = 0.7),axis.text.x = element_text(size=14, family="Arial"),axis.text.y = element_text(size=14, vjust = -1.7, family="Arial"), axis.ticks.y=element_blank(),text = element_text(size=14, family="Arial"), legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(), plot.margin=unit(c(1,1,1,1),"cm"),  plot.title = element_text(hjust = 0.2, vjust = 2)) + labs(title = "Ovex vs. Sham (Week 4)", y = "Log2 Fold Change", x = "") +   geom_hline(aes(yintercept = 0), size = 0.7)   +  ylim(-25, 25) +

annotate("text", x=3.5, y=-3.7, label =  "-2.0", size = 4.5) +
annotate("text", x=2.5, y=-4.0, label =  "-2.3", size = 4.5) +
annotate("text", x=4.7, y=5.3, label =  "3.7", size = 4.5) +
annotate("text", x=5.8, y=9.0, label =  "7.4", size = 4.5) +
annotate("text", x=1.4, y=8.9, label =  "7.3", size = 4.5) +
annotate("text", x=.2, y=-24.6, label =  "-22.3", size = 4.5) 

week4


ggsave(file="Log2_week4.svg", plot=week4, width=8, height=6)





```


