03 Prophage Burden
Jupyter notebook from the Pan-Bacterial Anti-Phage Defense Arsenal project.
NB03 — Prophage Burden per Species¶
Purpose: Re-derive per-species prophage burden using the classifier from
projects/prophage_ecology/src/prophage_utils.py (7 operationally defined
modules A-G: packaging, head morphogenesis, tail, lysis, integration, lysogenic
regulation, anti-defense). This gives us the "phage-side" data to correlate
against species-level defense-system counts (NB04, arms race).
Approach:
- Import
prophage_utils.build_spark_where_clause()andclassify_gene_to_module()from the sibling project. - Query
eggnog_mapper_annotationsfor all prophage-candidate gene clusters (single large OR-chain). - Classify each hit into one or more modules in Python.
- Aggregate to per-species presence of each module → per-species prophage burden = number of modules present out of 7.
Output: data/species_prophage_burden.tsv.gz.
In [1]:
import os
import sys
import pandas as pd
from berdl_notebook_utils.setup_spark_session import get_spark_session
# Import prophage_utils from the sibling prophage_ecology project
sys.path.insert(0, "../../prophage_ecology/src")
import prophage_utils
spark = get_spark_session()
DATA_DIR = "../data"
1. Sanity check: module list¶
In [2]:
print(prophage_utils.get_module_summary())
### Module A_packaging: Packaging Module - **Presence rule**: TerL alone OR Portal + TerS - **Key markers**: terminase large subunit - **Description keywords**: 9 - **PFam keywords**: 7 - **KEGG KOs**: 2 ### Module B_head_morphogenesis: Head Morphogenesis Module - **Presence rule**: MCP (major capsid protein) required - **Key markers**: major capsid protein - **Description keywords**: 8 - **PFam keywords**: 5 - **KEGG KOs**: 0 ### Module C_tail: Tail Module - **Presence rule**: >=1 structural tail protein - **Key markers**: tail tube protein, tail sheath protein, tape measure protein - **Description keywords**: 10 - **PFam keywords**: 10 - **KEGG KOs**: 0 ### Module D_lysis: Lysis Module - **Presence rule**: Holin OR endolysin - **Key markers**: holin, endolysin - **Description keywords**: 10 - **PFam keywords**: 12 - **KEGG KOs**: 0 ### Module E_integration: Integration Module - **Presence rule**: Integrase required - **Key markers**: integrase - **Description keywords**: 6 - **PFam keywords**: 6 - **KEGG KOs**: 1 ### Module F_lysogenic_regulation: Lysogenic Regulation Module - **Presence rule**: CI-like repressor required - **Key markers**: CI-like repressor - **Description keywords**: 12 - **PFam keywords**: 5 - **KEGG KOs**: 0 ### Module G_anti_defense: Anti-Defense Module - **Presence rule**: Any recognized anti-defense gene - **Key markers**: anti-CRISPR, anti-restriction - **Description keywords**: 7 - **PFam keywords**: 2 - **KEGG KOs**: 0
2. Build and run the prophage-candidate query¶
The WHERE clause matches any gene cluster whose eggNOG Description/PFAMs/
KEGG_ko field mentions any of the 7 modules' markers. It's a large OR-chain
— full-scans the 93M-row eggnog_mapper_annotations table, so budget a
couple of minutes.
In [3]:
where_clause = prophage_utils.build_spark_where_clause()
print(f"WHERE clause length: {len(where_clause):,} chars")
print("First 3 conditions:")
for line in where_clause.split("\n")[:3]:
print(" ", line.strip())
WHERE clause length: 5,377 chars First 3 conditions: LOWER(ann.Description) LIKE '%anti-cbass%' OR LOWER(ann.Description) LIKE '%anti-crispr%' OR LOWER(ann.Description) LIKE '%anti-defense%'
In [4]:
prophage_query = f"""
SELECT
ann.query_name AS gene_cluster_id,
gc.gtdb_species_clade_id,
ann.Description,
ann.PFAMs,
ann.KEGG_ko,
ann.COG_category
FROM kbase_ke_pangenome.eggnog_mapper_annotations ann
JOIN kbase_ke_pangenome.gene_cluster gc
ON ann.query_name = gc.gene_cluster_id
WHERE {where_clause}
"""
prophage_hits_df = spark.sql(prophage_query).toPandas()
print(f"Prophage-candidate hits: {len(prophage_hits_df):,} rows")
print(f"Unique gene clusters: {prophage_hits_df['gene_cluster_id'].nunique():,}")
print(f"Unique species: {prophage_hits_df['gtdb_species_clade_id'].nunique():,}")
Prophage-candidate hits: 4,005,537 rows
Unique gene clusters: 4,005,537 Unique species: 27,702
3. Classify hits into modules¶
Each gene cluster can match ≥1 module. classify_gene_to_module returns a
list of module IDs.
In [5]:
prophage_hits_df["modules"] = prophage_hits_df.apply(
lambda r: prophage_utils.classify_gene_to_module(
r["Description"], r["PFAMs"], r["KEGG_ko"], r["COG_category"]
),
axis=1,
)
# How many hits classified into at least one module?
n_classified = prophage_hits_df["modules"].apply(bool).sum()
print(f"Hits classified into at least one module: {n_classified:,} / {len(prophage_hits_df):,}")
# Explode into long-form: one row per (gene_cluster_id, module)
long = prophage_hits_df[["gene_cluster_id", "gtdb_species_clade_id", "modules"]].explode("modules")
long = long.dropna(subset=["modules"])
print(f"Exploded rows: {len(long):,}")
# Per-module hit counts
module_counts = long["modules"].value_counts()
print("\nHits per module:")
print(module_counts)
Hits classified into at least one module: 4,004,730 / 4,005,537
Exploded rows: 4,228,150 Hits per module: modules F_lysogenic_regulation 1682902 E_integration 773417 D_lysis 724499 A_packaging 707217 C_tail 189271 B_head_morphogenesis 87336 G_anti_defense 63508 Name: count, dtype: int64
4. Aggregate to per-species module presence + total burden¶
In [6]:
# Per-species per-module presence: 1 if any gene cluster in that species matches
species_module = (
long.groupby(["gtdb_species_clade_id", "modules"])
.size()
.rename("n_clusters")
.reset_index()
)
# Wide matrix: species x module → n_clusters
prophage_wide_counts = species_module.pivot_table(
index="gtdb_species_clade_id",
columns="modules",
values="n_clusters",
fill_value=0,
).astype(int)
prophage_wide_counts.columns.name = None
# Binary presence
prophage_wide_bin = (prophage_wide_counts > 0).astype(int)
prophage_wide_bin.columns = [f"prophage_{c}" for c in prophage_wide_bin.columns]
# Total prophage burden = number of modules present (0-7)
prophage_wide_bin["n_prophage_modules"] = prophage_wide_bin.sum(axis=1)
# Total prophage cluster count = sum across modules (double-counts multi-module clusters)
prophage_wide_bin["n_prophage_clusters"] = (
prophage_hits_df.groupby("gtdb_species_clade_id")["gene_cluster_id"].nunique()
).reindex(prophage_wide_bin.index, fill_value=0)
burden = prophage_wide_bin.reset_index()
print(f"Species with prophage data: {len(burden):,}")
print(f"Species with n_prophage_modules >= 4 (heavy prophage load): {(burden['n_prophage_modules'] >= 4).sum():,}")
burden.head(5)
Species with prophage data: 27,702 Species with n_prophage_modules >= 4 (heavy prophage load): 27,521
Out[6]:
| gtdb_species_clade_id | prophage_A_packaging | prophage_B_head_morphogenesis | prophage_C_tail | prophage_D_lysis | prophage_E_integration | prophage_F_lysogenic_regulation | prophage_G_anti_defense | n_prophage_modules | n_prophage_clusters | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | s__0-14-0-80-60-11_sp018897875--GB_GCA_0188978... | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 6 | 148 |
| 1 | s__0-14-3-00-41-53_sp002780895--GB_GCA_0027808... | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 4 | 95 |
| 2 | s__01-FULL-36-15b_sp001782035--GB_GCA_001782035.1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 4 | 19 |
| 3 | s__01-FULL-44-24b_sp001793235--GB_GCA_001793235.1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 4 | 29 |
| 4 | s__01-FULL-45-10b_sp001804205--GB_GCA_001804205.1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 4 | 55 |
5. Distribution of prophage burden¶
In [7]:
burden_dist = burden["n_prophage_modules"].value_counts().sort_index()
print("Distribution of prophage burden (# modules present out of 7):")
print(burden_dist)
Distribution of prophage burden (# modules present out of 7): n_prophage_modules 2 3 3 178 4 4577 5 6897 6 6389 7 9658 Name: count, dtype: int64
6. Save¶
In [8]:
out = os.path.join(DATA_DIR, "species_prophage_burden.tsv.gz")
burden.to_csv(out, sep="\t", index=False, compression="gzip")
print(f"Wrote: {out}")
print(f"Rows: {len(burden):,}")
print(f"Columns: {list(burden.columns)}")
Wrote: ../data/species_prophage_burden.tsv.gz Rows: 27,702 Columns: ['gtdb_species_clade_id', 'prophage_A_packaging', 'prophage_B_head_morphogenesis', 'prophage_C_tail', 'prophage_D_lysis', 'prophage_E_integration', 'prophage_F_lysogenic_regulation', 'prophage_G_anti_defense', 'n_prophage_modules', 'n_prophage_clusters']
Summary¶
- Per-species prophage burden derived from the
prophage_ecologyclassifier. - Wide table: 7 module presence flags +
n_prophage_modules(0-7) +n_prophage_clusters. - Species with high prophage burden (≥4 modules) are candidates for high phage pressure.
Next: 04_arms_race.ipynb — join defense matrix (NB02) with prophage
burden (NB03) and test whether defense-system count scales with prophage load.