Annotations and enrichment tab

database/Snakemake

@@ -14,9 +14,11 @@ Fill in database/config.yaml's `sources` URLs before running. The pipeline:
   6. load_mibig_compounds  } turn results into "compound" db entries (linked to MIBiG's URL)
   7. parse_mibig_gbks      - antiSMASH GBKs -> linear module readouts (PARAS-annotated)
   8. load_mibig_bgcs       - turn readouts into "bgc" db entries
+  9. annotate_npclassifier - chemical-class annotation (compounds only, NPClassifier API)
+  10. annotate_chebi       - bioactivity annotation (compounds only, local ChEBI flat files)
 
-Steps 4, 6, and 8 all mutate the same DuckDB file, so they're chained through marker
-files (rather than each declaring the database itself as `output`) to force
+Steps 4, 6, 8, 9, and 10 all mutate the same DuckDB file, so they're chained through
+marker files (rather than each declaring the database itself as `output`) to force
 Snakemake to serialize them -- DuckDB doesn't support concurrent writers.
 """
 
@@ -39,10 +41,45 @@ PARAS_TRAINING_DATA_PATH = config["paras"].get("training_data_path")
 PARSE_COMPOUNDS_WORKERS = config["compute"]["parse_compounds_workers"]
 PARSE_GBKS_WORKERS = config["compute"]["parse_gbks_workers"]
 
+TAXONOMY_ENABLED = config.get("taxonomy", {}).get("enabled", False)
+TAXDUMP_DIR = WORKDIR / "taxdump"
+
+NPCLASSIFIER_REQUESTS_PER_SECOND = config.get("npclassifier", {}).get("requests_per_second", 2.0)
+NPCLASSIFIER_WORKERS = config.get("npclassifier", {}).get("workers", 8)
+
+CHEBI_DIR = WORKDIR / "chebi"
+
+# Toggle whole branches of the pipeline off (see config.yaml's `enabled` comment) --
+# disabling npatlas/mibig skips their download+parse rules entirely, not just db
+# loading, since the "real" inputs below are only declared when enabled: nothing else
+# in the DAG needs npatlas/results.jsonl or mibig_gbk/readouts.jsonl, so Snakemake
+# never schedules the rules that would produce them. Same idea for chebi: disabling it
+# skips its bulk download entirely.
+ENABLED = config.get("enabled", {})
+NPATLAS_ENABLED = ENABLED.get("npatlas", True)
+MIBIG_ENABLED = ENABLED.get("mibig", True)
+NPCLASSIFIER_ENABLED = ENABLED.get("npclassifier", True)
+CHEBI_ENABLED = ENABLED.get("chebi", True)
+
 
 rule all:
     input:
-        MARKERS / "bgcs_loaded.done"
+        MARKERS / "bgcs_loaded.done",
+        MARKERS / "npclassifier_annotated.done",
+        MARKERS / "chebi_annotated.done"
+
+
+# ---------------------------------------------------------------------------
+# NCBI taxonomy dump (used to standardize phylogeny genus/species/type to taxids)
+# ---------------------------------------------------------------------------
+
+rule download_taxdump:
+    output:
+        names=TAXDUMP_DIR / "names.dmp",
+        nodes=TAXDUMP_DIR / "nodes.dmp"
+    run:
+        import taxonomy
+        taxonomy.download_taxdump(TAXDUMP_DIR)
 
 
 # ---------------------------------------------------------------------------
@@ -131,19 +168,24 @@ rule parse_npatlas:
 
 rule load_npatlas_compounds:
     input:
-        results=WORKDIR / "npatlas" / "results.jsonl",
-        db_created=MARKERS / "db_created.done"
+        db_created=MARKERS / "db_created.done",
+        # Only declared when enabled -- nothing else in the DAG needs npatlas/results.jsonl,
+        # so download_npatlas/resolve_npatlas_sdf/parse_npatlas never run when disabled.
+        **({"results": WORKDIR / "npatlas" / "results.jsonl"} if NPATLAS_ENABLED else {}),
+        **({"taxdump_names": TAXDUMP_DIR / "names.dmp", "taxdump_nodes": TAXDUMP_DIR / "nodes.dmp"} if TAXONOMY_ENABLED and NPATLAS_ENABLED else {})
     output:
         marker=touch(MARKERS / "npatlas_loaded.done")
     run:
-        import load_compounds
-        load_compounds.run(
-            results_path=input.results,
-            db_path=DB_PATH,
-            source="npatlas",
-            reaction_rules_path=RXN_RULES,
-            matching_rules_path=MXN_RULES,
-        )
+        if NPATLAS_ENABLED:
+            import load_compounds
+            load_compounds.run(
+                results_path=input.results,
+                db_path=DB_PATH,
+                source="npatlas",
+                reaction_rules_path=RXN_RULES,
+                matching_rules_path=MXN_RULES,
+                taxdump_dir=TAXDUMP_DIR if TAXONOMY_ENABLED else None,
+            )
 
 
 # ---------------------------------------------------------------------------
@@ -159,13 +201,18 @@ rule extract_mibig_compounds:
         # module docstring): MIBiG 4.0's GBKs dropped the ACCESSION.VERSION suffix
         # parse_gbks.py used to read this from, but the JSON's own top-level "version"
         # field still has it.
-        versions=WORKDIR / "mibig_json" / "versions.json"
+        versions=WORKDIR / "mibig_json" / "versions.json",
+        # accession -> {organism_name, ncbi_tax_id, biosyn_class}, consumed by both
+        # load_mibig_compounds and load_mibig_bgcs to populate phylogeny/chemical_class
+        # annotations (see RetroMolDuckDB.add_phylogeny_annotation/add_flat_annotation).
+        annotations=WORKDIR / "mibig_json" / "annotations.json"
     run:
         import extract_mibig_compounds
         extract_mibig_compounds.run(
             mibig_json_dir=input.extract_dir,
             output_path=output.compounds,
             versions_output_path=output.versions,
+            annotations_output_path=output.annotations,
         )
 
 
@@ -189,23 +236,35 @@ rule parse_mibig_compounds:
 
 rule load_mibig_compounds:
     input:
-        results=WORKDIR / "mibig_json" / "results.jsonl",
-        # MIBiG URLs need an accession's version, from the JSON (see
+        prev=MARKERS / "npatlas_loaded.done",
+        # Only declared when enabled -- nothing else in the DAG needs these, so
+        # download_mibig_json/extract_mibig_compounds/parse_mibig_compounds never run
+        # when disabled. MIBiG URLs need an accession's version, from the JSON (see
         # extract_mibig_compounds rule above) -- not from the GBKs, which no longer carry it.
-        versions=WORKDIR / "mibig_json" / "versions.json",
-        prev=MARKERS / "npatlas_loaded.done"
+        **(
+            {
+                "results": WORKDIR / "mibig_json" / "results.jsonl",
+                "versions": WORKDIR / "mibig_json" / "versions.json",
+                "annotations": WORKDIR / "mibig_json" / "annotations.json",
+            }
+            if MIBIG_ENABLED else {}
+        ),
+        **({"taxdump_names": TAXDUMP_DIR / "names.dmp", "taxdump_nodes": TAXDUMP_DIR / "nodes.dmp"} if TAXONOMY_ENABLED and MIBIG_ENABLED else {})
     output:
         marker=touch(MARKERS / "mibig_compounds_loaded.done")
     run:
-        import load_compounds
-        load_compounds.run(
-            results_path=input.results,
-            db_path=DB_PATH,
-            source="mibig",
-            reaction_rules_path=RXN_RULES,
-            matching_rules_path=MXN_RULES,
-            mibig_versions_path=input.versions,
-        )
+        if MIBIG_ENABLED:
+            import load_compounds
+            load_compounds.run(
+                results_path=input.results,
+                db_path=DB_PATH,
+                source="mibig",
+                reaction_rules_path=RXN_RULES,
+                matching_rules_path=MXN_RULES,
+                mibig_versions_path=input.versions,
+                mibig_annotations_path=input.annotations,
+                taxdump_dir=TAXDUMP_DIR if TAXONOMY_ENABLED else None,
+            )
 
 
 # ---------------------------------------------------------------------------
@@ -232,17 +291,97 @@ rule parse_mibig_gbks:
 
 rule load_mibig_bgcs:
     input:
-        readouts=WORKDIR / "mibig_gbk" / "readouts.jsonl",
-        versions=WORKDIR / "mibig_json" / "versions.json",
-        prev=MARKERS / "mibig_compounds_loaded.done"
+        prev=MARKERS / "mibig_compounds_loaded.done",
+        # Only declared when enabled -- nothing else in the DAG needs these, so
+        # download_mibig_gbk/parse_mibig_gbks never run when disabled.
+        **(
+            {
+                "readouts": WORKDIR / "mibig_gbk" / "readouts.jsonl",
+                "versions": WORKDIR / "mibig_json" / "versions.json",
+                "annotations": WORKDIR / "mibig_json" / "annotations.json",
+            }
+            if MIBIG_ENABLED else {}
+        ),
+        **({"taxdump_names": TAXDUMP_DIR / "names.dmp", "taxdump_nodes": TAXDUMP_DIR / "nodes.dmp"} if TAXONOMY_ENABLED and MIBIG_ENABLED else {})
     output:
         marker=touch(MARKERS / "bgcs_loaded.done")
     run:
-        import load_bgcs
-        load_bgcs.run(
-            readouts_path=input.readouts,
-            db_path=DB_PATH,
-            reaction_rules_path=RXN_RULES,
-            matching_rules_path=MXN_RULES,
-            mibig_versions_path=input.versions,
+        if MIBIG_ENABLED:
+            import load_bgcs
+            load_bgcs.run(
+                readouts_path=input.readouts,
+                db_path=DB_PATH,
+                reaction_rules_path=RXN_RULES,
+                matching_rules_path=MXN_RULES,
+                mibig_versions_path=input.versions,
+                mibig_annotations_path=input.annotations,
+                taxdump_dir=TAXDUMP_DIR if TAXONOMY_ENABLED else None,
+            )
+
+
+# ---------------------------------------------------------------------------
+# Step 9: NPClassifier chemical-class annotation (compounds only)
+# ---------------------------------------------------------------------------
+
+rule annotate_npclassifier:
+    input:
+        # Chained after bgcs_loaded.done (not just the compound-loading steps) purely to
+        # serialize this write against load_mibig_bgcs's -- both mutate DB_PATH and DuckDB
+        # doesn't support concurrent writers (see module docstring at the top of this file).
+        prev=MARKERS / "bgcs_loaded.done"
+    output:
+        marker=touch(MARKERS / "npclassifier_annotated.done")
+    # Declared like parse_npatlas/parse_mibig_gbks's worker counts, even though this
+    # rule is I/O-bound rather than CPU-bound: without a `threads:` declaration,
+    # Snakemake assumes a 1-core job while it actually opens NPCLASSIFIER_WORKERS
+    # concurrent connections, so `--cores` wouldn't see or cap it. Declaring it here
+    # means `--cores` fewer than npclassifier.workers automatically caps this job's
+    # `threads` value at run time -- passed through below instead of the raw config
+    # number, so a run with e.g. `--cores 6` gets at most 6 workers even if
+    # npclassifier.workers says 8.
+    threads: NPCLASSIFIER_WORKERS
+    run:
+        if NPCLASSIFIER_ENABLED:
+            import annotate_npclassifier
+            annotate_npclassifier.run(
+                db_path=DB_PATH,
+                cache_path=WORKDIR / "npclassifier" / "cache.jsonl",
+                requests_per_second=NPCLASSIFIER_REQUESTS_PER_SECOND,
+                workers=threads,
+            )
+
+
+# ---------------------------------------------------------------------------
+# Step 10: ChEBI bioactivity annotation (compounds only)
+# ---------------------------------------------------------------------------
+
+rule download_chebi:
+    output:
+        compounds=CHEBI_DIR / "compounds.tsv.gz",
+        structures=CHEBI_DIR / "structures.tsv.gz",
+        relation=CHEBI_DIR / "relation.tsv.gz"
+    run:
+        import chebi
+        chebi.download_chebi_flat_files(CHEBI_DIR)
+
+
+rule annotate_chebi:
+    input:
+        # Chained after npclassifier_annotated.done purely to serialize this write
+        # against annotate_npclassifier's (see module docstring at the top of this file).
+        prev=MARKERS / "npclassifier_annotated.done",
+        # Only declared when enabled -- download_chebi never runs when disabled.
+        **(
+            {
+                "compounds": CHEBI_DIR / "compounds.tsv.gz",
+                "structures": CHEBI_DIR / "structures.tsv.gz",
+                "relation": CHEBI_DIR / "relation.tsv.gz",
+            }
+            if CHEBI_ENABLED else {}
         )
+    output:
+        marker=touch(MARKERS / "chebi_annotated.done")
+    run:
+        if CHEBI_ENABLED:
+            import annotate_chebi
+            annotate_chebi.run(db_path=DB_PATH, chebi_dir=CHEBI_DIR)

database/config.yaml

@@ -6,17 +6,48 @@ sources:
   mibig_json_url: "https://dl.secondarymetabolites.org/mibig/mibig_json_4.0.tar.gz"
   mibig_gbk_url: "https://dl.secondarymetabolites.org/mibig/mibig_gbk_4.0.tar.gz"
 
+enabled:
+  # Toggle whole branches of the pipeline off. Disabling a source skips its downloads
+  # and RetroMol parsing entirely (not just db loading) -- e.g. npatlas: false means
+  # download_npatlas/parse_npatlas never run. mibig covers both MIBiG compounds and
+  # BGCs (they share the same JSON/GBK downloads). npclassifier gates step 9
+  # (chemical-class annotation) independently of which source(s) loaded the compounds
+  # it classifies.
+  npatlas: true
+  mibig: true
+  npclassifier: true
+  chebi: true
+
 paths:
   # Final DuckDB database produced by the pipeline.
-  database: "/Users/davidmeijer/Downloads/retromol.duckdb"
+  database: "/Users/davidmeijer/Desktop/retromol.duckdb"
 
   # Scratch space for downloads and intermediate per-step results.
-  workdir: "/Users/davidmeijer/retromol_tmp"
+  workdir: "/Users/davidmeijer/Desktop/retromol_tmp"
 
   # null -> RuleSet.load_default()'s bundled reaction/matching rules.
   reaction_rules: null
   matching_rules: null
 
+npclassifier:
+  # Rate limit for the free, GNPS2-hosted NPClassifier API (no published limit --
+  # kept conservative). This is the *combined* rate across all workers below, not
+  # per-worker. Classifications are cached in workdir/npclassifier/cache.jsonl, so
+  # reruns only pay for compounds not already classified.
+  requests_per_second: 50.0
+
+  # Concurrent requests. Classification is I/O-bound (waiting on the API), not
+  # CPU-bound, so this is safe to raise well past your core count -- it's bounded by
+  # requests_per_second above either way.
+  workers: 8
+
+taxonomy:
+  # NCBI taxdump (names.dmp/nodes.dmp), downloaded once into workdir/taxdump and reused
+  # across pipeline runs -- used to standardize phylogeny genus/species/type to NCBI
+  # taxids (see database/scripts/taxonomy.py). Set to null to skip taxid resolution
+  # entirely (phylogeny is then stored as raw, unstandardized text/no taxids).
+  enabled: true
+
 paras:
   threshold: 0.1
   keep_top: 3

database/envs/retromol.yaml

@@ -0,0 +1,34 @@
+# Conda env for running database/Snakemake (the database-construction pipeline only --
+# this has no bearing on the GUI/webapp).
+#
+# Every rule in that Snakefile uses Snakemake's `run:` directive (inline Python that
+# imports and calls database/scripts/*.py functions directly), not `shell:`. Snakemake
+# executes `run:` blocks in its own process rather than a subprocess, so a per-rule
+# `conda:` env (the --use-conda mechanism) would NOT actually reach any of this
+# pipeline's own code -- only `shell()` calls made from *within* a `run:` block, which
+# this pipeline never makes. There is deliberately no per-rule conda: directive in
+# database/Snakemake, and no --use-conda flag in the invocation below, because of this.
+#
+# This env is instead the one Snakemake itself runs in, so it needs both Snakemake and
+# every runtime dependency database/scripts/*.py imports (rdkit, duckdb, etc. -- all
+# already listed in the repo's own pyproject.toml, pulled in by the editable install
+# below rather than duplicated here).
+#
+# Usage:
+#
+#   conda env create -f database/envs/retromol.yaml
+#   conda activate retromol
+#   pip install -e .
+#   snakemake -p -s database/Snakemake --configfile database/config.yaml \
+#       --workflow-profile database/profiles/slurm
+name: retromol
+
+channels:
+  - conda-forge
+
+dependencies:
+  - python=3.11
+  - pip
+  - snakemake
+  - snakemake-executor-plugin-slurm
+  - mamba