Nb04B Paradox Ranking
Jupyter notebook from the AlphaFold MSA Depth as a Lens on the Bacterial Annotation Gap project.
NB04b — Paradox Protein Ranking: Core Genes with Low MSA Depth¶
Requires: BERDL JupyterHub (Spark)
Purpose: Extract and rank all core gene clusters with msa_depth < 10 (paradox proteins). Joins gene_cluster × bakta_annotations × alphafold_msa_depths to get product names, functional annotation flags, and MSA depth for each paradox cluster.
Output:
data/paradox_top1000.csv— top-1000 ranked paradox proteins (lowest msa_depth first)data/paradox_summary.csv— annotation stats for the full paradox set
In [1]:
from berdl_notebook_utils.setup_spark_session import get_spark_session
spark = get_spark_session()
print("Spark session ready:", spark.version)
Spark session ready: 4.0.1
In [2]:
from pyspark.sql.functions import col, regexp_replace, count, avg, sum as spark_sum
# --- gene_cluster table: filter to core only ---
gc_core = spark.table("kbase_ke_pangenome.gene_cluster").filter(
"is_core = true"
).select("gene_cluster_id", "gtdb_species_clade_id")
print(f"Core gene clusters: {gc_core.count():,}")
Core gene clusters: 62,062,686
In [3]:
# --- bakta_annotations: keep product name + annotation flags ---
bakta = spark.table("kbase_ke_pangenome.bakta_annotations").select(
"gene_cluster_id",
"hypothetical",
"ec",
"kegg_orthology_id",
"product",
"uniref100"
).filter(
"uniref100 IS NOT NULL AND uniref100 NOT LIKE 'UniRef100_UPI%'"
).withColumn(
"uniprot_accession",
regexp_replace(col("uniref100"), "UniRef100_", "")
)
print(f"Bakta rows with real UniProt accession: {bakta.count():,}")
Bakta rows with real UniProt accession: 38,804,903
In [4]:
# --- AlphaFold MSA depths ---
msa = spark.table("kescience_alphafold.alphafold_msa_depths").select(
"uniprot_accession", col("msa_depth").cast("double")
).filter("msa_depth < 10") # push filter early — only paradox-depth proteins
print(f"AlphaFold entries with msa_depth < 10: {msa.count():,}")
AlphaFold entries with msa_depth < 10: 12,748,556
In [5]:
# --- Join: core × bakta × msa (msa_depth < 10 already filtered) ---
paradox = (
gc_core
.join(bakta, on="gene_cluster_id", how="inner")
.join(msa, on="uniprot_accession", how="inner")
.select(
"gene_cluster_id",
"gtdb_species_clade_id",
"uniprot_accession",
"msa_depth",
"hypothetical",
"ec",
"kegg_orthology_id",
"product"
)
)
print(f"Paradox protein rows (core + msa_depth < 10): {paradox.count():,}")
paradox.printSchema()
paradox.show(5, truncate=True)
Paradox protein rows (core + msa_depth < 10): 415,733
root |-- gene_cluster_id: string (nullable = true) |-- gtdb_species_clade_id: string (nullable = true) |-- uniprot_accession: string (nullable = true) |-- msa_depth: double (nullable = true) |-- hypothetical: boolean (nullable = true) |-- ec: string (nullable = true) |-- kegg_orthology_id: string (nullable = true) |-- product: string (nullable = true)
+--------------------+---------------------+-----------------+---------+------------+----+-----------------+--------------------+ | gene_cluster_id|gtdb_species_clade_id|uniprot_accession|msa_depth|hypothetical| ec|kegg_orthology_id| product| +--------------------+---------------------+-----------------+---------+------------+----+-----------------+--------------------+ | DJIV01000140.1_6| s__Accumulibacter...| A0A011NL59| 9.0| true|NULL| NULL|hypothetical protein| | DEID01000170.1_27| s__Accumulibacter...| A0A011P2X1| 3.0| true|NULL| NULL|hypothetical protein| | JEMX01000078.1_72| s__Accumulibacter...| A0A011PMW2| 1.0| true|NULL| NULL|hypothetical protein| | DEID01000034.1_8| s__Accumulibacter...| A0A011PPW0| 5.0| true|NULL| NULL|hypothetical protein| |NZ_JEOB01000004.1...| s__Ruminococcus_D...| A0A011VVL2| 7.0| false|NULL| NULL|YcxB-like protein...| +--------------------+---------------------+-----------------+---------+------------+----+-----------------+--------------------+ only showing top 5 rows
In [6]:
# --- Summary stats for the full paradox set ---
from pyspark.sql.functions import countDistinct, percentile_approx
summary = paradox.agg(
count("*").alias("n_paradox_rows"),
countDistinct("gene_cluster_id").alias("n_distinct_clusters"),
countDistinct("gtdb_species_clade_id").alias("n_distinct_species_clades"),
avg("msa_depth").alias("mean_msa_depth"),
percentile_approx("msa_depth", 0.5).alias("median_msa_depth"),
spark_sum(col("hypothetical").cast("int")).alias("n_hypothetical"),
spark_sum((col("ec").isNotNull()).cast("int")).alias("n_has_ec"),
spark_sum((col("kegg_orthology_id").isNotNull()).cast("int")).alias("n_has_kegg")
)
summary.show(truncate=False)
summary_pd = summary.toPandas()
summary_pd.to_csv("/tmp/paradox_summary.csv", index=False)
print("Saved → /tmp/paradox_summary.csv")
+--------------+-------------------+-------------------------+-----------------+----------------+--------------+--------+----------+ |n_paradox_rows|n_distinct_clusters|n_distinct_species_clades|mean_msa_depth |median_msa_depth|n_hypothetical|n_has_ec|n_has_kegg| +--------------+-------------------+-------------------------+-----------------+----------------+--------------+--------+----------+ |415733 |415603 |14768 |4.569935511494156|4.0 |286439 |137 |346 | +--------------+-------------------+-------------------------+-----------------+----------------+--------------+--------+----------+
Saved → /tmp/paradox_summary.csv
In [7]:
# --- Top-1000 paradox proteins ranked by msa_depth ASC ---
# (msa_depth=1 is the most structurally unprecedented; ties broken by gene_cluster_id)
from pyspark.sql.functions import row_number
from pyspark.sql.window import Window
# Deduplicate: one row per gene_cluster_id (take the record with lowest msa_depth)
w = Window.partitionBy("gene_cluster_id").orderBy(col("msa_depth").asc())
paradox_dedup = (
paradox
.withColumn("rn", row_number().over(w))
.filter(col("rn") == 1)
.drop("rn")
)
print(f"Distinct paradox clusters: {paradox_dedup.count():,}")
top1000 = (
paradox_dedup
.orderBy(col("msa_depth").asc(), "gene_cluster_id")
.limit(1000)
.toPandas()
)
top1000.to_csv("/tmp/paradox_top1000.csv", index=False)
print(f"Saved top-1000 ({len(top1000)} rows) → /tmp/paradox_top1000.csv")
print("\nTop-10 paradox proteins:")
print(top1000[["gene_cluster_id","gtdb_species_clade_id","msa_depth","product","hypothetical","ec"]].head(10).to_string(index=False))
Distinct paradox clusters: 415,603
Saved top-1000 (1000 rows) → /tmp/paradox_top1000.csv Top-10 paradox proteins: gene_cluster_id gtdb_species_clade_id msa_depth product hypothetical ec AAVT01000002.1_305 s__Oceanicoccus_sp000169075--GB_GCA_000169075.1 1.0 hypothetical protein True NaN AAVT01000002.1_354 s__Oceanicoccus_sp000169075--GB_GCA_000169075.1 1.0 hypothetical protein True NaN AAVT01000003.1_307 s__Oceanicoccus_sp000169075--GB_GCA_000169075.1 1.0 hypothetical protein True NaN AAVT01000009.1_83 s__Oceanicoccus_sp000169075--GB_GCA_000169075.1 1.0 hypothetical protein True NaN AAVT01000010.1_63 s__Oceanicoccus_sp000169075--GB_GCA_000169075.1 1.0 hypothetical protein True NaN AEPR01000347.1_1 s__CAILRJ01_sp010365165--GB_GCA_000195205.2 1.0 DNA-directed RNA polymerase, omega subunit family protein False NaN AEPR01000360.1_2 s__CAILRJ01_sp010365165--GB_GCA_000195205.2 1.0 hypothetical protein True NaN AEPR01000527.1_4 s__CAILRJ01_sp010365165--GB_GCA_000195205.2 1.0 FXSXX-COOH protein False NaN AFXR01000179.1_4 s__Dwaynesavagella_sp000270205--RS_GCF_000270205.1 1.0 hypothetical protein True NaN AFXU01000074.1_3 s__Dwaynesavagella_sp000270205--RS_GCF_000270205.1 1.0 hypothetical protein True NaN
In [8]:
print("\n=== NB04b complete. Copy outputs: ===")
print(" cp /tmp/paradox_top1000.csv projects/alphafold_msa_annotation/data/")
print(" cp /tmp/paradox_summary.csv projects/alphafold_msa_annotation/data/")
=== NB04b complete. Copy outputs: === cp /tmp/paradox_top1000.csv projects/alphafold_msa_annotation/data/ cp /tmp/paradox_summary.csv projects/alphafold_msa_annotation/data/