Notebook

NB00 — Exploration & Feasibility¶

Purpose: Confirm Phase-A findings — the seven target anti-phage defense system families (CRISPR-Cas, R-M Type I/II, CBASS, Gabija, Retron, BREX, DISARM) are detectable at pangenome scale via interproscan_domains (primary) and eggnog_mapper_annotations (secondary), and to record baseline pangenome and gene_cluster sizes that inform the extraction strategy in NB01.

This notebook is intentionally lightweight — the substantive extraction work is in NB01. Detection Pfam accessions and per-system rules are locked in RESEARCH_PLAN.md.

In [1]:

import os
import pandas as pd
from berdl_notebook_utils.setup_spark_session import get_spark_session

spark = get_spark_session()
print("Spark version:", spark.version)

Spark version: 4.0.1

1. Load Phase-A detection feasibility table¶

Written during Phase A (see RESEARCH_PLAN.md §"Defense System Detection Rules"). This table records the raw hit counts per marker in interproscan_domains and eggnog_mapper_annotations and is our reference for interpreting the extraction in NB01.

In [2]:

feas = pd.read_csv("../data/detection_feasibility.csv")
print(f"Rows: {len(feas)}")
feas

Rows: 18

Out[2]:

	system	marker_type	marker	table	n_hits	notes
0	CRISPR-Cas	description	lower(Description) LIKE '%crispr%' OR '%cas9%'	eggnog_mapper_annotations	1012329	broad; anchor on Cas1 for specificity
1	CRISPR-Cas	pfam_accession	PF01867 (Cas1)	interproscan_domains	24666	universal CRISPR-Cas marker
2	CRISPR-Cas	pfam_accession	PF09707 (Cas2)	interproscan_domains	2695	universal marker; lower recall than Cas1
3	CRISPR-Cas	pfam_accession	PF18019 (Cas3_HD)	interproscan_domains	8053	Type I subtype marker
4	CRISPR-Cas	pfam_accession	PF22702 (Cas9_RuvC)	interproscan_domains	2318	Type II subtype marker
5	CRISPR-Cas	pfam_accession	PF07282 (Cas12f1-like)	interproscan_domains	42404	Type V subtype (includes broad TnpB homologs)
6	R-M Type II	description	lower(Description) LIKE '%type ii restriction%'	eggnog_mapper_annotations	27663	well-detected via eggNOG
7	R-M Type I	description	lower(Description) LIKE '%type i restriction%'...	eggnog_mapper_annotations	245609	well-detected via eggNOG
8	CBASS	pfam_accession	PF18178 (CD-NTase)	interproscan_domains	283	highly specific
9	CBASS	pfam_accession	PF14090 (SAVED)	interproscan_domains	2558	highly specific effector
10	CBASS	pfam_accession	PF19918 (Cap6)	interproscan_domains	155	rare effector variant
11	Gabija	pfam_accession	PF20473 (OLD_TOPRIM_C)	interproscan_domains	12941	GajA-specific
12	Gabija	pfam_accession	PF13361 (UvrD)	interproscan_domains	145561	broad helicase — use only with GajA anchor
13	Retron	pfam_accession	PF00078 (RVT_1)	interproscan_domains	72238	broad RT; needs retron-specific msr/msd context
14	BREX	pfam_accession	PF08843 (PglZ)	interproscan_domains	43816	BREX-diagnostic
15	BREX	pfam_accession	PF13175 (BrxC)	interproscan_domains	55701	BREX-diagnostic
16	BREX	pfam_accession	PF13401 (BrxL AAA_11)	interproscan_domains	41545	BREX-diagnostic
17	DISARM	pfam_accession	PF13091 (DrmC PLD)	interproscan_domains	106742	broad PLD — needs DrmB co-occurrence

2. Per-system marker summary¶

Confirms Phase A findings: all seven systems are detectable. Broad Pfams (Gabija UvrD, Retron RVT_1, DISARM PLD) are flagged for anchor+context filtering in NB01.

In [3]:

summary = (
    feas.groupby("system", as_index=False)
    .agg(markers=("marker", "count"), max_hits=("n_hits", "max"))
    .sort_values("max_hits", ascending=False)
)
summary

Out[3]:

	system	markers	max_hits
2	CRISPR-Cas	6	1012329
5	R-M Type I	1	245609
4	Gabija	2	145561
3	DISARM	1	106742
7	Retron	1	72238
0	BREX	3	55701
6	R-M Type II	1	27663
1	CBASS	3	2558

3. Confirm pangenome-side table sizes¶

Sanity-check the sizes we plan to filter against in NB01 (gene_cluster, pangenome, genome) — informs whether broadcast joins are appropriate.

In [4]:

sizes = {}
for tbl in ["gene_cluster", "pangenome", "genome", "gtdb_taxonomy_r214v1", "gtdb_metadata", "gtdb_species_clade"]:
    n = spark.sql(f"SELECT COUNT(*) FROM kbase_ke_pangenome.{tbl}").collect()[0][0]
    sizes[tbl] = n
    print(f"  {tbl:<28s}  {n:>15,}")

  gene_cluster                      132,531,501

  pangenome                              27,702

  genome                                293,059

  gtdb_taxonomy_r214v1                  293,059

  gtdb_metadata                         293,059

  gtdb_species_clade                     27,690

4. Species with sufficient genome sampling¶

Arms-race and syndrome analyses require reliable core/accessory calls, which depend on species with enough sequenced genomes. Restrict downstream analyses to species with no_genomes >= 5.

In [5]:

species_counts = spark.sql("""
    SELECT
        SUM(CASE WHEN no_genomes >= 5 THEN 1 ELSE 0 END) AS species_ge5,
        SUM(CASE WHEN no_genomes >= 10 THEN 1 ELSE 0 END) AS species_ge10,
        SUM(CASE WHEN no_genomes >= 50 THEN 1 ELSE 0 END) AS species_ge50,
        COUNT(*) AS species_total
    FROM kbase_ke_pangenome.pangenome
""").toPandas()
species_counts

Out[5]:

	species_ge5	species_ge10	species_ge50	species_total
0	7334	2812	457	27702

5. Confirm the Pfam-version-suffix issue is table-specific¶

In Phase A we discovered that bakta_pfam_domains.pfam_id includes a version suffix (e.g., PF01867.29) while interproscan_domains.signature_acc is version-free (PF01867). NB01's extraction query relies on the version-free form of signature_acc — confirm here.

In [6]:

sample = spark.sql("""
    SELECT DISTINCT signature_acc
    FROM kbase_ke_pangenome.interproscan_domains
    WHERE analysis = 'Pfam'
      AND signature_acc IN ('PF01867','PF08843','PF14090','PF20473','PF13091','PF00078')
""").toPandas()
print("Version-free Pfam accessions confirmed in interproscan_domains.signature_acc:")
sample

Version-free Pfam accessions confirmed in interproscan_domains.signature_acc:

Out[6]:

	signature_acc
0	PF20473
1	PF13091
2	PF01867
3	PF00078
4	PF08843
5	PF14090

Summary¶

All 7 defense system families are detectable.
interproscan_domains (833M rows) uses version-free Pfam accessions.
bakta_pfam_domains (18.8M rows) uses versioned accessions (PFXXXXX.Y) and is a narrower cross-check.
eggnog_mapper_annotations (93M rows) description-based classification is best for R-M and CRISPR confirmation.
Species with no_genomes >= 5 will be the analysis set for arms-race and syndromes tests.

Next: 01_extract_defense_clusters.ipynb runs the primary Pfam-based extraction and caches to data/defense_gene_clusters.tsv.gz.