9a49afb16653363a70d8e4d205513008b7b08df5 max Wed May 13 06:18:42 2026 -0700 varFreqs: add UK Biobank subtrack from Neale Lab Round 2 imputed-v3 variant manifest (13.7M variants, 361k white British samples). TSV → VCF conversion + CrossMap hg19→hg38, refs #36642 Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> diff --git src/hg/makeDb/doc/hg38/varFreqs.txt src/hg/makeDb/doc/hg38/varFreqs.txt index ab2669beda4..9fc3d9b5834 100644 --- src/hg/makeDb/doc/hg38/varFreqs.txt +++ src/hg/makeDb/doc/hg38/varFreqs.txt @@ -501,15 +501,62 @@ # trackDb filter UI. The individual `npm` subtrack still has # tableBrowser off (license), but its data is folded into varFreqsAll # same as for finngen / kova / mgrb / swefreq / tishkoff180 (precedent). cd /hive/data/genomes/hg38/bed/varFreqs/all/ # Force re-merge + re-csq (existing per-VCF normalize cache is reused). rm -f merged.vcf.gz merged.vcf.gz.tbi merged.annotated.vcf.gz \ merged.annotated.vcf.gz.tbi rm -rf extracted beds varFreqsAll.bed ./mergeAndAnnotate.sh python3 vcfToBigBed.py --output-prefix varFreqsAll --threads 8 # Build is in-place: /gbdb/hg38/varFreqs/varFreqsAll.bb stays symlinked # at /hive/.../all/varFreqsAll.bb, which is overwritten at the end of # the bedToBigBed step. +########## +# 2026-05-12 Claude max +# UK Biobank Neale Lab Round 2 variant manifest, refs #36642 (varFreqs). +# +# Source: variants.tsv.bgz (~916 MB, 13,791,467 rows), the Neale Lab +# Round 2 UKB GWAS variant table with imputation INFO score, post-QC +# allele counts and minor-allele frequencies from 361,194 imputed UK +# Biobank samples restricted to white British ancestry. The TSV ships +# bare chromosome names ("1".."22","X") in GRCh37 coordinates, so we +# convert to a sites-only VCF, then lift to hg38. +cd /hive/data/genomes/hg38/bed/varFreqs/ukbb +# Step 1: TSV -> VCF (hg19, chr-prefixed). Drops AC=0 rows (variants +# present in the imputation panel but not observed in the cohort), sets +# AN = 2 * n_called, and carries the imputation INFO score, HWE p-value, +# hom-ref / het / hom-alt sample counts, and the VEP consequence / +# Neale Lab consequence_category as INFO fields. +python3 ~/kent/src/hg/makeDb/scripts/varFreqs/ukbbToVcf.py +# total=13,791,467 written=13,751,808 skipped_AC0=39,659 +# Step 2: lift to hg38 with CrossMap. +CrossMap.py vcf /gbdb/hg19/liftOver/hg19ToHg38.over.chain.gz \ + ukbb.hg19.vcf \ + /hive/data/genomes/hg38/bed/varFreqs/all/hg38.fa \ + ukbb.hg38.unsorted.vcf +# Total entries: 13,751,808, Failed to map: 6,889 -> 13,744,919 lifted. +# Step 3: filter to canonical chromosomes (chr1-22/X/Y/M) in both the +# data lines and the contig header (CrossMap emits every hg38 contig). +awk ' + /^##contig=/ { + if (match($0, /ID=([^,>]+)/, m) && m[1] ~ /^chr([1-9]|1[0-9]|2[0-2]|X|Y|M)$/) print + next + } + /^#/ {print} +' ukbb.hg38.unsorted.vcf > header.vcf +awk -F'\t' 'BEGIN{OFS="\t"} !/^#/ && $1 ~ /^chr([1-9]|1[0-9]|2[0-2]|X|Y|M)$/ {print}' \ + ukbb.hg38.unsorted.vcf > body.vcf +cat header.vcf body.vcf > ukbb.hg38.canon.vcf +# 13,743,085 variants survive (1,834 of the lifted variants landed on +# alt/random/Un/fix contigs and were dropped). +# Step 4: sort, bgzip, tabix. +bcftools sort ukbb.hg38.canon.vcf -Oz -m 16G -T /data/tmp/ -o ukbb.vcf.gz +tabix -p vcf ukbb.vcf.gz +rm -f ukbb.hg19.vcf ukbb.hg38.unsorted.vcf ukbb.hg38.unsorted.vcf.unmap \ + ukbb.hg38.canon.vcf header.vcf body.vcf +# Final: 551 MB bgzip + 1.7 MB tabix index, 13,743,085 variants. +# 48,382 of the 13,791,467 source rows (0.35%) are not represented in +# hg38: 39,659 had AC=0 (panel-only), 6,889 failed CrossMap, 1,834 +# lifted to non-canonical contigs.