In krisrs1128/mvarVis: Visualize multivariate analysis results.

Loading / reshaping data

First we load the data from the phyloseq package.

library("phyloseq")
library("ggplot2")
library("mvarVis")

data("GlobalPatterns")
GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A = 0.5 * nsamples(GP))
GP1 = prune_taxa(wh0, GP)
# Pick low depth of the samples = 1000
GP1 = transform_sample_counts(GP1, function(x) 1e+03 * x/sum(x))
phylum.sum = tapply(taxa_sums(GP1), tax_table(GP1)[, "Phylum"], sum, na.rm = TRUE)
top5phyla = names(sort(phylum.sum, TRUE))[1:5]
GP1 = prune_taxa((tax_table(GP1)[, "Phylum"] %in% top5phyla), GP1)
human = get_variable(GP1, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
sample_data(GP1)$human <- factor(human)
rm(GlobalPatterns, GP, top5phyla, wh0, phylum.sum)

Boostrap Ordination

Define the distance we want to use, in this case a weighted UniFrac

distWUnifrac <- function (D) {
  if (!all((dim(D) == dim(GP1@otu_table)))) GP1 <- t(GP1)
  rownames(D) <- taxa_names(GP1)
  colnames(D) <- sample_names(GP1)
  otu_table(GP1) <- otu_table(D, TRUE)
  phyloseq::UniFrac(GP1, weighted = TRUE, normalized = TRUE)
}

bootOrd <- boot_ordination(GP1@otu_table, n = 20, method = "metaMDS", 
                          dist_method = distWUnifrac, common_depth = TRUE, 
                          replace_zero = FALSE, round = FALSE, scannf = F,
                          rows_annot = GP1@sam_data) 

plot_mvar(bootOrd, col = "SampleType", shape = "human", 
          layers_list = "point-and-contour", bins = 5) + theme_bw()
plot_mvar(bootOrd, col = "SampleType", shape = "human", 
          layers_list = "point-and-density")

bootOrd <- boot_ordination(t(GP1@otu_table), n = 20, method = "pco", 
                          dist_method = "euclidean", 
                          common_depth = TRUE, replace_zero = FALSE, 
                          round = FALSE, scannf = F, rows_annot = GP1@sam_data) 
plot_mvar(bootOrd, col = "SampleType",  shape = "human", layers_list = "point", bins = 5)
plot_mvar(bootOrd, col = "SampleType", shape = "human", 
          layers_list = "point-and-density")

bootOrd <- boot_ordination(GP1@otu_table, n = 20, method = "pco", 
                          dist_method = distWUnifrac, common_depth = TRUE, 
                          replace_zero = FALSE, round = FALSE, scannf = F,
                          rows_annot = GP1@sam_data) 

plot_mvar(bootOrd, col = "SampleType", shape = "human", layers_list = "point-and-contour")
plot_mvar(bootOrd, col = "SampleType", shape = "human", layers_list = "point-and-density")

bootOrd <- boot_ordination(GP1@otu_table, n = 20, method = "ade4_pca", 
                          dist_method = "bray", common_depth = TRUE, 
                          replace_zero = FALSE, round = FALSE, scannf = F,
                          rows_annot = GP1@tax_table, cols_annot = GP1@sam_data) 

plot_mvar(bootOrd, col = "Phylum", shape = "human", layers_list = "point-and-contour")
plot_mvar(bootOrd, col = "Phylum", shape = "human", layers_list = "point-and-density")

TO DO:

1. What is happening in "pco" from ade4 ?

############################ Check what is happening ########################

origDist <- dist_method(GP1@otu_table)
orig_ord <- ade4::dudi.pco(origDist, scannf = F)

# if only distance object is supplied then what is orig_ord$li and what is orig_ord$co????
# shouldn't there be only one set of cooredinates???

# Warning message:
#   In ordi_method(X, ...) : Non euclidean distance

krisrs1128/mvarVis documentation built on Oct. 13, 2019, 11:14 p.m.