Input: Metabolomics Data (Peak Table or Spectra)
    |
    v
Phase 1: Data Import & Metabolite Identification
    |-- Load peak table or process raw spectra
    |-- Match features to metabolite databases (HMDB, KEGG)
    |-- Annotate with chemical IDs, formulas, pathways
    |-- Confidence scoring for IDs
    |
    v
Phase 2: Quality Control & Filtering
    |-- Assess peak quality (CV, blank ratios)
    |-- Remove background peaks
    |-- Filter low-quality metabolites
    |-- Internal standard check
    |
    v
Phase 3: Normalization
    |-- Sample-wise normalization (TIC, PQN, internal standards)
    |-- Batch effect correction
    |-- Log-transform or scaling
    |
    v
Phase 4: Exploratory Analysis
    |-- PCA for sample clustering
    |-- Quality assessment plots
    |-- Outlier detection
    |-- Sample correlation
    |
    v
Phase 5: Differential Analysis
    |-- Statistical testing (t-test, ANOVA, Wilcoxon)
    |-- Fold change calculation
    |-- Multiple testing correction
    |-- Volcano plots, heatmaps
    |
    v
Phase 6: Pathway Analysis
    |-- Metabolite set enrichment (MSEA)
    |-- Pathway topology analysis
    |-- KEGG/Reactome pathway mapping
    |-- Identify dysregulated pathways
    |
    v
Phase 7: Multi-Omics Integration
    |-- Correlate with enzyme expression (RNA/protein)
    |-- Metabolite-gene associations
    |-- Pathway-level integration
    |-- Metabolic flux inference
    |
    v
Phase 8: Generate Report
    |-- Summary statistics
    |-- Differential metabolites
    |-- Pathway diagrams
    |-- Multi-omics integration plots
    |-- Biomarker panel

Sample_ID | Glucose | Lactate | Glutamine | ... | Cholesterol
----------|---------|---------|-----------|-----|------------
Control_1 | 125000  | 45000   | 78000     | ... | 23000
Control_2 | 130000  | 43000   | 82000     | ... | 25000
Disease_1 | 85000   | 92000   | 45000     | ... | 45000

def load_metabolomics_data(file_path, file_type='peak_table'):
    """
    Load metabolomics data.

    file_type options:
    - 'peak_table': CSV/TSV with metabolites as columns
    - 'mzml': Raw LC-MS data (requires processing)
    - 'nmr': NMR spectra
    """
    import pandas as pd

    if file_type == 'peak_table':
        # Load peak table
        data = pd.read_csv(file_path, index_col=0)
        # Rows = samples, Columns = metabolites
        return data

    elif file_type == 'mzml':
        # Process raw MS data (requires pymzml or similar)
        # Peak detection, alignment, quantification
        pass

def identify_metabolites(feature_data, mass_list, rt_list=None):
    """
    Match features to metabolite databases.

    Uses accurate mass and retention time (if available).
    Queries: HMDB, KEGG Compound, PubChem
    """
    from tooluniverse import ToolUniverse
    tu = ToolUniverse()

    identified_metabolites = []

    for i, mass in enumerate(mass_list):
        # Query HMDB by mass (±5 ppm tolerance)
        hmdb_result = tu.run_one_function({
            "name": "hmdb_search_by_mass",
            "arguments": {
                "mass": mass,
                "mass_tolerance": 0.005  # 5 ppm
            }
        })

        if hmdb_result and 'data' in hmdb_result:
            matches = hmdb_result['data']
            # Rank by confidence
            # (exact mass match, pathway context, etc.)
            identified_metabolites.append({
                'feature_id': i,
                'metabolite_name': matches[0]['name'],
                'hmdb_id': matches[0]['accession'],
                'formula': matches[0]['chemical_formula'],
                'confidence': calculate_confidence(matches[0])
            })

    return identified_metabolites

Level 1: Confirmed with authentic standard (MS + RT match)
Level 2: Probable structure (accurate mass + MS/MS)
Level 3: Tentative match (accurate mass only)
Level 4: Unknown metabolite

