README.md

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

xviz package demo

yeguanhua

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:

| SampleID | diagnosis | | :-------- | :---------------- | | s17118657 | healthy | | s17118661 | healthy | | s17118667 | healthy | | s17118714 | healthy | | s17118646 | intestinal cancer | | s17118664 | intestinal cancer | | s17118669 | intestinal cancer | | s17118686 | intestinal cancer | | s17118647 | liver cancer | | s17118684 | liver cancer | | s17118715 | liver cancer | | s17118730 | liver cancer | | s17118650 | lung cancer | | s17118680 | lung cancer | | s17118691 | lung cancer | | s17118703 | lung cancer |

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()

| diagnosis | samples | Observed | Chao1 | ACE | Shannon | Simpson | InvSimpson | Fisher | | :---------------- | :-------- | -------: | -------: | -------: | -------: | --------: | ---------: | -------: | | healthy | s17118657 | 256 | 259.5000 | 259.7298 | 3.886926 | 0.9528851 | 21.224687 | 32.61626 | | healthy | s17118661 | 236 | 238.8125 | 241.4922 | 3.583856 | 0.9452327 | 18.259087 | 30.78986 | | healthy | s17118667 | 214 | 225.0000 | 220.5107 | 2.379044 | 0.6687893 | 3.019226 | 27.28225 | | healthy | s17118714 | 265 | 274.0000 | 269.1680 | 3.405226 | 0.8781308 | 8.205520 | 34.21672 | | intestinal cancer | s17118646 | 317 | 322.0000 | 321.9582 | 3.840584 | 0.9298637 | 14.257945 | 41.17106 | | intestinal cancer | s17118664 | 383 | 408.0000 | 398.4252 | 4.615068 | 0.9716829 | 35.314377 | 53.15833 | | intestinal cancer | s17118669 | 280 | 314.0000 | 288.6257 | 3.642227 | 0.9299103 | 14.267425 | 36.75522 | | intestinal cancer | s17118686 | 463 | 499.9091 | 479.0042 | 4.765762 | 0.9826985 | 57.798311 | 66.58539 | | liver cancer | s17118647 | 246 | 252.8750 | 252.3801 | 3.885810 | 0.9543788 | 21.919626 | 32.03147 | | liver cancer | s17118684 | 416 | 426.9091 | 422.5903 | 4.762457 | 0.9827857 | 58.091385 | 60.36914 | | liver cancer | s17118715 | 325 | 377.5000 | 335.9106 | 4.119560 | 0.9571293 | 23.325945 | 43.60859 | | liver cancer | s17118730 | 265 | 268.7500 | 269.9306 | 3.901370 | 0.9587435 | 24.238613 | 35.69064 | | lung cancer | s17118650 | 258 | 266.2727 | 266.3957 | 1.928436 | 0.4850172 | 1.941812 | 33.23692 | | lung cancer | s17118680 | 344 | 359.1111 | 353.7075 | 4.479324 | 0.9733006 | 37.453969 | 46.98003 | | lung cancer | s17118691 | 197 | 213.5000 | 204.0969 | 3.331857 | 0.9208330 | 12.631526 | 24.27305 | | lung cancer | s17118703 | 286 | 293.2000 | 289.6606 | 4.050606 | 0.9521931 | 20.917464 | 37.23602 |

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,])

| | Healthy | IntestinalCancer | LiverCancer | LungCancer | | :---- | -------: | ---------------: | ----------: | ---------: | | OTU1 | 3.344196 | 3.398114 | 3.685652 | 4.761228 | | OTU2 | 4.459920 | 3.800442 | 4.330170 | 3.850462 | | OTU3 | 4.607680 | 3.105510 | 3.398808 | 0.000000 | | OTU4 | 3.562293 | 3.815246 | 3.679882 | 4.287667 | | OTU5 | 3.841422 | 3.063709 | 3.291369 | 4.328970 | | OTU6 | 3.335257 | 4.333326 | 3.739018 | 3.129368 | | OTU7 | 3.749118 | 3.952453 | 3.788946 | 3.879383 | | OTU8 | 3.247728 | 4.188141 | 3.401917 | 3.894980 | | OTU9 | 3.950413 | 4.148788 | 3.278754 | 3.049218 | | OTU10 | 3.778368 | 3.604982 | 3.967922 | 3.808481 |

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/visualization416S documentation built on March 22, 2022, 9:03 p.m.