High-Throughput Structure Prediction Pipelines

# Download required databases
# UniRef90 (~60 GB)
wget https://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz

# MGnify (~120 GB)
wget https://ftp.ebi.ac.uk/pub/databases/metagenomics/peptide_database/current_release/

# BFD (~1.8 TB) - optional, more comprehensive
wget https://bfd.mmseqs.com/

# PDB for templates (~80 GB)
rsync -rlpt -v -z --delete rsync.rcsb.org::ftp_data/structures/divided/pdb/ pdb/

# Build search indices
# This takes 1-2 days!
mmseqs createdb uniref90.fasta uniref90_db
mmseqs createindex uniref90_db tmp

#!/usr/bin/env python3
"""
High-throughput AlphaFold2 pipeline
Processes multiple proteins in parallel
"""

import multiprocessing as mp
from pathlib import Path
import subprocess

def run_alphafold(fasta_path, output_dir, gpu_id):
    """Run AlphaFold2 for single protein"""
    cmd = [
        'python', 'run_alphafold.py',
        f'--fasta_paths={fasta_path}',
        f'--output_dir={output_dir}',
        f'--model_preset=monomer',
        f'--max_template_date=2023-01-01',
        f'--db_preset=full_dbs',
        f'--use_gpu_relax=True',
        '--num_multimer_predictions_per_model=5'
    ]

    # Set GPU device
    env = os.environ.copy()
    env['CUDA_VISIBLE_DEVICES'] = str(gpu_id)

    subprocess.run(cmd, env=env, check=True)

def process_batch(fasta_files, output_dir, num_gpus=4):
    """Process batch of proteins across multiple GPUs"""

    with mp.Pool(num_gpus) as pool:
        results = []
        for i, fasta in enumerate(fasta_files):
            gpu_id = i % num_gpus  # Round-robin GPU assignment
            result = pool.apply_async(
                run_alphafold,
                args=(fasta, output_dir, gpu_id)
            )
            results.append(result)

        # Wait for all jobs
        for result in results:
            result.get()

if __name__ == '__main__':
    fasta_files = Path('input_fastas').glob('*.fasta')
    process_batch(list(fasta_files), 'output', num_gpus=4)

import json
import numpy as np
from pathlib import Path

def extract_qc_metrics(prediction_dir):
    """Extract QC metrics from AlphaFold2 output"""

    metrics = {}

    # Load ranking JSON
    ranking_file = prediction_dir / 'ranking_debug.json'
    with open(ranking_file) as f:
        ranking = json.load(f)

    # Get best model metrics
    best_model = ranking['order'][0]
    metrics['ptm'] = ranking['ptm'][best_model]
    metrics['iptm'] = ranking.get('iptm', {}).get(best_model, None)

    # Load pLDDT scores
    plddt_file = prediction_dir / f'{best_model}_plddt.json'
    with open(plddt_file) as f:
        plddt = json.load(f)

    metrics['mean_plddt'] = np.mean(plddt)
    metrics['min_plddt'] = np.min(plddt)
    metrics['fraction_confident'] = sum(p > 80 for p in plddt) / len(plddt)

    # MSA depth
    msa_file = prediction_dir / 'msas' / 'uniclust30_hits.sto'
    if msa_file.exists():
        num_seqs = count_msa_sequences(msa_file)
        metrics['msa_depth'] = num_seqs
    else:
        metrics['msa_depth'] = 0

    return metrics

def flag_low_quality(metrics, thresholds):
    """Flag predictions that don't meet QC thresholds"""
    flags = []

    if metrics['ptm'] < thresholds['ptm']:
        flags.append('LOW_PTM')
    if metrics['mean_plddt'] < thresholds['plddt']:
        flags.append('LOW_PLDDT')
    if metrics['msa_depth'] < thresholds['msa']:
        flags.append('POOR_MSA')

    return flags if flags else ['PASS']

# Process all predictions
all_results = []
for pred_dir in Path('output').glob('*/'):
    metrics = extract_qc_metrics(pred_dir)
    flags = flag_low_quality(metrics, {
        'ptm': 0.6,
        'plddt': 70,
        'msa': 30
    })

    all_results.append({
        'protein': pred_dir.name,
        'metrics': metrics,
        'qc_flags': flags
    })