def metabolomics_qc(data, sample_metadata):
    """
    Quality control for metabolomics data.

    QC metrics:
    - Coefficient of variation (CV) in QC samples
    - Blank ratios (signal in samples vs blanks)
    - Missing values per metabolite
    - Total ion current per sample
    """
    # 1. CV in QC samples (should be < 30%)
    qc_samples = sample_metadata['sample_type'] == 'QC'
    qc_data = data[qc_samples]

    cv_per_metabolite = qc_data.std() / qc_data.mean()
    high_cv = cv_per_metabolite > 0.3

    print(f"Metabolites with CV > 30%: {high_cv.sum()}")

    # 2. Blank subtraction
    blank_samples = sample_metadata['sample_type'] == 'Blank'
    blank_data = data[blank_samples]

    # Calculate blank ratios
    blank_means = blank_data.mean()
    sample_means = data[~blank_samples & ~qc_samples].mean()
    blank_ratio = sample_means / blank_means

    # Filter: keep metabolites with sample/blank > 3
    keep_metabolites = blank_ratio > 3

    # 3. Missing values (remove if >50% missing)
    missing_per_metabolite = (data == 0).sum() / data.shape[0]
    keep_metabolites &= (missing_per_metabolite < 0.5)

    # Filter data
    filtered_data = data.loc[:, keep_metabolites]

    return filtered_data

def normalize_tic(data):
    """
    Normalize by total ion current.
    Assumes total metabolite abundance is similar across samples.
    """
    tic = data.sum(axis=1)
    median_tic = tic.median()
    norm_factors = median_tic / tic
    normalized = data.multiply(norm_factors, axis=0)
    return normalized

def normalize_pqn(data, reference_sample=None):
    """
    Probabilistic quotient normalization.
    More robust than TIC to large metabolite changes.
    """
    import numpy as np

    # Use median sample as reference
    if reference_sample is None:
        reference = data.median(axis=0)
    else:
        reference = data.loc[reference_sample]

    # Calculate quotients
    quotients = data.div(reference, axis=1)

    # Median quotient per sample
    norm_factors = quotients.median(axis=1)

    # Normalize
    normalized = data.div(norm_factors, axis=0)

    return normalized

def normalize_internal_standard(data, is_metabolite):
    """
    Normalize by spiked-in internal standard.
    Most accurate if added before sample processing.
    """
    is_abundance = data[is_metabolite]
    norm_factors = is_abundance.median() / is_abundance
    normalized = data.multiply(norm_factors, axis=0)

    # Remove internal standard from data
    normalized = normalized.drop(columns=[is_metabolite])

    return normalized

def transform_data(data, method='log'):
    """
    Transform metabolite abundances.

    Methods:
    - 'log': log2 transform (stabilize variance)
    - 'pareto': Pareto scaling (mean-center, divide by sqrt(std))
    - 'auto': Auto-scaling (mean-center, divide by std)
    """
    import numpy as np

    if method == 'log':
        # Add small constant to avoid log(0)
        transformed = np.log2(data + 1)

    elif method == 'pareto':
        # Pareto scaling (common in metabolomics)
        mean = data.mean(axis=0)
        std = data.std(axis=0)
        transformed = (data - mean) / np.sqrt(std)

    elif method == 'auto':
        # Auto-scaling (z-score)
        mean = data.mean(axis=0)
        std = data.std(axis=0)
        transformed = (data - mean) / std

    return transformed

def perform_pca_metabolomics(data, sample_groups):
    """
    Principal component analysis for sample clustering.
    """
    from sklearn.decomposition import PCA
    import matplotlib.pyplot as plt

    # PCA
    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(data)

    # Plot
    plt.figure(figsize=(8, 6))
    for group in sample_groups.unique():
        mask = sample_groups == group
        plt.scatter(pca_result[mask, 0], pca_result[mask, 1], label=group)

    plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%})')
    plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%})')
    plt.legend()
    plt.title('PCA - Metabolomics Data')

def plsda_analysis(X, y, n_components=2):
    """
    PLS-DA for supervised dimensionality reduction.
    Better separation than PCA for classification tasks.
    """
    from sklearn.cross_decomposition import PLSRegression
    from sklearn.preprocessing import LabelEncoder

    # Encode labels
    le = LabelEncoder()
    y_encoded = le.fit_transform(y)

    # PLS
    pls = PLSRegression(n_components=n_components)
    X_pls = pls.fit_transform(X, y_encoded)[0]

    # Plot
    plt.scatter(X_pls[:, 0], X_pls[:, 1], c=y_encoded)
    plt.xlabel('PLS Component 1')
    plt.ylabel('PLS Component 2')
    plt.title('PLS-DA')

    return X_pls

