beta_diversity_template.Rmd

---
title: "MiSeq-${miseq_run_number}"
output:
  html_notebook:
    df_print: paged
    number_sections: yes
    theme: lumen
    toc: yes
    code_folding: hide
---

<style>
  html {
    font-size: 16pt;
  }
  
  body {
    font-size: 16pt;
  }
  
  h1 {
    font-size: 2.2rem;
  }
  
  h2 {
    font-size: 2rem;
  }
  
  h3 {
    font-size: 1.8rem;
  }
  
  h4 {
    font-size: 1.4rem;
  }
  
</style>

# {.tabset}

## [Close]

## Setup

### Start Conda ENV
```{r}
startCondaEnv('betadiv', lib='~/R35')
```

### Set Knitr Options
```{r}
library(rmarkdown)
library(knitr)

opts_chunk$set(
  echo=TRUE,
  dpi=300,
  fig.width=12
  )
```

### Load Libraries
```{r, echo=F}
library(tidyverse)
library(magrittr)
library(ggplot2)
library(ggbeeswarm)
library(openxlsx)
library(DT)
library(pheatmap)
library(kableExtra)
library(vegan)
library(phyloseq)
library(plotly)

### Custom libraries that can be loaded from GitHub
source('~/utilities/analysis-utilities/general_asv_data_analysis_utilities.R')
source('~/utilities/analysis-utilities/beta_diversity_utilities.R')
```
### Set Paths
```{r}
project_dir = "~/miseq-239"
dada2_table_dir = file.path(project_dir, "processed_data")
metadata_dir = file.path(project_dir, "metadata")
```

### Load Data
```{r}
asv_table = 
  read.delim(
    file.path(dada2_table_dir, 'miseq-239_tall_asv_table.tsv'), 
    stringsAsFactors=F
    )

taxonomy_table = 
  read.delim(
    file.path(dada2_table_dir, 'miseq-239_taxonomy_table.tsv'), 
    stringsAsFactors=F
    ) %>%
  mutate_all(as.character())

sample_metadata = 
  read.delim(
    file.path(metadata_dir, 'miseq-239_metadata_CNV_GA.tsv'), 
    stringsAsFactors=F
    ) %>%
  filter(SampleID != 'Undetermined') %>%
  mutate(AgeGroup = ifelse(Age<90, 'Age_lt_90', 'Age_ge_90')) %>%
  mutate(
    CaseString = factor(CaseString, levels=c('Control', 'AMD')),
    Gender = factor(Gender, levels=c('F', 'M')), 
    AREDS = factor(AREDS, levels=c('N','Y'))
  )

sample_names = sample_metadata$SampleName
# 
# asv_colanmes = 
#   data.frame(SampleID=colnames(asv_table)) %>%
#   inner_join(sample_metadata, by='SampleID') %>%
#   pull(SampleName) %>%
#   c('ASVs', .)
# 
# colnames(asv_table) = asv_colanmes
```

### Get Filtered Taxa Abundance
```{r}
taxa_abundance = getFilteredTaxaCounts(
  asv_table,
  taxonomy_table,
  sample_metadata, 
  lowest_rank='Genus',
  relative_abundance_cutoff=0.005,
  prevalence_cutoff=0.1,
  id_col="SampleName",
  normalize=T
  )


filtered_taxa = taxa_abundance$short_glommed_taxa

taxa_abundance %>%
  datatable(
    extensions = c('Buttons'), 
    options = list(
      paging = F,
      searching = T,
      autoWidth = TRUE,
      ordering = TRUE,
      dom = 'tB',
      buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
      scrollX=T,
      scrollY=600
      ),
    class="display",
    caption="Taxa Abundance"
  )

```

### Build Master Table
```{r}
master_table = getMasterTable(
  taxa_abundance,
  sample_metadata,
  'SampleName',
  'short_glommed_taxa'
)
```


## MDS Plots {.tabset}
### Scree Plot
```{r}
# plotNMDSScree <- function(x) 
# {
#   plot(
#     rep(1, 10), 
#     replicate(10, metaMDS(x, autotransform = F, k = 1)$stress), 
#     xlim = c(1, 10),
#     ylim = c(0, 0.30), 
#     xlab = "# of Dimensions", 
#     ylab = "Stress", 
#     main = "NMDS stress plot"
#     )
#   for (i in 1:10)
#   {
#     points(
#       rep(i + 1,10),
#       replicate(10, metaMDS(x, autotransform = F, k = i + 1)$stress)
#       )
#   }
# }
# 
# plotNMDSScree(master_table %>% select(filtered_taxa))
```


