--- title: "MeLSI Tutorial: Basic Usage and Examples" author: "Nathan Bresette" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true toc_float: true vignette: > %\VignetteIndexEntry{MeLSI Tutorial} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) ``` # MeLSI Tutorial This tutorial provides a comprehensive introduction to MeLSI (Metric Learning for Statistical Inference) for microbiome data analysis. ## Abstract MeLSI is a novel machine learning method that learns optimal distance metrics to improve statistical power in detecting group differences in microbiome data. Unlike traditional distance metrics (Bray-Curtis, Euclidean, Jaccard), MeLSI adapts to the specific characteristics of your dataset. This package provides a statistically rigorous approach to beta diversity analysis that automatically identifies which microbial features drive group differences, making it particularly valuable for biomarker discovery and hypothesis generation in microbiome research. ## Introduction MeLSI addresses a fundamental limitation in current microbiome beta diversity analysis: traditional distance metrics treat all microbial taxa equally, missing subtle but biologically important differences. MeLSI learns adaptive distance metrics optimized for your specific dataset while maintaining rigorous statistical validity through permutation testing. ### Comparison with Existing Bioconductor Packages MeLSI complements existing Bioconductor microbiome analysis packages: - **phyloseq** and **microbiome**: These packages provide comprehensive data structures and analysis pipelines for microbiome data. MeLSI integrates seamlessly with data from these packages (see "Working with Bioconductor Packages" section below) and adds a novel metric learning approach for beta diversity analysis. - **vegan**: MeLSI uses `vegan::adonis2()` for PERMANOVA calculations but extends beyond fixed distance metrics by learning optimal metrics from data. - **Existing distance-based methods**: Unlike packages that use fixed distance metrics (Bray-Curtis, Jaccard, UniFrac), MeLSI learns which features matter most for each specific comparison, providing both improved statistical power and biological interpretability through feature importance weights. ## Installation and Setup ```{r install, eval=FALSE} if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("MeLSI") ``` ## Basic Usage ### Step 1: Prepare Your Data ```{r basic_setup} # Load required packages library(MeLSI) library(vegan) library(ggplot2) # Generate synthetic microbiome data for demonstration # Using smaller dataset for faster vignette execution set.seed(42) test_data <- generate_test_data(n_samples = 40, n_taxa = 50, n_signal_taxa = 5) X <- test_data$counts y <- test_data$metadata$Group # Display data dimensions cat("Data dimensions:", nrow(X), "samples x", ncol(X), "taxa\n") cat("Groups:", paste(unique(y), collapse = ", "), "\n") ``` ### Step 2: Data Preprocessing ```{r preprocessing} # CLR transformation (recommended for microbiome data) # Use the built-in helper function for easy CLR transformation X_clr <- clr_transform(X) # Display preprocessing results cat("CLR transformation completed\n") cat("Data range:", round(range(X_clr), 3), "\n") ``` ### Step 3: Run MeLSI Analysis ```{r run_melsi} # Run MeLSI with default parameters # Note: For real analyses, use n_perms = 200+ for reliable p-values # Using reduced parameters here for faster vignette execution (<15 min for Bioconductor) cat("Running MeLSI analysis...\n") melsi_results <- melsi( X_clr, y, n_perms = 50, # Number of permutations (use 200+ for real analysis) B = 20, # Number of weak learners (default 30, reduced for speed) m_frac = 0.8, # Fraction of features per learner show_progress = TRUE ) # Display results cat("\nMeLSI Results:\n") cat(sprintf("F-statistic value: %.4f\n", melsi_results$F_observed)) cat(sprintf("P-value: %.4f\n", melsi_results$p_value)) cat(sprintf("Significant: %s\n", ifelse(melsi_results$p_value < 0.05, "Yes", "No"))) # VIP plot is automatically generated, but you can also create it manually: # Using top_n = 10 for faster execution plot_vip(melsi_results, top_n = 10) # PCoA plot using the learned MeLSI distance plot_pcoa(melsi_results, X_clr, y) ``` ## Comparison with Standard Methods ```{r comparison} # Compare with Euclidean distance dist_euclidean <- dist(X_clr) adonis_euclidean <- adonis2(dist_euclidean ~ y, permutations = 50) # Compare with Bray-Curtis distance dist_bray <- vegdist(X, method = "bray") adonis_bray <- adonis2(dist_bray ~ y, permutations = 50) # Display comparison results cat("\nMethod Comparison:\n") cat(sprintf("MeLSI: F-statistic = %.4f, p = %.4f\n", melsi_results$F_observed, melsi_results$p_value)) cat(sprintf("Euclidean: F-statistic = %.4f, p = %.4f\n", adonis_euclidean[["F"]][1], adonis_euclidean[["Pr(>F)"]][1])) cat(sprintf("Bray-Curtis: F-statistic = %.4f, p = %.4f\n", adonis_bray[["F"]][1], adonis_bray[["Pr(>F)"]][1])) # Calculate improvement improvement <- (melsi_results$F_observed - adonis_euclidean[["F"]][1]) / adonis_euclidean[["F"]][1] * 100 cat(sprintf("\nMeLSI improvement over Euclidean: %.1f%%\n", improvement)) ``` ## Visualizing Results ### Variable Importance Plot (VIP) MeLSI automatically generates a VIP plot showing which features drive group differences. You can customize it using the `plot_vip()` function: ```{r vip_plot} # Plot VIP with directionality (default) # Using top_n = 10 for faster execution plot_vip(melsi_results, top_n = 10, title = "Top 10 Features by Importance") # Plot VIP without directionality coloring plot_vip(melsi_results, top_n = 10, directionality = FALSE) # Extract and examine feature weights programmatically feature_weights <- melsi_results$feature_weights top_features <- head(sort(feature_weights, decreasing = TRUE), 10) cat("\nTop 10 Most Weighted Features:\n") for (i in 1:length(top_features)) { cat(sprintf("%2d. %s (weight: %.4f)\n", i, names(top_features)[i], top_features[i])) } ``` ### Principal Coordinates Analysis (PCoA) Visualize the learned MeLSI distance using PCoA: ```{r pcoa_plot} # Create PCoA plot using the learned MeLSI distance matrix plot_pcoa(melsi_results, X_clr, y, title = "PCoA: MeLSI Distance") ``` ## Advanced Usage: Custom Parameters You can customize MeLSI parameters for your specific needs: ```{r advanced_usage, eval=FALSE} # Run MeLSI with custom parameters # Note: This example is set to eval=FALSE for faster vignette execution # Uncomment and run locally for full demonstration cat("Running MeLSI with custom parameters...\n") custom_results <- melsi( X_clr, y, n_perms = 200, # More permutations for higher precision (recommended: 200+) B = 50, # Larger ensemble for more stable results m_frac = 0.8, # Fraction of features per learner show_progress = TRUE ) cat(sprintf("Custom MeLSI F-statistic value: %.4f (p = %.4f)\n", custom_results$F_observed, custom_results$p_value)) ``` **Note:** The advanced usage example above is disabled (`eval=FALSE`) to keep vignette execution time under 15 minutes for Bioconductor requirements. For real analyses, use the parameters shown above (n_perms = 200+, B = 30-50). **Parameter Guidelines:** - `n_perms`: Use 200+ for reliable p-values in real analyses (vignette uses 50 for speed) - `B`: Default 30 is usually sufficient; increase to 50-100 for more stability (vignette uses 20 for speed) - `m_frac`: Default 0.8 works well; lower values (0.6-0.7) can help with very high-dimensional data ## Working with Bioconductor Packages ### Integration with phyloseq MeLSI integrates seamlessly with the `phyloseq` package, a core Bioconductor package for microbiome data: ```{r phyloseq_integration, eval=FALSE} # Load phyloseq data library(phyloseq) data(GlobalPatterns) # Example phyloseq object # Extract OTU table and sample data X <- as(otu_table(GlobalPatterns), "matrix") y <- sample_data(GlobalPatterns)$SampleType # CLR transformation X_clr <- clr_transform(X) # Run MeLSI analysis results <- melsi(X_clr, y, n_perms = 200, B = 30, show_progress = FALSE) # Note: For vignette execution, use reduced parameters (n_perms=50, B=20) to stay under 15 minutes # Results can be used with other phyloseq functions # For example, you can add the learned distance matrix to a phyloseq object ``` This interoperability ensures MeLSI works within the existing Bioconductor ecosystem and can be integrated into comprehensive microbiome analysis workflows. ## Summary This tutorial demonstrated: 1. **Basic MeLSI usage** with synthetic microbiome data 2. **Easy data preprocessing** using `clr_transform()` 3. **Automatic visualization** with `plot_vip()` and `plot_pcoa()` 4. **Comparison with standard methods** (Euclidean, Bray-Curtis) 5. **Custom parameters** for advanced usage ### Key Takeaways - MeLSI consistently outperforms standard distance metrics - Helper functions (`clr_transform`, `plot_vip`, `plot_pcoa`) simplify the workflow - Learned metric weights provide interpretable feature importance - Default parameters work well for most analyses - Use `n_perms = 200+` for reliable p-values in real analyses ### Next Steps - Try MeLSI on your own microbiome datasets - Use `clr_transform()` for easy CLR preprocessing - Visualize results with `plot_vip()` and `plot_pcoa()` - Adjust parameters (`B`, `m_frac`) if needed for your specific data For more information, see the main README.md file and the comprehensive documentation. ## Session Info ```{r session_info} sessionInfo() ```