Calculate correlations between NMDS ordination axes (from beta diversity distances) and metadata, and plot a ggplot2 heatmap with numeric labels and significance stars.
For continuous variables, Pearson correlation is used. P values follow this convention: p < 0.001, p < 0.01, p < 0.05 = three, two, one stars, respectively.
Two-level categorical variables are coded by factor level order, so the sign of the correlation follows the level order. Categorical variables with more than two levels are summarized by the multiple correlation R (the square root of R^2 from a one-way linear model, score ~ variable). R is unsigned (0 to 1) and its p-value comes from the model F-test, so results do not depend on how levels are labelled. Variables that cannot be tested (a single level, constant values, or no residual degrees of freedom, e.g. a sample ID column) are dropped.
Default num_mds is 10, but function will internally set max number of NMDS to be n-1 where n is your number of samples.
Value
A list with elements: cor_matrix, pval_matrix, stars, plot (ggplot), and variance explained (squared correlations between ordination distances and observed dissimilarities.)
Examples
data(Snail1KO)
counts <- as.matrix(Snail1KO$rawCounts)
counts <- counts[rowSums(counts >= 10) >= 3, ]
meta <- Snail1KO$metadata
rownames(meta) <- meta$Sample
norm <- log2(t(t(counts) / colSums(counts)) * 1e6 + 1)
obj <- GenEDA(normalized = norm, metadata = meta, counts = counts)
set.seed(42)
oc <- PlotOrdCorr(obj, num_mds = 3, meta_cols = "Genotype")
#> Warning: stress is (nearly) zero: you may have insufficient data
oc$plot
