Complete prokaryotic stalling peptide discovery pipeline

[jayhesselberth]

Overview

Complete the implementation of the prokaryotic stalling peptide discovery pipeline. The workflow rules are defined in workflow/rules/prokaryotic.smk but only 8 out of 10 required scripts are implemented.

Background

The prokaryotic pipeline discovers ribosomal stalling peptides in bacteria, archaea, and phages using three complementary approaches:

1. Seed-based discovery - Start with 47 known stalling peptides (PMID 38565864), build HMMs, search proteomes
2. Domain-guided discovery - Extract inter-domain GP motifs (hypothesis: stalling occurs at domain boundaries)
3. Comprehensive GP discovery - Extract all GP motifs as unbiased control

Known stalling peptides include:

• SecM (E. coli): FSTPVWISQAQGIRAGP - monitors secretion
• TnaC (E. coli): WDP...IGP - regulates tryptophanase
• MifM (B. subtilis): ...GPXXGIAGP - monitors membrane insertion
• CydA and 43 others from systematic screens

See workflow/PROKARYOTIC-DISCOVERY.md for complete documentation.

Current Status

Implemented Scripts (8/10) ✅

• ✅ cluster_gp_motifs.py - Cluster GP motifs by sequence context
• ✅ expand_hits.py - Expand HMM hits (general utility)
• ✅ extract_gp_motifs.py - Extract GP-containing sequences from proteomes
• ✅ filter_alignment.py - Filter alignments (general utility)
• ✅ filter_interdomain_gp.py - Classify GP positions relative to domains
• ✅ merge_alignments.py - Merge alignments (general utility)
• ✅ parse_stalling_peptides.py - Parse known stalling peptides
• ✅ split_peptides_by_motif.py - Split peptides by motif type (RAGP, QAPP, etc.)

Missing Scripts (2/10) ❌

• ❌ annotate_gp_domains.py - Run Pfam/InterProScan on GP-containing proteins
• ❌ analyze_gp_conservation.py - Analyze conservation within clusters, generate logos
• ❌ identify_consensus_patterns.py - Identify consensus patterns from clusters
• ❌ validate_stalling_peptides.py - Validate discovered motifs against known peptides
• ❌ compare_gp_approaches.py - Compare all-GP vs inter-domain approaches
• ❌ prokaryotic_discovery_report.qmd - Generate comprehensive Quarto report

Implementation Plan

Phase 1: Domain Annotation Script

File: workflow/scripts/annotate_gp_domains.py

Purpose: Annotate domains in GP-containing proteins using Pfam

Implementation:

# Use hmmscan to search sequences against Pfam database
# Input: FASTA of GP-containing sequences
# Output: TSV with domain annotations (name, start, end, evalue)

Key features:

• Run hmmscan with Pfam database
• Parse domain hits (gathering thresholds)
• Record domain boundaries for each sequence
• Link to GP positions for inter-domain classification

Dependencies:

• HMMER (hmmscan) via conda environment
• Pfam database (download via rule)
• Biopython for parsing

Resources:

• 8 CPUs (parallel hmmscan)
• 16GB RAM
• 4-8 hours runtime for bacterial proteomes

---

Phase 2: Conservation Analysis Script

File: workflow/scripts/analyze_gp_conservation.py

Purpose: Analyze conservation patterns within each cluster

Implementation:

# For each cluster:
#   1. Align sequences (MUSCLE/MAFFT)
#   2. Calculate per-position conservation
#   3. Identify conserved motifs
#   4. Generate sequence logo via Skylign API (see issue #12)

Inputs:

• gp_clusters.tsv.gz - Cluster assignments
• interdomain_gp_motifs.tsv.gz - GP sequences

Outputs:

• cluster_conservation.tsv.gz - Conservation scores per position per cluster
• logos/cluster_{1..20}.png - Sequence logos for top 20 clusters

Key features:

• Per-cluster multiple sequence alignment
• Shannon entropy for conservation scoring
• Identification of conserved positions
• Integration with Skylign API for logo generation
• Statistics on cluster size, conservation, consensus sequence

Dependencies:

• MUSCLE or MAFFT for alignment
• Biopython for alignment I/O
• requests for Skylign API calls
• Logomaker or weblogo as fallback

---

Phase 3: Consensus Pattern Identification

File: workflow/scripts/identify_consensus_patterns.py

Purpose: Identify consensus patterns for major clusters

