156d289f517d82c4d5a8c59981dfec21a715af8d lrnassar Mon Sep 21 09:22:45 2026 -0700 Fiber-seq QA fixes: take sample class from the lab's own table, pin the overlay draw order, and correct the description pages. refs #36210 Sample class had been guessed from the free-text cell type with a two-name exception list, which filed five lymphoblastoid lines as HPRC that are not. It now comes from a fifth column in fiberSeqSamples.tsv carrying the classification the lab supplies, and sampleClass()/NOT_HPRC are gone. Facet counts go from 25/16 to 20/21. The CpG difference track is a solidOverlay whose four files hold the same value at a shared base, so the tier drawn last is the colour the reader sees. Its children inherited a single priority from the parent, trackPriCmp ties on that, and slSort is not stable, so the paint order was arbitrary and p<0.01 was covering p<0.0001. The four levels and the two haplotype children now carry explicit priorities. At chr20:29,300,000-29,305,000 red goes from 2 image columns to 88. CpG haplotype children take their shortLabel prefix from their container, so they read "<sample> CpG Hap1" rather than colliding with the accessibility children's "<sample> Hap1". 82 labels were duplicated. Description pages: fiberSeqAcc.html attached the FIRE score's "fewer than four elements" cutoff to the percent accessible signal, which is a different quantity; two pages overstated how many of the lymphoblastoid lines come from HPRC; both pages told readers to query the API with container track names, which it refuses by design. Also documents GM12878's trio phasing and the FDR ceiling at 100, corrects the difference track's stated range, replaces a non-ASCII author name with numeric entities, and adds db= to hgTrackUi links. makeDoc: corrects a basesCovered figure that was out by a factor of a hundred, refreshes the peak file sizes after the PM00001 reissue, rewrites the sampleClass rationale, and records why multiWig children need their own priority. diff --git src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py index 978de73183c..90d715505e6 100755 --- src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py +++ src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py @@ -65,64 +65,50 @@ ] # Shown behind an info icon on the Sample class column heading. SAMPLE_CLASS_DESCRIPTION = ( "HPRC = Lymphoblastoid (B-lymphocyte, EBV) cell lines from the NHGRI " "Human Pangenome Reference Consortium") # Sample class swatches, shown next to that facet's checkboxes. # Okabe-Ito colors for the swatches. SAMPLE_CLASS_COLORS = { "HPRC": "#0072B2", "Common Cell Line": "#D55E00", } -# Lymphoblastoid lines that are not from the Human Pangenome Reference -# Consortium: GM12878 is the ENCODE line and HG002 is Genome in a Bottle. Both -# are B-lymphocyte EBV lines like the HPRC samples, so cell type alone cannot -# tell them apart and they have to be named. -NOT_HPRC = {"GM12878", "HG002"} - - -def sampleClass(sample, cellType): - """Split a sample two ways for the Sample class column. - - HPRC is a lymphoblastoid (B-lymphocyte, EBV) line from the consortium; - everything else, including the two lymphoblastoid lines listed in NOT_HPRC, - is a common cell line.""" - if "lymphoblastoid" in cellType.lower() and sample not in NOT_HPRC: - return "HPRC" - return "Common Cell Line" - - def readSamples(path): - """Read fiberSeqSamples.tsv into a list of dicts, in file order.""" + """Read fiberSeqSamples.tsv into a list of dicts, in file order. + + sampleClass comes from the lab's own sample sheet, not from the cell type: + five of the lymphoblastoid lines are common cell lines rather than HPRC + samples, so there is nothing in the cell type that tells the two apart.""" samples = [] with open(path) as f: for line in f: if line.startswith("#") or not line.strip(): continue fields = line.rstrip("\n").split("\t") - if len(fields) < 4: - sys.exit("bad sample line, want 4 fields: %s" % line.rstrip()) - acc, sample, cellType, _hash = fields[:4] + if len(fields) < 5: + sys.exit("bad sample line, want 5 fields: %s" % line.rstrip()) + acc, sample, cellType, _hash, sampleClass = fields[:5] samples.append({ "accession": acc, "sample": sample, "cellType": cellType, - "sampleClass": sampleClass(sample, cellType), + "sampleClass": sampleClass, }) if not samples: sys.exit("no samples read from %s" % path) return samples def writeMetadata(path, samples): """The sample table shown on the track UI page. A plain column name gets facet checkboxes, a leading underscore means searchable and sortable but not faceted. Accession is the primaryKey but sits last, since it is the least interesting thing about a sample. Nothing requires the primaryKey to come first: facetedComposite.js only checks that the column exists, and every use of it is by name. It is still the default sort, because accession order keeps the @@ -373,73 +359,90 @@ "aggregate solidOverlay", "showSubtrackColorOnUi on", "type bigWig -100 100", "viewLimits -100:100", "autoScale off", "windowingFunction mean", "maxHeightPixels 100:50:8", "shortLabel %s CpG diff" % name, "longLabel %s CpG methylation difference between haplotypes, " "by significance threshold" % name, "onlyVisibility full", pri(), ]) # Least significant first, so the more significant levels draw on top. # Not "i": that is the sample index pri() builds its priority from. + # + # The priority is what makes that happen and is not decoration. This is + # a solid overlay and all four files hold the same value at a shared + # base, so the level that draws last is the colour the user sees. Order + # of declaration does not survive: makeContainerTrack() in + # hg/hgTracks/container.c sorts the children with trackPriCmp, which + # compares priority alone, and slSort is not a stable sort, so children + # left on the parent's inherited priority draw in an arbitrary order. for level, (fname, label, color) in enumerate(DIFF_LEVELS): out += stanza(12, [ "track fiberSeqCompendium_%s_cpgDiff_l%d" % (acc, level), "parent fiberSeqCompendium_%s_cpgDiff" % acc, "type bigWig", "bigDataUrl %s/%s/%s" % (gbdb, acc, fname), "color %s" % color, "shortLabel %s %s" % (name, label), "longLabel %s CpG haplotype difference, %s" % (name, label), + "priority %d" % (level + 1), ]) return out def hapOverlay(gbdb, acc, name, dataType, file1, file2, shortLabel, longLabel, childLongLabel, windowing, priority): """A haplotype 1 / haplotype 2 transparent overlay, used for both the accessibility and the CpG haplotype data types.""" out = stanza(8, [ "track fiberSeqCompendium_%s_%s" % (acc, dataType), "parent fiberSeqCompendium off", "container multiWig", "aggregate transparentOverlay", "showSubtrackColorOnUi on", "type bigWig 0 100", "viewLimits 0:100", "autoScale off", "alwaysZero on", "windowingFunction %s" % windowing, "maxHeightPixels 100:40:8", "shortLabel %s" % shortLabel, "longLabel %s" % longLabel, "onlyVisibility full", priority, ]) + # Haplotype 1 then haplotype 2. This overlay is transparent, so the order + # barely shows, but a child without its own priority inherits the parent's + # and then draws in whatever order slSort happens to leave it in - see the + # longer note on the cpgDiff children, where the same thing is visible. for hap, fname, color in (("h1", file1, HAP1_COLOR), ("h2", file2, HAP2_COLOR)): n = hap[1] out += stanza(12, [ "track fiberSeqCompendium_%s_%s_%s" % (acc, dataType, hap), "parent fiberSeqCompendium_%s_%s" % (acc, dataType), "type bigWig", "bigDataUrl %s/%s/%s" % (gbdb, acc, fname), "color %s" % color, - "shortLabel %s Hap%s" % (name, n), + "priority %s" % n, + # Take the prefix from the container's own shortLabel, so the CpG + # children come out "<sample> CpG Hap1" rather than colliding with + # the accessibility children's "<sample> Hap1". + "shortLabel %s%s" % (shortLabel.removesuffix("1/2"), n), "longLabel %s %s" % (childLongLabel, n), ]) return out def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--data-dir", default="/hive/data/genomes/hg38/bed/fiberSeq", help="where the metadata and color files are written") ap.add_argument("--gbdb-dir", default="/gbdb/hg38/fiberSeq", help="bigDataUrl prefix, i.e. the symlink directory") ap.add_argument("--trackdb-dir", default=os.path.expanduser( "~/kent/src/hg/makeDb/trackDb/human/hg38"),