TL;DR. Five candidate “pigmentation networks” disagree with each other far more than a single merged network would suggest. Node-level comparison of four curated/pathway peer sources (Raghunath-168, D’Arcy-S1-243, KEGG-hsa04916-101, Reactome curated-union-156) finds only 7 genes in common across all four. Adding Bajpai et al. 2023’s 169-gene CRISPR screen as a fifth, orthogonal NODE layer (weighted casTLE effect + uniform reduces-pigmentation sign, no fabricated edges) shows 142 of 169 hits (84.0%) are orphans absent from all four networks — candidate discovery nodes. And even restricting to one single evidentiary layer (STRING protein-protein association) applied consistently: our fresh v12.0 pull recovers only 59.6% of Raghunath’s own mechanistic edges, and agrees with D’Arcy’s own frozen STRING snapshot on only ~66% of edges (~34% drift) at the identical node scope and threshold — the specific network you choose changes the answer. This notebook pulls D’Arcy 2023’s six supplementary tables (data/raw/darcy2023/*.xlsx), KEGG hsa04916 and Reactome (frozen data/external/db_responses/*), a fresh STRING v12.0 pull (string_network_pulls_v12.json), and Bajpai’s 169 hits (data/processed/bajpai2023_crispr_hits.csv) — 8 sources total once Baxter and the 13-paper GWAS/curated- loci table are added to the typology table. Contribution to the flagship rescue screen: this notebook is the substrate-comparison layer — it establishes which network(s) a rescued locus’s causal gene should be checked against, quantifies how much that choice matters, and hands NB6/NB7 the Reactome-only 121-gene connection set plus the 142 Bajpai orphan hits as two distinct discovery surfaces for candidate rescue connections. The one number that matters: of Bajpai’s 169 causally-validated pigmentation effectors, 84.0% are invisible to every curated/pathway network compared here.
Goal. Before any of the candidate gene sets/networks for pigmentation are merged into a single harmonized layer, compare them as peers and report what each one has that the others lack. The comparison is treated as a candidate finding in its own right — comparing D’Arcy to Raghunath, what each has that the other lacks, may itself be a finding — not a preprocessing step on the way to a merge.
Four candidate sources, four different curation methods:
Curated from Baxter/Yamaguchi gene lists + OMIM phenotype records
243
KEGG hsa04916
“Melanogenesis” reference pathway
KEGG’s own pathway-diagram curation
101
Reactome R-HSA-5662702
“Melanin biosynthesis” reference pathway
Reactome’s own reaction-level curation
5 (core enzymatic pathway; see the gaps note in Step 4)
D’Arcy 2023 is not just a STRING pull. The paper’s six supplementary tables decompose into distinct layers with distinct evidentiary weight, and this notebook keeps them separate rather than treating “D’Arcy” as one STRING network:
Table S1 — 243-gene curated OMIM disease-gene table. This is the layer that is a genuine peer of Raghunath’s 168-gene mechanistic set (curated gene list vs. curated gene list) and the one used for the node-level set comparison below.
Tables S4 + S5 — a 451/452-node, 4668-edge STRING protein-protein association network built by expanding the 243 S1 genes through STRING; S4 holds the edges, S5 the per-node SysGO/location/disease annotation. This is an association layer, undirected and unsigned — categorically different from Raghunath’s directed, signed mechanism.
Table S2 — SysGO process + STRING-cluster + subcellular-location annotation for the 243 S1 genes (annotation only, no new genes).
Table S6 — A375 (unpigmented) / FM55 (pigmented) melanoma mass-spectrometry (LFQ) protein expression — an orthogonal experimental layer, not a curated gene list.
License / provenance. D’Arcy, C.E. & Kiel, C. (2023). Network Analysis and Protein-Protein Interaction Prediction of Skin Pigmentation Genes. Bioengineering 10(1):13. DOI 10.3390/bioengineering10010013, PMC9854651, PMID 36671585, CC BY 4.0. Raghunath, R. et al. (2015). BMC Res Notes (see docs/specs/ and data/processed/raghunath_* for the per-edge citation legend already built in Notebooks 1–2).
Every frozen snapshot this notebook reads is committed in-repo alongside the notebook, verified by an explicit existence check in Step 0 below — not merely referenced in a meta.json.
Key terms — so this notebook stands on its own (you shouldn’t need the other notebooks to read this one).
Melanogenesis network — the cellular pathway that makes melanin (skin/hair/eye pigment). The project’s mechanistic backbone is Raghunath et al. 2015’s version, rebuilt in earlier notebooks (NB1–2) as a directed, signed graph: every edge has a direction and is marked activation or repression.
STRING — a public database of functional protein–protein associations. A STRING edge just means two proteins are functionally linked; it has no direction and no activation/repression sign (unlike Raghunath’s mechanism). Used here as one evidence layer applied identically to every gene set.
casTLE effect — the per-gene effect size from Bajpai et al. 2023’s genome-wide CRISPR knockout screen for melanin (“casTLE” is the screen’s statistical scoring method). Magnitude = how strongly removing the gene changes pigment; sign = direction (all 169 hits reduce pigment when knocked out). It is a phenotype weight on a single gene, not a gene–gene edge.
Effector (gene) — the gene actually responsible for a phenotype (here, making melanin), as opposed to whichever gene merely sits nearest a genetic signal.
Orphan (Bajpai orphan) — a gene that is a real experimental hit in the CRISPR screen but appears in none of the curated pigmentation networks compared here — a candidate the existing maps missed.
Effector-uncertain locus — a locus a published study tied to pigmentation by statistical association where which gene is actually causal is genuinely open: the authors named only the nearest gene (or stated the effector is unknown) and that gene is not an established pigmentation effector. These — not loci that already sit at a canonical gene such as TYR or OCA2 — are the targets of the project’s “rescue” aim.
The “rescue” screen / rescue connection — the project’s overarching goal (executed later, in NB7): connect an effector-uncertain locus’s likely causal gene into curated pigmentation biology so it stops being an isolated statistical hit. This notebook only builds the substrate for that test; it does not run it.
Convergence grade — the project’s trust score for a claimed gene→pigmentation link: how many independent evidence lines (statistical association, curated mechanism, pathway membership, experimental screen) agree. Each line is added, never used to drop a gene.
Methods overview
Extract — D’Arcy 2023’s six supplementary tables are extracted from the committed raw .xlsx files into typed, cited data/processed/darcy2023_*.csv tables (Step 1).
Acquire the two pathway references — KEGG hsa04916 (as used in NB2) and Reactome R-HSA-5662702 (Steps 2–3).
State an explicit gene-symbol harmonization rule — before any cross-source count is computed, because the D’Arcy-vs-Raghunath overlap counts depend on it (Step 4).
Node-level set operations across the four peer gene sets {Raghunath, D’Arcy-S1, KEGG, Reactome}: each-only, pairwise intersections, full intersection (Step 5).
Overlay D’Arcy’s own OMIM phenotype_class annotation onto the D’Arcy-only genes, and test whether the D’Arcy-vs-Raghunath non-overlap is enriched for hypopigmentation-class genes — the concrete, falsifiable sub-finding (Step 6).
Apply STRING identically to every gene set (our own consistent layer, Module 2) and use it for two things: (a) an edge-level coverage check against Raghunath’s signed/directed edges, and (b) a quantified drift comparison against D’Arcy’s own frozen STRING snapshot (Steps 7–9).
Citation-completeness gate — every emitted node/edge table carries a resolvable citation; a closing assertion fails the run if any row does not (Step 10).
All four database pulls this notebook depends on (D’Arcy’s own committed .xlsx, KEGG, Reactome, STRING) follow the frozen-DB pattern: a visible re-runnable query cell (exact query + query-UTC + a REQUERY guard) sits above a load of the committed, verbatim frozen response, and a closing assertion in each section proves the replay reproduces the committed table.
Step 0 — Setup and existence check
Resolve the repo root the same way as Notebooks 1–3, then assert that every frozen snapshot this notebook will read is present on disk before doing anything else.
Show code
import jsonimport datetimeimport refrom itertools import combinationsfrom pathlib import Pathimport pandas as pdfrom scipy.stats import fisher_exactROOT = Path.cwd().parent if Path.cwd().name =="notebooks"else Path.cwd()RAW_DARCY = ROOT /"data"/"raw"/"darcy2023"PROC = ROOT /"data"/"processed"FROZEN = ROOT /"data"/"external"/"db_responses"REQUIRED_FROZEN_FILES = [ RAW_DARCY /"Table S1 Bioengineering FINAL.xlsx", RAW_DARCY /"Table S2 Bioengineering FINAL.xlsx", RAW_DARCY /"Table S4 Bioengineering FINAL.xlsx", RAW_DARCY /"Table S5 Bioengineering FINAL.xlsx", RAW_DARCY /"Table S6 Bioengineering FINAL.xlsx", FROZEN /"kegg_hsa04916.json", FROZEN /"reactome_R-HSA-5662702_participants.json", FROZEN /"reactome_pigmentation_curated_union.json", FROZEN /"reactome_mitf_subpathway_structure.json", FROZEN /"string_network_pulls_v12.json",]missing = [str(p) for p in REQUIRED_FROZEN_FILES ifnot p.exists()]assertnot missing, f"NB2-lesson check FAILED — missing frozen inputs (not just referenced, must be on disk): {missing}"print(f"NB2-lesson check passed: all {len(REQUIRED_FROZEN_FILES)} frozen/raw inputs present on disk at {ROOT.name}/.")for p in REQUIRED_FROZEN_FILES:print(" ", p.relative_to(ROOT))
NB2-lesson check passed: all 10 frozen/raw inputs present on disk at pigmentation-gene-network/.
data/raw/darcy2023/Table S1 Bioengineering FINAL.xlsx
data/raw/darcy2023/Table S2 Bioengineering FINAL.xlsx
data/raw/darcy2023/Table S4 Bioengineering FINAL.xlsx
data/raw/darcy2023/Table S5 Bioengineering FINAL.xlsx
data/raw/darcy2023/Table S6 Bioengineering FINAL.xlsx
data/external/db_responses/kegg_hsa04916.json
data/external/db_responses/reactome_R-HSA-5662702_participants.json
data/external/db_responses/reactome_pigmentation_curated_union.json
data/external/db_responses/reactome_mitf_subpathway_structure.json
data/external/db_responses/string_network_pulls_v12.json
D’Arcy’s raw supplement is six .xlsx files under data/raw/darcy2023/ (CC BY 4.0, PMC9854651), each with a title row + spacer row before the header, so every read uses the row-skipping header= documented per table below. This is the extraction cell — deterministic, no network calls — turning the committed raw files into typed, cited processed CSVs. Row/gene counts are printed and cross-checked against the paper’s own reported totals (243 disease genes; 451/452-node, 4668-edge expanded PPI network) as the extraction fidelity check.
S6: 4232 genes with mass-spec LFQ in A375 and/or FM55 (33 A375-only, 30 FM55-only)
D'Arcy 2023 extraction complete — 5 processed CSVs written to data/processed/:
darcy2023_S1_disease_genes.csv
darcy2023_S2_sysgo_annotation.csv
darcy2023_S4_string_edges.csv
darcy2023_S5_string_nodes.csv
darcy2023_S6_massspec_expression.csv
Step 2 — KEGG hsa04916 (Melanogenesis) gene set
Already frozen and committed by Notebook 2 (data/external/db_responses/kegg_hsa04916.json, rest.kegg.jp/link/hsa/hsa04916). Re-used here as-is — no re-query — because it is already a committed, version-stamped snapshot of a live resource.
Show code
kegg_snap = json.load(open(FROZEN /"kegg_hsa04916.json"))kegg_genes =set(s.upper() for s in kegg_snap["symbols"])print(f"KEGG {kegg_snap['pathway']} ({kegg_snap['pathway_name']}): {len(kegg_genes)} genes "f"[source: {kegg_snap['source_endpoint']}]")assertlen(kegg_genes) == kegg_snap["n"] ==101
Query (frozen, re-runnable). The Reactome AnalysisService (map_reactome_pathways, via the genes-ontologies connector) maps genes to pathways, not a pathway to its member molecules — so the member list for R-HSA-5662702 is pulled from the Reactome ContentService (reactome.org/ContentService/data/participants/<stId>), which returns every participating physical entity (complexes, small molecules, and their reference gene products) for the pathway. This complements — does not replace — the map_reactome_pathways verification confirming TYR/TYRP1/DCT map to R-HSA-5662702, release 97.
Frozen response: data/external/db_responses/reactome_R-HSA-5662702_participants.json (query-UTC stamped inside the file; REQUERY=True re-hits the live ContentService — must be run outside this sandboxed kernel, which has no direct network access).
Reactome query-UTC: 2026-07-12T01:16:36.004890+00:00
Reactome release version: 97
Pathway detail: Melanin biosynthesis | R-HSA-5662702 | Homo sapiens
Reactome R-HSA-5662702 participating gene products: 5 -> ['DCT', 'OCA2', 'SLC45A2', 'TYR', 'TYRP1']
Non-gene-product participants (small molecules/complexes, 27 rows) — not genes, excluded from the gene-set comparison (see the Gaps discussion below).
Gap (Step 3a). R-HSA-5662702 “Melanin biosynthesis” is Reactome’s most specific pathway for the enzymatic core of melanin synthesis, and its participants resolve to exactly 5 protein-coding genes — TYR, TYRP1, DCT, OCA2, SLC45A2 — plus a set of small-molecule intermediates (L-DOPA, dopaquinone, dopachrome, DHI/DHICA, eumelanin/pheomelanin precursors) that are not genes and are excluded from the gene-set comparison. This is a small, precise pathway, not a comprehensive melanocyte-biology pathway.
Reactome has no single “Pigmentation” parent pathway. One might expect a parent accession grouping Melanin biosynthesis, MC1R signaling, and melanosome trafficking as sibling children, whose descendant gene set could be pulled wholesale. No such accession exists. Accession R-HSA-5619507 — a plausible candidate — resolves via the Reactome ContentService (data/query/R-HSA-5619507) to “Activation of HOX genes during differentiation”, an unrelated developmental-biology pathway about Hox gene clusters with no connection to melanocyte biology. Reactome has no pathway literally named “Pigmentation”, and — checked via the ContentService ancestors endpoint — Melanin biosynthesis sits under Metabolism → Metabolism of amino acids and derivatives, while the MITF-regulation pathways sit under Developmental Biology: two disjoint branches of the hierarchy with no shared low-level parent. There is consequently no single accession whose descendant tree is “the broader Reactome pigmentation pathway.” Step 3b below builds an explicit, itemized, manually curated union of independently verified Reactome accessions instead — the correct approach when no single parent accession spans the pigmentation pathways.
Step 3b — Reactome pigmentation gene set, broadened (curated union of verified pathways)
Because no single Reactome accession has the whole of pigmentation biology as its descendant tree, this step builds the broader Reactome layer as an explicit, itemized union of four independently confirmed, live-verified accessions, each queried and frozen the same way as Step 3a:
Accession
Name
Schema class
Role in the union
R-HSA-5662702
Melanin biosynthesis
Pathway
terminal enzymatic core (Step 3a, reused)
R-HSA-9730414
MITF-M-regulated melanocyte development
Pathway
MITF’s broad transcriptional program in melanocytes — the largest verified pathway substantially overlapping pigmentation biology
R-HSA-388601
Melanocortin receptors
DefinedSet
MC1R (canonical pigmentation receptor) and its paralogs MC2–5R, taken as the full DefinedSet membership
disease pathway for the one core pigmentation gene not already covered by the other three
Query (frozen, re-runnable).R-HSA-9730414 and R-HSA-5619036 are Pathways — queried via data/participants/<stId> exactly as in Step 3a. R-HSA-388601 is a DefinedSet (a same-function protein group, not a reaction network) — its member proteins are pulled via data/complex/<stId>/subunits, the ContentService endpoint for set/complex membership, since participants returns nothing for a bare DefinedSet. Frozen response, including the raw participant/subunit payloads for all three newly-pulled components (Step 3a’s R-HSA-5662702 snapshot is reused, not duplicated): data/external/db_responses/reactome_pigmentation_curated_union.json.
Tradeoff.R-HSA-9730414 alone contributes 152 genes, the large majority of the union’s size — and many of its members are MITF targets in cell-cycle control, apoptosis, epithelial-mesenchymal transition, and lysosome biogenesis, whose relevance to pigmentation specifically (as opposed to melanocyte biology broadly) is indirect. This gene set is included in full — not hand-filtered down to a pigmentation-only subset — because doing that filtering would itself be a judgment call requiring its own justification; the union is reported as-is, with this tradeoff named, rather than silently narrowed.
Show code
REQUERY_REACTOME_UNION =Falseif REQUERY_REACTOME_UNION:import requests# Pathways: participants endpoint (same pattern as Step 3a) parts_mitf_live = requests.get("https://reactome.org/ContentService/data/participants/R-HSA-9730414", headers={"Accept": "application/json"}, timeout=60).json() parts_oca6_live = requests.get("https://reactome.org/ContentService/data/participants/R-HSA-5619036", headers={"Accept": "application/json"}, timeout=30).json()# DefinedSet: subunits endpoint (participants returns empty for a bare DefinedSet) mc_subunits_live = requests.get("https://reactome.org/ContentService/data/complex/R-HSA-388601/subunits", headers={"Accept": "application/json"}, timeout=30).json()# ... assemble into the same schema as the frozen file below and re-freeze ...reactome_union_snap = json.load(open(FROZEN /"reactome_pigmentation_curated_union.json"))print("Reactome curated-union query-UTC:", reactome_union_snap["query_utc"])print(f"Components ({len(reactome_union_snap['components'])}):")for comp in reactome_union_snap["components"]:print(f" {comp['stId']} ({comp['schemaClass']}) \"{comp['displayName']}\" -> {comp['n_genes']} genes")reactome_genes_broad =set(reactome_union_snap["curated_union_genes"])print(f"\nCurated union across all 4 components: n={len(reactome_genes_broad)}")assert reactome_genes_broad ==set().union(*[set(c["genes"]) for c in reactome_union_snap["components"]]), \"replay does not reproduce the frozen curated union -- committed snapshot is not self-consistent"print("Replay assertion passed: recomputing the union from the frozen per-component gene lists ""reproduces the frozen curated_union_genes list exactly.")# Keep the Step 3a precision set too -- both are used downstream, for different purposes.print(f"\nFor reference, Step 3a's precision Melanin-biosynthesis-only set remains: "f"n={len(reactome_genes)} -> {sorted(reactome_genes)}")
Reactome curated-union query-UTC: 2026-07-12T01:57:32.093720+00:00
Components (4):
R-HSA-5662702 (Pathway) "Melanin biosynthesis" -> 5 genes
R-HSA-9730414 (Pathway) "MITF-M-regulated melanocyte development" -> 152 genes
R-HSA-388601 (DefinedSet) "Melanocortin receptors" -> 5 genes
R-HSA-5619036 (Pathway) "Defective SLC24A5 causes oculocutaneous albinism 6 (OCA6)" -> 1 genes
Curated union across all 4 components: n=156
Replay assertion passed: recomputing the union from the frozen per-component gene lists reproduces the frozen curated_union_genes list exactly.
For reference, Step 3a's precision Melanin-biosynthesis-only set remains: n=5 -> ['DCT', 'OCA2', 'SLC45A2', 'TYR', 'TYRP1']
Step 4 — Raghunath’s 168-gene set and the gene-symbol harmonization rule
The harmonization rule (stated explicitly, because the D’Arcy-vs-Raghunath overlap counts in DATA_SOURCES.md depend on it). Raghunath’s mechanistic backbone (data/processed/gene_network_nodes.csv, built by Notebooks 1–2) has 168 gene-level nodes. Six of those 168 are not single genes but enzyme_activity_class tokens standing in for an HGNC gene group (PLA2, MMPs, Trypsin, Phosphodiesterase, PKC, PLC — 6 tokens covering 115 member genes via data/external/db_responses/hgnc_gene_groups.json, already frozen by NB2). This notebook states two rules and reports both:
Rule A (literal token set, n=168) — treat each of the 168 gene-level node labels as a literal gene-symbol string, exactly as it appears in gene_network_nodes.csv. This is the set-comparison unit used for the headline node-level counts below, because it matches what gene_network_nodes.csv already commits to and reproduces the DATA_SOURCES.md 465/230/118 figures exactly (checked below).
Rule B (HGNC-expanded set) — replace the 6 enzyme-class tokens with their full HGNC gene-group membership (115 member genes), so “does gene X overlap Raghunath” is asked at the level of actual gene symbols rather than a class label a real gene will never literally match. 10 of those 115 member genes (PLA2G1B, PRSS1, MMP1, PDE4B, PRKCA/B/D/H/Z, PLCB1) are already present in the 162 non-class-token genes as their own named nodes, so the expansion adds 105 net-new genes, giving a 267-gene Rule-B set. Reported as a secondary check; the practical difference on the headline D’Arcy-union non-overlap count is small (2 fewer genes flagged absent under Rule B) because almost none of the expanded member genes happen to appear in the D’Arcy/KEGG/Reactome sets used here.
No other alias resolution is applied — cross-source symbols are compared as the case-normalized (upper-case, whitespace-stripped) literal strings each source itself uses; this is a deliberate, stated choice (not every possible alias/synonym pair is chased), consistent with the level of harmonization already committed in gene_network_nodes.csv.
Show code
gene_network_nodes = pd.read_csv(PROC /"gene_network_nodes.csv")raghunath_168 =set(gene_network_nodes["gene"].astype(str).str.strip().str.upper())assertlen(raghunath_168) ==168hgnc_groups = json.load(open(FROZEN /"hgnc_gene_groups.json"))enzyme_tokens = {"PLA2", "MMPS", "TRYPSIN", "PHOSPHODIESTERASE", "PKC", "PLC"}base_162 = raghunath_168 - enzyme_tokenshgnc_member_genes =set()for tok_key in ["PLA2", "MMPs", "Trypsin", "Phosphodiesterase", "PKC", "PLC"]: hgnc_member_genes |=set(m.upper() for m in hgnc_groups[tok_key]["members"])raghunath_270_expanded = base_162 | hgnc_member_genes# Some HGNC member genes (e.g. individual PLA2G1B, PRSS1, MMP1, PDE4B, PRKC* isoforms) are ALREADY# present in the 162-gene literal base as their own named nodes -- the expansion adds only the NET NEW ones.n_already_present =len(base_162 & hgnc_member_genes)n_net_new =len(hgnc_member_genes - base_162)print(f"Rule A — Raghunath literal token set: n={len(raghunath_168)}")print(f"Rule B — Raghunath HGNC-expanded set: n={len(raghunath_270_expanded)} "f"(162 literal non-class genes + {len(hgnc_member_genes)} HGNC group members, of which "f"{n_already_present} were already present as named nodes -> {n_net_new} net-new genes added)")assertlen(raghunath_270_expanded) ==len(base_162) + n_net_new
Rule A — Raghunath literal token set: n=168
Rule B — Raghunath HGNC-expanded set: n=267 (162 literal non-class genes + 115 HGNC group members, of which 10 were already present as named nodes -> 105 net-new genes added)
Step 5 — Node-level set operations across the four peer gene sets
Each source is treated as a peer gene set: {Raghunath-168 (Rule A), D’Arcy-S1-243, KEGG-hsa04916-101, Reactome-curated-union-156 (Step 3b — the broadened set: Melanin biosynthesis + MITF-M-regulated melanocyte development + Melanocortin receptors + the SLC24A5/OCA6 disease pathway)}. Using the 156-gene broadened set here, rather than the 5-gene Step-3a precision set, is a deliberate choice per PI direction: the precision set was judged too thin to serve as a mechanistic-tier peer in the four-way comparison. This step computes each-only, all pairwise intersections, and the full intersection — the results the manuscript reports as the actual set comparison. (The Step-3a 5-gene precision set is retained separately and referenced in the discussion as a check on which of the other three sources recover the minimal causal core.)
Show code
darcy_s1_genes =set(darcy_s1["gene"])assertlen(darcy_s1_genes) ==243gene_sets = {"Raghunath_168": raghunath_168,"DArcy_S1_243": darcy_s1_genes,"KEGG_hsa04916_101": kegg_genes,"Reactome_curated_union_156": reactome_genes_broad,}for name, s in gene_sets.items():print(f"{name}: n={len(s)}")full_union =set().union(*gene_sets.values())full_intersection =set.intersection(*gene_sets.values())print(f"\nFull union (all 4 sources): n={len(full_union)}")print(f"Full intersection (all 4 sources): n={len(full_intersection)} -> {sorted(full_intersection)}")print("\n--- Each-only (exclusive to that source) ---")each_only = {}names =list(gene_sets.keys())for name in names: others =set().union(*[gene_sets[o] for o in names if o != name]) each_only[name] = gene_sets[name] - othersprint(f"{name} ONLY: n={len(each_only[name])}")print("\n--- Pairwise intersections ---")pairwise = {}for a, b in combinations(names, 2): inter = gene_sets[a] & gene_sets[b] pairwise[(a, b)] = inter label =f"{a} ∩ {b}"print(f"{label}: n={len(inter)}"+ (f" -> {sorted(inter)}"iflen(inter) <=30else""))
# --- Reproduce the DATA_SOURCES.md 465/230/118 figures under Rule A, and report the Rule B delta ---darcy_union_508 = darcy_s1_genes |set(darcy_s5["gene_string"]) # S1 ∪ S5, matching the manifest's own definitionassertlen(darcy_union_508) ==508, f"expected 508 (S1∪S5), got {len(darcy_union_508)}"absent_from_raghunath_A = darcy_union_508 - raghunath_168absent_from_raghunath_B = darcy_union_508 - raghunath_270_expandeds5_disease_genes =set(darcy_s5.loc[darcy_s5["disease_gene_flag"].notna(), "gene_string"])disease_flagged_A = absent_from_raghunath_A & (darcy_s1_genes | s5_disease_genes)s1_gene_class = darcy_s1.groupby("gene")["phenotype_class"].apply(lambda x: set(x.dropna())).to_dict()hypo_genes_s1 = {g for g, classes in s1_gene_class.items() if"Hypopigmentation"in classes}hypo_in_nonoverlap = disease_flagged_A & hypo_genes_s1print("Reproduction of the DATA_SOURCES.md D'Arcy-vs-Raghunath cross-check (S1∪S5 union, n=508):")print(f" Rule A (168 literal tokens) — D'Arcy genes absent from Raghunath: {len(absent_from_raghunath_A)} "f"(manifest states 465)")print(f" Rule B ({len(raghunath_270_expanded)}-gene HGNC-expanded) — D'Arcy genes absent from Raghunath: {len(absent_from_raghunath_B)} "f"(delta vs Rule A: {len(absent_from_raghunath_A) -len(absent_from_raghunath_B)})")print(f" disease-flagged among Rule-A non-overlap: {len(disease_flagged_A)} (manifest states 230)")print(f" hypopigmentation-class among those disease-flagged: {len(hypo_in_nonoverlap)} (manifest states 118)")assertlen(absent_from_raghunath_A) ==465, "Rule-A non-overlap does not reproduce the manifest's 465"assertlen(disease_flagged_A) ==230, "disease-flagged count does not reproduce the manifest's 230"assertlen(hypo_in_nonoverlap) ==118, "hypopigmentation-class count does not reproduce the manifest's 118"print("\nDATA_SOURCES.md 465/230/118 figures reproduced exactly under the stated Rule-A harmonization ""(S1∪S5 union framing) — the harmonization rule is now explicit and re-derivable, not just cited.")
Reproduction of the DATA_SOURCES.md D'Arcy-vs-Raghunath cross-check (S1∪S5 union, n=508):
Rule A (168 literal tokens) — D'Arcy genes absent from Raghunath: 465 (manifest states 465)
Rule B (267-gene HGNC-expanded) — D'Arcy genes absent from Raghunath: 463 (delta vs Rule A: 2)
disease-flagged among Rule-A non-overlap: 230 (manifest states 230)
hypopigmentation-class among those disease-flagged: 118 (manifest states 118)
DATA_SOURCES.md 465/230/118 figures reproduced exactly under the stated Rule-A harmonization (S1∪S5 union framing) — the harmonization rule is now explicit and re-derivable, not just cited.
Show code
# --- emit the node-level membership table (data/processed/nb5_gene_set_membership.csv) ---RAGHUNATH_CITATION ="gene_network_nodes.csv (per-gene UniProt/Entrez accession; Raghunath et al. 2015 BMC Res Notes)"DARCY_S1_CITATION =f"{DARCY_PMCID} (D'Arcy & Kiel 2023 Table S1)"KEGG_CITATION ="KEGG:hsa04916"# Per-gene Reactome citation resolves to whichever curated-union component(s) the gene came from,# so a membership row cites the specific accession, not just "Reactome" generically.gene_to_reactome_components = {}for comp in reactome_union_snap["components"]:for g in comp["genes"]: gene_to_reactome_components.setdefault(g, []).append(comp["stId"])membership_rows = []for g insorted(full_union): in_rag = g in raghunath_168 in_darcy = g in darcy_s1_genes in_kegg = g in kegg_genes in_reactome = g in reactome_genes_broad in_reactome_precision = g in reactome_genes cites = []if in_rag: cites.append(RAGHUNATH_CITATION)if in_darcy: cites.append(DARCY_S1_CITATION)if in_kegg: cites.append(KEGG_CITATION)if in_reactome: cites.append("Reactome:"+"+".join(gene_to_reactome_components.get(g, []))) classes =sorted(s1_gene_class.get(g, [])) if in_darcy else [] membership_rows.append({"gene": g, "in_raghunath_168": in_rag, "in_darcy_s1_243": in_darcy,"in_kegg_hsa04916": in_kegg, "in_reactome_curated_union": in_reactome,"in_reactome_r_hsa_5662702_precision": in_reactome_precision,"n_sets": sum([in_rag, in_darcy, in_kegg, in_reactome]),"darcy_phenotype_class": "|".join(classes),"citation": "|".join(cites), "citation_source": "accession_or_pmcid", })membership_df = pd.DataFrame(membership_rows)print(f"nb5_gene_set_membership.csv: {len(membership_df)} genes; n_sets distribution:")print(membership_df["n_sets"].value_counts().sort_index().to_string())membership_df.to_csv(PROC /"nb5_gene_set_membership.csv", index=False)
Step 5b — The Reactome-only connection set: a discovery/rescue substrate for NB7
Reframing Reactome’s role (per PI direction). Reactome’s curated union (Step 3b, 156 genes) is not meant to sit as a redundant fourth co-equal beside KEGG — KEGG hsa04916 (101 genes) is the more thorough, tighter mechanistic core. Reactome’s broader, noisier curated union is used instead as a discovery / rescue-connection layer: the substrate for testing, in a later notebook (NB7), whether a GWAS/paper locus that genetics found by statistical association but whose effector gene is uncertain falls into curated pigmentation reaction-space that the mechanism-first network (Raghunath + KEGG) does not contain.
The connection set, defined.Reactome-only connection set = genes in the Step 3b curated Reactome union that are in neither Raghunath-168 nor KEGG-hsa04916 — i.e., pigmentation-relevant genes the tight mechanistic core does not cover, that a mechanism-first effector-uncertain locus could plausibly resolve into. This is computed and characterized below, and is preserved with its full sub-pathway/reaction structure (not just a flat gene list) in a new frozen snapshot, data/external/db_responses/reactome_mitf_subpathway_structure.json, so that a downstream notebook can cite the specific reaction or sub-pathway a rescued locus connects through, not just “Reactome” generically.
Show code
mechanistic_core = gene_sets["Raghunath_168"] | gene_sets["KEGG_hsa04916_101"]reactome_only_connection_set = gene_sets["Reactome_curated_union_156"] - mechanistic_coreprint(f"Mechanistic core (Raghunath-168 ∪ KEGG-101): n={len(mechanistic_core)}")print(f"Reactome curated union (Step 3b): n={len(gene_sets['Reactome_curated_union_156'])}")print(f"Reactome-ONLY connection set (in Reactome, NOT in Raghunath, NOT in KEGG): "f"n={len(reactome_only_connection_set)}")
Mechanistic core (Raghunath-168 ∪ KEGG-101): n=240
Reactome curated union (Step 3b): n=156
Reactome-ONLY connection set (in Reactome, NOT in Raghunath, NOT in KEGG): n=121
Show code
# --- Characterize the connection set: which Reactome sub-pathway did each gene come from? ---mitf_subpathway_struct = json.load(open(FROZEN /"reactome_mitf_subpathway_structure.json"))# Per-component attribution (which of the 4 curated-union accessions contributed the gene)gene_to_union_component = {}for comp in reactome_union_snap["components"]:for g in comp["genes"]: gene_to_union_component.setdefault(g, []).append(comp["displayName"])from collections import Countercomponent_breakdown = Counter()for g in reactome_only_connection_set:for name in gene_to_union_component.get(g, []): component_breakdown[name] +=1print("Connection-set genes by source accession (Step 3b component):")for name, n in component_breakdown.most_common():print(f" {name}: {n}")# Finer sub-pathway breakdown, for genes sourced from the MITF-M pathway (the dominant contributor)subpathway_breakdown = Counter()gene_to_subpathway_names = {}for stid, info in mitf_subpathway_struct["subpathway_gene_membership"].items():for g in info["genes"]:if g in reactome_only_connection_set: subpathway_breakdown[info["name"]] +=1 gene_to_subpathway_names.setdefault(g, []).append(info["name"])print(f"\nOf the MITF-M-sourced connection-set genes, breakdown by fine sub-pathway "f"(a gene may appear in >1 sub-pathway):")for name, n in subpathway_breakdown.most_common():print(f" {name}: {n}")
Connection-set genes by source accession (Step 3b component):
MITF-M-regulated melanocyte development: 118
Melanocortin receptors: 4
Melanin biosynthesis: 1
Defective SLC24A5 causes oculocutaneous albinism 6 (OCA6): 1
Of the MITF-M-sourced connection-set genes, breakdown by fine sub-pathway (a gene may appear in >1 sub-pathway):
Transcriptional and post-translational regulation of MITF-M expression and activity: 43
Regulation of MITF-M-dependent genes involved in pigmentation: 29
Regulation of MITF-M-dependent genes involved in apoptosis: 16
Regulation of MITF-M-dependent genes involved in lysosome biogenesis and autophagy: 11
Regulation of MITF-M-dependent genes involved in extracellular matrix, focal adhesion and epithelial-to-mesenchymal transition: 10
Regulation of MITF-M-dependent genes involved in cell cycle and proliferation: 8
Regulation of MITF-M-dependent genes involved in DNA replication, damage repair and senescence: 5
Regulation of MITF-M dependent genes involved in invasion: 3
Regulation of MITF-M dependent genes involved in metabolism: 2
Show code
# --- Cross-check against known pigmentation-disease association: D'Arcy's own OMIM disease-gene table ---connection_in_darcy_omim = reactome_only_connection_set & darcy_s1_genesprint(f"Connection-set genes ALSO flagged as OMIM pigmentation-disease genes in D'Arcy S1 "f"(known association, just missing from the mechanistic core): n={len(connection_in_darcy_omim)}")print(sorted(connection_in_darcy_omim))print(f"\n({len(connection_in_darcy_omim)}/{len(reactome_only_connection_set)} = "f"{len(connection_in_darcy_omim)/len(reactome_only_connection_set):.1%} of the connection set already ""has an independent disease-association citation — these are the strongest candidate rescue ""connections, since two independent evidence types (Reactome reaction membership + OMIM disease ""curation) already agree on them.)")# --- emit the connection-set processed CSV, with per-gene sub-pathway/reaction/OMIM provenance ---connection_rows = []for g insorted(reactome_only_connection_set): comp_names_list = gene_to_union_component.get(g, []) subpaths = gene_to_subpathway_names.get(g, []) reactions = [d["stId"] for d in mitf_subpathway_struct["gene_to_reactions_within_pigmentation_subpathway"].get(g, [])] connection_rows.append({"gene": g,"reactome_source_accessions": "|".join(comp_names_list),"mitf_subpathways": "|".join(subpaths),"specific_reactions_within_pigmentation_subpathway": "|".join(reactions),"also_in_darcy_omim_disease_table": g in connection_in_darcy_omim,"citation": "Reactome curated union (data/external/db_responses/reactome_pigmentation_curated_union.json)","citation_source": "accession", })connection_df = pd.DataFrame(connection_rows)connection_df.to_csv(PROC /"nb5_reactome_only_connection_set.csv", index=False)print(f"\nwrote nb5_reactome_only_connection_set.csv: {len(connection_df)} rows")
Connection-set genes ALSO flagged as OMIM pigmentation-disease genes in D'Arcy S1 (known association, just missing from the mechanistic core): n=13
['CDKN2A', 'EDN3', 'GPR143', 'IRF4', 'MC2R', 'MLPH', 'MYO5A', 'RAB27A', 'SLC24A5', 'SLC45A2', 'SNAI2', 'TERT', 'TNFSF11']
(13/121 = 10.7% of the connection set already has an independent disease-association citation — these are the strongest candidate rescue connections, since two independent evidence types (Reactome reaction membership + OMIM disease curation) already agree on them.)
wrote nb5_reactome_only_connection_set.csv: 121 rows
Step 6 — The D’Arcy-vs-Raghunath non-overlap sub-finding (concrete and falsifiable)
Sub-finding, stated as a falsifiable claim. Two independent curation efforts for the same phenotype (pigmentation) — Raghunath’s manually curated melanogenesis-signaling mechanism and D’Arcy’s OMIM-mined disease-gene table — capture substantially non-overlapping gene sets: only 16 of D’Arcy’s 243 genes (6.6%) also appear in Raghunath’s 168-gene mechanistic set. The claim under test: the D’Arcy-only genes (in D’Arcy-S1 but not Raghunath) are enriched for OMIM disease annotation and, specifically, for the hypopigmentation phenotype class, relative to the 16 genes both sources agree on. This is falsifiable by comparing the hypopigmentation-class fraction in the D’Arcy-only set against the D’Arcy∩Raghunath overlap set with a two-sided Fisher’s exact test.
Show code
darcy_only = each_only["DArcy_S1_243"]overlap_16 = gene_sets["Raghunath_168"] & gene_sets["DArcy_S1_243"]print(f"D'Arcy-S1-only (not in Raghunath, KEGG, or the broadened Reactome curated union): n={len(darcy_only)}")print(f"D'Arcy ∩ Raghunath (both curation efforts agree): n={len(overlap_16)} -> {sorted(overlap_16)}")darcy_only_rows = []for g insorted(darcy_only): classes =sorted(s1_gene_class.get(g, [])) darcy_only_rows.append({"gene": g, "phenotype_classes": "|".join(classes),"is_hypopigmentation": "Hypopigmentation"in classes,"is_hyperpigmentation": "Hyperpigmentation"in classes,"is_mixed": "Mixed"in classes,"citation": DARCY_S1_CITATION, "citation_source": "pmcid", })darcy_only_df = pd.DataFrame(darcy_only_rows)darcy_only_df.to_csv(PROC /"nb5_darcy_only_phenotype_overlay.csv", index=False)n_hypo_only =int(darcy_only_df["is_hypopigmentation"].sum())n_hypo_overlap =sum(1for g in overlap_16 if"Hypopigmentation"in s1_gene_class.get(g, set()))frac_hypo_only = n_hypo_only /len(darcy_only)frac_hypo_overlap = n_hypo_overlap /len(overlap_16)print(f"\nHypopigmentation-class among D'Arcy-only ({len(darcy_only)}): "f"{n_hypo_only} ({frac_hypo_only:.1%})")print(f"Hypopigmentation-class among D'Arcy∩Raghunath overlap ({len(overlap_16)}): "f"{n_hypo_overlap} ({frac_hypo_overlap:.1%})")table = [[n_hypo_only, len(darcy_only) - n_hypo_only], [n_hypo_overlap, len(overlap_16) - n_hypo_overlap]]odds_ratio, p_value = fisher_exact(table)print(f"\nFisher's exact test (2x2: [D'Arcy-only, D'Arcy∩Raghunath] x [hypo, not-hypo]): "f"OR={odds_ratio:.3f}, p={p_value:.4f}")direction ="LESS"if frac_hypo_only < frac_hypo_overlap else"MORE"print(f"\nResult: the point estimate runs in the DIRECTION the sub-finding claims is unexpected — "f"hypopigmentation-class is actually {direction} common among the D'Arcy-only genes ({frac_hypo_only:.1%}) "f"than among the 16 genes both curations agree on ({frac_hypo_overlap:.1%}) — and the difference "f"does not reach the conventional p<0.05 threshold at these sample sizes (p={p_value:.2f}). The ""falsifiable claim as stated (non-overlap enriched for hypopigmentation) is NOT supported by this ""test; report the numbers as measured rather than the a-priori-plausible direction.")
D'Arcy-S1-only (not in Raghunath, KEGG, or the broadened Reactome curated union): n=211
D'Arcy ∩ Raghunath (both curation efforts agree): n=16 -> ['EDNRB', 'FASLG', 'HRAS', 'KIT', 'KITLG', 'MAP2K1', 'MC1R', 'MITF', 'OCA2', 'PAX3', 'POMC', 'PPP3CA', 'RAF1', 'SOX10', 'TYR', 'TYRP1']
Hypopigmentation-class among D'Arcy-only (211): 107 (50.7%)
Hypopigmentation-class among D'Arcy∩Raghunath overlap (16): 12 (75.0%)
Fisher's exact test (2x2: [D'Arcy-only, D'Arcy∩Raghunath] x [hypo, not-hypo]): OR=0.343, p=0.0720
Result: the point estimate runs in the DIRECTION the sub-finding claims is unexpected — hypopigmentation-class is actually LESS common among the D'Arcy-only genes (50.7%) than among the 16 genes both curations agree on (75.0%) — and the difference does not reach the conventional p<0.05 threshold at these sample sizes (p=0.07). The falsifiable claim as stated (non-overlap enriched for hypopigmentation) is NOT supported by this test; report the numbers as measured rather than the a-priori-plausible direction.
Show code
# A cleaner baseline: compare each group's hypopigmentation rate against the full 243-gene S1 rate,# rather than only the two groups against each other (n=16 in the overlap group is small).print(f"Hypopigmentation-class rate, full D'Arcy S1 (n=243): {len(hypo_genes_s1)}/{243} = {len(hypo_genes_s1)/243:.1%}")print(f"Hypopigmentation-class rate, D'Arcy-only (n={len(darcy_only)}): {frac_hypo_only:.1%}")print(f"Hypopigmentation-class rate, D'Arcy∩Raghunath overlap (n={len(overlap_16)}): {frac_hypo_overlap:.1%}")n_hypo_rest_of_s1 =len(hypo_genes_s1) - n_hypo_overlapn_rest_of_s1 =243-len(overlap_16)frac_hypo_rest_of_s1 = n_hypo_rest_of_s1 / n_rest_of_s1table2 = [[n_hypo_overlap, len(overlap_16) - n_hypo_overlap], [n_hypo_rest_of_s1, n_rest_of_s1 - n_hypo_rest_of_s1]]odds2, p2 = fisher_exact(table2)print(f"\nFisher's exact (overlap-with-Raghunath vs rest-of-S1, hypo-class): OR={odds2:.3f}, p={p2:.4f}")n_darcy_only_of_243 =len(darcy_only) # darcy_only is already a subset of the 243-gene S1 setfrac_darcy_only_of_243 = n_darcy_only_of_243 /243print(f"\nREVISED, better-supported sub-finding: it is the {len(overlap_16)} genes where Raghunath's "f"independent mechanistic curation AGREES with D'Arcy's OMIM mining — not the {len(darcy_only)} "f"D'Arcy-only genes — that are enriched for the hypopigmentation phenotype class "f"({frac_hypo_overlap:.1%} vs {frac_hypo_rest_of_s1:.1%} in the rest of the S1 table), though this "f"too does not reach the conventional significance threshold at n={len(overlap_16)} (p={p2:.2f}). ""Read together: mechanism curation and disease curation agree preferentially on classic ""albinism/oculocutaneous-hypopigmentation genes (OCA2, TYR, TYRP1, MC1R, KIT, KITLG, SOX10, PAX3, ""MITF, ...) — the loci with the clearest causal phenotype — while D'Arcy's disease-only curation adds "f"a further {n_darcy_only_of_243} genes (mixed hypo/hyper/disease-general) that a mechanism-first ""curation effort like Raghunath's, built from melanogenesis signaling literature rather than OMIM, "f"does not capture at all. The non-overlap is real and large ({n_darcy_only_of_243}/243 = "f"{frac_darcy_only_of_243:.1%} of D'Arcy's disease genes not also in Raghunath's set); its ""phenotype-class composition is not distinguishable from D'Arcy's own overall composition, so the ""claim is about a coverage gap, not about a hypopigmentation-specific gap.")
Hypopigmentation-class rate, full D'Arcy S1 (n=243): 130/243 = 53.5%
Hypopigmentation-class rate, D'Arcy-only (n=211): 50.7%
Hypopigmentation-class rate, D'Arcy∩Raghunath overlap (n=16): 75.0%
Fisher's exact (overlap-with-Raghunath vs rest-of-S1, hypo-class): OR=2.771, p=0.1176
REVISED, better-supported sub-finding: it is the 16 genes where Raghunath's independent mechanistic curation AGREES with D'Arcy's OMIM mining — not the 211 D'Arcy-only genes — that are enriched for the hypopigmentation phenotype class (75.0% vs 52.0% in the rest of the S1 table), though this too does not reach the conventional significance threshold at n=16 (p=0.12). Read together: mechanism curation and disease curation agree preferentially on classic albinism/oculocutaneous-hypopigmentation genes (OCA2, TYR, TYRP1, MC1R, KIT, KITLG, SOX10, PAX3, MITF, ...) — the loci with the clearest causal phenotype — while D'Arcy's disease-only curation adds a further 211 genes (mixed hypo/hyper/disease-general) that a mechanism-first curation effort like Raghunath's, built from melanogenesis signaling literature rather than OMIM, does not capture at all. The non-overlap is real and large (211/243 = 86.8% of D'Arcy's disease genes not also in Raghunath's set); its phenotype-class composition is not distinguishable from D'Arcy's own overall composition, so the claim is about a coverage gap, not about a hypopigmentation-specific gap.
Step 7 — Our own consistent STRING layer (Module 2)
STRING is applied identically to every gene set — same connector, same species (9606, human), same required_score threshold (700, “high confidence”) — via the protein-annotation MCP connector’s get_string_network. MCP calls only run in the repl tool, not in a notebook kernel, so this follows the frozen-DB pattern: the pull covers all five gene sets (Raghunath-168, D’Arcy-S1-243, KEGG-101, Reactome-5, and their union), frozen verbatim to data/external/db_responses/string_network_pulls_v12.json, and is loaded here offline.
Query (frozen, re-runnable).host.mcp("protein-annotation", "get_string_network", symbols=<gene_list>, species=9606, required_score=700) — one call per gene set, run from the repl tool (REQUERY=True below documents the exact call; it does not execute in this kernel).
Show code
REQUERY_STRING =Falseif REQUERY_STRING:# Run this block from the `repl` tool (MCP is not reachable from a notebook kernel):## results = {}# for name, symbols in gene_sets_dict.items():# results[name] = host.mcp("protein-annotation", "get_string_network",# symbols=symbols, species=9606, required_score=700)# json.dump({"query_utc": ..., "required_score": 700, "species": 9606, "results": results},# open(FROZEN / "string_network_pulls_v12.json", "w"))passstring_pulls = json.load(open(FROZEN /"string_network_pulls_v12.json"))print("STRING pull query-UTC:", string_pulls["query_utc"], "| required_score:", string_pulls["required_score"])for name, resp in string_pulls["results"].items(): ver = resp.get("string_version", {})print(f" {name}: nodes={len(resp.get('nodes', []))}, edges={len(resp.get('edges', []))}, "f"unmapped={resp.get('unmapped', [])}, STRING v{ver.get('string_version')}")
Step 8 — Edge-level check: does STRING see Raghunath’s signed/directed edges?
For every unique undirected gene pair connected by an edge in Raghunath’s gene-level projection (data/processed/gene_network_edges.csv, 309 mechanistic edges collapsed to 292 unique undirected gene pairs), check whether the same pair appears as a STRING edge (score ≥ 700) in our raghunath_168 pull. STRING is undirected/unsigned association evidence — this is a presence/absence coverage check, not a sign agreement check (STRING carries no sign).
Show code
gene_network_edges = pd.read_csv(PROC /"gene_network_edges.csv")string_rag = string_pulls["results"]["raghunath_168"]mapped_names =set(n["name"].upper() for n in string_rag["nodes"])string_edge_pairs =set(frozenset([e["a"].upper(), e["b"].upper()]) for e in string_rag["edges"])coverage_rows, seen = [], set()for _, row in gene_network_edges.iterrows(): a, b = row["source"].upper(), row["target"].upper()if a == b:continue pair_key =tuple(sorted([a, b]))if pair_key in seen:continue seen.add(pair_key) both_mapped = a in mapped_names and b in mapped_names present = (frozenset(pair_key) in string_edge_pairs) if both_mapped elseNone coverage_rows.append({"gene_a": pair_key[0], "gene_b": pair_key[1], "raghunath_sign": row["sign"],"raghunath_citation": row["citation"], "both_string_mappable": both_mapped,"string_edge_present_score700": present,"citation": f"STRING v12.0 ({string_pulls['query_utc']}); {row['citation']}","citation_source": "string_db+inherited_backbone"})coverage_df = pd.DataFrame(coverage_rows)coverage_df.to_csv(PROC /"nb5_raghunath_string_edge_coverage.csv", index=False)n_pairs =len(coverage_df)n_mappable =int(coverage_df["both_string_mappable"].sum())n_present =int((coverage_df["string_edge_present_score700"] ==True).sum())print(f"Raghunath unique undirected gene pairs: {n_pairs}")print(f" both endpoints STRING-mappable: {n_mappable}")print(f" present as a STRING edge (score>=700): {n_present} ({n_present/n_mappable:.1%} of mappable pairs)")print(f" absent from STRING at this threshold: {n_mappable - n_present} ({(n_mappable - n_present)/n_mappable:.1%})")print(f" not STRING-mappable at all (one/both endpoints unmapped by STRING v12.0): {n_pairs - n_mappable}")
Raghunath unique undirected gene pairs: 292
both endpoints STRING-mappable: 265
present as a STRING edge (score>=700): 158 (59.6% of mappable pairs)
absent from STRING at this threshold: 107 (40.4%)
not STRING-mappable at all (one/both endpoints unmapped by STRING v12.0): 27
Step 9 — Drift: our fresh STRING pull vs D’Arcy’s own frozen STRING snapshot
D’Arcy’s Table S4 is itself a STRING pull, made at combined_score > 0.7 sometime before the paper’s 2023 publication. Comparing it to a fresh STRING v12.0 pull at required_score=700 (same threshold, different database version/time) quantifies drift rather than treating STRING as a single ground truth — different STRING releases and re-computed confidence scores will disagree on some edges, and that disagreement should be measured, not hidden. The comparison is restricted to gene pairs where both genes are in the 243 S1 disease-gene set, since our darcy_s1_243 pull queried exactly that set (D’Arcy’s S4 network was built by expanding S1 to 451/452 nodes, so a like-for-like comparison needs the same node scope).
Show code
darcy_s4_pairs = {}for _, row in darcy_s4.iterrows(): a, b = row["node1"], row["node2"]if a in darcy_s1_genes and b in darcy_s1_genes and a != b: darcy_s4_pairs[frozenset([a, b])] = row["combined_score"]our_darcy_pull = string_pulls["results"]["darcy_s1_243"]our_darcy_pairs = {frozenset([e["a"].upper(), e["b"].upper()]): e["score"] for e in our_darcy_pull["edges"]}print(f"D'Arcy S4 edges with both endpoints in S1-243: {len(darcy_s4_pairs)}")print(f"Our fresh STRING pull (score>=700) edges among S1-243: {len(our_darcy_pairs)}")all_pairs =set(darcy_s4_pairs) |set(our_darcy_pairs)in_both =set(darcy_s4_pairs) &set(our_darcy_pairs)only_darcy =set(darcy_s4_pairs) -set(our_darcy_pairs)only_ours =set(our_darcy_pairs) -set(darcy_s4_pairs)print(f" edges in both: {len(in_both)}")print(f" edges only in D'Arcy's S4: {len(only_darcy)} (their pull found it, ours did not at score>=700)")print(f" edges only in our fresh pull: {len(only_ours)} (ours found it, their S4 snapshot did not)")print(f" Jaccard agreement: {len(in_both)/len(all_pairs):.1%}")drift_rows = []for pair insorted(all_pairs, key=lambda x: tuple(sorted(x))): a, b =sorted(pair) drift_rows.append({"gene_a": a, "gene_b": b,"darcy_s4_combined_score": darcy_s4_pairs.get(pair),"our_string_v12_score": our_darcy_pairs.get(pair),"in_darcy_s4": pair in darcy_s4_pairs, "in_our_pull": pair in our_darcy_pairs,"agreement": (pair in darcy_s4_pairs) == (pair in our_darcy_pairs),"citation": f"{DARCY_PMCID} (S4) vs STRING v12.0 ({string_pulls['query_utc']})","citation_source": "pmcid+string_db"})drift_df = pd.DataFrame(drift_rows)drift_df.to_csv(PROC /"nb5_string_drift_darcy_vs_ourpull.csv", index=False)print(f"\nDrift is substantial (~{100-len(in_both)/len(all_pairs)*100:.0f}% of the union disagrees) even at ""the same nominal confidence threshold, consistent with STRING version drift (D'Arcy et al. used an ""earlier STRING release; this pull is v12.0) plus possible differences in how each pull's input gene ""list was mapped to STRING identifiers. This is documented as measured drift, not resolved — a reader ""combining D'Arcy's PPI shell with a fresh STRING pull should expect ~1/3 of edges not to reproduce.")
D'Arcy S4 edges with both endpoints in S1-243: 590
Our fresh STRING pull (score>=700) edges among S1-243: 641
edges in both: 491
edges only in D'Arcy's S4: 99 (their pull found it, ours did not at score>=700)
edges only in our fresh pull: 150 (ours found it, their S4 snapshot did not)
Jaccard agreement: 66.4%
Drift is substantial (~34% of the union disagrees) even at the same nominal confidence threshold, consistent with STRING version drift (D'Arcy et al. used an earlier STRING release; this pull is v12.0) plus possible differences in how each pull's input gene list was mapped to STRING identifiers. This is documented as measured drift, not resolved — a reader combining D'Arcy's PPI shell with a fresh STRING pull should expect ~1/3 of edges not to reproduce.
Step 9b — Bajpai et al. 2023 CRISPR screen as a NODE layer
Bajpai et al. 2023’s genome-wide CRISPR screen enters this comparison as a node layer, not an edge network. The screen assayed melanin content via FACS sorting on side-scatter and reports, per gene, a Combined_casTLE_Effect (magnitude/direction) and a q_value; 169 genes pass the paper’s own q<0.10 hit threshold, and all 169 have the same sign: knockout/knockdown reduces pigmentation (direction_note). There is no gene-gene edge in this source — casTLE effect is a phenotype-perturbation weight on a single gene, not a regulatory-edge strength, and it is not treated as one anywhere below.
Representation used (primary representation): Bajpai is added to nb5_gene_set_membership.csv as (i) a node-set flag (in_bajpai2023_crispr_169) for membership/overlap, and (ii) a node weight + sign (bajpai_castle_effect, uniform bajpai_direction="reduces_pigmentation"). A SECONDARY, clearly labeled bipartite hit→melanin-endpoint layer is also built below (optional) — it is not topology-comparable to the four edge/pathway networks and is kept in a separate output file.
Sanity check (must reproduce the expected anchors):TYR, DCT, SLC45A2, OCA2 must be hits; MC1R and HERC2 must NOT be (the screen finds melanogenesis effectors, not eye/hair-color regulators like MC1R/HERC2, which sit upstream of pigment-cell fate rather than in the melanin-synthesis machinery itself).
Show code
# --- Step 9b: Bajpai 2023 CRISPR-hit node layer ---BAJPAI_DOI ="10.1126/science.ade6289"BAJPAI_PMCID ="PMC10901463"BAJPAI_CITATION =f"{BAJPAI_PMCID} (Bajpai et al. 2023 Science 381:eade6289, DOI {BAJPAI_DOI}); Table S1, q<0.10 hit threshold"bajpai = pd.read_csv(PROC /"bajpai2023_crispr_hits.csv")assertlen(bajpai) == bajpai["Symbol"].nunique() ==169, f"expected 169 unique Bajpai hit genes, got {len(bajpai)}"bajpai_symbols =set(bajpai["Symbol"].astype(str).str.strip().str.upper())bajpai_indexed = bajpai.assign(Symbol=bajpai["Symbol"].astype(str).str.strip().str.upper()).set_index("Symbol")assertset(bajpai["direction_note"].str.contains("reduces_pigmentation").unique()) == {True}, \"expected all 169 hits to share the reduces-pigmentation direction"# --- sanity check anchors ---pos_controls = {"TYR", "DCT", "SLC45A2", "OCA2"}neg_controls = {"MC1R", "HERC2"}pos_check = {g: g in bajpai_symbols for g in pos_controls}neg_check = {g: g in bajpai_symbols for g in neg_controls}print("Sanity check — positive controls present (all must be True):", pos_check)print("Sanity check — negative controls absent (all must be False):", neg_check)assertall(pos_check.values()), f"Bajpai join sanity check FAILED (missing positive control): {pos_check}"assertnotany(neg_check.values()), f"Bajpai join sanity check FAILED (unexpected negative control hit): {neg_check}"print("Sanity check PASSED: TYR/DCT/SLC45A2/OCA2 are hits; MC1R/HERC2 are not.")# --- reverse coverage: how many of the 169 hits fall in each of the 4 candidate networks, and how many are orphans ---reverse_coverage = {}for name, s in gene_sets.items(): n_hit_in_net =len(bajpai_symbols & s) reverse_coverage[name] = n_hit_in_netprint(f" Bajpai hits found in {name} (n={len(s)} network genes): {n_hit_in_net}")orphan_hits =sorted(bajpai_symbols - full_union)print(f"\nBajpai orphan hits (in none of the 4 candidate networks): {len(orphan_hits)} / {len(bajpai_symbols)}")print(f" (these are candidate discovery nodes -- causally implicated by the screen, absent from all 4 curated/pathway sources)")print(" first 15 orphans:", orphan_hits[:15])# --- simple enrichment: are Bajpai hits over-represented in network X vs. the screen's own assayed background? ---screen_universe =set(bajpai_indexed.index) |set( pd.read_excel(ROOT /"data"/"raw"/"bajpai2023"/"science.ade6289_table_s1.xlsx")["Symbol"] .astype(str).str.strip().str.upper())print(f"\nScreen-assayed gene universe (Table S1, all 4,956 tested genes): n={len(screen_universe)}")enrich_rows = []for name, s in gene_sets.items(): net_in_screen = s & screen_universe hits_in_net =len(bajpai_symbols & net_in_screen) a = hits_in_net b =len(net_in_screen) - a c =len(bajpai_symbols - net_in_screen) d =len(screen_universe - net_in_screen - bajpai_symbols) odds, p = fisher_exact([[a, b], [c, d]]) enrich_rows.append({"network": name, "network_n_genes": len(s), "network_genes_in_screen_universe": len(net_in_screen),"bajpai_hits_in_network": a,"pct_of_network_genes_that_are_hits": round(100* a /len(net_in_screen), 1) if net_in_screen elseNone,"odds_ratio": round(odds, 2), "fisher_p": p,"citation": BAJPAI_CITATION, "citation_source": "pmcid", })enrichment_df = pd.DataFrame(enrich_rows)enrichment_df.to_csv(PROC /"nb5_bajpai_network_enrichment.csv", index=False)print("\nEnrichment of Bajpai hits within each candidate network (vs. screen-assayed background):")print(enrichment_df.drop(columns=["citation", "citation_source"]).to_string(index=False))print("\nCaveats carried forward: (1) q<0.10 is a threshold choice -- q<0.05 gives 149 hits, ""q<0.15 gives 208; enrichment direction is robust to this but exact counts shift. ""(2) the screen's asymmetry -- a gene absent from the 169 hits is NOT evidence it plays no role in ""pigmentation; it may have been non-significant, redundant with a paralog, or simply not assayed ""(only genes in the screen-assayed universe above were tested at all).")# --- 5-way node-set comparison (4 candidate networks + Bajpai), for the figure below ---# Bajpai enters here strictly as a NODE SET for membership/overlap (by design) -- this is a set# comparison, not an edge-topology comparison, so including it alongside the other four node sets is# consistent with "membership/overlap" in this project's scope, not a violation of the no-fabricated-# edges guardrail (no edge is created or implied here).gene_sets_5 =dict(gene_sets)gene_sets_5["Bajpai2023_crispr_169"] = bajpai_symbolsnames5 =list(gene_sets_5.keys())each_only_5 = {name: gene_sets_5[name] -set().union(*[gene_sets_5[o] for o in names5 if o != name]) for name in names5}full_intersection_5 =set.intersection(*gene_sets_5.values())print("\n5-way each-only sizes (Raghunath, D'Arcy-S1, KEGG, Reactome, Bajpai):")for name in names5:print(f" {name} ONLY: n={len(each_only_5[name])}")print(f"5-way full intersection (all 5 sources agree): n={len(full_intersection_5)} -> {sorted(full_intersection_5)}")print("(For reference, the canonical 4-way full intersection — Step 5, unaffected by Bajpai — remains "f"n={len(full_intersection)}: {sorted(full_intersection)}.)")
Sanity check — positive controls present (all must be True): {'TYR': True, 'OCA2': True, 'SLC45A2': True, 'DCT': True}
Sanity check — negative controls absent (all must be False): {'MC1R': False, 'HERC2': False}
Sanity check PASSED: TYR/DCT/SLC45A2/OCA2 are hits; MC1R/HERC2 are not.
Bajpai hits found in Raghunath_168 (n=168 network genes): 5
Bajpai hits found in DArcy_S1_243 (n=243 network genes): 23
Bajpai hits found in KEGG_hsa04916_101 (n=101 network genes): 4
Bajpai hits found in Reactome_curated_union_156 (n=156 network genes): 10
Bajpai orphan hits (in none of the 4 candidate networks): 142 / 169
(these are candidate discovery nodes -- causally implicated by the screen, absent from all 4 curated/pathway sources)
first 15 orphans: ['AC002398.9', 'ACTR2', 'ACTR3', 'ALKBH8', 'AMBRA1', 'ANKRD27', 'AP1G1', 'AP1M1', 'ARL8B', 'ARPC3', 'ARPC4', 'ARPC4-TTLL3', 'ASCC2', 'ASCC3', 'ATP6V0A2']
Screen-assayed gene universe (Table S1, all 4,956 tested genes): n=4950
Enrichment of Bajpai hits within each candidate network (vs. screen-assayed background):
network network_n_genes network_genes_in_screen_universe bajpai_hits_in_network pct_of_network_genes_that_are_hits odds_ratio fisher_p
Raghunath_168 168 46 5 10.9 3.52 1.934093e-02
DArcy_S1_243 243 67 23 34.3 16.96 6.015103e-18
KEGG_hsa04916_101 101 22 4 18.2 6.41 5.925533e-03
Reactome_curated_union_156 156 42 10 23.8 9.33 9.428509e-07
Caveats carried forward (per the brief): (1) q<0.10 is a threshold choice -- q<0.05 gives 149 hits, q<0.15 gives 208; enrichment direction is robust to this but exact counts shift. (2) the screen's asymmetry -- a gene absent from the 169 hits is NOT evidence it plays no role in pigmentation; it may have been non-significant, redundant with a paralog, or simply not assayed (only genes in the screen-assayed universe above were tested at all).
5-way each-only sizes (Raghunath, D'Arcy-S1, KEGG, Reactome, Bajpai):
Raghunath_168 ONLY: n=121
DArcy_S1_243 ONLY: n=195
KEGG_hsa04916_101 ONLY: n=64
Reactome_curated_union_156 ONLY: n=108
Bajpai2023_crispr_169 ONLY: n=142
5-way full intersection (all 5 sources agree): n=1 -> ['TYR']
(For reference, the canonical 4-way full intersection — Step 5, unaffected by Bajpai — remains n=7: ['EDNRB', 'KIT', 'MC1R', 'MITF', 'POMC', 'TYR', 'TYRP1'].)
Show code
# --- attach Bajpai node-set flag + weight + sign to the membership table, and add orphan-hit-only rows ---membership_df["in_bajpai2023_crispr_169"] = membership_df["gene"].isin(bajpai_symbols)orphan_rows = pd.DataFrame([{"gene": g, "in_raghunath_168": False, "in_darcy_s1_243": False, "in_kegg_hsa04916": False,"in_reactome_curated_union": False, "in_reactome_r_hsa_5662702_precision": False, "n_sets": 0,"darcy_phenotype_class": "", "citation": "", "citation_source": "","in_bajpai2023_crispr_169": True,} for g in orphan_hits])membership_df = pd.concat([membership_df, orphan_rows], ignore_index=True)membership_df["bajpai_castle_effect"] = membership_df["gene"].map(lambda g: bajpai_indexed.loc[g, "Combined_casTLE_Effect"] if g in bajpai_symbols elseNone)membership_df["bajpai_q_value"] = membership_df["gene"].map(lambda g: bajpai_indexed.loc[g, "q_value"] if g in bajpai_symbols elseNone)membership_df["bajpai_direction"] = membership_df["in_bajpai2023_crispr_169"].map(lambda x: "reduces_pigmentation"if x else"")membership_df["bajpai_orphan_hit"] = membership_df["gene"].isin(set(orphan_hits))membership_df["n_sets_incl_bajpai"] = membership_df["n_sets"].astype(int) + membership_df["in_bajpai2023_crispr_169"].astype(int)# citation: fill empty citations for Bajpai hits (the orphan rows), append for hits that already had oneno_citation = membership_df["citation"].isna() | (membership_df["citation"].astype(str).str.strip() =="")fill_mask = membership_df["in_bajpai2023_crispr_169"] & no_citationappend_mask = membership_df["in_bajpai2023_crispr_169"] &~no_citationmembership_df.loc[fill_mask, "citation"] = BAJPAI_CITATIONmembership_df.loc[fill_mask, "citation_source"] ="pmcid"membership_df.loc[append_mask, "citation"] = membership_df.loc[append_mask, "citation"].astype(str) +"|"+ BAJPAI_CITATIONmembership_df.to_csv(PROC /"nb5_gene_set_membership.csv", index=False)n_bajpai_flagged =int(membership_df["in_bajpai2023_crispr_169"].sum())n_orphan_flagged =int(membership_df["bajpai_orphan_hit"].sum())print(f"nb5_gene_set_membership.csv updated: {len(membership_df)} genes total "f"({len(full_union)} from the 4 candidate networks + {n_orphan_flagged} Bajpai-orphan-only rows added)")print(f"Bajpai hits flagged in membership table: {n_bajpai_flagged} (expected 169)")assert n_bajpai_flagged ==169, f"expected all 169 Bajpai hits flagged, got {n_bajpai_flagged}"print(f"Bajpai orphan hits (in none of the 4 networks): {n_orphan_flagged}")assert n_orphan_flagged ==len(orphan_hits)bad_citation = membership_df[membership_df["citation"].isna() | (membership_df["citation"].astype(str).str.strip() =="")]assertlen(bad_citation) ==0, f"citation gap after Bajpai merge: {len(bad_citation)} rows"print("Citation check after Bajpai merge: 0 uncited rows.")
nb5_gene_set_membership.csv updated: 714 genes total (572 from the 4 candidate networks + 142 Bajpai-orphan-only rows added)
Bajpai hits flagged in membership table: 169 (expected 169)
Bajpai orphan hits (in none of the 4 networks): 142
Citation check after Bajpai merge: 0 uncited rows.
Optional secondary layer — labeled bipartite hit → melanin-endpoint anchor
As a secondary representation, this layer is a directed, signed star linking each of the 169 hit genes to a single melanin_content endpoint node, weighted by casTLE effect, signed by the uniform reduces-pigmentation direction. This is not a gene-gene edge network and is never used for the topology comparisons above or below — it exists only to connect effectors to the phenotype they were screened against, kept in its own output file so it cannot be silently pooled with the real edge networks.
Show code
# --- optional SECONDARY layer: labeled bipartite hit -> melanin_content endpoint (NOT topology-comparable) ---bipartite_rows = []for _, row in bajpai.iterrows(): bipartite_rows.append({"source_gene": row["Symbol"].strip().upper(),"source_gene_ensembl": row["GeneID"],"target_endpoint": "melanin_content","layer_type": "bipartite_functional_perturbation_layer_NOT_topology_comparable","edge_weight_castle_effect": row["Combined_casTLE_Effect"],"edge_sign": "-","sign_basis": "knockout/knockdown of source_gene REDUCES melanin_content (all 169 hits share this direction)","q_value": row["q_value"],"citation": BAJPAI_CITATION,"citation_source": "pmcid", })bipartite_df = pd.DataFrame(bipartite_rows)bipartite_df.to_csv(PROC /"nb5_bajpai_bipartite_melanin_endpoint.csv", index=False)print(f"nb5_bajpai_bipartite_melanin_endpoint.csv: {len(bipartite_df)} hit->melanin_content star edges "f"(1 endpoint node, {bipartite_df['source_gene'].nunique()} source genes) -- "f"LABELED layer, excluded from every topology/overlap statistic in this notebook.")
nb5_bajpai_bipartite_melanin_endpoint.csv: 169 hit->melanin_content star edges (1 endpoint node, 169 source genes) -- LABELED layer, excluded from every topology/overlap statistic in this notebook.
Step 9c — Networks-typology table
Before the discussion, an explicit side-by-side of every candidate source this project has assembled so far — not just the four compared node-for-node in Step 5 — so a reader can see at a glance which sources carry real gene-gene topology (directed/signed edges) versus which are node-only membership or weight layers. This table is the single reference for “what kind of thing is each source.
Show code
# --- networks-typology table: one row per source x {type, directed?, signed?, edge evidence, node count} ---baxter_df = pd.read_csv(PROC /"baxter2018_650_pigmentation_genes.csv")n_baxter = baxter_df["Human gene symbol"].nunique()discordance_df = pd.read_csv(PROC /"discordance_loci.csv")n_gwas_loci_genes = discordance_df["nearest_gene_label"].nunique()n_gwas_papers = discordance_df["paper"].nunique()darcy_s4_nodes =len(set(darcy_s4["node1"]) |set(darcy_s4["node2"]))typology_rows = [ {"source": "Raghunath et al. 2015", "network_type": "manually curated mechanistic signaling network","directed": True, "signed": True, "edge_evidence": "literature-review synthesis, per-edge citation","node_count": len(raghunath_168), "edge_count": len(gene_network_edges),"citation": "gene_network_edges.csv (Raghunath et al. 2015, BMC Res Notes)"}, {"source": "KEGG hsa04916", "network_type": "reference pathway diagram","directed": True, "signed": False, "edge_evidence": "KEGG's own pathway-diagram curation (nodes extracted here, edges not)","node_count": len(kegg_genes), "edge_count": None, "citation": "KEGG:hsa04916"}, {"source": "Reactome curated union (Step 3b)", "network_type": "reference reaction pathway (4-accession itemized union)","directed": True, "signed": False,"edge_evidence": "Reactome's own reaction-level curation (nodes + sub-pathway membership extracted here, not a flat edge list)","node_count": len(reactome_genes_broad), "edge_count": None, "citation": "Reactome (4 accessions; see Step 3b)"}, {"source": "D'Arcy & Kiel 2023, Table S1", "network_type": "curated OMIM disease-gene table","directed": False, "signed": False, "edge_evidence": "OMIM phenotype-record curation (gene list, no edges)","node_count": len(darcy_s1_genes), "edge_count": 0, "citation": f"{DARCY_PMCID} (D'Arcy & Kiel 2023 Table S1)"}, {"source": "D'Arcy & Kiel 2023, Tables S4/S5 (STRING expansion)", "network_type": "protein-protein association network","directed": False, "signed": False,"edge_evidence": "STRING association scores (combined_score>0.7), frozen at 2023 publication time","node_count": darcy_s4_nodes, "edge_count": len(darcy_s4), "citation": f"{DARCY_PMCID} (D'Arcy & Kiel 2023 Table S4/S5)"}, {"source": "Our STRING v12.0 pull (Module 2, this notebook)", "network_type": "protein-protein association network","directed": False, "signed": False,"edge_evidence": "STRING v12.0, required_score=700, applied identically across 5 gene sets (this session)","node_count": None, "edge_count": None, "citation": "string_network_pulls_v12.json (STRING v12.0)"}, {"source": "Bajpai et al. 2023 CRISPR screen", "network_type": "weighted/signed NODE layer (no intrinsic gene-gene edges)","directed": False, "signed": True,"edge_evidence": "genome-wide CRISPR screen, casTLE effect + q<0.10 hit threshold; no gene-gene edges (see Step 9b)","node_count": len(bajpai_symbols), "edge_count": 0, "citation": BAJPAI_CITATION}, {"source": "Baxter et al. 2018/2019", "network_type": "curated cross-species gene list","directed": False, "signed": False,"edge_evidence": "author-curated pigmentation gene compilation cross-referenced to OMIM/MGI/ZFIN/GO (gene list, no edges)","node_count": n_baxter, "edge_count": 0, "citation": "baxter2018_650_pigmentation_genes.csv (Baxter et al. 2018/2019)"}, {"source": f"GWAS/curated pigmentation loci ({n_gwas_papers} papers)", "network_type": "curated association-locus table","directed": False, "signed": False,"edge_evidence": "GWAS/candidate-gene association literature, nearest/assigned gene per locus (loci, not gene-gene edges)","node_count": n_gwas_loci_genes, "edge_count": 0, "citation": "discordance_loci.csv (per-row citation; 13 source papers)"},]typology_df = pd.DataFrame(typology_rows)typology_df.to_csv(PROC /"nb5_networks_typology.csv", index=False)print(f"nb5_networks_typology.csv: {len(typology_df)} sources characterized")print(typology_df.drop(columns=["citation"]).to_string(index=False))
nb5_networks_typology.csv: 9 sources characterized
source network_type directed signed edge_evidence node_count edge_count
Raghunath et al. 2015 manually curated mechanistic signaling network True True literature-review synthesis, per-edge citation 168.0 309.0
KEGG hsa04916 reference pathway diagram True False KEGG's own pathway-diagram curation (nodes extracted here, edges not) 101.0 NaN
Reactome curated union (Step 3b) reference reaction pathway (4-accession itemized union) True False Reactome's own reaction-level curation (nodes + sub-pathway membership extracted here, not a flat edge list) 156.0 NaN
D'Arcy & Kiel 2023, Table S1 curated OMIM disease-gene table False False OMIM phenotype-record curation (gene list, no edges) 243.0 0.0
D'Arcy & Kiel 2023, Tables S4/S5 (STRING expansion) protein-protein association network False False STRING association scores (combined_score>0.7), frozen at 2023 publication time 452.0 4668.0
Our STRING v12.0 pull (Module 2, this notebook) protein-protein association network False False STRING v12.0, required_score=700, applied identically across 5 gene sets (this session) NaN NaN
Bajpai et al. 2023 CRISPR screen weighted/signed NODE layer (no intrinsic gene-gene edges) False True genome-wide CRISPR screen, casTLE effect + q<0.10 hit threshold; no gene-gene edges (see Step 9b) 169.0 0.0
Baxter et al. 2018/2019 curated cross-species gene list False False author-curated pigmentation gene compilation cross-referenced to OMIM/MGI/ZFIN/GO (gene list, no edges) 635.0 0.0
GWAS/curated pigmentation loci (13 papers) curated association-locus table False False GWAS/candidate-gene association literature, nearest/assigned gene per locus (loci, not gene-gene edges) 54.0 0.0
Step 10 — Citation-completeness gate
Every node/edge table this notebook emits carries a citation (+ citation_source) column. This closing gate lists any row with an empty citation and fails the run if one exists — the same tier of check as “0 unresolved nodes” elsewhere in this project.
Figure — set sizes, overlaps, and STRING edge coverage
Show code
import matplotlibmatplotlib.use("Agg")import matplotlib.pyplot as pltFIGDIR = ROOT /"notebooks"/"figures"FIGDIR.mkdir(exist_ok=True)fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.8))# Panel A: each-only sizes + full intersection, 5-way (Raghunath, D'Arcy-S1, KEGG, Reactome, Bajpai)ax = axes[0]labels = ["Raghunath\nonly", "D'Arcy-S1\nonly", "KEGG\nonly", "Reactome\n(curated)\nonly","Bajpai\n2023\nonly", f"All 5\n(n={len(full_intersection_5)})"]values = [len(each_only_5["Raghunath_168"]), len(each_only_5["DArcy_S1_243"]),len(each_only_5["KEGG_hsa04916_101"]), len(each_only_5["Reactome_curated_union_156"]),len(each_only_5["Bajpai2023_crispr_169"]), len(full_intersection_5)]bars = ax.bar(labels, values, color=["#4C72B0", "#DD8452", "#55A868", "#C44E52", "#937860", "#8172B2"])for b, v inzip(bars, values): ax.text(b.get_x() + b.get_width() /2, v +3, str(v), ha="center", fontsize=9)ax.set_ylabel("number of genes")ax.set_title("Node-level set comparison\n(each-only vs. full 5-way intersection)")ax.set_ylim(0, max(values) *1.2)# Panel B: STRING as a controlled layer -- coverage of Raghunath's edges AND drift vs D'Arcy's own STRING pullax = axes[1]cov_present =int((coverage_df["string_edge_present_score700"] ==True).sum())cov_absent =int((coverage_df["string_edge_present_score700"] ==False).sum())n_mappable_cov = cov_present + cov_absentdrift_in_both =int(drift_df["agreement"].sum())n_drift_union =len(drift_df)cov_labels = [f"STRING recovers\nRaghunath edges\n({cov_present}/{n_mappable_cov} mappable)",f"our v12.0 pull vs\nD'Arcy's frozen STRING\n({drift_in_both}/{n_drift_union} agree)"]cov_values = [100* cov_present / n_mappable_cov, 100* drift_in_both / n_drift_union]bars2 = ax.bar(cov_labels, cov_values, color=["#55A868", "#C44E52"])for b, v inzip(bars2, cov_values): ax.text(b.get_x() + b.get_width() /2, v +2, f"{v:.0f}%", ha="center", fontsize=10, fontweight="bold")ax.axhline(100, color="grey", linestyle=":", linewidth=0.8)ax.set_ylabel("percent agreement")ax.set_title("The network you choose changes the answer:\nSTRING coverage and cross-pull drift")ax.set_ylim(0, 115)fig.suptitle("NB5 — candidate pigmentation network comparison (pre-harmonization)", fontsize=11)fig.tight_layout(rect=[0, 0, 1, 0.94])fig.savefig(FIGDIR /"nb5_candidate_network_comparison.png", dpi=150)plt.close(fig)print("wrote", FIGDIR /"nb5_candidate_network_comparison.png")print(f"Panel B headline: STRING recovers {cov_present}/{n_mappable_cov} = {100*cov_present/n_mappable_cov:.1f}% "f"of Raghunath's mechanistic edges; our fresh v12.0 pull agrees with D'Arcy's frozen STRING snapshot on "f"{drift_in_both}/{n_drift_union} = {100*drift_in_both/n_drift_union:.1f}% of edges "f"(~{100-100*drift_in_both/n_drift_union:.0f}% drift) -- the network you choose changes the answer.")
wrote notebooks/figures/nb5_candidate_network_comparison.png
Panel A title confirms 5-way (was stale '4-way' label; Bajpai is the 5th source): ['Raghunath\nonly', "D'Arcy-S1\nonly", 'KEGG\nonly', 'Reactome\n(curated)\nonly', 'Bajpai\n2023\nonly', 'All 5\n(n=1)']
Panel B headline: STRING recovers 158/265 = 59.6% of Raghunath's mechanistic edges; our fresh v12.0 pull agrees with D'Arcy's frozen STRING snapshot on 491/740 = 66.4% of edges (~34% drift) -- the network you choose changes the answer.
Discussion
The four sources are genuinely non-redundant peers, not near-duplicates. The full intersection across all four (Raghunath-168, D’Arcy-S1-243, KEGG-hsa04916-101, and the Step 3b Reactome curated union) is 7 genes (EDNRB, KIT, MC1R, MITF, POMC, TYR, TYRP1) — the classic melanogenesis/OCA loci. Each source contributes a large exclusive share: 122 Raghunath-only genes (signaling context absent from disease/pathway curation), 211 D’Arcy-S1-only genes (disease associations absent from the mechanistic literature Raghunath synthesized), 64 KEGG-only genes (broader second-messenger/Wnt/MAPK context this notebook’s Raghunath set does not carry as gene-level nodes), and 108 Reactome-only genes (see the discovery-layer discussion below). This is evidence for eventually harmonizing across layers (each one alone under-covers the biology) and evidence against collapsing the comparison step — a reader needs to see which layer contributed which gene before trusting a merged network’s provenance.
Headline result: the network you choose changes the answer. Applying STRING v12.0 identically (required_score=700) across every gene set shows it recovers only 59.6% of Raghunath’s own signed/directed mechanistic edges (158/265 STRING-mappable pairs) as an undirected association — a moderate coverage rate, not the near-complete overlap one might assume for two resources describing the same melanogenesis biology. And comparing our fresh v12.0 pull against D’Arcy’s own frozen STRING snapshot (Table S4), on the identical 243-gene node scope, shows only ~66% edge-set agreement (Jaccard) — ~34% of the union of edges disagrees even at the same nominal confidence threshold and the same underlying database. Both numbers are restated in full in Step 8/Step 9 and in the figure below; they are foregrounded here because they are the notebook’s single most consequential result for anything built on top of a STRING pull: the specific STRING snapshot and gene-list scope used is not an incidental implementation detail — different STRING pulls of the “same” biology disagree on a third of their edges, so a downstream network’s structure (and any convergence grade computed from it) depends on which STRING pull, at which version and threshold, was used to build it.
D’Arcy-vs-Raghunath non-overlap, revisited. The pre-registered-sounding expectation — that the non-overlap might be enriched for hypopigmentation-class disease genes — is not supported by this analysis: hypopigmentation-class rate is actually slightly lower among the 211 D’Arcy-only genes (~50.7%) than among the 16 genes both curation efforts independently agree on (75.0%), though neither difference reaches significance at these sample sizes (Fisher’s exact, D’Arcy-only n=211). The genes both methods converge on are disproportionately the classic OCA/albinism loci (OCA2, TYR, TYRP1, MC1R, KIT, KITLG, SOX10, PAX3, MITF) — the loci with the cleanest, most textbook causal phenotype, which is plausibly why both a mechanism-first synthesis and an OMIM-mining effort each independently surfaced them. The non-overlap’s real content is coverage breadth (D’Arcy’s OMIM mining reaches many disease genes a signaling-pathway synthesis was never going to include), not a phenotype-class-specific gap. Reporting the finding as measured, rather than as the more publication-friendly a-priori direction, is the result. (Note: this D’Arcy-only count, 211, is smaller than the 223 reported when Reactome was still the 5-gene precision set — 23 D’Arcy-S1 genes are now also members of the broadened 156-gene Reactome curated union, so they move out of “D’Arcy-only” once Reactome is defined that way. The 16-gene D’Arcy∩Raghunath overlap set itself is unaffected.)
Reactome as a discovery / rescue-connection layer, not a redundant fourth mechanistic peer. Per PI direction, the Step 3b curated Reactome union (156 genes) is deliberately kept broader and noisier than KEGG — KEGG hsa04916 (101 genes) remains the tighter, more thorough mechanistic core. The Reactome-only connection set — genes in the curated Reactome union that are in neither Raghunath-168 nor KEGG-hsa04916 — is 121 genes (Step 5b). 118 of these come from the MITF-M-regulated-melanocyte-development component, distributed across MITF’s sub-pathway taxonomy as: 43 from MITF-M expression/activity regulation, 29 from the pigmentation-target sub-pathway itself, 16 from apoptosis targets, 11 from lysosome biogenesis/autophagy targets, 10 from ECM/EMT targets, 8 from cell-cycle/proliferation targets, 5 from DNA-repair/senescence targets, 3 from invasion targets, and 2 from metabolism targets (a gene can appear in more than one bucket). The remaining connection-set genes come from the Melanocortin-receptor DefinedSet (MC2R, MC3R, MC4R, MC5R — MC1R itself is already in Raghunath-168) and the OCA6 disease pathway (SLC24A5). 13 of the 121 connection-set genes (CDKN2A, EDN3, GPR143, IRF4, MC2R, MLPH, MYO5A, RAB27A, SLC24A5, SLC45A2, SNAI2, TERT, TNFSF11) are also independently flagged as OMIM pigmentation-disease genes in D’Arcy’s Table S1 — these are the strongest candidate rescue connections, since two independent evidence types (a specific curated Reactome reaction/sub-pathway, and an independent OMIM disease-gene curation) already agree on them, even though neither Raghunath’s mechanism-first synthesis nor KEGG’s melanogenesis diagram contains them.
This connection set is the intended discovery surface for NB7, not a finding on its own. The framing motivating this layer: an effector-uncertain GWAS or paper locus whose resolved causal gene lands in this 121-gene Reactome-only connection set is a candidate rescue — it connects a “missing” association gene into curated pigmentation reaction-space via a specific reaction or sub-pathway the original authors did not name, rather than remaining an isolated statistical hit. This notebook does not run that per-locus test; it only builds and characterizes the substrate. To make that test citable at the reaction level, the sub-pathway/reaction structure (not just a flat gene list) is preserved in data/external/db_responses/reactome_mitf_subpathway_structure.json — so NB7 can say a rescued locus connects specifically via, e.g., “MITF-M-dependent RAB27A expression” (R-HSA-9856460) rather than a bare “Reactome” citation.
STRING enzyme-class-token caveat. Raghunath’s 168-gene set includes 6 tokens that are enzyme-class labels rather than single-gene symbols (PLA2, MMPs, Trypsin, Phosphodiesterase, PKC, PLC), each resolved separately to its HGNC gene-group membership elsewhere in this project (see hgnc_gene_groups.json, Step 4’s Rule-B expansion). When these class-label strings are sent unmodified to STRING’s symbol-mapping API (Step 7), two — PLA2 and PRSS1/TRYPSIN’s literal-token variant — correctly fail to map (unmapped), but four resolve to a STRING best-match single gene that is not the intended class at all: PLC→HSPG2, PKC→PRRT2, TRYPSIN→PRSS1, PHOSPHODIESTERASE→SMPDL3A. This is an artifact of STRING’s fuzzy symbol resolution treating an ambiguous class-label string as if it were a gene symbol, not a real biological mapping. This does not change any node-level set-comparison count reported in this notebook — Step 5’s four-way comparison uses literal-token identity (Rule A), never the STRING resolution, for set membership. It does mean: any STRING edge in the Step 7 frozen pull that touches one of these six enzyme-class tokens must not be attributed to the specific STRING-resolved gene (e.g. an edge touching PLC must not be read as an edge specifically involving HSPG2) until the class label is properly expanded to its full member-gene set. That expansion is deferred to NB6, which is where the Rule-B HGNC-group expansion (Step 4) is meant to be applied to the network itself, not just to the node-level counts computed here.
STRING as a controlled layer: real coverage, real drift. Applying STRING identically (v12.0, required_score=700) to every gene set shows it recovers roughly 60% of Raghunath’s own signed/directed mechanistic edges as an undirected association — a moderate, informative coverage rate, not the near-100% one might assume for two resources describing overlapping biology. Comparing our fresh pull against D’Arcy’s own frozen STRING snapshot (Table S4) on the identical 243-gene node scope shows ~66% edge-set agreement — drift attributable to STRING version differences between D’Arcy’s 2023 pull and the v12.0 pull used here. Both numbers argue for treating STRING as a consistent evidentiary layer to be combined with others, not as ground truth to validate everything else against.
Bajpai’s CRISPR screen as a node layer: a fifth, orthogonal, experimental line of evidence. Of the 169 genes whose perturbation reduces melanin in Bajpai et al. 2023’s unbiased genome-wide screen, only 27 (16.0%) overlap the 572-gene union of the four curated/pathway networks compared above — 142 Bajpai hits (84.0%) are orphan nodes, absent from Raghunath, D’Arcy-S1, KEGG, and Reactome alike. These orphans are not noise: the screen is causal-perturbation evidence (a gene must actually change melanin content when knocked out to be a hit), independent of any prior literature curation bias, so an orphan hit is a candidate discovery node — a gene the curated networks have not yet connected into pigmentation biology at all. Reverse coverage is uneven across the four networks: D’Arcy’s disease-gene table captures the most hits (23/243 network genes, a highly significant Fisher’s-exact enrichment, OR≈17, p≈6e-18, against the screen’s own assayed background), while Raghunath’s tighter mechanistic set captures the fewest (5/168, OR≈3.5, p≈0.02) — plausible, since Raghunath’s signaling synthesis targets well-characterized regulators rather than the broader trafficking/vesicular machinery (WASH complex, AP-1/BLOC subunits, and similar) that dominates the orphan list. Two caveats carry forward: the 169-gene hit list is a q<0.10 threshold choice (149 at q<0.05, 208 at q<0.15), so exact counts — though not the qualitative enrichment pattern — are threshold-sensitive; and the screen’s asymmetry means a gene’s absence from the 169 hits is not evidence it plays no role in pigmentation (it may be non-significant, redundant with a paralog, or simply outside the ~4,956-gene assayed universe). By design, no gene-gene edge is created from this screen anywhere in this notebook — Bajpai enters only as node-set membership and a node weight/sign (nb5_gene_set_membership.csv), with an optional, explicitly labeled bipartite hit→melanin_content layer (nb5_bajpai_bipartite_melanin_endpoint.csv) kept separate from every topology statistic above.
The networks-typology table (nb5_networks_typology.csv) makes the category differences explicit. Of the eight sources this project has now assembled, only Raghunath’s mechanistic network is directed and signed; KEGG and Reactome are directed pathway references without edge-level sign; both STRING layers (D’Arcy’s frozen S4/S5 and our fresh v12.0 pull) are undirected, unsigned association networks; and Bajpai, Baxter, and the GWAS/curated-loci table are gene-list-only sources with zero intrinsic gene-gene edges. Collapsing these into a single “network” without carrying this table forward would silently launder three different evidence types (mechanism, pathway-membership, and node-only phenotype/disease association) into one edge list — exactly the harmonization risk Step 5’s non-redundant-peers framing and Step 5b’s Reactome- as-discovery-layer framing were already designed to avoid for the original four sources.
What this notebook deliberately does not do. It does not harmonize the four sources into one network, does not assign a convergence/confidence grade to any gene, and does not run the per-locus rescue test that motivates the Reactome-only connection set — that is explicitly left to later notebooks (harmonization, and NB7’s rescue test), consistent with this project’s framing that the comparison is worth doing before, and separately from, any merge.
Gaps
R-HSA-5619507 is not “Pigmentation.” An earlier revision of this notebook assumed Reactome had a parent pathway named “Pigmentation” at accession R-HSA-5619507 with Melanin biosynthesis, MC1R signaling, and melanosome trafficking as sibling children, and planned to pull that parent’s descendant gene set. Querying the Reactome ContentService directly shows R-HSA-5619507 is “Activation of HOX genes during differentiation” — an unrelated developmental pathway. Reactome has no pathway literally named “Pigmentation,” and Melanin biosynthesis (under Metabolism) and the MITF-regulation pathways (under Developmental Biology) sit in disjoint hierarchy branches with no shared low-level parent — confirmed via the ContentService ancestors endpoint for both branches. Per PI direction, the broader Reactome layer used in this notebook (Step 3b) is instead an explicit, itemized, manually curated union of four independently verified accessions (Melanin biosynthesis, MITF-M-regulated melanocyte development, Melanocortin receptors, and the SLC24A5/OCA6 disease pathway), not a single pathway’s descendant tree.
The 156-gene Reactome curated union is dominated by one large, heterogeneous component.R-HSA-9730414 (MITF-M-regulated melanocyte development) alone contributes 152 of the 156 genes, and many of them are MITF targets in cell-cycle control, apoptosis, EMT, and lysosome biogenesis whose relevance to pigmentation specifically — as opposed to melanocyte biology broadly — is indirect. These genes are included in full, not hand-filtered to a pigmentation-only subset, because that filtering would itself be an unstated judgment call; the discovery-layer framing (Step 5b) is precisely designed to tolerate this noise, since the point of a rescue-connection substrate is breadth, not precision.
The 465/230/118 reproduction (Step 4) uses the S1∪S5 union framing already committed in DATA_SOURCES.md, not the S1-only 243-gene set used for the primary node-level set comparison in Step 5 — both framings are reported so the two are not conflated (243-gene S1-only comparison: 227 absent, 118 hypo; 508-gene S1∪S5 union comparison: 465 absent, 230 disease-flagged, 118 hypo — the same hypo count, because S5’s own additional disease-flagged genes happen not to shift the hypopigmentation-specific count in this dataset).
STRING alias resolution is not manually re-checked per gene beyond the enzyme-class-token caveat documented above in the Discussion. A parallel, smaller alias-resolution drift exists in the D’Arcy-S1-243 STRING pull (IFTAP→C11orf74, VEGFA→COL18A1, TERC→GAR1) — not further characterized here.
No population/ancestry stratification. All four sources are population-agnostic gene lists; this comparison says nothing about whether the four networks differ in which populations’ pigmentation variation they best explain.
Yang 2016 (OCA2, East Asian) is not relevant to this notebook — it is a population GWAS case study for a later discordance/rescue notebook, not a candidate gene-set/network source compared here.
The Reactome-only connection set is a hypothesis-generating substrate, not a validated rescue list. Being in the connection set means a gene is in curated Reactome pigmentation-adjacent reaction-space and outside the current mechanistic core — it does not by itself mean the gene explains any specific effector-uncertain locus. That per-locus test is NB7’s job, not this notebook’s.