Implementation:

# 1. Load conservation data for all clusters
# 2. Filter clusters by size and conservation
# 3. Generate consensus sequences
# 4. Create Stockholm alignments for HMM building
# 5. Classify by motif type (RAGP, QAPP, etc.)

Inputs:

• cluster_conservation.tsv.gz - Conservation data
• gp_clusters.tsv.gz - Cluster membership

Outputs:

• consensus_patterns.tsv - Consensus sequences and statistics
• cluster_{1..20}.sto - Stockholm alignments for top clusters

Key features:

• Consensus calling (>60% conservation)
• Pattern classification (RAGP, RAPG, QAPP, HGPP, etc.)
• Quality metrics (cluster size, conservation, entropy)
• Stockholm format for downstream HMM building

Filtering criteria:

• Minimum cluster size: 5 sequences
• Minimum conservation: 60% at GP position
• Focus on top 20 largest clusters

---

Phase 4: Validation Script

File: workflow/scripts/validate_stalling_peptides.py

Purpose: Validate discovered motifs against known stalling peptides

Implementation:

# 1. Load all cluster HMMs
# 2. Search known stalling peptides (resources/stalling-peptides/)
# 3. Record hits to each known peptide
# 4. Generate validation summary

Inputs:

• cluster_{1..20}.hmm - Discovered HMMs
• known_stalling_peptides.fasta - 47 known peptides

Outputs:

• known_peptide_hits.tsv - Hits to each known peptide
• validation_summary.txt - Summary of recovery rate

Key features:

• Run hmmsearch for each HMM against known peptides
• Track which known peptides are recovered
• Identify closest matching cluster for each known peptide
• Calculate sensitivity (% known peptides recovered)

Success criteria:

• Should recover SecM, TnaC, MifM, CydA
• Should identify motif families (RAGP, QAPP, etc.)
• High-quality clusters should match known patterns

---

Phase 5: Comparative Analysis Script

File: workflow/scripts/compare_gp_approaches.py

Purpose: Compare all-GP vs inter-domain GP approaches

Implementation:

# 1. Load all GP motifs
# 2. Load inter-domain GP motifs  
# 3. Load cluster assignments
# 4. Load validation results
# 5. Statistical comparison
# 6. Generate comparison plots

Inputs:

• all_gp_motifs.tsv.gz - All GP motifs
• interdomain_gp_motifs.tsv.gz - Inter-domain only
• gp_clusters.tsv.gz - Clustering results
• known_peptide_hits.tsv - Validation results

Outputs:

• approach_comparison.tsv - Statistics for each approach
• comparison_plots/ - Visualizations

Analyses:

1. Count statistics:

   • Total GP motifs found
   • Inter-domain vs intra-domain vs unstructured
   • Cluster size distributions

2. Enrichment analysis:

   • Are known stalling peptides enriched in inter-domain?
   • Conservation differences between categories

3. Visualizations:

   • Venn diagram of approaches
   • Position distribution (N-term, middle, C-term)
   • Motif type frequencies
   • Conservation heatmaps

---

Phase 6: Quarto Report

File: workflow/scripts/prokaryotic_discovery_report.qmd

Purpose: Generate comprehensive HTML report

Sections:

1. Executive Summary

   • Total sequences analyzed
   • GP motifs found
   • Clusters identified
   • Known peptides recovered

2. Approach Comparison

   • All-GP vs inter-domain statistics
   • Domain annotation summary
   • Position distributions

3. Pattern Discovery

   • Top 20 clusters with logos
   • Consensus sequences
   • Conservation plots
   • Motif classification

4. Validation Results

   • Known peptide recovery
   • Comparison to literature
   • Novel vs known patterns

5. Candidate Stalling Peptides

   • High-confidence candidates (large clusters, high conservation)
   • Novel patterns not matching known peptides
   • Recommendations for experimental validation

6. Methods

   • Databases searched
   • Parameters used
   • Clustering algorithm
   • HMM building approach

Visualizations:

• Sequence logos for top clusters (from Skylign API)
• Conservation heatmaps
• Cluster size distributions
• Domain junction enrichment plots
• Phylogenetic/taxonomic distribution

---

Dependencies

Conda Environments

Update workflow/envs/python.yaml:

dependencies:
  - python=3.11
  - biopython
  - pandas
  - click
  - logomaker  # For sequence logos
  - requests   # For Skylign API
  - matplotlib
  - seaborn

