54766047e2e79d545162f894d20d631553175d3d mspeir Tue Aug 4 15:05:35 2026 -0700 singleCellSignalsPeaks: cite the SEA-AD resource paper, fail loudly on a hub-format change The Hawrylycz 2024 reference was listed but never cited. It is the SEA-AD resource paper and the source for the two brain regions the SEA-AD tracks cover, so it now appears on that dataset's bullet alongside Gabitto. BrainVar and the Allen basal-ganglia set have no paper to cite, so those bullets stay as they are. makeSingleCellSignalsPeaksRa.py matched the parent line by exact string, so if the hub ever emitted "parent <composite> off" or changed its spacing, every stanza would be skipped and the script would write a header-only .ra and exit 0 -- the next trackDb load would then quietly drop all 925/587 subtracks. It now matches on the parent's first token and refuses to write when nothing matched, or when the subtrack count disagrees with the facet metadata that build_stanzas emits alongside it. Verified: a renamed composite and a truncated stanza file both exit non-zero without writing a file, a "parent <composite> off" line still yields the full 587, and the regenerated .ra files are byte-identical to the previous commit. refs #37914 Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> diff --git src/hg/makeDb/scripts/singleCellSignalsPeaks/makeSingleCellSignalsPeaksRa.py src/hg/makeDb/scripts/singleCellSignalsPeaks/makeSingleCellSignalsPeaksRa.py index 9eefca0a207..74a09b868e6 100644 --- src/hg/makeDb/scripts/singleCellSignalsPeaks/makeSingleCellSignalsPeaksRa.py +++ src/hg/makeDb/scripts/singleCellSignalsPeaks/makeSingleCellSignalsPeaksRa.py @@ -1,132 +1,163 @@ #!/usr/bin/env python3 """ Regenerate the native singleCellSignalsPeaks trackDb .ra for a genome assembly from the Cell Browser all-tracks-hub build (Redmine #37820 for hg38, #37914 for mm10). The native "singleCellSignalsPeaks" faceted composite is the Genome Browser version of the hub's main signal-&-peaks composite for that assembly (cellBrowserHg38 for hg38, cellBrowserMm10 for mm10). This script takes that composite's stanzas from the hub build and rewrites them into a native track: - the composite is renamed cellBrowser<Asm> -> singleCellSignalsPeaks - each subtrack bigDataUrl is repointed to the local /gbdb copy - subtrack colors / labels / types are carried through unchanged (so the harmonized cell-type labels and any per-track colors come along for free) The data files themselves are copied into /hive/data/genomes/<asm>/bed/singleCellSignalsPeaks/<served-relpath> (see copySingleCellSignalsPeaksFiles.py) and served via the /gbdb/<asm>/bbi/singleCellSignalsPeaks symlink; this script only (re)writes the trackDb .ra. Usage: makeSingleCellSignalsPeaksRa.py [--assembly hg38|mm10] [--stanzas STANZAS] [--out OUT] """ -import re, os, argparse +import re, os, sys, argparse from urllib.parse import urlparse # Where the hub build (build_manifest.py / build_stanzas.py) writes its stanzas and # metadata. That machinery builds the whole Cell Browser super hub, not just this track, # so it lives outside the kent tree; override with HUB_BUILD when it moves. HUB_BUILD = os.environ.get( "HUB_BUILD", "/hive/users/mspeir/claude/cell-browser/all-tracks-hub-build") TRACK = "singleCellSignalsPeaks" GROUP = "regulation" # ATAC-seq signal/peaks live with the ENCODE # regulatory tracks, not under singleCell ORG = {"hg38": "human", "mm10": "mouse"} # trackDb organism subdir per assembly def main(): ap = argparse.ArgumentParser() ap.add_argument("--assembly", default="hg38", choices=sorted(ORG)) ap.add_argument("--stanzas") ap.add_argument("--out") args = ap.parse_args() asm = args.assembly hub_composite = "cellBrowser" + asm.capitalize() # cellBrowserHg38 / cellBrowserMm10 gbdb = "/gbdb/%s/bbi/%s" % (asm, TRACK) stanzas = args.stanzas or os.path.join(HUB_BUILD, "stanzas/%s.trackDb.txt" % asm) out = args.out or os.path.join( os.path.dirname(os.path.abspath(__file__)), "../../trackDb/%s/%s/%s.ra" % (ORG[asm], asm, TRACK)) header = "\n".join([ "track " + TRACK, "compositeTrack faceted", "group " + GROUP, "visibility hide", "type bigBed 3", "shortLabel Single-cell ATAC-seq", "longLabel Single-cell ATAC-seq Peaks and Signals for UCSC Cell Browser datasets", "metaDataUrl %s/%s_metadata.tsv" % (gbdb, TRACK), "primaryKey Track", "subtrackUrls Dataset=https://cells.ucsc.edu/?ds=$$", "defaultSortField Cell_class", "maxCheckboxes 200", ]) # class ordering for subtrack priority: palette line order (neurons, glia, # vascular, immune, ...) so same-class tracks group together in the display, # with the source (hub/dataset) order preserved within a class. The subtrack's # broad class is recovered from its color (palette is 1:1 class<->color). color_rank = {} # prefer the palette archived alongside this script (the copy of record, written by # build_celltype_crosswalks.py); fall back to the hub build dir _palf = os.path.join(os.path.dirname(os.path.abspath(__file__)), "celltype-crosswalks", "celltype-palette.tsv") if not os.path.isfile(_palf): _palf = os.path.join(HUB_BUILD, "celltype-crosswalks", "celltype-palette.tsv") for _i, _l in enumerate(open(_palf)): _pp = _l.rstrip("\n").split("\t") if len(_pp) >= 2: color_rank[_pp[1]] = _i class_seq = {} # rank -> running counter within that class # a source path segment of "old" / "*.old" / "*_old" marks deprecated data # (e.g. cortex-atac/hub/interact.old/); the hub may keep it, but the native # track must not carry it. OLD_SEG = re.compile(r"(^|/)[^/]*(\.old|_old|\bold)($|/)", re.I) out_stanzas = [header] - n = skipped_old = 0 + n = skipped_old = n_stanzas = n_parented = 0 for s in re.split(r"\n\s*\n", open(stanzas).read().strip()): lines = s.splitlines() - parent = next((l for l in lines if l.startswith("parent ")), "") - if parent.strip() != "parent " + hub_composite: + n_stanzas += 1 + # Match on the parent's first token rather than the whole line. The exact-string + # compare this replaces would have silently skipped every stanza if the hub ever + # emitted "parent <composite> off" or changed its spacing, leaving a header-only + # .ra and a zero exit status. + parent = next((l for l in lines if l.strip().startswith("parent ")), "") + ptoks = parent.split() + if len(ptoks) >= 2: + n_parented += 1 + if len(ptoks) < 2 or ptoks[1] != hub_composite: continue bdu = next((l for l in lines if l.strip().startswith("bigDataUrl ")), "") rel_check = urlparse(bdu.split(None, 1)[1].strip()).path if bdu else "" if OLD_SEG.search(rel_check): skipped_old += 1 continue n += 1 # priority groups tracks by broad class (via color), source order within color = next((l.split(None, 1)[1].strip() for l in lines if l.strip().startswith("color ")), "") rank = color_rank.get(color, len(color_rank)) # unknown/uncolored last seq = class_seq.get(rank, 0); class_seq[rank] = seq + 1 priority = rank * 100000 + seq newl = [] for l in lines: if l.startswith("track "): suffix = l.split(None, 1)[1][len(hub_composite) + 1:] newl.append("track %s_%s" % (TRACK, suffix)) - elif l.strip() == "parent " + hub_composite: + elif l.strip().startswith("parent ") and l.split()[1] == hub_composite: # "off" so every subtrack is unchecked by default; the user turns # on individual tracks via the faceted selector newl.append("parent " + TRACK + " off") newl.append("priority " + str(priority)) elif l.strip().startswith("bigDataUrl "): rel = urlparse(l.split(None, 1)[1].strip()).path.lstrip("/") newl.append("bigDataUrl %s/%s" % (gbdb, rel)) else: newl.append(l) out_stanzas.append("\n".join(newl)) + # Sanity checks: fail loudly rather than write a truncated .ra. A hub-format change + # that stops the parent line matching would otherwise produce a header-only file and + # exit 0, and the next trackDb load would quietly drop every subtrack. + if n == 0: + sys.exit("ERROR: no subtracks matched composite '%s' in %s " + "(%d stanzas, %d with a parent line). Has the hub stanza format " + "changed?" % (hub_composite, stanzas, n_stanzas, n_parented)) + # The facet metadata is the parallel artifact: build_stanzas writes one row per + # subtrack of this composite, so the counts must agree once the old-dir skips are + # added back. A mismatch means the .ra and the metadata disagree, which shows up in + # the browser as subtracks with no facet row (or facet rows with no track). + meta = os.path.join(HUB_BUILD, "meta", "%s.metadata.tsv" % asm) + if os.path.isfile(meta): + with open(meta) as fh: + meta_rows = sum(1 for _ in fh) - 1 # minus the header + if meta_rows != n + skipped_old: + sys.exit("ERROR: %s has %d rows but %d subtracks were kept (+%d old-dir " + "skipped); the .ra and the facet metadata must match 1:1" + % (meta, meta_rows, n, skipped_old)) + else: + sys.stderr.write("WARNING: no facet metadata at %s, skipping the 1:1 check\n" + % meta) + with open(os.path.abspath(out), "w") as fh: fh.write("\n\n".join(out_stanzas) + "\n") print("wrote %s: %d subtracks (assembly=%s, composite=%s, group=%s; skipped %d old-dir)" % ( out, n, asm, hub_composite, GROUP, skipped_old)) if __name__ == "__main__": main()