README.Rmd

---
title: "visual16S package demo"
author: "yeguanhua"
output: github_document
---

```{r setup, echo=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.width=10, fig.height=6)
```

# Install dependencies

Before install the visual16S, there are some packages needed to be manually install from Bioconductor.

```{r dependancy, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("dada2", "phyloseq", "DESeq2", "EnhancedVolcano"))
```


# Install and load package

Note: 

You need to run install from the parent working directory that contains the visual16S folder.

```{r install, eval=FALSE}
devtools::install('visual16S')
```

#### If you can't install package in R Console, try this: open visual16S.Rproj first, then use RStudio → Build → Install and Restart.

# Load package and demo data.

```{r load}
library(visual16S)
demo_phyloseq_object <- visual16S::demo_phyloseq_object
demo_dada2_result <- visual16S::demo_dada2_result
```

# Data status

Primer: 

CCTAYGGGRBGCASCAG ; GGACTACNNGGGTATCTAAT

DADA2 filter parameters: 

dada2::filterAndTrim(truncLen=c(0,0), maxEE=c(2,2))

DADA2 taxonomy database: 

silva_nr_v132

Metadata: 

```{r matadata, echo=FALSE}
extract_metadata_phyloseq(demo_phyloseq_object) %>% arrange(diagnosis) %>% knitr::kable()
```

# DADA2 workflow reads track

First use the track_reads_dada2 function to check reads drop associated with every step in DADA2.

Note:

You can plot reads track in either absolute abundance or relative abundance.

If the sample size is too large to show on legend, set legend_position to "none".

```{r track reads}
track_reads_dada2(demo_dada2_result$reads_track, 
                  single_end = FALSE, 
                  relative_abundance = TRUE, 
                  legend_position = "top")
```

# Plot sparsity - plot_sparsity()

```{r sparsity}
plot_sparsity(demo_dada2_result$seq_tab, binwidth = 10)
```

# Stacked barplot of phylogenetic composition - plot_stacked_bar()

Use plot_stacked_bar function to plot the taxonomy level abundance in every sample. You can change the 'level' argument to plot abundance in different level, or change the 'feature' to choose different feature you want to show in x-axis.

### Minimum usage

```{r Stacked barplot 1, fig.height=12}
plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Class")
```

### Plot in relative abundance

```{r Stacked barplot 5, fig.height=12}
plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Class", 
                 relative_abundance = TRUE)
```

### Set "feature" parameter to show feature information

```{r Stacked barplot 3, fig.height=12}
plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Class", 
                 relative_abundance = TRUE, 
                 feature = "diagnosis")
```

### Pass ordered sample names to "order" parameter to plot in specific order

```{r Stacked barplot 4, fig.height=12}
sample_order <- extract_metadata_phyloseq(demo_phyloseq_object) %>% arrange(diagnosis) %>% .$SampleID
plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Class", 
                 relative_abundance = TRUE, 
                 feature = "diagnosis", 
                 order = sample_order)
```

# Alpha diversity - plot_alpha_diversity()

Use plot_alpha_diversity to plot alpha diversity. 

### All alpha diversity

```{r All alpha}
plot_alpha_diversity(phyloseq = demo_phyloseq_object, feature = "diagnosis") %>% 
  # Show result in a table in markdown
  knitr::kable()
```

Or you can change 'measures' argument to use different measurement.

### Chao1

```{r Chao1}
plot_alpha_diversity(phyloseq = demo_phyloseq_object, 
                     feature = "diagnosis", 
                     measures = "Chao1", 
                     p_test = "kruskal")
```

# Beta diversity - plot_beta_diversity()

Use plot_beta_diversity to plot beta diversity. Change 'method' to draw different beta diversity plot. You can locate specific sample in beta diversity plot by the table printed to the screen.

### Bray-Curtis

```{r Bray-Curtis}
plot_beta_diversity(phyloseq = demo_phyloseq_object, 
                    feature = "diagnosis", 
                    method = "bray")
```

# Log2 fold change - log2fc()

Use log2fc function to show differential analysis result.

### Minimum usage

```{r log2fc 1}
log2fc(phyloseq = demo_phyloseq_object, 
       feature = "diagnosis", 
       p_value = 0.05)
```

### Choose a taxonomy level to calculate log2fc.

```{r log2fc 3}
log2fc(phyloseq = demo_phyloseq_object, 
       feature = "diagnosis", 
       p_value = 0.05, 
       level = "Genus")
```

### Set "reference" and "treatment" parameters to change log2fc treatment vs reference. Both "reference" and "treatment" should be one of the levels in "feature".

```{r log2fc 2}
log2fc(phyloseq = demo_phyloseq_object, 
       feature = "diagnosis", 
       p_value = 0.05, 
       level = "Genus", 
       reference = 'healthy', 
       treatment = 'liver cancer')
```

# Correlation - plot_correlation()

```{r correlation table}
# Construct correlation table
cor_tab <- demo_dada2_result$seq_tab %>% t() %>% as.data.frame() %>% 
  rownames_to_column("OTU") %>% 
  mutate(Healthy = s17118657 + s17118661 + s17118667 + s17118714) %>% 
  mutate(IntestinalCancer = s17118646 + s17118664 + s17118669 + s17118686) %>% 
  mutate(LiverCancer = s17118647 + s17118684 + s17118715 + s17118730) %>% 
  mutate(LungCancer = s17118650 + s17118680 + s17118691 + s17118703) %>% 
  select(OTU, Healthy, IntestinalCancer, LiverCancer, LungCancer) %>% 
  column_to_rownames("OTU")
# Normalization
cor_tab <- (cor_tab + 1) %>% log10()
knitr::kable(cor_tab[1:10,])
```

### Minimum usage

```{r correlation 1}
plot_correlation(cor_tab, x = "IntestinalCancer", y = "Healthy")
```

### Multiple correlation in one plot

```{r correlation 2}
plot_correlation(cor_tab, x = c("IntestinalCancer", "LiverCancer", "LungCancer"), y = "Healthy")
```