--- name: bio-microbiome-diversity-analysis description: Alpha and beta diversity analysis for microbiome data. Calculate within-sample richness, evenness, and between-sample dissimilarity with phyloseq and vegan. Use when comparing community composition across samples or testing for group differences in microbiome structure. tool_type: r primary_tool: phyloseq --- ## Version Compatibility Reference examples tested with: R stats (base), ggplot2 3.5+, phyloseq 1.46+, scanpy 1.10+, vegan 2.6+ Before using code patterns, verify installed versions match. If versions differ: - R: `packageVersion('')` then `?function_name` to verify parameters If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. # Diversity Analysis **"Compare microbial diversity across my samples"** → Calculate alpha diversity (within-sample richness/evenness) and beta diversity (between-sample dissimilarity) to test for community composition differences across groups. - R: `phyloseq::estimate_richness()` for alpha, `phyloseq::ordinate()` for beta - R: `vegan::adonis2()` for PERMANOVA testing ## Create phyloseq Object ```r library(phyloseq) library(vegan) library(ggplot2) seqtab <- readRDS('seqtab_nochim.rds') taxa <- readRDS('taxa.rds') metadata <- read.csv('sample_metadata.csv', row.names = 1) ps <- phyloseq(otu_table(seqtab, taxa_are_rows = FALSE), tax_table(taxa), sample_data(metadata)) taxa_names(ps) <- paste0('ASV', seq(ntaxa(ps))) ``` ## Alpha Diversity ```r # Calculate multiple metrics alpha_div <- estimate_richness(ps, measures = c('Observed', 'Chao1', 'Shannon', 'Simpson')) alpha_div$SampleID <- rownames(alpha_div) alpha_div <- merge(alpha_div, sample_data(ps), by = 'row.names') # Statistical test kruskal.test(Shannon ~ Group, data = alpha_div) # Pairwise comparisons pairwise.wilcox.test(alpha_div$Shannon, alpha_div$Group, p.adjust.method = 'BH') ``` ## Alpha Diversity Plots ```r plot_richness(ps, x = 'Group', measures = c('Observed', 'Shannon')) + geom_boxplot() + theme_minimal() # Custom plot ggplot(alpha_div, aes(x = Group, y = Shannon, fill = Group)) + geom_boxplot() + geom_jitter(width = 0.2, alpha = 0.5) + theme_minimal() + labs(y = 'Shannon Diversity Index') ``` ## Faith's Phylogenetic Diversity **Goal:** Calculate phylogenetic alpha diversity (Faith's PD) from ASV data by building a de novo phylogeny and summing branch lengths. **Approach:** Align ASV sequences with DECIPHER, construct a neighbor-joining tree with phangorn, root at midpoint, and compute PD using picante. ```r library(picante) # Requires phylogenetic tree in phyloseq object # Build tree from ASV sequences library(DECIPHER) library(phangorn) seqs <- refseq(ps) alignment <- AlignSeqs(seqs, anchor = NA) phang_align <- phyDat(as(alignment, 'matrix'), type = 'DNA') dm <- dist.ml(phang_align) tree <- NJ(dm) tree <- midpoint(tree) phy_tree(ps) <- tree # Calculate Faith's PD otu_mat <- as.matrix(t(otu_table(ps))) faith_pd <- pd(otu_mat, phy_tree(ps), include.root = TRUE) alpha_div$PD <- faith_pd$PD ``` ## Rarefaction Curves ```r # Check if sequencing depth is adequate rarecurve_data <- vegan::rarecurve(t(otu_table(ps)), step = 100, sample = min(sample_sums(ps))) # ggplot version with ggrare (install from GitHub) # devtools::install_github('gauravsk/ranacapa') library(ranacapa) p_rare <- ggrare(ps, step = 100, color = 'Group', se = FALSE) p_rare + theme_minimal() + labs(title = 'Rarefaction Curves') ``` ## Rarefaction ```r # Check sequencing depth sample_sums(ps) # Rarefy to minimum depth ps_rarefied <- rarefy_even_depth(ps, sample.size = min(sample_sums(ps)), rngseed = 42, replace = FALSE) ``` ## Beta Diversity ```r # Calculate distance matrices bray <- phyloseq::distance(ps, method = 'bray') # Bray-Curtis jaccard <- phyloseq::distance(ps, method = 'jaccard') # Jaccard unifrac <- UniFrac(ps, weighted = TRUE) # Weighted UniFrac (requires tree) # Ordination ord_bray <- ordinate(ps, method = 'PCoA', distance = bray) # Plot plot_ordination(ps, ord_bray, color = 'Group') + stat_ellipse(level = 0.95) + theme_minimal() ``` ## PERMANOVA ```r # Test for group differences metadata <- data.frame(sample_data(ps)) permanova_result <- adonis2(bray ~ Group, data = metadata, permutations = 999) permanova_result # With covariates adonis2(bray ~ Group + Age + Sex, data = metadata, permutations = 999) ``` ## Beta Dispersion ```r # Test homogeneity of dispersions (assumption of PERMANOVA) beta_disp <- betadisper(bray, metadata$Group) permutest(beta_disp) plot(beta_disp) ``` ## NMDS Ordination ```r ord_nmds <- ordinate(ps, method = 'NMDS', distance = bray) # Check stress ord_nmds$stress # Should be < 0.2 plot_ordination(ps, ord_nmds, color = 'Group') + theme_minimal() ``` ## Distance Metrics Comparison | Metric | Type | Considers Abundance | Phylogeny | |--------|------|---------------------|-----------| | Bray-Curtis | Quantitative | Yes | No | | Jaccard | Binary | No | No | | UniFrac (unweighted) | Binary | No | Yes | | UniFrac (weighted) | Quantitative | Yes | Yes | ## Related Skills - amplicon-processing - Generate ASV table - differential-abundance - Identify taxa driving differences - data-visualization/ggplot2-fundamentals - Custom diversity plots - ecological-genomics/biodiversity-metrics - Hill number coverage-based rarefaction for ecological data - ecological-genomics/community-ecology - Constrained ordination and indicator species

bio-microbiome-diversity-analysis

Install this skill

Comments