freedomintelligence/openclaw-medical-skills/tooluniverse-proteomics-analysis

Input: MS Proteomics Data
    |
    v
Phase 1: Data Import & QC
    |-- Load MaxQuant/Spectronaut/DIA-NN output
    |-- Parse protein groups, intensities, modifications
    |-- Quality control plots (missing values, intensity distributions)
    |-- Sample correlation and PCA
    |
    v
Phase 2: Preprocessing
    |-- Filter low-confidence proteins
    |-- Handle missing values (imputation or filtering)
    |-- Log-transform intensities
    |-- Normalize across samples
    |
    v
Phase 3: Differential Expression Analysis
    |-- Statistical testing (limma, t-test, ANOVA)
    |-- Multiple testing correction (BH, Bonferroni)
    |-- Fold change calculation
    |-- Significance thresholds (p < 0.05, |log2FC| > 1)
    |
    v
Phase 4: PTM Analysis (if applicable)
    |-- Identify modified peptides
    |-- Localization probability filtering
    |-- PTM site quantification
    |-- Kinase-substrate prediction
    |-- PTM enrichment analysis
    |
    v
Phase 5: Functional Enrichment
    |-- Gene Ontology enrichment
    |-- KEGG/Reactome pathway enrichment
    |-- Protein complex enrichment (CORUM)
    |-- Tissue-specific enrichment
    |
    v
Phase 6: Protein-Protein Interactions
    |-- Query STRING for interaction networks
    |-- Identify protein complexes
    |-- Network clustering and modules
    |-- Hub protein identification
    |
    v
Phase 7: Multi-Omics Integration (optional)
    |-- Correlate with RNA-seq data
    |-- Identify translation-regulated proteins
    |-- Compare with variant/CNV data
    |-- Integrate with metabolomics
    |
    v
Phase 8: Generate Report
    |-- Summary statistics
    |-- Volcano plots and heatmaps
    |-- Pathway diagrams
    |-- Protein network visualizations
    |-- Multi-omics integration plots

def load_maxquant_proteins(protein_groups_file):
    """
    Load MaxQuant proteinGroups.txt file.
    Returns:
    - DataFrame with proteins as rows, samples as columns
    - Metadata (protein names, gene names, sequence coverage)
    """
    import pandas as pd
    # Read file
    df = pd.read_csv(protein_groups_file, sep='\t')
    # Extract intensity columns (LFQ or raw)
    intensity_cols = [col for col in df.columns if 'LFQ intensity' in col or 'Intensity ' in col]
    # Create intensity matrix
    intensity_matrix = df[intensity_cols].copy()
    intensity_matrix.columns = [col.replace('LFQ intensity ', '').replace('Intensity ', '')
                                 for col in intensity_cols]
    # Add protein metadata
    metadata = df[['Protein IDs', 'Gene names', 'Fasta headers',
                   'Peptides', 'Sequence coverage [%]']].copy()
    return intensity_matrix, metadata

def assess_missing_values(intensity_matrix):
    """
    Calculate percentage of missing values per protein and sample.
    """
    # Per protein
    missing_per_protein = (intensity_matrix == 0).sum(axis=1) / intensity_matrix.shape[1]
    # Per sample
    missing_per_sample = (intensity_matrix == 0).sum(axis=0) / intensity_matrix.shape[0]
    # Visualize
    plot_missing_value_heatmap(intensity_matrix)
    return missing_per_protein, missing_per_sample

def plot_intensity_distributions(intensity_matrix):
    """
    Plot log10 intensity distributions per sample.
    Check for consistent distributions across samples.
    """
    import matplotlib.pyplot as plt
    import numpy as np
    log_intensities = np.log10(intensity_matrix.replace(0, np.nan))
    # Boxplot per sample
    log_intensities.plot(kind='box')
    plt.ylabel('log10 Intensity')
    plt.title('Intensity Distribution per Sample')
    # Should see similar median and spread across samples

def plot_sample_correlation(intensity_matrix):
    """
    Calculate and visualize sample-sample correlation.
    Expect: High correlation within replicates, lower between conditions.
    """
    # Log-transform and remove zeros
    log_data = np.log2(intensity_matrix.replace(0, np.nan))
    # Correlation matrix
    corr_matrix = log_data.corr(method='pearson')
    # Heatmap
    import seaborn as sns
    sns.heatmap(corr_matrix, annot=True, cmap='RdYlBu_r', vmin=0.8, vmax=1.0)

