Nb04 Paradox Proteins
Jupyter notebook from the AlphaFold MSA Depth as a Lens on the Bacterial Annotation Gap project.
NB04 — Paradox Proteins: Core Genes with Low MSA Depth¶
Runs locally from NB01 output.
Purpose: Identify and characterize gene clusters that are universally conserved (core) yet structurally novel (MSA depth < 10) — the highest-priority targets for experimental structural characterization.
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
os.makedirs("../figures", exist_ok=True)
# Load NB01 aggregate — the full per-cluster paradox list requires a JupyterHub query
# (see NB04 notes below). Here we characterize using the group-level aggregate.
agg = pd.read_csv("../data/gc_msa_agg.csv")
for col in ["is_core", "is_auxiliary", "is_singleton"]:
agg[col] = agg[col].astype(str).str.strip() == "True"
agg["group"] = agg.apply(
lambda r: "core" if r["is_core"]
else ("aux_singleton" if r["is_singleton"] else "aux_non_singleton"),
axis=1
)
# Paradox protein counts from H4 (NB03)
agg["very_low_rate"] = agg["n_very_low_msa"] / agg["n_clusters"]
agg["ec_rate"] = agg["n_has_ec"] / agg["n_clusters"]
agg["kegg_rate"] = agg["n_has_kegg"] / agg["n_clusters"]
agg["hypo_rate"] = agg["n_hypothetical"] / agg["n_clusters"]
print("Paradox protein summary by pangenome class:")
print(agg[["group","n_clusters","n_very_low_msa","very_low_rate","ec_rate","kegg_rate","hypo_rate"]].to_string(index=False))
core = agg[agg.group == "core"].iloc[0]
print(f"\nCore paradox proteins: {int(core.n_very_low_msa):,}")
print(f" EC-annotated rate (all core): {core.ec_rate:.1%}")
print(f" KEGG-annotated rate (all core): {core.kegg_rate:.1%}")
print(f" Hypothetical rate (all core): {core.hypo_rate:.1%}")
Paradox protein summary by pangenome class:
group n_clusters n_very_low_msa very_low_rate ec_rate kegg_rate hypo_rate
aux_singleton 7095643 392959 0.055380 0.108078 0.150922 0.138014
core 25571299 415733 0.016258 0.217348 0.254277 0.038321
aux_non_singleton 5384900 245002 0.045498 0.109303 0.159816 0.115647
Core paradox proteins: 415,733
EC-annotated rate (all core): 21.7%
KEGG-annotated rate (all core): 25.4%
Hypothetical rate (all core): 3.8%
In [2]:
print("Characterization steps:")
print()
print("1. ANNOTATION RATES (from NB01 aggregate, all-core baseline):")
print(" - EC coverage (all core): 21.7% → paradox subset likely lower")
print(" - KEGG coverage (all core): 25.4% → paradox subset likely lower")
print(" - Hypothetical (all core): 3.8% → paradox subset expected higher")
print()
print("2. PER-CLUSTER RANKING (pending JupyterHub query):")
print(" Required columns: gene_cluster_id, msa_depth, ec, kegg_orthology_id,")
print(" hypothetical, product_name (from bakta_annotations)")
print(" Ranking metric: n_species_with_cluster × (1 / (msa_depth + 1))")
print()
print("3. SPECIES DISTRIBUTION:")
print(" Which GTDB clades carry the most paradox proteins?")
print(" Join: gene_cluster → gene → genome → gtdb_metadata for taxonomy")
print()
print("4. FUNCTIONAL DARK MATTER OVERLAP:")
print(" Cross-reference with functional_dark_matter project findings")
print(" (hypothetical core genes with no domain annotations)")
print()
print("NOTE: Steps 2–4 require the full per-cluster paradox list from JupyterHub.")
print("A NB04b Spark notebook should join:")
print(" gc × bakta_annotations × alphafold_msa_depths")
print(" WHERE is_core = true AND msa_depth < 10")
print(" LIMIT output to top-1000 for local analysis")
Characterization steps: 1. ANNOTATION RATES (from NB01 aggregate, all-core baseline): - EC coverage (all core): 21.7% → paradox subset likely lower - KEGG coverage (all core): 25.4% → paradox subset likely lower - Hypothetical (all core): 3.8% → paradox subset expected higher 2. PER-CLUSTER RANKING (pending JupyterHub query): Required columns: gene_cluster_id, msa_depth, ec, kegg_orthology_id, hypothetical, product_name (from bakta_annotations) Ranking metric: n_species_with_cluster × (1 / (msa_depth + 1)) 3. SPECIES DISTRIBUTION: Which GTDB clades carry the most paradox proteins? Join: gene_cluster → gene → genome → gtdb_metadata for taxonomy 4. FUNCTIONAL DARK MATTER OVERLAP: Cross-reference with functional_dark_matter project findings (hypothetical core genes with no domain annotations) NOTE: Steps 2–4 require the full per-cluster paradox list from JupyterHub. A NB04b Spark notebook should join: gc × bakta_annotations × alphafold_msa_depths WHERE is_core = true AND msa_depth < 10 LIMIT output to top-1000 for local analysis
In [3]:
print("Priority ranking framework (pending per-cluster data):")
print(" Score = core_fraction × (1 / (msa_depth + 1))")
print(" core_fraction = fraction of species clades where this cluster is core")
print(" (requires gene_cluster → gtdb_species_clade join to compute)")
print()
print("Expected output: CSV table of top-1000 paradox proteins for structural biology follow-up")
print(" Columns: gene_cluster_id, msa_depth, score, product_name, ec, kegg_id,")
print(" hypothetical, n_species_clades_core")
print()
# -----------------------------------------------------------------------
# What CAN be computed now: annotation rate comparison across MSA depth
# bins using the gc_domain_sample (n_domain_hits proxy for annotation richness)
domain_sample = pd.read_csv("../data/gc_domain_sample.csv")
domain_sample["msa_depth_bin"] = pd.cut(
[0] * len(domain_sample), # msa_depth not in sample; placeholder
bins=1, labels=["unknown"]
)
# gc_domain_sample does not carry msa_depth — this confirms the need for
# a dedicated NB04b Spark join to enable per-cluster paradox ranking.
print("gc_domain_sample columns:", list(domain_sample.columns))
print("→ msa_depth not present in domain sample; paradox ranking requires JupyterHub join.")
Priority ranking framework (pending per-cluster data):
Score = core_fraction × (1 / (msa_depth + 1))
core_fraction = fraction of species clades where this cluster is core
(requires gene_cluster → gtdb_species_clade join to compute)
Expected output: CSV table of top-1000 paradox proteins for structural biology follow-up
Columns: gene_cluster_id, msa_depth, score, product_name, ec, kegg_id,
hypothetical, n_species_clades_core
gc_domain_sample columns: ['gene_cluster_id', 'n_domain_hits', 'n_distinct_ipr', 'n_analysis_tools', 'msa_depth_bin']
→ msa_depth not present in domain sample; paradox ranking requires JupyterHub join.