AustriaCadaverOverview.rmd

---
title: "AustriaCadaverOverview"
author: "JReceveur"
output: 
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    code_folding: hide
---
  
  ```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width=14, fig.height=10)
knitr::opts_chunk$set(fig.align="center")
```
36 Samples remained after filtering and rarefaction (4,000 reads). Case A and B refer to the different donors A male, B female. 
#Import
``` {r import, message=FALSE, warning=FALSE,echo=TRUE}
set.seed(3246)
#Parse Silva function
parse_taxonomy_silva_128 <- function(char.vec){
  # Use default to assign names to elements in case problem with greengenes prefix
  char.vec = parse_taxonomy_default(char.vec)
  # Check for unassigned taxa
  if (char.vec["Rank1"] == "Unassigned") {
    char.vec <- c(Rank1="D_0__Unassigned", Rank2="D_1__Unassigned", Rank3="D_2__Unassigned", Rank4="D_3__Unassigned",
                  Rank5="D_4__Unassigned", Rank6="D_5__Unassigned", Rank7="D_6__Unassigned")
  }
  # Define the meaning of each prefix according to SILVA taxonomy
  Tranks = c(D_0="Kingdom", D_1="Phylum", D_2="Class", D_3="Order", D_4="Family", D_5="Genus", D_6="Species")
  # Check for prefix using regexp, warn if there were none. trim indices, ti
  ti = grep("[[:alpha:]]\\_[[:digit:]]{1}\\_\\_", char.vec)
  if( length(ti) == 0L ){
    warning(
      "No silva prefixes were found. \n",
      "Consider using parse_taxonomy_delfault() instead if true for all OTUs. \n",
      "Dummy ranks may be included among taxonomic ranks now."
    )
    # Will want to return without further modifying char.vec
    taxvec = char.vec
    # Replace names of taxvec according to prefix, if any present...
  } else {
    # Format character vectors for Ambiguous taxa
    if( length(ti) < 7 ){
      for (key in names(char.vec)) {
        if ( char.vec[key] == "Ambiguous_taxa" ) {
          tax_no <- (as.numeric(substr(key, 5, 5)) - 1)
          char.vec[key] = sprintf("D_%s__Ambiguous_taxa", tax_no)
        }
      }
      # Reset the trimmed indicies if Ambiguous taxa
      ti = grep("[[:alpha:]]\\_[[:digit:]]{1}\\_\\_", char.vec)
    }
    # Remove prefix using sub-"" regexp, call result taxvec
    taxvec = gsub("[[:alpha:]]\\_[[:digit:]]{1}\\_\\_", "", char.vec)
    # Define the ranks that will be replaced
    repranks = Tranks[substr(char.vec[ti], 1, 3)]
    # Replace, being sure to avoid prefixes notK present in Tranks
    names(taxvec)[ti[!is.na(repranks)]] = repranks[!is.na(repranks)]
  }
  return(taxvec)
}

library(vegan)
library(MASS)
library(ggplot2)
library(plyr)
library(dplyr)
library(magrittr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(randomForest)
library(knitr)
library(ape)
library(ggpubr)
library(viridis)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#000000","#CC79A7")
theme_set(theme_bw(base_size = 18)+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))
biom=import_biom("AustriaCadaverJan2020.biom",parseFunction=parse_taxonomy_silva_128)
#biom
#tax_table(biom)
#tax_table(biom) <- tax_table(biom)[,-c(5:10,14)]#remove dummy ranks

metadata=read.table("AustriaBodyFarmMetadataWDivJan2020.tsv",header = TRUE)

metadata$Timepoint<-factor(metadata$Timepoint, levels = c("8/4","9/1","11/17"))
metadata$Timepoint2<-factor(metadata$Timepoint2, levels = c("Aug","Sep","Nov"))


ddply(metadata, c("Case","SampleType"), summarise, N    = length(id))


tree=read_tree("AustriaCadaverTree.nwk")

sampdat=sample_data(metadata)
sample_names(sampdat)=metadata$id
physeq=merge_phyloseq(biom,sampdat,tree)
physeq

```

#Alpha Diversity{.tabset}
Samples from the Mouth Nose Rectum and Ears are show with the rest of the body sites as well as seperated out since they were the only sites sampled at multiple timepoints.

##All Sample Types
###Observed Species
```{r}
plot_richness(physeq, x="Timepoint2",color="Case",shape="Case", measures=c("Observed"))+ylab("Observed Species")+facet_wrap(~SampleType)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)
```

###Inv Simpson
```{r}
plot_richness(physeq, x="Timepoint2",color="Case",shape="Case", measures=c("InvSimpson"))+ylab("Inverse Simpson")+facet_wrap(~SampleType)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)


print("Inv Simpson without soil")
physeqWOSoil<-subset_samples(physeq,SampleType!="soil")
plot_richness(physeqWOSoil, x="Timepoint2",color="Case",shape="Case", measures=c("InvSimpson"))+ylab("Inverse Simpson")+facet_wrap(~SampleType)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)

```