def perform_pca(intensity_matrix, sample_groups):
    """
    Principal component analysis for sample clustering.
    """
    from sklearn.decomposition import PCA
    # Prepare data (log, impute, scale)
    log_data = np.log2(intensity_matrix.replace(0, np.nan))
    # Simple imputation with minimum value
    imputed = log_data.fillna(log_data.min().min())
    # PCA
    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(imputed.T)
    # Plot with group colors
    plt.scatter(pca_result[:, 0], pca_result[:, 1], c=sample_groups)
    plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%})')
    plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%})')

def filter_proteins(intensity_matrix, metadata, min_valid=3):
    """
    Filter out low-confidence proteins.
    Criteria:
    - At least 2 unique peptides (from metadata)
    - At least min_valid samples with detected intensity
    - Remove contaminants and reverse sequences
    """
    # Filter by peptide count
    valid_proteins = metadata['Peptides'] >= 2
    # Filter by detection in samples
    n_detected = (intensity_matrix > 0).sum(axis=1)
    valid_detection = n_detected >= min_valid
    # Remove contaminants (from MaxQuant)
    is_contaminant = metadata['Protein IDs'].str.contains('CON__', na=False)
    is_reverse = metadata['Protein IDs'].str.contains('REV__', na=False)
    # Combined filter
    keep = valid_proteins & valid_detection & ~is_contaminant & ~is_reverse
    return intensity_matrix[keep], metadata[keep]

def impute_missing_values(intensity_matrix, method='MinProb'):
    """
    Impute missing protein intensities.
    Methods:
    - MinProb: Random from minimum observed + normal noise (for MNAR)
    - KNN: K-nearest neighbors imputation
    - QRILC: Quantile regression-based imputation
    """
    if method == 'MinProb':
        # Assume missing = low abundance (MNAR assumption)
        min_val = intensity_matrix[intensity_matrix > 0].min().min()
        width = 0.3  # Standard deviation of noise
        shift = 1.8  # Downshift from minimum
        # Replace zeros with random low values
        imputed = intensity_matrix.copy()
        missing_mask = imputed == 0
        n_missing = missing_mask.sum().sum()
        random_vals = np.random.normal(
            loc=min_val - shift,
            scale=width,
            size=n_missing
        )
        imputed.values[missing_mask.values] = random_vals
        return imputed
    elif method == 'KNN':
        from sklearn.impute import KNNImputer
        imputer = KNNImputer(n_neighbors=5)
        imputed = pd.DataFrame(
            imputer.fit_transform(intensity_matrix.replace(0, np.nan)),
            index=intensity_matrix.index,
            columns=intensity_matrix.columns
        )
        return imputed

def normalize_intensities(intensity_matrix, method='median'):
    """
    Normalize protein intensities across samples.
    Methods:
    - median: Divide by median intensity per sample
    - quantile: Quantile normalization (same distribution)
    - TMM: Trimmed mean of M-values (from edgeR)
    - VSN: Variance-stabilizing normalization
    """
    if method == 'median':
        # Median normalization
        medians = intensity_matrix.median(axis=0)
        global_median = medians.median()
        norm_factors = global_median / medians
        normalized = intensity_matrix * norm_factors
        return normalized
    elif method == 'quantile':
        # Quantile normalization
        from sklearn.preprocessing import quantile_transform
        normalized = pd.DataFrame(
            quantile_transform(intensity_matrix, axis=1),
            index=intensity_matrix.index,
            columns=intensity_matrix.columns
        )
        return normalized

def differential_expression_limma(log2_intensities, group1_samples, group2_samples):
    """
    Perform differential expression using limma-like approach.
    Returns:
    - log2 fold changes
    - p-values
    - adjusted p-values (BH)
    """
    from scipy import stats
    results = []
    for protein in log2_intensities.index:
        # Extract intensities for each group
        group1 = log2_intensities.loc[protein, group1_samples]
        group2 = log2_intensities.loc[protein, group2_samples]
        # Calculate statistics
        mean1 = group1.mean()
        mean2 = group2.mean()
        log2fc = mean2 - mean1
        # t-test
        t_stat, p_value = stats.ttest_ind(group1, group2, equal_var=False)