In yeguanhuav/visual16S: Visualization functions for 16S and Metagenomics analysis

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

Install xviz

# Before installing the xviz, there are some packages needed to be manually installed from Bioconductor.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("dada2", "phyloseq", "DESeq2", "EnhancedVolcano", "ComplexHeatmap"))
# Note: You need to run devtools::install() from parent working directory that contains the xviz folder.
setwd("../")
devtools::install('xviz')
# If you can't install package in R Console, try this: open xviz.Rproj first, then use RStudio → Build → Install and Restart.

Load package and demo data.

library(xviz)
demo_phyloseq_object <- xviz::demo_phyloseq_object
demo_dada2_result <- xviz::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:

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".

track_reads_dada2(demo_dada2_result$reads_track, 
                  single_end = FALSE, 
                  relative_abundance = TRUE, 
                  legend_position = "top")

Plot sparsity - plot_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.

Minimum usage: plot in relative abundance

plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Family")

Set feature parameter to show feature information

plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Family", 
                 relative_abundance = TRUE, 
                 feature = "diagnosis")

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

sample_order <- extract_metadata_phyloseq(demo_phyloseq_object) %>% arrange(diagnosis) %>% .$SampleID
plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Family", 
                 relative_abundance = TRUE, 
                 feature = "diagnosis", 
                 order = sample_order)

Use facet_wrap(vars, scale="free") funciton to facet stacked barplot

When faceting, feature parameter will not work.

plot_stacked_bar(phyloseq = demo_phyloseq_object, 
                 level = "Family", 
                 relative_abundance = TRUE) + 
  facet_wrap("diagnosis", scale="free")

Alpha diversity - plot_alpha_diversity()

Use plot_alpha_diversity to plot alpha diversity.

All alpha diversity

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.

Shannon Index

plot_alpha_diversity(phyloseq = demo_phyloseq_object, 
                     feature = "diagnosis", 
                     measures = "Shannon", 
                     label = TRUE)

Shannon Index (Genus Level)

plot_alpha_diversity(phyloseq = demo_phyloseq_object, 
                     level = "Genus",
                     feature = "diagnosis", 
                     measures = "Shannon", 
                     label = TRUE)

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

plot_beta_diversity(phyloseq = demo_phyloseq_object, 
                    feature = "diagnosis", 
                    method = "bray", 
                    label = TRUE)

Bray-Curtis with PC1 and PC2 info

plot_beta_diversity(phyloseq = demo_phyloseq_object, 
                    feature = "diagnosis", 
                    method = "bray", 
                    label = TRUE, 
                    add_pc1 = TRUE, add_pc1_otu_level = "Genus", add_pc1_cor_method = "spearman", 
                    add_pc2 = TRUE, add_pc2_otu_level = "Genus", add_pc2_cor_method = "spearman")

Bray-Curtis (Genus Level)

plot_beta_diversity(phyloseq = demo_phyloseq_object, 
                    level = "Genus",
                    feature = "diagnosis", 
                    method = "bray", 
                    label = TRUE)

Bray-Curtis (Genus Level) with PC1 and PC2 info

plot_beta_diversity(phyloseq = demo_phyloseq_object, 
                    level = "Genus",
                    feature = "diagnosis", 
                    method = "bray", 
                    label = TRUE, 
                    add_pc1 = TRUE, add_pc1_cor_method = "spearman", 
                    add_pc2 = TRUE, add_pc2_cor_method = "spearman")

Log2 fold change - log2fc()

Use log2fc function to show differential analysis result.

Minimum usage

log2fc(phyloseq = demo_phyloseq_object, 
       feature = "diagnosis", 
       p_value = 0.05)

Choose a taxonomy level to calculate log2 fold change.

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.

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

Correlation - plot_correlation()

# 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

plot_correlation(cor_tab, x = "IntestinalCancer", y = "Healthy")

Multiple correlation in one plot

plot_correlation(cor_tab, x = c("IntestinalCancer", "LiverCancer", "LungCancer"), y = "Healthy")

Correlation heatmap

plot_correlation(cor_tab, heatmap = TRUE)

Correlation heatmap in given order

plot_correlation(cor_tab, heatmap = TRUE, 
                 row_order = c("Healthy", "LungCancer", "LiverCancer", "IntestinalCancer"), 
                 col_order = c("Healthy", "LungCancer", "LiverCancer", "IntestinalCancer"))

yeguanhuav/visual16S documentation built on Feb. 19, 2022, 10:32 a.m.