c656fab2dcb67f0789af54ea1d04afcf6ab3d07b max Wed Jul 15 23:38:55 2026 -0700 varFreqs: expand SPARK phenotype-split documentation in makeDoc #Preview2 week - bugs introduced now will need a build patch to fix Rework the SFARI SPARK section of the varFreqs makeDoc to describe, step by step, what sparkMergeVcfAddCounts.sh does per chromosome (fill-tags overall + optional per-group AC/AN/AF, drop genotypes, concat) and how the individuals_registration TSV drives the ASD / non-ASD phenotype split. diff --git src/hg/makeDb/doc/hg38/varFreqs.txt src/hg/makeDb/doc/hg38/varFreqs.txt index dffe798c1b3..9c6457c46b0 100644 --- src/hg/makeDb/doc/hg38/varFreqs.txt +++ src/hg/makeDb/doc/hg38/varFreqs.txt @@ -83,30 +83,34 @@ # NPM Singapore cd /hive/data/genomes/hg38/bed/varFreqs/npm/ # downloaded data manually from chorus website, https://chorus.grids-platform.io/vcfdl bcftools concat --threads 10 -Oz -o SG10K_Health_r5.3.2.sites.vcf.bgz SG10K_Health_r5.3.2.sites.chr{1..22}.vcf.bgz SG10K_Health_r5.3.2.sites.chrX.vcf.bgz SG10K_Health_r5.3.2.sites.chrY.vcf.bgz tabiv -p vcf SG10K_Health_r5.3.2.sites.vcf.bgz # Saudi 300 genomes cd /hive/data/genomes/hg38/bed/varFreqs/saudi wget https://figshare.com/ndownloader/files/51297884 -O 51297884.tsv.gz python3 ~/kent/src/hg/makeDb/scripts/varFreqs/saudiToVcf.py bgzip saudi.vcf tabix -p vcf saudi.vcf.gz # SFARI SPARK +# Initial build (overall AC/AN/AF only). The current build adds ASD/non-ASD +# phenotype-stratified counts and is documented under the 2026-05-28 entry +# below ("Phenotype-stratified SPARK counts"), which also describes step by +# step what sparkMergeVcfAddCounts.sh and the bcftools csq annotation do. cd /hive/data/genomes/hg38/bed/varFreqs/sparkExomes/ # used globus to download into vcf/ sh ~/kent/src/hg/makeDb/scripts/varFreqs/sparkMergeVcfAddCounts.sh vcf/SPARK.iWES_v3.2024_08.deepvariant 8 bcftools norm -m- SPARK.iWES_v3.2024_08.deepvariant.sites.vcf.gz -Oz > SPARK.iWES_v3.2024_08.deepvariant.norm.vcf.gz && tabix -p vcf SPARK.iWES_v3.2024_08.deepvariant.norm.vcf.gz cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs/ # used globus to download into vcf/ sh ~/kent/src/hg/makeDb/scripts/varFreqs/sparkMergeVcfAddCounts.sh vcf/wgs_12519_genome.deepvariant 8 bcftools norm -m- wgs_12519_genome.deepvariant.sites.vcf.gz -Oz > wgs_12519_genome.deepvariant.norm.vcf.gz tabix -p vcf wgs_12519_genome.deepvariant.norm.vcf.gz # NCBI ALFA bigBed to VCF, Max Jan 26 2026 # Source: ALFA R4 bigBed files, 904M variants, output 163M with non-zero AF cd /hive/data/genomes/hg38/bed/varFreqs/alfa python3 ~/kent/src/hg/makeDb/scripts/varFreqs/alfa_to_vcf.py --out ALFA.vcf --zero-af-file ALFA_zero.txt @@ -876,37 +880,52 @@ # sites split, 3,559,402 indels left-realigned, 0 dup/mismatch). # GoNL added to databases.tsv as well, so the next mergeAndAnnotate.sh # run will pick it up into varFreqsAll.bb. The current /gbdb _all # bigBed predates GoNL. ########## # 2026-05-28 Claude max # varFreqs reorg: split SFARI SPARK into ASD/non-ASD counts, add SCHEMA # case/control counts, drop genotyping-array cohorts (TPMI + MexBB) from the # WGS/WES varFreqsAll track, and add two new combined tracks: # varFreqsDisease (SPARK, SFARI WGS, TOPMed, SCHEMA, GREGoR, GA4K) # varFreqsArray (TPMI, MexBB, UKBB) # refs #36642 # --- 1. Phenotype-stratified SPARK counts ----------------------------------- -# sparkMergeVcfAddCounts.sh now takes an optional 3rd arg: the SPARK -# individuals_registration TSV. With it, bcftools +fill-tags is given a -# sample-group file (col1=sample_sp_id, col2=AUT|NON_AUT) so the output -# sites VCF carries AC_AUT/AN_AUT/AF_AUT and AC_NON_AUT/AN_NON_AUT/AF_NON_AUT -# in addition to the overall AC/AN/AF. The same registration file (col 8 -# "asd": TRUE/FALSE) covers both WES (142,357 samples) and WGS (12,519 -# samples); samples with a blank asd value are not assigned to either group. +# What sparkMergeVcfAddCounts.sh does, per chromosome and in parallel: +# 1. bcftools +fill-tags <chr.vcf.gz> -- -t AC,AN,AF [-S groups.txt] +# computes the overall AC/AN/AF from the genotypes. The optional -S +# sample-group file additionally emits per-group AC_<grp>/AN_<grp>/AF_<grp>. +# 2. bcftools view -G drops the genotype columns, leaving an anonymous, +# sites-only VCF (required: the genotype-level pVCFs cannot be shared). +# 3. bcftools concat joins the per-chromosome sites VCFs back into one file. +# After the script, bcftools norm -m- splits multi-allelic records into one +# row per ALT allele. There is no allele-frequency cutoff. +# +# The script now takes an optional 3rd arg: the SPARK individuals_registration +# TSV. With it, the -S sample-group file (col1=sample_sp_id, col2=AUT|NON_AUT) +# makes the output sites VCF carry AC_AUT/AN_AUT/AF_AUT and +# AC_NON_AUT/AN_NON_AUT/AF_NON_AUT in addition to the overall AC/AN/AF. The +# same registration file (col 8 "asd": TRUE/FALSE) covers both WES (142,357 +# samples) and WGS (12,519 samples); samples with a blank asd value are not +# assigned to either group. +# +# These .norm.vcf.gz files are what the per-cohort sfariSparkExomes / +# sfariSparkWgs tracks display. Predicted protein consequences (BCSQ) are +# added later, in the combined-track build (section 3 below), by the +# bcftools csq step inside mergeAndAnnotate.sh. REG=/hive/data/genomes/hg38/bed/varFreqs/sparkExomes/SPARKDataRelease_2025-12-15/individuals_registration-2025-12-15.tsv cd /hive/data/genomes/hg38/bed/varFreqs/sparkExomes/ sh ~/kent/src/hg/makeDb/scripts/varFreqs/sparkMergeVcfAddCounts.sh \ vcf/SPARK.iWES_v3.2024_08.deepvariant 8 "$REG" bcftools norm -m- SPARK.iWES_v3.2024_08.deepvariant.sites.vcf.gz \ -Oz -o SPARK.iWES_v3.2024_08.deepvariant.norm.vcf.gz tabix -p vcf SPARK.iWES_v3.2024_08.deepvariant.norm.vcf.gz # 13,299,149 variants out, identical position list to the pre-pheno build. cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs/ sh ~/kent/src/hg/makeDb/scripts/varFreqs/sparkMergeVcfAddCounts.sh \ vcf/wgs_12519_genome.deepvariant 8 "$REG" bcftools norm -m- wgs_12519_genome.deepvariant.sites.vcf.gz \ -Oz -o wgs_12519_genome.deepvariant.norm.vcf.gz @@ -929,30 +948,49 @@ mv SCHEMA_variant_results_hg38_sorted.vcf.gz SCHEMA_variant_results_hg38_sorted.vcf.gz.preCC.bak mv SCHEMA_variant_results_hg38_sorted.vcf.gz.tbi SCHEMA_variant_results_hg38_sorted.vcf.gz.tbi.preCC.bak mv SCHEMA_variant_results_hg38_sorted.withCC.vcf.gz SCHEMA_variant_results_hg38_sorted.vcf.gz mv SCHEMA_variant_results_hg38_sorted.withCC.vcf.gz.tbi SCHEMA_variant_results_hg38_sorted.vcf.gz.tbi # 8,864,488 variants, AC_CASE+AC_CTRL = AC, sums verified. # --- 3. Three combined-track builds ----------------------------------------- # mergeAndAnnotate.sh and vcfToBigBed.py now accept --databases / # --databases-file flags, so the same code drives three builds from three # config TSVs: databases.tsv, databases_disease.tsv, databases_array.tsv. # The normalized/ cache under all/normalized is shared across builds # (positions are identical regardless of which combined track they feed). # mergeAndAnnotate.sh pins /cluster/software/src/bcftools-1.22/ in front of # PATH because bcftools csq's --unify-chr-names was added in 1.22; the conda # default (1.14) silently fails the annotate step otherwise. +# +# What mergeAndAnnotate.sh does, given a databases.tsv list of cohort VCFs: +# - strip + normalize each input VCF (bcftools annotate -x to drop ID/QUAL/ +# FILTER/INFO, fix chr naming, then bcftools norm -m-any against hg38.fa). +# - bcftools merge -m none all normalized VCFs into one sites-only VCF. +# - bcftools csq annotates predicted protein consequences against the +# Ensembl GRCh38 release-115 GFF3 gene models. The exact call is: +# bcftools csq -p a -l -n 64 --unify-chr-names chr,-,chr --force \ +# -f hg38.fa -g Homo_sapiens.GRCh38.115.chr.gff3.gz \ +# merged.vcf.gz -Oz -o merged.annotated.vcf.gz +# The GFF3 is from +# https://ftp.ensembl.org/pub/release-115/gff3/homo_sapiens/Homo_sapiens.GRCh38.115.chr.gff3.gz +# (Ensembl, not GENCODE: csq needs the Parent=transcript: attribute format), +# sorted, bgzipped and tabix-indexed. --unify-chr-names reconciles the +# chr1-style VCF/fasta names with the bare "1" names in the Ensembl GFF3. +# The consequence + amino-acid change land in INFO/BCSQ; vcfToBigBed.py +# parses BCSQ into the bigBed for the combined tracks. +# The same csq recipe can be run directly on a single cohort's .norm.vcf.gz +# (e.g. the SFARI files above) to annotate it with protein effects. # Reorg removed TPMI and MexBB from databases.tsv (moved to # databases_array.tsv) and added UKBB to databases_array.tsv. Both # databases_disease.tsv and databases_array.tsv reuse populations.tsv: # load_config only attaches a population row to a database that is in the # active config, so the same populations file works for all three builds. # varFreqsAll (28 cohorts, all WGS/WES; ~1.34 billion variants, 48 GB): bash ~/kent/src/hg/makeDb/scripts/varFreqs/mergeAndAnnotate.sh python3 ~/kent/src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py \ --annotated-vcf /hive/data/genomes/hg38/bed/varFreqs/all/merged.annotated.vcf.gz \ --output-prefix varFreqsAll \ --threads 6 \ --work-dir /hive/data/genomes/hg38/bed/varFreqs/all ln -sfn /hive/data/genomes/hg38/bed/varFreqs/all/varFreqsAll.bb \