# Save QC report
import pandas as pd
df = pd.DataFrame(all_results)
df.to_csv('qc_report.csv', index=False)

print(f"Processed {len(df)} predictions")
print(f"Passed QC: {sum(df['qc_flags'] == 'PASS')}")
print(f"Low quality: {len(df) - sum(df['qc_flags'] == 'PASS')}")

# Fast MSA for screening (reduced database)
jackhmmer --cpu 8 -N 3 \
  -E 0.001 \
  --incE 0.001 \
  query.fasta \
  uniref90_reduced.fasta \
  > msa_fast.sto

# Full MSA for high-quality prediction
jackhmmer --cpu 16 -N 5 \
  -E 0.0001 \
  --incE 0.0001 \
  query.fasta \
  uniref90_full.fasta \
  > msa_full.sto

#!/usr/bin/env nextflow

// Nextflow pipeline for high-throughput AlphaFold2
// Handles dependencies, retries, and resource allocation

params.input_fastas = 'input/*.fasta'
params.output_dir = 'results'
params.num_models = 5

process generateMSA {
    cpus 16
    memory '32 GB'
    time '4h'

    input:
    path fasta

    output:
    path 'msa.a3m'

    """
    jackhmmer --cpu ${task.cpus} \\
        -N 5 \\
        ${fasta} \\
        ${params.uniref90_db} \\
        > msa.a3m
    """
}

process predictStructure {
    cpus 4
    memory '64 GB'
    time '2h'
    accelerator 1, type: 'nvidia-tesla-a100'

    input:
    path fasta
    path msa

    output:
    path 'prediction/*'

    """
    python run_alphafold.py \\
        --fasta=${fasta} \\
        --msa=${msa} \\
        --output_dir=prediction \\
        --num_models=${params.num_models}
    """
}

process qualityControl {
    cpus 2
    memory '8 GB'

    input:
    path prediction_dir

    output:
    path 'qc_report.json'

    """
    python qc_analysis.py \\
        --input=${prediction_dir} \\
        --output=qc_report.json
    """
}

workflow {
    fasta_ch = Channel.fromPath(params.input_fastas)

    msa_ch = generateMSA(fasta_ch)
    pred_ch = predictStructure(fasta_ch, msa_ch)
    qc_ch = qualityControl(pred_ch)

    qc_ch.collectFile(name: 'all_qc.json', storeDir: params.output_dir)
}

import time
import logging
from tenacity import retry, stop_after_attempt, wait_exponential

logger = logging.getLogger(__name__)

@retry(
    stop=stop_after_attempt(3),
    wait=wait_exponential(multiplier=1, min=4, max=60),
    reraise=True
)
def run_alphafold_with_retry(fasta_path, output_dir):
    """Run AlphaFold2 with automatic retries on failure"""
    try:
        run_alphafold(fasta_path, output_dir)
        return True
    except Exception as e:
        logger.error(f"Failed to predict {fasta_path}: {e}")
        # Clean up partial results
        if output_dir.exists():
            shutil.rmtree(output_dir)
        raise

def process_proteome(fasta_files, output_base):
    """Process entire proteome with error recovery"""
    results = {
        'success': [],
        'failed': [],
        'retried': []
    }

    for fasta in fasta_files:
        protein_id = fasta.stem
        output_dir = output_base / protein_id

        try:
            run_alphafold_with_retry(fasta, output_dir)
            results['success'].append(protein_id)
        except Exception as e:
            logger.error(f"Permanent failure for {protein_id}: {e}")
            results['failed'].append(protein_id)

            # Save failed protein for manual review
            failed_dir = output_base / 'failed'
            failed_dir.mkdir(exist_ok=True)
            shutil.copy(fasta, failed_dir / fasta.name)

    return results

def archive_prediction(prediction_dir, archive_base):
    """Archive completed prediction to save space"""

    # Keep only essential files
    essential_files = [
        'ranked_0.pdb',  # Best model
        'ranking_debug.json',  # Confidence metrics
        'pae_matrix.json',  # PAE data
    ]

    archive_dir = archive_base / prediction_dir.name
    archive_dir.mkdir(parents=True, exist_ok=True)

    for filename in essential_files:
        src = prediction_dir / filename
        if src.exists():
            shutil.copy(src, archive_dir)

    # Compress MSAs and features
    subprocess.run([
        'tar', 'czf',
        archive_dir / 'msas_features.tar.gz',
        '-C', prediction_dir,
        'msas/', 'features.pkl'
    ])

    # Delete original large files
    shutil.rmtree(prediction_dir / 'msas')
    (prediction_dir / 'features.pkl').unlink()

    logger.info(f"Archived {prediction_dir.name}")

