a5c699a7301156154700f51a20ae571bc6987051 max Sat Sep 26 17:53:41 2026 -0700 varFreqs: remove the AF-table-based sfariSparkWgs45k subtrack, superseded by the genotype-based sfariSparkWgs45kAsd; drop its scripts, the DSCAM demo script and their makeDoc sections, refs #38424 diff --git src/hg/makeDb/doc/hg38/varFreqs.txt src/hg/makeDb/doc/hg38/varFreqs.txt index ba0e8fbb278..0a2720c7d48 100644 --- src/hg/makeDb/doc/hg38/varFreqs.txt +++ src/hg/makeDb/doc/hg38/varFreqs.txt @@ -1248,137 +1248,75 @@ --split-affected \ --threads 8 \ --work-dir /hive/data/genomes/hg38/bed/varFreqs/all # Merged variants: 1,374,129,993 (HostSeq added ~38M sites not seen in any # other cohort). Outputs: varFreqsAffected.bb 18.7 GB / 133,290,997 items, # varFreqsBackground.bb 60.0 GB / 1,278,354,588 items; both 185 fields # (up from 165: HostSeq adds overall AC/AF + 9 ancestry-group AC/AF = 20 cols). # Spot-check APOE rs429358 in the background bb matches the standalone HostSeq # VCF exactly: HostSeqAC=2315, AF=0.1258, afr AF=0.234, nfe 0.130, oth 0.079. # The /gbdb _affected/_background symlinks point at the in-place .bb (unchanged). # varFreqs.ra: added HostSeq|HostSeq Canada to filterValues.backgroundSources in # both stanzas (HostSeq is NOT in affectedCohorts) plus commented per-DB # HostSeqAF/AC filter blocks. Length and AC/AN filter ranges were unchanged by # HostSeq (confirmed against the auto-generated all/varFreqs.trackDb.ra). -########## -# 2026-09-23 Claude max -# SFARI SPARK WGS August 2026 release (DS0000135), 45,178 genomes -> new -# subtrack sfariSparkWgs45k (doc page shared with sfariSparkExomes.html). -# Downloaded resources/ with the Globus UI. There is no sites VCF in the release, -# only an AF table made from the GLnexus pVCFs with -# bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' -# AF is rounded to 6 decimals, there is no AC or AN. All low AFs in the file, -# including chrX/chrY, are exact multiples of 1/90356 = 1/(2*45178): the 179M -# AF=1 alleles are all 1.1e-05, never 1.2e-05, so missingness is ~0 and AN is set -# to 90356 everywhere, AC=round(AF*AN). Multiallelics are split, singletons -# (AC=1) are dropped, then bcftools norm left-aligns against hg38. -cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs45k -bash ~/kent/src/hg/makeDb/scripts/varFreqs/sparkWgs45kToVcf.sh \ - resources/SPARK.WGS.2026_08.gatk.pvcf_variant_frequencies.tsv \ - /hive/data/genomes/hg38/bed/varFreqs/sparkWgs45k/sparkWgs45k.vcf.gz 24 > build.log 2>&1 -# Accounting: 436,319,161 input records (chr1-22, X, Y; no chrM, no alts) with -# 521,466,014 ALT alleles -> 179,411,092 singleton alleles dropped -> -# 342,054,922 VCF records. No '*' alleles, 0 REF mismatches in bcftools norm, -# 24.1M records realigned. Output 2.2 GB + tabix index. -# Spot check APOE: rs429358 chr19:44908684 T>C AF=0.1455, rs7412 -# chr19:44908822 C>T AF=0.0755. -rm -rf work # keep only the per-chrom logs, in logs/ - ########## # 2026-09-24 Claude max -# SFARI SPARK WGS 2026_08 from the genotype pVCFs -> new subtrack -# sfariSparkWgs45kAsd, with exact AC/AN/AF and ASD/non-ASD splits, no AC cutoff. -# For now only a demo: the DSCAM locus on chr21, while the ~80 TB of pVCFs download. -# Downloaded the pVCFs with Globus into pvcf/<chrom>/, 2.5 Mb chunks named -# SPARK.WGS.2026_08.gatk.<chrom>_<start>_<end>.vcf.gz, 60-85 GB each, no index. -# INFO has only AF/AQ; the samples are the sample_id column of the metadata. -# Chunks compress ~10:1 (1.8 MB of text per record), so bcftools parses ~150 MB/s -# of text, ~1.5 h per chunk. sparkWgs45kPvcfSlice.py bisects over the BGZF -# blocks to find a position without an index, so a range can be split into -# pieces that run in parallel. +# SFARI SPARK WGS August 2026 release (DS0000135), 45,178 genomes -> subtrack +# sfariSparkWgs45kAsd (doc page shared with sfariSparkExomes.html), with exact +# AC/AN/AF counted from the genotypes, ASD/non-ASD splits, no AC cutoff, refs #38424 +# Downloaded the GLnexus pVCFs with Globus into pvcf/<chrom>/: 1,248 chunks of +# 2.5 Mb (chr1-22, X, Y), SPARK.WGS.2026_08.gatk.<chrom>_<start>_<end>.vcf.gz, +# 75 TB, 60-85 GB each. INFO has only AF/AQ; the samples are the sample_id +# column of resources/SPARK.WGS.2026_08.sample_metadata.tsv. +# The chunks compress ~10:1 (1.8 MB of text per record), so each chunk is split +# into 500 kb pieces, one parasol job each. sparkWgs45kPvcfSlice.py bisects over +# the BGZF blocks to find a position without needing an index. +# Per piece, sparkWgs45kPvcfToSites.sh strips every sample column to GT (the +# other FORMAT fields are 90% of the text, and PL has 528 values per sample at a +# 31-allele STR such as chr1:10992130, which ran a job out of memory), then +# bcftools +fill-tags -S (AC/AN/AF overall and per group), view -G, norm -m- +# -f hg38.fa, drops AC=0 alleles. MONOALLELIC records (GLnexus alleles it +# could not merge into an overlapping site) have genotypes only for carriers, +# so fill-tags gives AN=1-2 and AF=1; for these AN is set to 2 x group size and +# AF recomputed (the old 12k WGS track has the same problem). VARLEN = +# len(ALT)-len(REF) is added, and the per-allele INFO fields are declared +# Number=1 (one ALT per record after norm -m-), as the VCF track filters need. +# The scripts use /cluster/software/src/bcftools-1.22: the conda bcftools in +# ~max/software links libopenblas, which segfaults under the parasol -ram +# address-space limit. cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs45k -mkdir -p pvcfSites/dscam +mkdir -p pvcfSites/run # sample groups from the asd column: 20,858 AUT, 24,320 NON_AUT, none unassigned awk -F'\t' 'NR==1{for(i=1;i<=NF;i++)c[$i]=i; next} {a=$c["asd"]; print $c["sample_id"]"\t"(a=="True"?"AUT":(a=="False"?"NON_AUT":"NA"))}' \ resources/SPARK.WGS.2026_08.sample_metadata.tsv | grep -v 'NA$' > pvcfSites/groups.txt -# DSCAM, knownGene ENST00000400454.6 chr21:40010999-40847158 (1-based), 64 pieces -bash ~/kent/src/hg/makeDb/scripts/varFreqs/sparkWgs45kPvcfRange.sh \ - pvcf/chr21/SPARK.WGS.2026_08.gatk.chr21_40000001_42500000.vcf.gz pvcfSites/groups.txt \ - 40010999 40847159 pvcfSites/dscam/sparkWgs45kAsd.dscam.vcf.gz 64 48 > pvcfSites/dscam/build.log 2>&1 -# 80 seconds on hgwdev. 127,615 input records -> norm split 12,625, realigned -# 11,261 -> 151,578 output records, of which 49,818 singletons (AC=1). No -# duplicate CHROM/POS/REF/ALT. -# For the SFARI demo the DSCAM file was also run through bcftools csq (same -# options and Ensembl 115 GFF3 as mergeAndAnnotate.sh): 464 missense, 78 -# frameshift, 6 stop_gained, 1,057 synonymous, 149,151 intronic. The genome-wide -# build below leaves BCSQ out, as for the other cohorts; it is only added when -# the cohorts are merged by mergeAndAnnotate.sh. -# MONOALLELIC records (GLnexus alleles it could not merge into an overlapping -# site): only carriers are called, so fill-tags gave AN=1-2 and AF=1 for them -# (4,117 records in DSCAM). The old 12k WGS track has the same problem. GLnexus's -# own AF uses all samples, so sparkWgs45kPvcfToSites.sh sets AN to 2 x group size -# for these and recomputes the AFs. 238 records still have AN < 45,178: real -# missing genotypes, e.g. in homopolymer runs. -# Compared with the AF-table track sfariSparkWgs45k in DSCAM: 99,147 alleles -# are in both, 6,985 of them differ in AF by more than 1e-4 (1,233 at sites with -# AN < 88,000). 2,835 alleles with AC >= 2 here are missing from the table track. -# The table's AF is GLnexus's INFO/AF, and it does not always match the -# genotypes: e.g. chr21:40011372 G>A has INFO AF=1.1e-05 (one allele) but two -# 0/1 calls, and 40011375 G>A has AF=0.000111 (10 alleles) but 22 0/1 calls. -# Counting the GT fields by hand confirms the fill-tags numbers. - -# Whole genome on parasol. The download finished: 1,248 chunks (chr1-22, X, Y), -# 75 TB, each ending in a valid BGZF EOF block, now also with .tbi files. Each -# chunk is split into 500 kb pieces, one job per piece (a 25 kb test piece: 43 s, -# ~2 cores, 0.6 GB peak RSS). parasol applies -ram as an address-space limit -# (ulimit -v and -d) to every process. The first try used the conda bcftools in -# ~max/software, which links libopenblas; OpenBLAS reserves virtual memory per -# node core, so every bcftools in the pipeline segfaulted at startup, even -# "view -h", at both -ram=2g and 4g. The scripts now use the plain -# /cluster/software/src/bcftools-1.22 build, which runs under 4g; its output for -# chr21:40300001-40300400 is identical to the demo build (records at the same -# position come out in a different order, which the merge sort fixes). -# The second try still crashed one job after 14 min at chr1:10992130, a -# 31-allele STR record with 95 MB of text: PL has 528 values per sample. Only -# GT is needed, so sparkWgs45kPvcfToSites.sh now strips every sample column to -# GT with perl before bcftools; that cuts the text ~10-fold, runs the STR -# region under a 4 GB limit, and gives identical output (25 kb test piece and -# DSCAM test region, same md5 of the sorted records). -cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs45k/pvcfSites -mkdir -p run +cd pvcfSites bash ~/kent/src/hg/makeDb/scripts/varFreqs/sparkWgs45kPvcfJobs.sh \ ../pvcf groups.txt pieces > run/jobList # 6,188 jobs cd run +# a 25 kb piece: 43 s, ~2 cores, 0.6 GB RSS; bcftools needs ~3 GB of address space para create -cpu=2 -ram=4g jobList para try para push -maxJob=200 # capped to spare /hive, 75 TB read in total # All 6,188 jobs ran OK (average 26 min, longest 5.6 h while reading /hive, # 29 h from submission to the last job with at most ~130 running). -# Pieces hold 436,319,161 input records, the same as the AF table: norm split -# 47,268,026, realigned 32,593,833, no REF mismatches, no duplicates. -cd /hive/data/genomes/hg38/bed/varFreqs/sparkWgs45k/pvcfSites +# The pieces hold 436,319,161 input records, the same number as SFARI's AF table +# (SPARK.WGS.2026_08.gatk.pvcf_variant_frequencies.tsv): norm split 47,268,026, +# realigned 32,593,833, no REF mismatches, no duplicates. +cd .. mkdir -p chroms ls pieces | sort -V | parallel -j 8 'bash ~/kent/src/hg/makeDb/scripts/varFreqs/sparkWgs45kPvcfMerge.sh chroms/{}.vcf.gz pieces/{}/*.bcf 2> chroms/{}.merge.log' /cluster/software/src/bcftools-1.22/bcftools concat --threads 8 -Oz -o sparkWgs45kAsd.vcf.gz \ $(ls chroms/*.vcf.gz | sort -V) tabix -p vcf sparkWgs45kAsd.vcf.gz # 13 GB +rm -rf pieces chroms # 518,290,049 records, 166,084,635 singletons (AC=1), 12,030,738 MONOALLELIC # records with corrected AN, 7,665,197 with AN < 45,178 (1,508,111 of them on # chrY, where only the 20,367 males are called; chrX is called diploid in males, # as in the older SPARK tracks). VARLEN -356 to 1,066. No adjacent duplicate -# CHROM/POS/REF/ALT. The DSCAM region is identical to the demo build without its -# BCSQ. APOE rs429358 AF=0.1456 (AUT 0.1472, NON_AUT 0.1442), rs7412 AF=0.0755. -# /gbdb/hg38/varFreqs/_sfari/sparkWgs45kAsd.vcf.gz now points to this file. -# The VCF track filters (Christopher's code, refs #37617) only accept INFO -# fields with Number=1, and fill-tags writes AC/AF and the per-group fields as -# Number=A. Every record has one ALT after norm -m-, so Number=1 is correct: -# sparkWgs45kPvcfToSites.sh now writes it, and the finished file was fixed with -B=/cluster/software/src/bcftools-1.22/bcftools -$B view -h sparkWgs45kAsd.vcf.gz | sed -E '/^##INFO=<ID=[^,]*,Number=A,/s/,Number=A,/,Number=1,/' > newHeader.txt -$B reheader -h newHeader.txt -o sparkWgs45kAsd.n1.vcf.gz sparkWgs45kAsd.vcf.gz -tabix -p vcf sparkWgs45kAsd.n1.vcf.gz -mv -f sparkWgs45kAsd.n1.vcf.gz sparkWgs45kAsd.vcf.gz -mv -f sparkWgs45kAsd.n1.vcf.gz.tbi sparkWgs45kAsd.vcf.gz.tbi -rm newHeader.txt -# Check: AF_AUT >= 0.001 at chr19:44908600-44908900 hides 122 of 125 variants, -# the same as bcftools view -i 'AF_AUT>=0.001'. +# CHROM/POS/REF/ALT. APOE rs429358 AF=0.1456 (AUT 0.1472, NON_AUT 0.1442), +# rs7412 AF=0.0755. AF_AUT >= 0.001 at chr19:44908600-44908900 hides 122 of 125 +# variants in the browser, the same as bcftools view -i 'AF_AUT>=0.001'. +# Note: GLnexus's INFO/AF, which SFARI's AF table is made from, does not always +# match the genotypes, e.g. chr21:40011372 G>A has AF=1.1e-05 (one allele) but +# two 0/1 calls. The counts here come from the genotypes.