def differential_metabolites(data, group1_samples, group2_samples):
    """
    Identify differential metabolites between two groups.
    """
    from scipy import stats
    import numpy as np

    results = []

    for metabolite in data.columns:
        # Extract abundances
        group1 = data.loc[group1_samples, metabolite]
        group2 = data.loc[group2_samples, metabolite]

        # Statistics
        mean1 = group1.mean()
        mean2 = group2.mean()
        fold_change = mean2 / mean1
        log2fc = np.log2(fold_change)

        # t-test
        t_stat, p_value = stats.ttest_ind(group1, group2, equal_var=False)

        results.append({
            'metabolite': metabolite,
            'fold_change': fold_change,
            'log2FC': log2fc,
            'mean_group1': mean1,
            'mean_group2': mean2,
            'p_value': p_value,
            't_statistic': t_stat
        })

    results_df = pd.DataFrame(results)

    # Multiple testing correction
    from statsmodels.stats.multitest import multipletests
    results_df['adj_p_value'] = multipletests(results_df['p_value'], method='fdr_bh')[1]

    # Significance
    results_df['significant'] = (
        (results_df['adj_p_value'] < 0.05) &
        (np.abs(results_df['log2FC']) > 1.0)
    )

    return results_df

def plot_metabolite_volcano(de_results):
    """Visualize differential metabolite results."""
    plt.figure(figsize=(8, 6))

    # Non-significant
    non_sig = de_results[~de_results['significant']]
    plt.scatter(non_sig['log2FC'], -np.log10(non_sig['p_value']),
                c='gray', alpha=0.5, s=20)

    # Significant
    sig = de_results[de_results['significant']]
    plt.scatter(sig['log2FC'], -np.log10(sig['p_value']),
                c='red', alpha=0.7, s=30)

    plt.axhline(-np.log10(0.05), color='blue', linestyle='--')
    plt.axvline(-1, color='blue', linestyle='--')
    plt.axvline(1, color='blue', linestyle='--')

    plt.xlabel('log2 Fold Change')
    plt.ylabel('-log10(p-value)')
    plt.title('Differential Metabolites')

def pathway_enrichment_metabolites(metabolite_list, organism='human'):
    """
    Perform pathway enrichment for differential metabolites.

    Uses KEGG metabolic pathways.
    """
    from tooluniverse import ToolUniverse
    tu = ToolUniverse()

    # Get KEGG compound IDs for metabolites
    kegg_ids = []
    for metabolite in metabolite_list:
        # Convert HMDB ID to KEGG Compound ID
        result = tu.run_one_function({
            "name": "kegg_find_compound",
            "arguments": {"query": metabolite}
        })
        if result and 'data' in result:
            kegg_ids.append(result['data'][0]['entry_id'])

    # Pathway enrichment
    enrichment = tu.run_one_function({
        "name": "kegg_enrich_pathway",
        "arguments": {
            "compound_list": ",".join(kegg_ids),
            "organism": organism
        }
    })

    return enrichment

def pathway_topology_analysis(metabolites, pathway_id):
    """
    Analyze pathway dysregulation considering topology.

    Metabolites at key pathway positions (hubs, bottlenecks)
    have more impact than peripheral metabolites.
    """
    # Load pathway structure from KEGG
    # Calculate impact score based on:
    # - Betweenness centrality (bottleneck metabolites)
    # - Degree centrality (hub metabolites)
    # - Pathway position (early vs late)

    pass

def correlate_metabolite_enzyme(metabolite_data, enzyme_expression):
    """
    Correlate metabolite levels with enzyme expression.

    Expected correlations:
    - Substrate + enzyme → negative correlation (consumption)
    - Product + enzyme → positive correlation (production)
    """
    from scipy.stats import spearmanr

    # For each metabolite-enzyme pair
    correlations = {}

    for metabolite in metabolite_data.columns:
        # Find enzymes that produce/consume this metabolite
        enzymes = find_metabolite_enzymes(metabolite)

        for enzyme in enzymes:
            if enzyme in enzyme_expression.index:
                met_levels = metabolite_data[metabolite]
                enz_expr = enzyme_expression.loc[enzyme]

                r, p = spearmanr(met_levels, enz_expr)

                correlations[f'{metabolite}_{enzyme}'] = {
                    'r': r,
                    'p': p,
                    'relationship': 'product' if r > 0 else 'substrate'
                }

    return correlations

def integrate_omics_pathway(metabolite_fc, gene_fc, pathway_id):
    """
    Integrate metabolite and gene fold changes at pathway level.

    For each reaction:
    - Check if metabolites are changed
    - Check if enzymes are changed
    - Score pathway dysregulation (combined evidence)
    """
    # Load pathway reactions
    # For each reaction:
    #   - Metabolite substrate/product changes
    #   - Enzyme expression changes
    #   - Concordance score

    pathway_score = calculate_pathway_dysregulation(
        metabolite_fc, gene_fc, pathway_id
    )

    return pathway_score

undefined

---

undefined