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/scripts/clinvarMapped/clinvarMappedParalogPairs.sh src/hg/makeDb/scripts/clinvarMapped/clinvarMappedParalogPairs.sh new file mode 100755 index 00000000000..b7f4cc35c88 --- /dev/null +++ src/hg/makeDb/scripts/clinvarMapped/clinvarMappedParalogPairs.sh @@ -0,0 +1,58 @@ +#!/bin/bash +# clinvarMappedParalogPairs.sh - fetch genome-wide within-species paralog pairs from +# Ensembl BioMart (release chosen by www.ensembl.org, currently 116). +# +# A single genome-wide martservice request does not close its socket cleanly and +# curl always hits --max-time, so we cannot tell a complete payload from a +# truncated one. Instead we query one chromosome at a time (chromosome_name +# filter): each chunk completes and closes normally, and the concatenation is +# provably complete. +# +# Output columns (TSV): +# geneId geneSym paralogId paralogSym percIdQ percIdT +# percIdQ = % of the query gene's protein identical to the paralog +# percIdT = % of the paralog's protein identical to the query gene +# Usage: clinvarMappedParalogPairs.sh +set -beEu -o pipefail + +outDir="${1:?usage: clinvarMappedParalogPairs.sh }" +out="$outDir/paralogPairs.tsv" +chunkDir="$outDir/paralogPairsChunks" +mkdir -p "$chunkDir" + +chroms=$(echo {1..22} X Y MT) + +fetchChrom() { + local chrom="$1" dest="$2" + local query + query=' + + + + + + + + + + + +' + curl -sS --retry 3 --retry-delay 10 --max-time 600 \ + -o "$dest" --data-urlencode "query=$query" \ + 'https://www.ensembl.org/biomart/martservice' + if grep -q "Query ERROR\|Exception" "$dest"; then + echo "ERROR: BioMart error for chrom $chrom:" 1>&2; head -3 "$dest" 1>&2; return 1 + fi +} + +for chrom in $chroms; do + dest="$chunkDir/chr${chrom}.tsv" + echo "$(date '+%F %T') fetching chromosome $chrom ..." 1>&2 + fetchChrom "$chrom" "$dest" + echo " chr${chrom}: $(wc -l < "$dest") rows" 1>&2 +done + +# Concatenate, drop rows with no paralog (empty paralog id). +cat "$chunkDir"/chr*.tsv | awk -F'\t' 'NF>=3 && $3!=""' | sort -u > "$out" +echo "$(date '+%F %T') wrote $(wc -l < "$out") paralog pair rows to $out" 1>&2