e8144b87c64b0e60e475886e3ced5353f139cf43 max Fri Sep 18 07:24:42 2026 -0700 hg38 episignatures: EpigenCentral CpG probes as a second subtrack #Preview2 week - bugs introduced now will need a build patch to fix The EpigenCentral group at the Centre for Computational Medicine and the Weksberg lab, Hospital for Sick Children, publishes its curated episignatures as a track hub and asked us to host it natively instead. 15,035 CpG sites from 24 episignatures for 23 rare disorders, alongside MethaDory in the same container. The two were compiled independently and neither contains the other: 14,236 sites are in both, at identical coordinates, and 799 are only here. Built from the lab's own bigBed by makeEpigenCentral.sh, so a refresh is one command rather than the hand edits it started as. Every feature in the hub is in the track; no coordinates or values were changed. Four things were rewritten on the way in: the OMIM column went from a full URL to the entry number so trackDb builds the link, the hub's pre-rendered mouse-over column became the direction alone with the text assembled from the fields, the disorder of the strongest signature was added as its own column, and 215 rows of the per-site comparison table on 137 sites were dropped as exact duplicates, which the site's own signature count had already collapsed. Two columns were renamed to what MethaDory calls the same numbers, height to maxAbsDelta and nSignatures to sigCount. Coordinates check out two ways: 5,000 sampled sites all land on a CG dinucleotide, and every probe shared with MethaDory, which was positioned from the Illumina manifests independently, is at the same base in both. The description page keeps the structure and wording written for the hub but not its markup. Its reference table had all 24 links pointing at one PMID while displaying another, so that table is now generated from the data joined to a checked-in PMID list, each one verified against PubMed. tableBrowser is off at the request of the data providers, so Data Access points at EpigenCentral's own portal and repository. refs #38112 diff --git src/hg/makeDb/doc/hg38/episignatures.txt src/hg/makeDb/doc/hg38/episignatures.txt index 0d3b7a55368..7ec1c97b09a 100644 --- src/hg/makeDb/doc/hg38/episignatures.txt +++ src/hg/makeDb/doc/hg38/episignatures.txt @@ -1,171 +1,292 @@ # 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 # "<!-- BEGIN generated ... -->" and "<!-- END generated -->" 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. # 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 "<!-- BEGIN generated ... -->" and "<!-- END generated -->" 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. +# +# 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. + +# 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 200000 as methaDory. Above 200 kb the track draws as a coverage +# graph, 81 px, and below it stays in pack with the probe IDs, which is the point of +# zooming in. The 110 kb view is still tall; a reader who does not want the labels can +# switch that one to squish, 121 px.