###Faith's PD
```{r}
Faiths<-ggplot(metadata, aes(x=Timepoint,y=faith_pd,shape=Case,group=Case,color=Case))+geom_point(size=5)+geom_line()+
  xlab("Sample Location")+ylab("Faith's Phylogenetic Diversity") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+
  theme(axis.title.x=element_blank(),axis.text.x = element_text(size = 7))+facet_wrap(~ SampleType)+  scale_fill_manual(values=cbPalette)
Faiths

```

###Shannon Diversity
```{r}
Shannon<-ggplot(metadata, aes(x=Timepoint,y=shannon,shape=Case,group=Case,color=Case))+geom_point(size=5)+geom_line()+
  facet_grid(~Phylum)+xlab("Sample Location")+ylab("Shannon Diversity") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+
  theme(axis.title.x=element_blank(),axis.text.x = element_text(size = 7))+facet_wrap(~ SampleType)+  scale_fill_manual(values=cbPalette)
Shannon

```



##Mouth Nose Rectum and Ears only (Sites collected at multiple timepoints)

###Observed Species
```{r}
PhyseqSubset<-subset_samples(physeq,SampleType=="mouth"|SampleType=="nose"|SampleType=="rectum"|SampleType=="ears")
plot_richness(PhyseqSubset, x="Timepoint2",color="Case",shape="Case", measures=c("Observed"))+ylab("Observed Species")+facet_wrap(~SampleType)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)
```

###Inv Simpson
```{r}
plot_richness(PhyseqSubset, x="Timepoint2",color="Case",shape="Case", measures=c("InvSimpson"))+ylab("Inverse Simpson")+facet_wrap(~SampleType)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)

```

###Faith's PD
```{r}
Subset<-subset(metadata,SampleType=="mouth"|SampleType=="nose"|SampleType=="rectum"|SampleType=="ears")
Faiths<-ggplot(Subset, aes(x=Timepoint,y=faith_pd,shape=Case,group=Case,color=Case))+geom_point(size=5)+geom_line()+
  facet_grid(~Phylum)+xlab("Sample Location")+ylab("Faith's Phylogenetic Diversity") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+
  theme(axis.title.x=element_blank(),axis.text.x = element_text(size = 7))+facet_wrap(~ SampleType)+  scale_fill_manual(values=cbPalette)+ theme(legend.position = "none")+scale_x_discrete(labels=c("8/4" = "Aug 4","9/1"= "Sep 1","11/17"="Nov 17"))
Faiths
```



###Shannon Diversity
```{r}

Subset<-subset(metadata,SampleType=="mouth"|SampleType=="nose"|SampleType=="rectum"|SampleType=="ears")
Shannon<-ggplot(Subset, aes(x=Timepoint,y=shannon,shape=Case,group=Case,color=Case))+geom_point(size=5)+geom_line()+
  facet_grid(~Phylum)+xlab("Sample Location")+ylab("Shannon Diversity") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+
  theme(axis.title.x=element_blank(),axis.text.x = element_text(size = 15))+facet_wrap(~SampleType)+  scale_fill_manual(values=cbPalette)+scale_x_discrete(labels=c("8/4" = "Aug 4","9/1"= "Sep 1","11/17"="Nov 17"))
Shannon

```


#Taxa Plots

```{r filteringForFamily2}


GPr  = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPrPhylum=tax_glom(GPr, "Phylum")
PhylumLevel = filter_taxa(GPrPhylum, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower tha 0.1%
GPrFamily=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(GPrFamily, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%
GPrGenus=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GPrGenus, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%

```

##Total Abundances
Relative Abundances of All bacterial phyla (35 present total)
```{r TotalABUs,warning=F}
GPr  = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
PhylumAll=tax_glom(GPr, "Phylum")

PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 3e-3, TRUE) #filter out any taxa lower tha 1%
FamilyAll=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(FamilyAll, function(x) mean(x) > 3e-2, TRUE) #filter out any taxa lower tha 3%
GenusAll=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GenusAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower tha 1%

df <- psmelt(PhylumAll)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
Trtdata
print("Family level abundances across all samples (> 3% of total relative abundances")
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
Trtdata


```


##Phylum Level





```{r TreatmentPlot5}
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum","SampleType","Case","Timepoint2"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance)
)
#Trtdata


PhylumPlot=ggplot(Trtdata, aes(x=SampleType,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+
  facet_grid(~Phylum)+xlab("Sample Location")+ylab("Relative Abundance (%, > 0.3%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(axis.title.x=element_blank())+facet_grid(Case~Timepoint2)+  scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
PhylumPlot

```


