c9446936b026a803d97e7f5128404a327829fe2e max Sun Jul 19 00:01:22 2026 -0700 ClinVar Mapped: add Pfam-domain mapping method Maps ClinVar coding variants through Pfam protein domains as a second method under the clinvarMapped container. Each Pfam family's profile HMM gives a shared coordinate, so a variant inside a domain is projected to the residue in the same HMM match-state column in every other MANE Select gene carrying that domain, then mapped back to that gene's genomic codon. This links genes that share a domain without being paralogs (e.g. the SCN sodium channels and the CACNA1 calcium channels via Ion_trans). Each variant is shown at its 25 most conserved equivalents. Uses Pfam-A 38.2. Adds two alpha subtracks (Pfam Domain Variants, Pfam Domains), the hmmsearch/hmmalign build scripts, autoSql, makeDoc and description page, refs #37883 diff --git src/hg/makeDb/doc/hg38/clinvarMapped.txt src/hg/makeDb/doc/hg38/clinvarMapped.txt index f37a28a21a4..be39b0689bc 100644 --- src/hg/makeDb/doc/hg38/clinvarMapped.txt +++ src/hg/makeDb/doc/hg38/clinvarMapped.txt @@ -1,89 +1,155 @@ # 2026-07-13 Claude (max): ClinVar variants mapped to related loci (clinvarMapped) # The clinvarMapped track container maps protein-changing ClinVar variants onto # related loci by sequence homology. The first (and so far only) method is # paralog mapping: every variant is projected onto the equivalent (aligned) # residue of each of its gene's paralogs, with the pairwise protein alignments # provided as evidence. Two bigBed outputs: # clinvarMappedParalog projected variants (bigBed 12+) # clinvarMappedParalogAln bigPsl alignments (with query CDS sequence) # (Future methods, e.g. mapping through PFAM domains, will add sibling subtracks.) # # Inputs, all already on hg38: # /gbdb/hg38/bbi/clinvar/clinvarMain.bb ClinVar short variants # /gbdb/hg38/mane/mane.bb MANE Select + Plus Clinical # /gbdb/hg38/hg38.2bit genome # Ensembl BioMart release 116 within-species paralog pairs # # All scripts are in kent/src/hg/makeDb/scripts/clinvarMapped . Python steps # run in the micromamba "claude" env (needs parasail, biopython, py2bit). # Scripts named clinvarMapped* are reusable infrastructure; clinvarMappedParalog* # are specific to the paralog mapping method. dir=/hive/data/genomes/hg38/bed/clinvarMapped scr=~/kent/src/hg/makeDb/scripts/clinvarMapped mkdir -p $dir && cd $dir # 1. Paralog pairs from Ensembl BioMart, one chromosome per request (a single # genome-wide request never closes its socket and always times out). $scr/clinvarMappedParalogPairs.sh $dir # -> paralogPairs.tsv (3,552,265 pairs; geneId geneSym paralogId paralogSym percIdQ percIdT) # 2. MANE Select transcript set: one representative transcript+protein per gene, # translated straight from the genome so protein and coordinate mapping agree. $scr/clinvarMappedMane.sh $dir # -> maneSelect.gp, maneSelect.faa, maneMeta.tsv (19,293 protein-coding genes) # 3. Coverage stats (genes with/without a paralog, percent-id distribution). $scr/clinvarMappedParalogStats.sh $dir # -> paralogStats.log (15,278 of 19,293 genes have a MANE-mappable paralog) # 4. Restrict to unordered pairs whose proteins share >=20% identity, and make a # per-gene protein FASTA keyed by Ensembl gene id. awk -F'\t' 'BEGIN{OFS="\t"} { q=$5+0; t=$6+0; m=(q>t?q:t); if(m<20) next; a=$1; b=$3; if(a<b){k=a"\t"b; sa=$2; sb=$4; p1=$5; p2=$6} else {k=b"\t"a; sa=$4; sb=$2; p1=$6; p2=$5} if(!(k in seen) || m>best[k]){seen[k]=sa"\t"sb; best[k]=m; pid[k]=p1"\t"p2} } END{ for(k in seen) print k, seen[k], pid[k], best[k] }' \ paralogPairs.mane.tsv | sort -k1,1 -k2,2 > pairs.ge20.tsv # -> pairs.ge20.tsv (94,471 pairs) awk -F'\t' '{print $1"\t"$2}' maneMeta.tsv > .enst2ensg awk 'BEGIN{while((getline l < ".enst2ensg")>0){split(l,a,"\t"); m[a[1]]=a[2]}} /^>/{ id=substr($1,2); g=m[id]; if(g=="") {skip=1; next} skip=0; print ">"g; next } !skip{print}' maneSelect.faa > geneProt.faa && rm -f .enst2ensg # 5. Pairwise global protein alignment (BLOSUM62) of every >=20% pair; output is # ungapped aligned blocks in protein coordinates (fast, ~3 s with parasail). python $scr/clinvarMappedParalogAlign.py geneProt.faa pairs.ge20.tsv -o alnBlocks.tsv -j 48 # -> alnBlocks.tsv (94,471 alignments) # 6. Assign every protein-changing ClinVar variant to its MANE codon. --selftest # verifies the genomic<->codon math; the result matches ClinVar's own protein # HGVS for 99.9% of missense variants. python $scr/clinvarMappedCodons.py /gbdb/hg38/bbi/clinvar/clinvarMain.bb \ maneSelect.gp geneProt.faa maneMeta.tsv -o clinvarCodons.tsv --selftest # -> clinvarCodons.tsv (2,623,549 variant-codon rows from 2,648,470 coding variants) # 7. Project each variant across its paralog alignments -> variant bigBed. Colors # match the ClinVar track (clinVarToBed). python $scr/clinvarMappedParalogProject.py $dir -o clinvarMappedParalog.bed -j 48 export LC_ALL=C sort -S8G --parallel=16 -k1,1 -k2,2n clinvarMappedParalog.bed -o clinvarMappedParalog.sorted.bed bedToBigBed -tab -type=bed12+14 -as=$scr/clinvarMappedParalog.as \ clinvarMappedParalog.sorted.bed /hive/data/genomes/hg38/chrom.sizes clinvarMappedParalog.bb # -> clinvarMappedParalog.bb (16,750,740 features) # 8. Turn the same alignments into PSL (both orientations) -> bigPsl evidence # track. --faOut writes a per-gene CDS FASTA (keyed by gene symbol = the PSL # qName) built from the same coordinates, so the embedded query sequence lets # the browser shade base matches/mismatches in the alignment. python $scr/clinvarMappedParalogAlnPsl.py $dir /hive/data/genomes/hg38/chrom.sizes \ -o clinvarMappedParalogAln.psl --faOut geneCds.fa --twoBit /gbdb/hg38/hg38.2bit pslToBigPsl clinvarMappedParalogAln.psl -fa=geneCds.fa stdout \ | sort -S4G --parallel=16 -k1,1 -k2,2n > clinvarMappedParalogAln.bigPslInput bedToBigBed -as=$HOME/kent/src/hg/lib/bigPsl.as -type=bed12+13 -tab \ clinvarMappedParalogAln.bigPslInput /hive/data/genomes/hg38/chrom.sizes clinvarMappedParalogAln.bb # -> clinvarMappedParalogAln.bb (188,942 alignment features, with query sequence) + +############################################################################# +# ClinVar Mapped - Pfam domain method (Claude, 2026-07) +# +# Second mapping method under the clinvarMapped container. Instead of pairwise +# paralog alignments, it uses each Pfam-A family's profile HMM as a shared +# coordinate: every instance of a domain is aligned to the same model, so a +# residue in gene A and a residue in gene B are the same position in the domain +# when they align to the same HMM match state (consensus column). This links all +# genes that carry a domain, not only annotated paralogs. Reuses the MANE Select +# proteins, CDS->genome maps and ClinVar codon assignments from the paralog steps +# above (geneProt.faa, maneSelect.gp, maneMeta.tsv, clinvarCodons.tsv). + +dir=/hive/data/genomes/hg38/bed/clinvarMapped +scr=$HOME/kent/src/hg/makeDb/scripts/clinvarMapped +hmmerBin=/hive/data/outside/pfam/Pfam29.0/PfamScan/hmmer-3.1b2-linux-intel-x86_64/binaries +mkdir -p $dir/pfam; cd $dir/pfam + +# 1. Pfam-A (Pfam 38.2), downloaded and pressed/indexed for hmmsearch/hmmfetch. +mkdir -p /hive/data/outside/pfam/Pfam38.2; cd /hive/data/outside/pfam/Pfam38.2 +wget https://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam38.2/Pfam-A.hmm.gz +gunzip Pfam-A.hmm.gz +cd $dir/pfam +ln -sf /hive/data/outside/pfam/Pfam38.2/Pfam-A.hmm Pfam-A.hmm +$hmmerBin/hmmpress Pfam-A.hmm +$hmmerBin/hmmfetch --index Pfam-A.hmm + +# 2. Domain membership: every Pfam family vs the 19,293 MANE proteins, at Pfam's +# curated per-family gathering thresholds (--cut_ga). +$hmmerBin/hmmsearch --cpu 15 --cut_ga --domtblout maneHits.domtbl -o hmmsearch.out \ + Pfam-A.hmm $dir/geneProt.faa +# -> 66,043 domain instances, 9,094 families, 18,750 of 19,293 genes + +# 3. Per-residue HMM column map. For each family present in >=2 genes, fetch the +# model and realign the matched subsequences with hmmalign; the #=GC RF line of +# the alignment gives each residue its match-state column. +python $scr/clinvarMappedPfamAlign.py $dir --hmmDir $dir/pfam --hmmerBin $hmmerBin -j 24 +# -> alnCols.tsv (5,826,029 residue-column rows), famStats.tsv +# 4,394 mappable families, 689,578 distinct (family,column) cells + +# 4. Project each ClinVar coding variant that sits on a domain column onto the +# equivalent residue of every other gene at that column, keeping the 25 most +# conserved targets per variant (identical > similar > different), then map +# each back to its genomic codon. Needs BLOSUM62 (biopython, env "claude"). +micromamba run -p $HOME/miniconda3/envs/claude python $scr/clinvarMappedPfamProject.py \ + $dir -o pfam/clinvarMappedPfam.bed --maxFanout 25 -j 40 +# -> 1,066,320 of 2,623,549 variants (40.6%) land on a domain column; +# 15,947,052 projected features (unbounded would be several hundred M) +# residue match: identical 9,513,991; similar 2,523,863; different 3,909,198 +export LC_ALL=C +sort -S8G --parallel=16 -k1,1 -k2,2n clinvarMappedPfam.bed -o clinvarMappedPfam.sorted.bed +bedToBigBed -tab -type=bed12+16 -as=$scr/clinvarMappedPfam.as \ + clinvarMappedPfam.sorted.bed /hive/data/genomes/hg38/chrom.sizes clinvarMappedPfam.bb + +# 5. Evidence track: each MANE domain instance placed on the genome (the +# MANE-restricted analog of ucscGenePfam), carrying its HMM column range. +python $scr/clinvarMappedPfamAln.py $dir -o pfam/clinvarMappedPfamAln.bed +# -> 66,043 domain footprints +sort -S4G --parallel=8 -k1,1 -k2,2n clinvarMappedPfamAln.bed -o clinvarMappedPfamAln.sorted.bed +bedToBigBed -tab -type=bed12+8 -as=$scr/clinvarMappedPfamAln.as \ + clinvarMappedPfamAln.sorted.bed /hive/data/genomes/hg38/chrom.sizes clinvarMappedPfamAln.bb + +# Validation: 400 random projected codons re-translated from hg38.2bit all matched +# the claimed residue (strand, split codons, blocks correct). SCN5A pathogenic +# Ion_trans variants map to the sodium-channel paralogs and, via the shared +# domain, to the calcium channels (CACNA1x, NALCN) that paralog annotation misses.