### 2D
```{r, echo=F}
makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,2),
  color='CaseString',
  title_extra='Case v Control'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(2,3),
  color='CaseString',
  title_extra='Case v Control'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,3),
  color='CaseString',
  title_extra='Case v Control'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,2),
  color='Gender',
  title_extra='Gender'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(2,3),
  color='Gender',
  title_extra='Gender'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,3),
  color='Gender',
  title_extra='Gender'
)


makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,2),
  color='Tissue_code',
  title_extra='Tissue_code'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(2,3),
  color='Gender',
  title_extra='Tissue_code'
)

makeNMDSPlot(
  master_table,
  features=filtered_taxa,
  axes=c(1,3),
  color='Gender',
  title_extra='Tissue_code'
)
```

### 3D
```{r}
pc = getNMDS(master_table, filtered_taxa, axes=c(1:3))
pcdf = data.frame(pc$points) %>% rownames_to_column('SampleName')
pcmt = inner_join(pcdf, master_table, by='SampleName')

axx <- list(
  gridcolor='rgb(255, 255, 255)',
  zerolinecolor='rgb(255, 255, 255)',
  showbackground=TRUE,
  backgroundcolor='rgb(230, 230,230)'
)

plot_ly(
  data=pcmt, 
  x = ~MDS1, 
  y = ~MDS2, 
  z = ~MDS3, 
  color=~CaseString
  ) %>%
  add_markers(size=3) %>%
  # layout(autosize = F, width = 1000, height = 1000)
  layout(
    # autosize=F,
    # width=700,
    # height=700,
    aspectmode='cube',
    title = 'MiSeq-239 Principal Components',
    scene = list(xaxis=axx, yaxis=axx, zaxis=axx),
    paper_bgcolor = 'rgb(243, 243, 243)',
    plot_bgcolor = 'rgb(243, 243, 243)'
    )
```


## PCoA Plots
### Scree Plot
```{r}
pc = getPrinceComps(master_table, filtered_taxa, 'SampleName')
screeplot(pc)
```

### 2D
```{r}
makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,2),
  color_group='CaseString',
  title_extra='Case v Control'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(2,3),
  color_group='CaseString',
  title_extra='Case v Control'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,3),
  color_group='CaseString',
  title_extra='Case v Control'
)

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,2),
  color_group='Gender',
  title_extra='Gender'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(2,3),
  color_group='Gender',
  title_extra='Gender'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,3),
  color_group='Gender',
  title_extra='Gender'
)

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,2),
  color_group='Tissue_code',
  title_extra='Tissue_code'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(2,3),
  color_group='Tissue_code',
  title_extra='Tissue_code'
) 

makePCOAPlot(
  master_table,
  method='bray',
  features=filtered_taxa,
  axes=c(1,3),
  color_group='Tissue_code',
  title_extra='Tissue_code'
)

# + aes(x=PC1, y=PC2, color=IGA, shape=CaseString)
```

### 3D
```{r}
pc = getPrinceComps(master_table, filtered_taxa, 'SampleName')
pcdf = data.frame(pc$x) %>% rownames_to_column('SampleName')
pcmt = inner_join(pcdf, master_table, by='SampleName')
axx <- list(
  gridcolor='rgb(255, 255, 255)',
  zerolinecolor='rgb(255, 255, 255)',
  showbackground=TRUE,
  backgroundcolor='rgb(230, 230,230)'
)

plot_ly(
  data=pcmt, 
  x = ~PC1, 
  y = ~PC2, 
  z = ~PC3, 
  color=~CaseString
  ) %>%
  add_markers(size=3) %>%
  # layout(autosize = F, width = 1000, height = 1000)
  layout(
    # autosize=F,
    # width=700,
    # height=700,
    aspectmode='cube',
    title = 'MiSeq-239 Principal Components',
    scene = list(xaxis=axx, yaxis=axx, zaxis=axx),
    paper_bgcolor = 'rgb(243, 243, 243)',
    plot_bgcolor = 'rgb(243, 243, 243)'
    )
```