7e87cadb469b4e0eb4fb7f973154357cfe7fc345 lrnassar Mon Sep 21 15:56:18 2026 -0700 QA fixes for the mei (Mobile Insertions) track collection. refs #37524 Fix two data bugs found during QA and rebuild the affected bigBeds. meiEul1dbToBed.py looked up samples and individuals by name, but euL1db joins on 1-based row numbers, so neither join ever matched and the individual count, tissues, clinical conditions and populations were empty on all 8,991 insertions while the contributing-samples table printed row numbers. Both loaders now key on the row number, the table prints the sample name, and the adjacent population filter is case-insensitive so it actually drops "unknown". meiHgsvc3CsvToBed.py took alt[1:] on every record, which dropped the first base of the element on the 96 GRCh38 and 111 T2T-CHM13 records where PALMER2 is the only caller and ALT carries no anchor base; it now prefers INFO SEQ, which always matches SVLEN. Correct seven statements on the description pages against their sources: the HGSVC3 single-caller split was attributed to PALMER rather than L1ME-AID, its orthogonal concordance was 90.8% rather than 92.5%, euL1db was credited with aligning the L1HS consensus when the paper says it was processed from our RepeatMasker track, DeepMEI's network was described as a classifier rather than a genotyper and given the wrong training set, euL1db listed two detection methods absent from the data, and HMEID contradicted itself on the MELT ASSESS cutoff. Also: the SweGen bigDataUrl now points at _swegen.bb so the restricted callset is kept off the download server; the container page no longer claims the whole collection is long-read, lists the two euL1db subtracks, scopes its display conventions to the subtracks they describe, and cites all six papers; dead and wrong track links are repointed and pinned to a db; $db replaces hardcoded hg38 in paths on pages that serve three assemblies; the euL1db labels no longer carry hg38 counts and a lift note that made no sense on hg19; all six subtracks gain a dataVersion; the euL1db filter ranges match the data; and five autoSql field descriptions match what the files contain. Document the gbdb symlinks and the QA changes in doc/hg38/mei.txt, correct the HMEID bedToBigBed type there, and add an hg19.txt pointer since hg19 carries the two euL1db subtracks. diff --git src/hg/makeDb/doc/hg38/mei.txt src/hg/makeDb/doc/hg38/mei.txt index 14a96541bec..51ad3ef85d1 100644 --- src/hg/makeDb/doc/hg38/mei.txt +++ src/hg/makeDb/doc/hg38/mei.txt @@ -1,274 +1,333 @@ # 2026-05-09 Claude (max) - Mobile Element Insertions track collection (mei) # Source: HGSVC3 (Logsdon et al. 2025, Nature, PMID 40702183) # https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC3/release/Mobile_Elements/1.0/ # README: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC3/release/Mobile_Elements/1.0/README.20241211.MEI.txt # This track collection holds polymorphic Mobile Element Insertions (MEIs). # The first subtrack, meiHgsvc3, is the HGSVC3 MEI callset: mobile element # insertions identified in 65 long-read assembled samples relative to the # reference assembly. Two parallel callsets are released, one against # GRCh38 and one against T2T-CHM13, and we build a bigBed for each. # Each item is drawn as a 1-bp anchor block at the insertion attachment # site; per-sample genotypes are summarised into alt-allele count, allele # number, alt-allele frequency, and a list of carrier samples. ############################################################ # GRCh38 / hg38 mkdir -p /hive/data/genomes/hg38/bed/mei cd /hive/data/genomes/hg38/bed/mei wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC3/release/Mobile_Elements/1.0/MEI_Callset_GRCh38.ALL.20241211.csv.gz # Convert CSV (VCF-like, 65 sample genotype columns + Caller_Count, # TE_Designation, L1ME-AID, PALMER, L1ME-AID_INFO, PALMER_INFO, # PAVMergedCalls) to bed9+15. The script tallies per-record alt-allele # counts and carrier sample lists, and colors items by mobile element # class. # Source: ~/kent/src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py \ MEI_Callset_GRCh38.ALL.20241211.csv.gz \ /hive/data/genomes/hg38/chrom.sizes \ meiHgsvc3.bed # -> Read 12642 records, wrote 12642, skipped 0 + 0 # Class distribution: Alu 10270, L1 1604, SVA 764, HERVK 3, snRNA 1. sort -k1,1 -k2,2n meiHgsvc3.bed > meiHgsvc3.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiHgsvc3.as \ -type=bed9+16 \ meiHgsvc3.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ meiHgsvc3.bb ############################################################ # T2T-CHM13 / hs1 mkdir -p /hive/data/genomes/hs1/bed/mei cd /hive/data/genomes/hs1/bed/mei wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC3/release/Mobile_Elements/1.0/MEI_Callset_T2T-CHM13.ALL.20241211.csv.gz python3 ~/kent/src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py \ MEI_Callset_T2T-CHM13.ALL.20241211.csv.gz \ /hive/data/genomes/hs1/chrom.sizes \ meiHgsvc3.bed # -> Read 12919 records, wrote 12919, skipped 0 + 0 # Class distribution: Alu 10458, L1 1664, SVA 791, HERVK 5, snRNA 1. sort -k1,1 -k2,2n meiHgsvc3.bed > meiHgsvc3.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiHgsvc3.as \ -type=bed9+16 \ meiHgsvc3.sorted.bed \ /hive/data/genomes/hs1/chrom.sizes \ meiHgsvc3.bb ############################################################ # DeepMEI 1000G callset (hg38 only) # Source: Xu et al. 2023, bioRxiv (10.1101/2023.03.07.531451) # https://github.com/xuxif/DeepMEI/tree/main/DeepMEI/1000g_high_callset # DeepMEI is a CNN MEI caller; the authors released a high-confidence # callset for the 3,202 high-coverage 1000 Genomes samples (NYGC). # The VCF uses symbolic ALTs (<INS:ME:ALU>, <INS:ME:LINE1>, <INS:ME:SVA>) # and does not report inserted sequence or insertion length, so the # resulting bigBed schema is a subset of the HGSVC3 one. mkdir -p /hive/data/genomes/hg38/bed/mei/deepmei cd /hive/data/genomes/hg38/bed/mei/deepmei wget https://github.com/xuxif/DeepMEI/raw/refs/heads/main/DeepMEI/1000g_high_callset/merge_1000g.latested.vcf.gz # Convert VCF (91617 MEIs, 3202 samples) to bed9+7. # Source: ~/kent/src/hg/makeDb/scripts/mei/meiDeepmei1kgVcfToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiDeepmei1kgVcfToBed.py \ merge_1000g.latested.vcf.gz \ /hive/data/genomes/hg38/chrom.sizes \ deepmei.bed # -> Read 91617 records, wrote 91617, skipped 0 + 0 + 0 # Class distribution: Alu 68282, L1 16891, SVA 6444. sort -k1,1 -k2,2n deepmei.bed > deepmei.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiDeepmei1kg.as \ -type=bed9+7 \ deepmei.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ deepmei1kg.bb ############################################################ # 2026-05-12 Claude (max) - HMEID v1.1 (hg38 only) # Source: Niu et al. 2022, Nucleic Acids Research, PMID 35212372 # http://bigdata.ibp.ac.cn/HMEID/ # HMEID is a site-level catalogue of 36,699 non-reference MEIs called # by MELT v2.1.5 on Illumina short-read WGS of 5,675 individuals: # 2,998 NyuWa (Chinese, ~26.2x) + 2,677 1000 Genomes (~7.4x), aligned # to GRCh38. The VCF carries per-cohort (NyuWa, 1KGP) and per-1KGP- # super-population (AFR, AMR, EAS, EUR, SAS) AC/AN/AF in INFO; there # are no per-sample genotype columns. SVTYPE is one of ALU/LINE1/SVA/ # HERVK, plus the MELT TSD and ASSESS fields. mkdir -p /hive/data/genomes/hg38/bed/mei/hmei cd /hive/data/genomes/hg38/bed/mei/hmei wget http://bigdata.ibp.ac.cn/HMEID/static/download/MEI.GRCh38.HMEIDv1.1.vcf.gz wget http://bigdata.ibp.ac.cn/HMEID/static/download/sample_info.HMEIDv1.1.txt.gz # Convert site-level VCF (36699 MEIs) to bed9+27. # Source: ~/kent/src/hg/makeDb/scripts/mei/meiHmeidVcfToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiHmeidVcfToBed.py \ MEI.GRCh38.HMEIDv1.1.vcf.gz \ /hive/data/genomes/hg38/chrom.sizes \ meiHmeid.bed # -> Read 36699 records, wrote 36699, skipped 0 + 0 + 0 # Class distribution: Alu 26553, HERVK 126, L1 7353, SVA 2667. sort -k1,1 -k2,2n meiHmeid.bed > meiHmeid.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiHmeid.as \ - -type=bed9+27 \ + -type=bed9+28 \ meiHmeid.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ meiHmeid.bb ############################################################ # 2026-05-13 Claude (max) - SweGen MELT MEI callset (hg38, lifted from hg19) # Source: Ameur et al. 2017, Eur J Hum Genet, PMID 28832569 (SweGen cohort) # Gardner et al. 2017, Genome Res, PMID 28855259 (MELT tool) # https://swefreq.nbis.se/dataset/SweGen # Site-level MELT v2.0.2 MEI callset on 1,000 Swedish WGS samples # (SweGen, Illumina HiSeq X, 150 bp PE, BWA-MEM v0.7.12 to GRCh37). # The VCF has no per-sample columns; INFO carries MELT_AN (allele # count, despite the name) and MELT_AF (allele frequency) plus the # usual MELT fields SVTYPE/SVLEN/TSD/ASSESS/MEIINFO/INTERNAL. # Coordinates are GRCh37 with contig names like '1' (no chr prefix); # we add chr in Python and then liftOver hg19 -> hg38. mkdir -p /hive/data/genomes/hg38/bed/mei/swegen cd /hive/data/genomes/hg38/bed/mei/swegen # Source VCF must be requested from https://swefreq.nbis.se/dataset/SweGen/download # (Swegen_MELT_16032018.zip); place MELT_SWEGEN.20180314.ALU_HERVK_LINE1_SVA.vcf # under Swegen_MELT_16032018/. # Parse the VCF to a bed9+9 file with GRCh37 coords (adds chr prefix). # Source: ~/kent/src/hg/makeDb/scripts/mei/meiSwegenVcfToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiSwegenVcfToBed.py \ Swegen_MELT_16032018/MELT_SWEGEN.20180314.ALU_HERVK_LINE1_SVA.vcf \ meiSwegen.hg19.bed # -> Read 18100 records, wrote 18100, skipped 0 (unknown SVTYPE) # Class distribution: Alu 14467, HERVK 73, L1 2429, SVA 1131. # Lift to hg38 (-tab -bedPlus=9 to keep extra fields intact). liftOver -tab -bedPlus=9 \ meiSwegen.hg19.bed \ /gbdb/hg19/liftOver/hg19ToHg38.over.chain.gz \ meiSwegen.hg38.bed \ meiSwegen.unmapped # -> 18090 mapped, 10 unmapped (all "Deleted in new"). sort -k1,1 -k2,2n meiSwegen.hg38.bed > meiSwegen.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiSwegen.as \ -type=bed9+9 \ meiSwegen.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ meiSwegen.bb ############################################################ # 2026-05-13 Claude (max) - euL1db (hg19 + hg38 by liftOver) # Source: Mir et al. 2015, Nucleic Acids Research, PMID 25352549 # http://eul1db.unice.fr/ (Download tab -> eul1db.zip) # euL1db is a curated database of L1-HS retrotransposon insertion # polymorphisms catalogued from 32 published studies (~140k sample-level # SRIPs aggregating into ~9k non-redundant MRIPs). The original # coordinates are hg19. We build the bigBed on hg19 from the MRIP table # joined with SRIP/Sample/Individual/Study/Methods, then liftOver to # hg38. A second small bigBed for the reference-genome L1HS catalogue # (ReferenceL1HS.txt) is built and lifted the same way. # Source files (released as eul1db.zip, version 1.00, 2014-10-14): # Family.txt Individuals.txt MRIP.txt Methods.txt ReferenceL1HS.txt # SRIP.txt Samples.txt Study.txt # We keep the source under hg38/bed/mei/eul1db (where it was first placed) # and the hg19 build outputs under hg19/bed/mei/eul1db. mkdir -p /hive/data/genomes/hg19/bed/mei/eul1db cd /hive/data/genomes/hg19/bed/mei/eul1db # 1) Build hg19 MRIP bigBed (8,991 MRIPs in the source; chr23/chr24 in # Helman2014 rows are renamed to chrX/chrY in the script). # Source: ~/kent/src/hg/makeDb/scripts/mei/meiEul1dbToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiEul1dbToBed.py \ --src /hive/data/genomes/hg38/bed/mei/eul1db \ --chrom-sizes /hive/data/genomes/hg19/chrom.sizes \ -o eul1db.hg19.bed # -> SRIPs read: 142,495; MRIPs read: 8,991; MRIPs written: 8,991 # -> chr23/24 renamed to chrX/chrY: 85; no chrom or range skips. sort -k1,1 -k2,2n eul1db.hg19.bed > eul1db.hg19.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiEul1db.as \ -type=bed9+19 \ eul1db.hg19.sorted.bed \ /hive/data/genomes/hg19/chrom.sizes \ eul1db.hg19.bb # 2) Build hg19 reference-L1HS bigBed (1,544 elements). # Source: ~/kent/src/hg/makeDb/scripts/mei/meiEul1dbRefToBed.py python3 ~/kent/src/hg/makeDb/scripts/mei/meiEul1dbRefToBed.py \ --src /hive/data/genomes/hg38/bed/mei/eul1db \ --chrom-sizes /hive/data/genomes/hg19/chrom.sizes \ -o eul1dbRef.hg19.bed # -> Reference L1HS read: 1,544; written: 1,544 (all chroms in hg19). sort -k1,1 -k2,2n eul1dbRef.hg19.bed > eul1dbRef.hg19.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiEul1dbRef.as \ -type=bed9+6 \ eul1dbRef.hg19.sorted.bed \ /hive/data/genomes/hg19/chrom.sizes \ eul1dbRef.hg19.bb # 3) liftOver to hg38 (-tab -bedPlus=9 to preserve extra fields). liftOver -tab -bedPlus=9 \ eul1db.hg19.sorted.bed \ /gbdb/hg19/liftOver/hg19ToHg38.over.chain.gz \ eul1db.hg38.bed eul1db.hg38.unmapped # -> 8,988 mapped, 3 unmapped (1 "Deleted in new", 1 "Partially deleted", # 1 "Deleted in new" on chr13/chrX). 99.97% mapped. liftOver -tab -bedPlus=9 \ eul1dbRef.hg19.sorted.bed \ /gbdb/hg19/liftOver/hg19ToHg38.over.chain.gz \ eul1dbRef.hg38.bed eul1dbRef.hg38.unmapped # -> 1,540 mapped, 4 unmapped (2 "Split in new", 2 "Partially deleted"). # 4) Build hg38 bigBeds (written to hg38/bed/mei/eul1db/ next to source). sort -k1,1 -k2,2n eul1db.hg38.bed \ > /hive/data/genomes/hg38/bed/mei/eul1db/eul1db.hg38.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiEul1db.as \ -type=bed9+19 \ /hive/data/genomes/hg38/bed/mei/eul1db/eul1db.hg38.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ /hive/data/genomes/hg38/bed/mei/eul1db/eul1db.hg38.bb sort -k1,1 -k2,2n eul1dbRef.hg38.bed \ > /hive/data/genomes/hg38/bed/mei/eul1db/eul1dbRef.hg38.sorted.bed bedToBigBed -tab \ -as=$HOME/kent/src/hg/makeDb/scripts/mei/meiEul1dbRef.as \ -type=bed9+6 \ /hive/data/genomes/hg38/bed/mei/eul1db/eul1dbRef.hg38.sorted.bed \ /hive/data/genomes/hg38/chrom.sizes \ /hive/data/genomes/hg38/bed/mei/eul1db/eul1dbRef.hg38.bb + +############################################################ +# /gbdb symlinks for the whole collection +# The trackDb stanza lives in trackDb/human/mei.ra and uses $D, so each +# assembly needs its own symlink directory. The SweGen file is +# underscore-prefixed because the callset cannot be redistributed (see +# below), which keeps it off the public download server; the stanza also +# carries "tableBrowser off", which blocks both hgTables and the REST API. + +mkdir -p /gbdb/hg38/mei /gbdb/hs1/mei /gbdb/hg19/mei + +ln -s /hive/data/genomes/hg38/bed/mei/meiHgsvc3.bb /gbdb/hg38/mei/hgsvc3.bb +ln -s /hive/data/genomes/hg38/bed/mei/deepmei/deepmei1kg.bb /gbdb/hg38/mei/deepmei1kg.bb +ln -s /hive/data/genomes/hg38/bed/mei/hmei/meiHmeid.bb /gbdb/hg38/mei/hmeid.bb +ln -s /hive/data/genomes/hg38/bed/mei/swegen/meiSwegen.bb /gbdb/hg38/mei/_swegen.bb +ln -s /hive/data/genomes/hg38/bed/mei/eul1db/eul1db.hg38.bb /gbdb/hg38/mei/eul1db.bb +ln -s /hive/data/genomes/hg38/bed/mei/eul1db/eul1dbRef.hg38.bb /gbdb/hg38/mei/eul1dbRef.bb +ln -s /hive/data/genomes/hs1/bed/mei/meiHgsvc3.bb /gbdb/hs1/mei/hgsvc3.bb +ln -s /hive/data/genomes/hg19/bed/mei/eul1db/eul1db.hg19.bb /gbdb/hg19/mei/eul1db.bb +ln -s /hive/data/genomes/hg19/bed/mei/eul1db/eul1dbRef.hg19.bb /gbdb/hg19/mei/eul1dbRef.bb + +############################################################ +# 2026-09-21 Lou - QA fixes (refs #37524) + +# 1) meiHgsvc3: inserted sequence on the PALMER-only records. +# The source CSV merges three different INFO schemas. Most records follow +# the VCF convention (ALT[0] == REF[0], len(ALT)-1 == SVLEN), but the 96 +# GRCh38 / 111 T2T-CHM13 records whose INFO carries SEQ/CALLERS do not: +# there ALT is the element itself (len(ALT) == SVLEN, and on 70 of the 96 +# ALT[0] differs from REF[0]), and 13 of them carry a truncated ALT. +# meiHgsvc3CsvToBed.py took alt[1:] unconditionally, which dropped the +# first base of the element on all 96 and left insertSeq disagreeing with +# svLen. The script now prefers INFO SEQ when present, which always +# matches SVLEN. Rebuilt both assemblies with the commands above; +# item counts and every other column are unchanged. + +# 2) meiEul1db: sample and individual joins. +# euL1db's foreign keys are 1-based row numbers, not names: SRIP.txt's +# sampleID column indexes Samples.txt (943 rows) and Samples.txt's +# Individual_id column indexes Individuals.txt (741 rows). Both totals +# match Table 1 of Mir et al. 2015. meiEul1dbToBed.py keyed both lookup +# tables by name, so neither join ever matched and individualCount, +# tissues, diseases and populations were empty on all 8,991 MRIPs, while +# the contributing-samples table printed row numbers instead of sample +# names. Both loaders now key on the row number and the sample table +# prints the name. The adjacent population filter compared against +# "Unknown" while the data says "unknown", so it never fired either; it is +# now case-insensitive. After the fix: individualCount set on 8,976 +# records, tissues 8,976, diseases 6,781, populations 7,589 (hg19; 7,586 +# after the lift to hg38, which loses 3 MRIPs). Rebuilt the +# hg19 bigBed and re-lifted to hg38 with the commands above; item counts +# (8,991 / 8,988) and all other columns are unchanged. + +# 3) meiSwegen: /gbdb/hg38/mei/swegen.bb renamed to _swegen.bb. The +# callset cannot be redistributed, and the underscore prefix is what keeps +# a file off the public download server. "tableBrowser off" already +# blocked hgTables and the REST API (verified: the API returns HTTP 403 +# "protected data"), but without the prefix the bigBed itself would have +# been mirrored to hgdownload. mei.ra was updated to match.