如何为基于phyloseq的微生物组排序图添加物种箭头/向量

First, let's break down the core issue: MDS/PCoA ordinations (the ones you’re creating with ordinate(..., "MDS")) are built from sample-to-sample distance matrices, not directly from raw species abundance data. Unlike methods like PCA or RDA, they don’t calculate inherent "species scores" because their sole focus is arranging samples based on pairwise similarity. That’s exactly why scores(MDS, display = "species") throws an error—there are no species-related scores stored in that MDS object to retrieve, no matter which distance metric (UniFrac or Bray-Curtis) you use.

Luckily, there are two reliable ways to add species-driven arrows to your ordination plot, depending on whether you want to stick with your MDS/PCoA workflow or switch to a different ordination method:

Option 1: Switch to a Species-Based Ordination (PCA/RDA) #

If you’re open to using an ordination method that natively calculates species loadings (scores), PCA or RDA are perfect choices. These methods use raw species abundance data directly, so you can easily extract species scores and plot them as arrows.

Example with PCA: #

# Run PCA on your phyloseq object
pca_ord <- ordinate(physeqobject, "PCA")

# Create the base sample plot
plot_pca <- plot_ordination(physeqobject, pca_ord, color = "variable1")

# Extract species scores from the PCA object and format as a data frame
species_scores <- scores(pca_ord, display = "species")
species_scores_df <- as.data.frame(species_scores)
species_scores_df$Species <- rownames(species_scores_df)

# Add species arrows and labels to the plot
plot_pca +
  geom_segment(data = species_scores_df,
               aes(x = 0, xend = PC1, y = 0, yend = PC2),
               arrow = arrow(length = unit(0.02, "npc")),
               color = "darkred", alpha = 0.7) +
  geom_text(data = species_scores_df,
            aes(x = PC1 * 1.1, y = PC2 * 1.1, label = Species),
            color = "darkred", size = 3)

If you have environmental variables to incorporate, swap PCA for RDA—its species scores work the exact same way for plotting arrows.

Option 2: Keep Using MDS/PCoA, Fit Species with envfit (vegan Package) #

If you need to stick with MDS/PCoA (e.g., because you rely on UniFrac distances), use the envfit() function from the vegan package to fit species abundances as "environmental variables" onto your existing ordination axes. This tests which species are significantly correlated with the MDS axes, then lets you plot those meaningful correlations as arrows.

Step-by-Step Code: #

# Load vegan (install first if needed: install.packages("vegan"))
library(vegan)

# Your existing MDS ordination (using precomputed UniFrac matrix)
MDS <- ordinate(physeqobject, "MDS", distance = unifracmatrix)
base_plot <- plot_ordination(physeqobject, MDS, color = "variable1")

# Extract species abundance matrix (transpose so samples are rows, species are columns)
otu_matrix <- t(as.matrix(otu_table(physeqobject)))

# Fit species abundances to MDS axes (permutations calculate significance)
species_fit <- envfit(MDS, otu_matrix, permutations = 999)

# Filter to only keep species with a significant correlation (p < 0.05, adjust as needed)
sig_species <- as.data.frame(species_fit$vectors[species_fit$vectors$pvals < 0.05, ])
sig_species$Species <- rownames(sig_species)

# Add significant species arrows and labels to your base plot
base_plot +
  geom_segment(data = sig_species,
               aes(x = 0, xend = MDS1, y = 0, yend = MDS2),
               arrow = arrow(length = unit(0.02, "npc")),
               color = "darkgreen", alpha = 0.8) +
  geom_text(data = sig_species,
            aes(x = MDS1 * 1.1, y = MDS2 * 1.1, label = Species),
            color = "darkgreen", size = 3)

Key Notes for envfit: #

• The permutations argument runs random permutations to calculate p-values—only plot species with low p-values to avoid cluttering the plot with noise.
• Arrow length corresponds to the strength of the correlation between the species and the ordination axes.

内容的提问来源于stack exchange，提问作者eok