8c809888b88a970d262577c31ad2d5a874659788 max Thu Sep 17 07:55:41 2026 -0700 hg38: episignatures container with the MethaDory CpG probes #Preview2 week - bugs introduced now will need a build patch to fix New alpha superTrack for the CpG positions that published DNA methylation episignatures are built from. Its first subtrack, MethaDory, holds 268,900 sites drawn from 189 episignatures for 105 rare developmental disorders in 74 studies, compiled by Federico Ferraro and Dmitrijs Rots at Erasmus MC and given to us for open release. One row per CpG, merging every episignature that reports it, since a site is shared by 3.1 signatures on average and by up to 41. The per-signature values line up as a real table on the details page via detailsDynamicTable; the same values are also kept as plain columns for the Table Browser and for the disorder, gene, study and direction filters. Probe IDs were placed from the Illumina manifests already on the browser, EPIC 850k first, then EPIC v2. Of 847,863 probe-by-signature records, 373 were dropped because the probe is in none of the manifests; 19,996 of 20,000 sampled sites land on a CG dinucleotide. Colour is the direction and size of the strongest effect at the site, split at a delta-beta of 0.10, plus a conflicting class. Quantile bins were avoided on purpose: the studies used different reporting cut-offs, so the low end of the distribution reflects what each paper chose to publish rather than biology. refs #38371 diff --git src/hg/makeDb/doc/hg38/episignatures.txt src/hg/makeDb/doc/hg38/episignatures.txt new file mode 100644 index 00000000000..0d3b7a55368 --- /dev/null +++ src/hg/makeDb/doc/hg38/episignatures.txt @@ -0,0 +1,171 @@ +# 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