Update workflow/envs/r-quarto.yaml:

dependencies:
  - r-base=4.3
  - r-tidyverse
  - r-knitr
  - quarto

External Tools

• Pfam database: Add download rule to prokaryotic.smk
• MUSCLE/MAFFT: Already in HMMER environment
• Skylign API: No installation needed (web service)

---

Testing Strategy

Unit Tests

Test each script independently:

# Test GP extraction
python workflow/scripts/extract_gp_motifs.py \
  --fasta resources/stalling-peptides/known_stalling_peptides.fasta \
  --motifs-out test_motifs.tsv \
  --sequences-out test_seqs.fasta

# Test clustering
python workflow/scripts/cluster_gp_motifs.py \
  --motifs test_motifs.tsv \
  --clusters test_clusters.tsv

Integration Tests

Test workflow segments:

# Test seed-based discovery (fast)
snakemake --configfile workflow/config-prokaryotic.yaml \
  --use-conda -n \
  search_with_comprehensive_hmm

# Test full workflow (dry-run)
snakemake --configfile workflow/config-prokaryotic.yaml \
  --use-conda -n \
  prokaryotic_discovery

Validation Tests

Known peptide recovery:

• Must identify SecM (RAGP family)
• Must identify TnaC (GP family)
• Should find MifM
• Target: >80% recovery of 47 known peptides

Pattern validation:

• RAGP motifs should cluster together
• QAPP motifs should cluster together
• Should identify 10-20 major clusters

---

Implementation Order

1. ✅ Phase 0: Fix wildcard mismatch in analyze_cluster_conservation rule (DONE - commit 4f2b897)

2. Phase 1: Implement annotate_gp_domains.py (enables domain-based filtering)

3. Phase 2: Implement analyze_gp_conservation.py (enables pattern analysis)

   • Initially without logos (generate text-based consensus)
   • Later integrate Skylign API (issue Add Skylign sequence logo generation to pipeline #12)

4. Phase 3: Implement identify_consensus_patterns.py (enables HMM building)

5. Phase 4: Implement validate_stalling_peptides.py (quality control)

6. Phase 5: Implement compare_gp_approaches.py (comparative analysis)

7. Phase 6: Implement prokaryotic_discovery_report.qmd (final reporting)

8. Phase 7: End-to-end testing with bacterial proteomes

9. Phase 8: Documentation and README updates

---

Expected Results

Pattern Classes

Based on known peptides, expect:

1. RAGP family (SecM-like)

   • Strong RAGP conservation
   • Aromatic upstream (W, F)
   • Largest cluster expected

2. RAPG family (broader)

   • R at -3 position
   • Moderately conserved

3. QAPP family (B. subtilis)

   • Q at -3 position
   • Multiple clusters expected

4. HGPP/QGPP families (rare)

   • Smaller clusters
   • Specific organisms

5. Novel patterns

   • GP variants not yet characterized
   • Validate experimentally

Success Metrics

• ✅ Recover >80% of 47 known stalling peptides
• ✅ Identify 10-20 major clusters
• ✅ Show inter-domain enrichment for known peptides
• ✅ Generate publication-quality sequence logos
• ✅ Comprehensive report with recommendations

---

Resources Required

Computational

• Local testing: 12 cores, 32GB RAM
• SLURM cluster:
  • Domain annotation: 8 cores, 16GB, 4-8 hours
  • Clustering: 8 cores, 16GB, 2-4 hours
  • Other steps: 2-4 cores, 8GB, <1 hour

Storage

• Bacterial proteomes: ~5GB
• Pfam database: ~5GB
• Results: ~2GB (clusters, alignments, logs)

Time Estimate

• Script development: 2-3 days per script
• Testing and debugging: 1-2 days
• Documentation: 1 day
• Total: 2-3 weeks for complete implementation

---

Related Issues

• Add Skylign sequence logo generation to pipeline #12 - Skylign sequence logo generation (will integrate with Phase 2)

---

References

• Wilson, D.N., et al. (2016). Ribosome-targeting antibiotics and mechanisms of bacterial resistance. Nat. Rev. Microbiol.
• Ito, K., et al. (2010). Nascentome analysis uncovers novel translational arrest sequences. Mol. Cell.
• Woolstenhulme, C.J., et al. (2013). Nascent peptides that block protein synthesis in bacteria. PNAS.
• PMID 38565864 - Source of 47 known stalling peptides