752cb1161a3dc006932fd3ba179c9b2fbe449c8d jnavarr5 Fri Sep 25 14:27:43 2026 -0700 Updating the episignatures makedoc for the semicolon filter labels and adding a recount command for the opposite-direction sites count, refs #38112 diff --git src/hg/makeDb/doc/hg38/episignatures.txt src/hg/makeDb/doc/hg38/episignatures.txt index 826e2a9c354..550d4a4aca8 100644 --- src/hg/makeDb/doc/hg38/episignatures.txt +++ src/hg/makeDb/doc/hg38/episignatures.txt @@ -1,306 +1,319 @@ # 2026-09-16 Claude max: MethaDory episignature CpG probes, hg38 # The episignatures superTrack collects the CpG positions that published DNA # methylation signatures of rare developmental disorders are built from. Its # first subtrack is MethaDory. ############################################################################## # MethaDory episignature loci (DONE 2026-09-17) ############################################################################## # Federico Ferraro and Dmitrijs Rots (Erasmus MC Rotterdam) sent two files by # email, for open release: # 20260916_episignatures_loci.tsv probe, gene/locus, deltaBeta, # study, pval, padj, label # 20260916_episignatures_loci_metadata.xlsx study -> PMID, signature, # disorder, notes # There is no public download URL for these; MethaDory itself, with its trained # classifiers, is at https://github.com/f-ferraro/MethaDory and the preprint is # doi:10.1101/2025.03.28.25324859. # # A first version dated 20260914 was built and then superseded; it is kept in # prev-20260914/ for reference. Two things changed in the second version, both # of which the build has to cope with: # - the three provisional CHD3 signatures (Santini20252A/2B/2C, alternative # versions still being tested for the original authors) were removed, and # the missing metadata for SPEN Radio2021 and the missing probes for # MORC2 Peymani2026 were added # - the 'signature' column of the metadata sheet, which is the key that joins # to the Label column of the loci table, switched from "CHD3_Santini2025" to # "CHD3 Santini2025". methaDoryToBed.py normalises spaces to underscores on # both sides so either spelling joins. # The loci tsv is also CRLF. Python's text mode strips the CR, but a shell # pipeline over it does not: counting labels with "zcat | cut -f7 | sort -u" # gives 190 rather than 189, because one line ends LF and the rest CRLF. mkdir -p /hive/data/genomes/hg38/bed/episignatures/methaDory cd /hive/data/genomes/hg38/bed/episignatures/methaDory # ... the two files above were copied in here, the tsv gzipped ... # The probe table gives Illumina probe IDs, not coordinates. Positions come from # the Illumina manifests already served as tracks on the browser, preferring # EPIC v1 850K, then EPIC v2, then 450K. probeCoords.sh does that and is called # by the build script below; makeMethaDory.sh runs the whole thing: bash ~/kent/src/hg/makeDb/scripts/episignatures/makeMethaDory.sh # Numbers it reports, for the record: # input rows: 847863 # rows with no probe position: 373 (293 distinct probes) # bed features written: 268900 # distinct disorders: 105 # distinct genes/loci: 113 # distinct studies: 74 # # The 373 dropped rows are probes that appear in none of the three manifests. # They are listed in methaDory.bed.unmapped. The remaining 847,490 probe-by- # signature rows merge into 268,900 CpG sites, since a site is usually reported # by more than one signature (up to 41, at cg05654765). 267,974 of the sites are # positioned from the EPIC 850k manifest and 926 from EPIC v2; 450K adds nothing # that the other two do not already have. # # Every signature label in the loci table now has a metadata row. The build # still warns about 5 metadata rows that name one study in the 'signature' # column and a different author in StudyID. Three are spelling variants, but two # are genuine crossings: the CREBBP-EP300 and NSD2 pairs swap PMID and disorder # between Levy2022 and ArefEshghi2020 / Kawai2024 respectively. Reported to the # authors, not yet resolved; the rows are used as the sheet has them, the script # does not guess. "Smith-Magenis Syndrome" and "Smith-Magenis syndrome" are also # still two separate disorder names and so two separate filter entries. # Sanity check that the coordinates are on the right base: every cg probe should # land on a CG dinucleotide. awk -F'\t' 'BEGIN{OFS="\t"} $4 ~ /^cg/ {print $1,$2,$3,$4}' methaDory.bed \ | shuf -n 20000 > /tmp/samp.bed twoBitToFa /gbdb/hg38/hg38.2bit stdout -bed=/tmp/samp.bed | grep -v '^>' \ | tr 'a-z' 'A-Z' | sort | uniq -c | sort -rn | head # 19996 CG, and 4 singletons at known polymorphic manifest positions. # Compare the position of a probe against the existing EPIC 850k track, which # is where it came from: bigBedToBed /gbdb/hg38/bbi/illumina/epic850K.bb stdout \ | awk -F'\t' '$4=="cg21870274"{print $1,$2,$3}' awk -F'\t' '$4=="cg21870274"{print $1,$2,$3}' methaDory.bed # The filter menus in the trackDb stanza (filterValues.disorders, # filterValues.loci, filterValues.studies) are generated by the build into # methaDoryFilters.ra. When the source data is refreshed, paste those three # pairs of lines back into human/hg38/episignatures.ra rather than editing them # by hand. # The two tables on the description page are generated the same way: python3 ~/kent/src/hg/makeDb/scripts/episignatures/makeHtmlTables.py \ studies studySummary.tsv > studyTable.html python3 ~/kent/src/hg/makeDb/scripts/episignatures/makeHtmlTables.py \ loci locusSummary.tsv > locusTable.html # and pasted into human/hg38/methaDory.html between the # "" and "" markers under the # "Studies" and "Genes and loci" headings. The counts quoted in the Description # and Methods paragraphs of that page have to be updated by hand at the same # time. # # Both tables link out to PubMed with a "Lastname Year" label rather than the bare # PMID, looked up from scripts/episignatures/pubAuthorYear.tsv (also used by # epigenCentralToBed.py for the EpigenCentral locus table). On a refresh, add an # entry for every new PMID before regenerating the tables; the scripts errAbort on # an id with no entry rather than falling back to the raw number. # Colour: five classes, from the strongest effect at the site. Warm = hyper, # cool = hypo, darker = |delta-beta| >= 0.10, purple = conflicting. Three things # decided this, all measured on the data rather than assumed: # # - A linear ramp on delta-beta would be useless. The distribution is heavily # right-skewed: median 0.110, p75 0.159, p90 0.232, p99 0.416, max 0.835, # and 43% of sites sit in 0.05-0.10. Nearly everything would land in the # bottom fifth of the ramp. # # - The low end is a reporting artefact, not biology. Per study, the minimum # |delta-beta| reveals the cut-off each paper applied before publishing its # site list: 4 studies cut at 0.20 (Velasco2021 among them, whose median is # therefore 0.254 against Levy2022's 0.076), 25 at ~0.10, 33 at ~0.05, and # 12 applied none. 98.1% of rows come from a study that cut at >= 0.05. # So quantile-derived bins would partly encode "which study reported this". # A fixed, interpretable line at 0.10 is used instead; it is also the modal # reporting cut-off. Check it on a refresh with: # python3 -c "..." # per-study min/median/p90 of abs(deltaBeta) # and see the caution paragraph on methaDory.html. # # - "Conflicting" needs a size test, not a unanimity test. Requiring all # signatures at a site to agree flags 39% of sites, and that class is almost # a proxy for how many signatures share the probe: 0% of single-signature # sites are mixed, 43.8% at two, 75.0% at 3-5, 95.5% at 6-10, 99.7% at 11+. # Requiring an opposing signature to reach 80% of the top effect # (CONFLICT_FRAC in methaDoryToBed.py) leaves 10.0%, which really are sites # with no dominant direction. At 50% it would be 27.3%, at 67% 17.2%. # # Resulting class sizes, all comfortably populated: # Strong hypermethylation 69879 26.0% | Weak hypermethylation 68907 25.6% # Weak hypomethylation 36045 13.4% | Strong hypomethylation 67250 25.0% # Conflicting 26819 10.0% # # The class names must not contain a comma: they go into filterValues.direction, # and comma is its separator. "Hypermethylated, strong" was silently chopped # into two menu entries before they were renamed to "Strong hypermethylation". # Details page: the per-signature values are stored twice. Once as parallel # comma-separated columns (signatures, disorders, loci, deltaBetas, pvalues, # studies), which is what the Table Browser and bigBedToBed hand back and what # the filters run on, and once as a JSON object in a _jsonSignatures field that # "detailsDynamicTable _jsonSignatures|..." renders as a real table, so a reader # can see which delta-beta belongs to which disorder. The comma columns are in # skipFields so the details page shows only the table. # # Two things to know before copying this pattern: # - detailsDynamicTable has two encodings. A field whose name starts with # "json" or "_json" is treated as JSON and handed to hgc.js; any other field # name is treated as the ";"-and-"|" encoding, which hgc expands inside a # fixed char[4096] and errAborts with "Error substituting" past that. One # probe here carries 41 signatures (3,243 bytes of JSON), so only the JSON # encoding is usable. # - hgc.js makeGenericTable() draws one row per key of the JSON object, and # when a key's value is an array it emits one cell per element and does not # draw the key. So an object of arrays gives a plain grid. Use string keys # that are not integer-like ("r000", "r001", ...): JavaScript reorders # integer-like keys numerically, other string keys keep insertion order. # The field is built by jsonTable() in methaDoryToBed.py, ASCII-escaped, since # kent's jsonParse passes \uXXXX through untouched. # Track search: the position box has to find a probe by its cg number. The # bigBed carries -extraIndex=name, and episignatures.ra has a matching # searchTable stanza. That stanza needs termRegex and semiShortCircuit, not just # searchPriority: hgFind splits specs into a short-circuiting list and a long # list (hgPositionsFind in hg/lib/hgFind.c), only specs with a termRegex land on # the short list, and as soon as one short spec matches the long list is never # run at all. Without the termRegex the spec loaded fine and searching a probe # silently returned only the Illumina array hits. semiShortCircuit then stops # our spec from suppressing those array hits in turn. hgsql hg38 -Ne "select searchName, shortCircuit, searchPriority from hgFindSpec_max where searchTable='methaDory'" # methaDory 1 50 ############################################################################## # EpigenCentral episignature loci (DONE 2026-09-18) ############################################################################## # 2026-09-18 Claude max: refs #38112 # # The second subtrack of the episignatures container. The EpigenCentral group at the # Centre for Computational Medicine / Weksberg lab (Hospital for Sick Children, # Toronto) publishes its episignatures as a track hub: # https://github.com/ccmbioinfo/EpigenCentral-UCSC-Genome-Browser # and asked us to host it as a native track instead. The portal itself is at # https://epigen.ccm.sickkids.ca/ and is described in PMID 32623772. # # Barali Kitiyakara and Eliza Alde reviewed the hub, cloned it and worked out the # trackDb settings; that work is on the ticket and at # ~bkitiyak/public_html/epigenCentralHub/ # What is here is the same set of changes redone as a repeatable build from the # upstream file, so a refresh is one command. # # The hub's episignatures.bb already has the probe coordinates on hg38, so there is # nothing to lift or map. The build downloads it and rewrites four things: # - the OMIM column holds a full https://omim.org/entry/NNNNNN URL; only the number # is kept so trackDb's "urls displayOmim=" can build the link # - the hub renders its own mouse-over into a column ("NSD1|Loss|-0.346") and points # mouseOverField at it. That column becomes the direction alone and trackDb builds # the mouse-over from the fields, which is what the rest of our tracks do # - 215 rows of the per-site comparison table, on 137 sites, are exact duplicates of # another row for the same signature. The site's own nSignatures/signatureList # columns already count each signature once, so the duplicates are dropped # - the disorder of the strongest signature is added as its own column, so the # mouse-over can name the disorder and not only the gene # Two columns are also renamed to what the MethaDory subtrack calls the same numbers: # height -> maxAbsDelta and nSignatures -> sigCount. mkdir -p /hive/data/genomes/hg38/bed/episignatures/epigenCentral cd /hive/data/genomes/hg38/bed/episignatures/epigenCentral bash ~/kent/src/hg/makeDb/scripts/episignatures/makeEpigenCentral.sh # Numbers it reports, for the record: # upstream md5: 1f2cd0ef372fefa3b42c796fb89f6f40 # features written: 15035 (same as the upstream itemCount) # episignatures: 24, over 23 disorders # displayed direction: Gain 5972, Loss 9063 # duplicate table rows dropped: 215, on 137 probes # table rows with no delta-beta: 1 # widest comparison table: 478 bytes (cg26004771), hgc limit 4096 # # No feature is dropped: every probe in the hub is in the track. The one row with no # delta-beta is the Dup7 signature at cg19457237, which the source has as NA in both # the direction and the delta-beta column while giving it an adjusted p-value; it is # shown as NA. The script errAborts rather than guessing if the displayed signature # itself were the one missing a value. # # The widest comparison table matters because that field uses the ';' rows / '|' cells # encoding of detailsDynamicTable, which printEmbeddedTable() in hg/hgc/hgc.c expands # inside a fixed char[4096]. At 478 bytes there is a wide margin, so this track does # not need the _json encoding that methaDory has to use. The script checks the limit # on every build and stops if a refresh ever crosses it. # Sanity check on the coordinates, the same one the MethaDory build does: every cg # probe should sit on a CG dinucleotide. The features are 1 bp, so take start..start+2. awk -F'\t' 'BEGIN{OFS="\t"} {print $1,$2,$2+2,$4}' epigenCentral.bed \ | shuf -n 5000 > /tmp/samp.bed twoBitToFa /gbdb/hg38/hg38.2bit stdout -bed=/tmp/samp.bed | grep -v '^>' \ | tr 'a-z' 'A-Z' | sort | uniq -c # 5000 CG, no exceptions. # Cross-check against the MethaDory subtrack, which got its positions from the Illumina # manifests rather than from EpigenCentral. 14,236 of the 15,035 probes are in both and # every one of them is at the same position; 799 are only here. bigBedToBed /gbdb/hg38/episignatures/methaDory.bb stdout \ | awk -F'\t' '{print $4"\t"$1":"$2}' | sort -u > /tmp/md.pos awk -F'\t' '{print $4"\t"$1":"$2}' epigenCentral.bed | sort -u > /tmp/ec.pos join -t$'\t' /tmp/md.pos /tmp/ec.pos | awk -F'\t' '$2!=$3' | wc -l # 0 # The signature filter menu and the "Included episignatures" table of the description # page are generated by the build, into epigenCentralFilters.ra and # epigenCentralTable.html. On a refresh, paste the two .ra lines back into # human/hg38/episignatures.ra, and the table back into human/hg38/epigenCentral.html # between the "" and "" markers, # rather than editing either by hand. The counts quoted in the Description and Methods # paragraphs of that page have to be updated at the same time. # +# One of those counts, in Display Conventions, is how many multi-signature sites also carry +# a signature in the opposite direction to the strongest one, which the color and the +# direction filter do not show (1,618 of 2,848 in this build). Recount it with: +python3 -c " +m = x = 0 +for l in open('epigenCentral.bed'): + rows = [r.split('|') for r in l.rstrip('\\n').split('\\t')[16].split(';')][1:] + if len(rows) > 1: + m += 1 + x += len({r[3] for r in rows if r[3] != 'NA'}) > 1 +print(x, 'of', m)" +# # The publication behind each signature comes from the hub's README and is kept in # scripts/episignatures/epigenCentralRefs.tsv, since the README is not machine # readable. All 18 distinct PMIDs were checked against PubMed esummary before being # written down; EHMT1 is the one signature with no PubMed record and is cited by DOI. # The hub's own HTML page has every one of these links pointing at PMID 31311581 while # displaying a different number, which is why they are generated here instead. # A comma in a filterValues entry is the entry separator and cannot be escaped when the -# filterType is one of the *List* kinds, so the commas inside four of the disorder names -# ("Dystonia 28, childhood-onset" and friends) are dropped in the menu labels by -# writeRa() in epigenCentralToBed.py. The bigBed keeps the names as the source has them. +# filterType is one of the *List* kinds, so the commas inside three of the disorder names +# ("Dystonia 28, childhood-onset" and friends) become semicolons in the menu labels, as in +# the methaDory menu, by writeRa() in epigenCentralToBed.py. The bigBed keeps the names as +# the source has them. (Changed from dropping the commas during QA, 2026-09-25.) # Track search: the position box finds a probe by its cg number. The bigBed carries # -extraIndex=name and episignatures.ra has a matching searchTable stanza. Both this # spec and the methaDory one carry termRegex and semiShortCircuit, so a probe ID that # is in both tracks returns a hit in both, plus the Illumina array tracks. hgsql hg38 -Ne "select searchName, shortCircuit, searchPriority from hgFindSpec_max where searchTable='epigenCentral'" # epigenCentral 1 51 # Lou asked on the ticket for downloads to be off, so the stanza has "tableBrowser off" # and the Data Access section points at EpigenCentral's own portal and repository. That # switch covers the Table Browser, the Data Integrator and the REST API; the bigBed # still sits in /gbdb and is reachable on hgdownload, so it is not a hard block, and QA # should confirm with the lab that this is what they wanted. # Visibility: the track is sparse almost everywhere, 15,035 sites over the genome, but it # has one hot spot at the HOXA cluster. Measured heights at 1100px wide, pack, with no # other track on: # chr7:27,140,000-27,250,000 (110 kb) 961 px # chr7:27,000,000-28,000,000 (1 Mb) 1937 px # chr7 (whole chromosome) 3345 px # There is no automatic fallback at any of those, so the stanza has the same # maxWindowCoverage as methaDory. It was 200000 at first and raised to 10000000 on # 2026-09-18, because 200 kb sent most gene-neighbourhood views to the coverage graph. # Above that the track draws as a coverage graph, 81 px, and below it stays in pack. # # Heights at 1100 px wide, HOXA hot spot unless noted, for the tall views this allows: # epigenCentral methaDory # 200 kb pack 1281 px 1889 px # 1 Mb pack 1937 px 5729 px # 5 Mb pack 2129 px 7777 px # 10 Mb pack 2193 px 8273 px # 10 Mb squish 553 px 825 px # 10 Mb, typical region pack 353 px 2769 px