5fe93cbef0d7692552e37e62bbf394cc7d222a28 max Thu Jul 16 00:32:55 2026 -0700 ClinVar Mapped track: map ClinVar coding variants through paralog alignments #Preview2 week - bugs introduced now will need a build patch to fix Adds the clinvarMapped composite (hg38) with two subtracks: - clinvarMappedParalog: every protein-changing ClinVar variant projected onto the equivalent (aligned) residue of each of its gene's paralogs - clinvarMappedParalogAln: the pairwise protein alignments used to do the mapping, as bigPsl, so the evidence for each projection can be inspected Pipeline: paralog pairs from Ensembl BioMart, one representative transcript per gene from MANE Select, pairwise global protein alignment (BLOSUM62) for pairs at >=20% identity, variants projected residue-to-residue and mapped back to the paralog's genomic codon. Colors and filters mirror the ClinVar track; the map track is filterable by source gene, classification, review stars, residue conservation and percent identity. refs #37883 diff --git src/hg/makeDb/doc/hg38/clinvarMapped.txt src/hg/makeDb/doc/hg38/clinvarMapped.txt new file mode 100644 index 00000000000..f37a28a21a4 --- /dev/null +++ src/hg/makeDb/doc/hg38/clinvarMapped.txt @@ -0,0 +1,89 @@ +# 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)