# Archive all completed predictions
for pred_dir in Path('predictions').glob('*/'):
    if is_complete(pred_dir):
        archive_prediction(pred_dir, Path('archive'))

-- SQLite database schema for tracking predictions

CREATE TABLE predictions (
    protein_id TEXT PRIMARY KEY,
    sequence TEXT NOT NULL,
    length INTEGER,
    prediction_date TIMESTAMP DEFAULT CURRENT_TIMESTAMP,
    status TEXT, -- 'RUNNING', 'COMPLETED', 'FAILED'
    ptm REAL,
    mean_plddt REAL,
    msa_depth INTEGER,
    model_path TEXT,
    notes TEXT
);

CREATE TABLE qc_flags (
    protein_id TEXT,
    flag TEXT, -- 'LOW_PTM', 'POOR_MSA', etc.
    severity TEXT, -- 'WARNING', 'ERROR'
    FOREIGN KEY (protein_id) REFERENCES predictions(protein_id)
);

CREATE INDEX idx_status ON predictions(status);
CREATE INDEX idx_quality ON predictions(ptm, mean_plddt);

from flask import Flask, render_template, jsonify
import sqlite3

app = Flask(__name__)

@app.route('/dashboard')
def dashboard():
    """Web dashboard for pipeline monitoring"""
    conn = sqlite3.connect('predictions.db')

    # Get summary statistics
    stats = {
        'total': conn.execute('SELECT COUNT(*) FROM predictions').fetchone()[0],
        'completed': conn.execute("SELECT COUNT(*) FROM predictions WHERE status='COMPLETED'").fetchone()[0],
        'running': conn.execute("SELECT COUNT(*) FROM predictions WHERE status='RUNNING'").fetchone()[0],
        'failed': conn.execute("SELECT COUNT(*) FROM predictions WHERE status='FAILED'").fetchone()[0],
    }

    # Recent completions
    recent = conn.execute('''
        SELECT protein_id, ptm, mean_plddt, prediction_date
        FROM predictions
        WHERE status='COMPLETED'
        ORDER BY prediction_date DESC
        LIMIT 10
    ''').fetchall()

    conn.close()

    return render_template('dashboard.html',
                          stats=stats,
                          recent=recent)

@app.route('/api/progress')
def progress():
    """API endpoint for progress tracking"""
    # Return JSON for plotting
    conn = sqlite3.connect('predictions.db')

    data = conn.execute('''
        SELECT DATE(prediction_date) as date,
               COUNT(*) as count
        FROM predictions
        WHERE status='COMPLETED'
        GROUP BY date
        ORDER BY date
    ''').fetchall()

    conn.close()

    return jsonify({
        'dates': [row[0] for row in data],
        'counts': [row[1] for row in data]
    })

if __name__ == '__main__':
    app.run(host='0.0.0.0', port=5000)

Project: E. coli K-12 proteome (4,300 proteins)
Infrastructure: 8x A100 GPUs, 128 CPU cores

Results:
- Total time: 5 days
- Average: 2.8 hours per protein
- 4,156 successful predictions (96.7%)
- 144 failed (poor MSA or too large)

Quality distribution:
- pTM > 0.8: 3,421 (82.3%)
- pTM 0.6-0.8: 587 (14.1%)
- pTM < 0.6: 148 (3.6%)

Bottlenecks:
- MSA generation: 60% of time
- GPU inference: 35% of time
- Post-processing: 5% of time

Project: 10,000 proteins with no PDB homologs
Goal: Identify druggable targets

Pipeline:
1. Fast screening with reduced MSA (1 week)
2. Full prediction for high-quality targets (2 weeks)
3. Binding site prediction and docking (1 week)

Results:
- 8,742 completed (87.4%)
- 2,341 high-confidence structures (pTM > 0.8)
- 456 with predicted druggable pockets
- 23 selected for experimental validation

Outcomes:
- 4 validated drug targets
- 2 lead series initiated
- Multiple publications in preparation