##Family Level Relative Abundance




```{r TreatmentPlot}
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100


Trtdata <- ddply(df, c("Family","SampleType","Case","Timepoint2"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance)
)
#Trtdata
write.csv(Trtdata,"FamilyLevelWOOther.csv")


FamilyPlot=ggplot(Trtdata, aes(x=SampleType,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+
  facet_grid(~Family)+xlab("Sample Location")+ylab("Relative Abundance (%, > 3%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(axis.title.x=element_blank())+facet_grid(Case~Timepoint2)#+  scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
FamilyPlot

```



##Genus level

```{r Summarizing3,warning=FALSE}
df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100


Trtdata <- ddply(df, c("Genus","SampleType","Case","Timepoint2"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance)
)
#Trtdata
write.csv(Trtdata,"GenusLevelWOOther.csv")


GenusPlot=ggplot(Trtdata, aes(x=SampleType,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+
  facet_grid(~Genus)+xlab("Sample Location")+ylab("Relative Abundance (%, > 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(axis.title.x=element_blank())+facet_grid(Case~Timepoint2)#+  scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
GenusPlot

```

#Ordinations
-All ordinations here were conducted using weighted unifrac (for other distince metrics see QIIME outputs)

-Ellipses with 95% CI for the mean of each group not shown 

##PCoA 
###
```{r PCoA, warning=FALSE}

ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="SampleType",shape="Timepoint2")+geom_point(size=5)+scale_colour_manual(values=cbPalette)+
  scale_fill_manual(values=cbPalette)
ordplot

```

##With Soil Removed
```{r}
StatSubset<-subset_samples(physeq,SampleType!="soil")

ord=ordinate(StatSubset,"PCoA", "wunifrac")
ordplot=plot_ordination(StatSubset, ord,"samples", color="SampleType",shape="Timepoint2")+geom_point(size=5)+scale_colour_manual(values=cbPalette)+
  scale_fill_manual(values=cbPalette)
ordplot

```



#PERMANOVAs

##Homogenieity of Multivariate Dispersions
```{r, warning=FALSE}
GPdist=phyloseq::distance(physeq, "wunifrac")
beta=betadisper(GPdist, sample_data(physeq)$Timepoint2)
permutest(beta)
boxplot(beta)
```

#RandomForest{.tabset}

##AllSamples 
```{r}
GenusAll
GenusAllSubset2<-subset_taxa(GenusAll, Genus != "NA")
GenusAllSubset2
ForestData=GenusAll#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Timepoint)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:20, ]#22

#write.csv(imp.20,"Top20GenusLevelPredictorsAllTimepoint.csv")

ImpWGenus<-read.csv("Top20GenusLevelPredictorsAllTimepoint.csv",header=T)
ImpWGenus$FamilyGenus<-paste0(ImpWGenus$Family,": ",ImpWGenus$Genus)
ImpWGenus<-arrange(ImpWGenus, desc(MeanDecreaseGini))
ImpWGenus$Genus <- factor(ImpWGenus$Genus, levels = ImpWGenus$Genus)
ggplot(ImpWGenus, aes(x = Genus, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip()
```

##Mouth Nose Ears subset

Only samples types present at all timepoints were included 
```{r}
GenusSubset<-subset_samples(GenusAll,SampleType=="mouth"|SampleType=="nose"|SampleType=="ears")

ForestData=GenusSubset#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Timepoint)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)

imp.20 <- imp.sort[1:20, ]
#imp.20
write.csv(imp.20,"Top20GenusLevelPredictorsSubsetTimepoint.csv")
ImpWGenus<-read.csv("Top20GenusLevelPredictorsSubsetTimepointWTaxa.csv",header=T)
ImpWGenus$FamilyGenus<-paste0(ImpWGenus$Family,": ",ImpWGenus$Genus)
ImpWGenus<-arrange(ImpWGenus, desc(MeanDecreaseGini))
ImpWGenus$Genus <- factor(ImpWGenus$Genus, levels = ImpWGenus$Genus)

ggplot(ImpWGenus, aes(x = Genus, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip()

#write.csv(tax_table(GenusAll),"GenusLevelTaxTable.csv")

```


Top 10 genus level indicator taxa for timepoint (MeanDecreaseGINI)
```{r}

df <- psmelt(GenusSubset)
df<- subset(df,Genus %in% ImpWGenus$Genus[1:10])
#df

df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus","SampleType","Case","Timepoint2"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance)
)
cdataplot=ggplot(Trtdata, aes(x=SampleType,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+ facet_grid(Case~Timepoint2)+xlab("Sample Type")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+scale_fill_viridis_d()
cdataplot

```