--- title: "Metabolomics Analysis Code Examples" task: "" lineage_type: import upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-metabolomics-analysis/code_examples.md upstream_sha: e2520a96 imported_at: 2026-06-26 prompt_class: prompt upstream_changes: accepted author: upstream validated: false --- # Metabolomics Analysis Code Examples ## Phase 1: Data Import & Metabolite Identification ### Data Loading ```python def load_metabolomics_data(file_path, file_type='peak_table'): """ Load metabolomics data. file_type: 'peak_table' (CSV/TSV), 'mzml' (raw LC-MS), 'nmr' """ import pandas as pd if file_type == '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) pass ``` ### Metabolite Identification ```python def identify_metabolites(feature_data, mass_list, rt_list=None): """Match features to metabolite databases (HMDB, KEGG, PubChem).""" from tooluniverse import ToolUniverse tu = ToolUniverse() identified_metabolites = [] for i, mass in enumerate(mass_list): # Mass-based annotation via Metabolomics Workbench. NOTE: `mass` is the # observed m/z and the result `data` is a TSV string (parse it); set # `adduct` to your ionization mode (M+H positive, M-H negative). hmdb_result = tu.run_one_function({ "name": "MetabolomicsWorkbench_search_by_mz", "arguments": {"mz_value": mass, "adduct": "M+H", "tolerance": 0.005} }) if hmdb_result and hmdb_result.get('status') == 'success': rows = _parse_tsv(hmdb_result['data']) # split TSV into list[dict] if rows: identified_metabolites.append({ 'feature_id': i, 'metabolite_name': rows[0].get('Name'), 'formula': rows[0].get('Formula'), 'matched_mz': rows[0].get('Matched m/z'), 'confidence': calculate_confidence(rows[0]) }) return identified_metabolites ``` **Confidence scoring**: ``` 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 ``` ## Phase 2: Quality Control & Filtering ```python def metabolomics_qc(data, sample_metadata): """ QC metrics: CV in QC samples (<30%), blank ratios (>3x), missing values (<50%), total ion current per sample. """ qc_samples = sample_metadata['sample_type'] == 'QC' qc_data = data[qc_samples] cv_per_metabolite = qc_data.std() / qc_data.mean() blank_samples = sample_metadata['sample_type'] == 'Blank' blank_data = data[blank_samples] blank_means = blank_data.mean() sample_means = data[~blank_samples & ~qc_samples].mean() blank_ratio = sample_means / blank_means keep_metabolites = blank_ratio > 3 missing_per_metabolite = (data == 0).sum() / data.shape[0] keep_metabolites &= (missing_per_metabolite < 0.5) return data.loc[:, keep_metabolites] ``` ## Phase 3: Normalization ### Total Ion Current (TIC) ```python def normalize_tic(data): """Assumes total metabolite abundance is similar across samples.""" tic = data.sum(axis=1) median_tic = tic.median() norm_factors = median_tic / tic return data.multiply(norm_factors, axis=0) ``` ### Probabilistic Quotient Normalization (PQN) ```python def normalize_pqn(data, reference_sample=None): """More robust than TIC to large metabolite changes.""" import numpy as np if reference_sample is None: reference = data.median(axis=0) else: reference = data.loc[reference_sample] quotients = data.div(reference, axis=1) norm_factors = quotients.median(axis=1) return data.div(norm_factors, axis=0) ``` ### Internal Standard Normalization ```python def normalize_internal_standard(data, is_metabolite): """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) return normalized.drop(columns=[is_metabolite]) ``` ### Transformation ```python def transform_data(data, method='log'): """Methods: 'log' (log2), 'pareto' (mean-center/sqrt(std)), 'auto' (z-score).""" import numpy as np if method == 'log': return np.log2(data + 1) elif method == 'pareto': mean, std = data.mean(axis=0), data.std(axis=0) return (data - mean) / np.sqrt(std) elif method == 'auto': mean, std = data.mean(axis=0), data.std(axis=0) return (data - mean) / std ``` ## Phase 4: Exploratory Analysis ```python def perform_pca_metabolomics(data, sample_groups): """PCA for sample clustering and outlier detection.""" from sklearn.decomposition import PCA import matplotlib.pyplot as plt pca = PCA(n_components=2) pca_result = pca.fit_transform(data) 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 separation (better than PCA for classification).""" from sklearn.cross_decomposition import PLSRegression from sklearn.preprocessing import LabelEncoder le = LabelEncoder() y_encoded = le.fit_transform(y) pls = PLSRegression(n_components=n_components) X_pls = pls.fit_transform(X, y_encoded)[0] return X_pls ``` ## Phase 5: Differential Metabolite Analysis ```python def differential_metabolites(data, group1_samples, group2_samples): """Identify differential metabolites with FDR correction.""" from scipy import stats from statsmodels.stats.multitest import multipletests import numpy as np import pandas as pd results = [] for metabolite in data.columns: group1 = data.loc[group1_samples, metabolite] group2 = data.loc[group2_samples, metabolite] fold_change = group2.mean() / group1.mean() t_stat, p_value = stats.ttest_ind(group1, group2, equal_var=False) results.append({ 'metabolite': metabolite, 'fold_change': fold_change, 'log2FC': np.log2(fold_change), 'mean_group1': group1.mean(), 'mean_group2': group2.mean(), 'p_value': p_value, 't_statistic': t_stat }) results_df = pd.DataFrame(results) results_df['adj_p_value'] = multipletests(results_df['p_value'], method='fdr_bh')[1] results_df['significant'] = ( (results_df['adj_p_value'] < 0.05) & (np.abs(results_df['log2FC']) > 1.0) ) return results_df ``` ## Phase 6: Metabolic Pathway Analysis ```python def pathway_enrichment_metabolites(metabolite_list, organism='human'): """Perform pathway enrichment using KEGG metabolic pathways.""" from tooluniverse import ToolUniverse tu = ToolUniverse() # Metabolite-set pathway enrichment takes metabolite NAMES directly (it # resolves them to KEGG internally) -- no separate ID lookup needed. enrichment = tu.run_one_function({ "name": "MetaboAnalyst_pathway_enrichment", "arguments": {"metabolites": metabolite_list, "organism": organism} }) return enrichment ``` ## Phase 7: Multi-Omics Integration ```python def correlate_metabolite_enzyme(metabolite_data, enzyme_expression): """Correlate metabolite levels with enzyme expression (Spearman).""" from scipy.stats import spearmanr correlations = {} for metabolite in metabolite_data.columns: enzymes = find_metabolite_enzymes(metabolite) for enzyme in enzymes: if enzyme in enzyme_expression.index: r, p = spearmanr( metabolite_data[metabolite], enzyme_expression.loc[enzyme] ) correlations[f'{metabolite}_{enzyme}'] = { 'r': r, 'p': p, 'relationship': 'product' if r > 0 else 'substrate' } return correlations ```