---
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 implements the PERMANOVA F-statistic directly (no runtime dependency on vegan) but extends beyond fixed distance metrics by learning optimal metrics from data. This vignette uses `vegan::adonis2()` only to benchmark MeLSI against standard fixed-metric PERMANOVA.

- **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(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, eval=requireNamespace("vegan", quietly=TRUE)}
# This comparison uses vegan (a Suggested package); it is skipped if vegan
# is not installed.
# Compare with Euclidean distance
dist_euclidean <- dist(X_clr)
adonis_euclidean <- vegan::adonis2(dist_euclidean ~ y, permutations = 50)

# Compare with Bray-Curtis distance
dist_bray <- vegan::vegdist(X, method = "bray")
adonis_bray <- vegan::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

## Parallel Permutation Testing

Permutation testing dominates the runtime of a MeLSI analysis, and the
permutations are independent of one another. To run them in parallel, pass a
`BiocParallelParam` object from the `BiocParallel` package via the `BPPARAM`
argument:

```{r parallel_usage, eval=FALSE}
library(BiocParallel)

# Distribute the permutation loop across 4 worker cores (macOS/Linux)
parallel_results <- melsi(
    X_clr, y,
    n_perms = 200,
    BPPARAM = MulticoreParam(workers = 4)
)

# On Windows, use SnowParam instead:
# parallel_results <- melsi(X_clr, y, n_perms = 200,
#                           BPPARAM = SnowParam(workers = 4))
```

When `BPPARAM` is left at its default (`NULL`), permutations run sequentially.
Parallel execution is supported for omnibus, pairwise, and multi-group analyses
and produces identical results to the sequential path. This example is disabled
(`eval=FALSE`) so the vignette stays within the Bioconductor build time limit.

## 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()
```
