dde908e5381a6da2e78aea5c8ae38c6dcff2cd46 chmalee Tue Feb 6 09:46:42 2024 -0800 Staging gnomAD v4, including reorganizing tracks so v4.1 bigBed are first after the v4 VCF, description pages changes and makedoc, refs #33746 diff --git src/hg/makeDb/doc/hg38/gnomad.txt src/hg/makeDb/doc/hg38/gnomad.txt index 4f33ae6..824880f 100644 --- src/hg/makeDb/doc/hg38/gnomad.txt +++ src/hg/makeDb/doc/hg38/gnomad.txt @@ -291,15 +291,204 @@ tail -n +2 $f \ | tawk '{print $1,$2,$24,$25,$27,$32,$37,$38,$40,$45,$12,$13,$17,$15,$42,$43,$29,$30,$20,$21}' \ | sort -k2 | join -t $'\t' -1 4 -2 2 transcripts.coords - \ | tawk '{for (i=1; i<=12; i++) {printf "%s\t", $i} printf "%s\t%s\t%s\t\t\t", $2, $3, $4; for (i=13; i <= NF; i++) {printf "%s", $i; if (i != NF) {printf "\t"}}; printf "\n"} ' \ | ~/kent/src/hg/makeDb/gnomad/combine.awk -v doTranscripts=true bedSort pliByTranscript.tab pliByTranscript.tab bedSort missenseByTranscript.tab missenseByTranscript.tab # Copy the old autosql file: cp ../../gnomAD.2/constraint/{missense,pli}Metrics.as . # Turn into a bigBed and link sizes=/hive/data/genomes/hg38/chrom.sizes bedToBigBed -type=bed12+6 -as=pliMetrics.as -tab -extraIndex=name,geneName pliByTranscript.tab $sizes pliByTranscript.bb bedToBigBed -type=bed12+5 -as=missenseMetrics.as -tab -extraIndex=name,geneName missenseByTranscript.tab $sizes missenseByTranscript.bb + +############################################################################## +# gnomAD v4.1 Update to bigBed - May 02, 2024 - ChrisL +############################################################################## +# start a screen as this will take a while +screen -S gnomADv4 +WORKDIR=/hive/data/inside/gnomAD/v4/v4.1 +cd $WORKDIR +db="hg38" +mkdir -p $db/{exomes,genomes} + +# get the headers: +bcftools view -h /hive/data/outside/gnomAD.4/v4.1/genomes/gnomad.genomes.v4.1.sites.chr1.vcf.bgz > gnomad.v4.1.genomes.header +bcftools view -h /hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.chr1.vcf.bgz > gnomad.v4.1.exomes.header + +# Looks like there are INFO field differences between the two files, for some reason +# there are different populations included for each? +# The Genomes data has Amish population information, while the exomes data has UK BioBank +# data and faf95_mid? +# For the full list of differences, use this command: +sdiff <(grep -Eo "^##.*=<[^,]+," gnomad.v4.1.genomes.header | sort) <(grep -Eo "^##.*=<[^,]+," gnomad.v4.1.exomes.header | sort ) | less +# List out the INFO fields for each set: +grep -Eo "^##INFO=<ID=[^,]+" gnomad.v4.1.exomes.header | cut -d'=' -f3 > gnomad.v4.1.exomes.infoFields +grep -Eo "^##INFO=<ID=[^,]+" gnomad.v4.1.genomes.header | cut -d'=' -f3 > gnomad.v4.1.genomes.infoFields + +# The VEP field has most of the important information, format: +bcftools view -h /hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.chr1.vcf.bgz | grep vep | grep -o "Format: .*$" | cut -d' ' -f2- | tr '|' '\t' | tl > gnomad.exomes.v4.1.vep.fields +bcftools view -h /hive/data/outside/gnomAD.4/v4.1/genomes/gnomad.genomes.v4.1.sites.chr1.vcf.bgz | grep vep | grep -o "Format: .*$" | cut -d' ' -f2- | tr '|' '\t' | tl > gnomad.genomes.v4.1.vep.fields +# at least the VEP fields are the same + +# Use the fields from v3.1.1 as a starting point: +grep "^fields" ../../v3.1.1/runVcfToBed.sh | cut -d'=' -f2 | tr ',' '\n' | tr -d '"' > v3.1.1.vcfToBed.fieldList +# check at least those fields are present in both sets: +for line in $(cat v3.1.1.vcfToBed.fieldList); do grep $line gnomad.v4.1.exomes.infoFields &> /dev/null; if [ $? != 0 ]; then echo "$line not in v4.1 exomes"; fi; done +popmax not in v4.1 exomes +AC_popmax not in v4.1 exomes +AN_popmax not in v4.1 exomes +AF_popmax not in v4.1 exomes +AC_ami not in v4.1 exomes +AN_ami not in v4.1 exomes +AF_ami not in v4.1 exomes +nhomalt_ami not in v4.1 exomes +AC_oth not in v4.1 exomes +AN_oth not in v4.1 exomes +AF_oth not in v4.1 exomes +nhomalt_oth not in v4.1 exomes +revel_score not in v4.1 exomes +splice_ai_max_ds not in v4.1 exomes +splice_ai_consequence not in v4.1 exomes +primate_ai_score not in v4.1 exomes + +for line in $(cat v3.1.1.vcfToBed.fieldList); do grep $line gnomad.v4.1.genomes.infoFields &> /dev/null; if [ $? != 0 ]; then echo "$line not in v4.1 genomes"; fi; done + +# So they changed popmax to grpmax, and spliceAI and revel scores: +grep -i "revel\|splice\|primate" gnomad.v4.1.exomes.infoFields +revel_max +spliceai_ds_max +grep -i "revel\|splice\|primate" gnomad.v4.1.genomes.infoFields +revel_max +spliceai_ds_max +cp v3.1.1.vcfToBed.fieldList v4.1.vcfToBed.exomes.fieldList +cp v3.1.1.vcfToBed.fieldList v4.1.vcfToBed.genomes.fieldList +# Make the fixups in vim, including the population information changes: remove ami from exomes, oth to remaining and the revel and splice_ai field changes, add variant_type +# I used this command fixing up until there was nothing only on the left +# and no '|' on the right +sdiff <(sort v4.1.vcfToBed.exomes.fieldList) <(sort gnomad.v4.1.exomes.infoFields) | grep -v "ukb" | less + +# Now make two control scripts for specifying fields, see +# /hive/data/inside/gnomAD/v4/v4.1/runVcfToBed.{exomes,genomes}.sh and +# /hive/data/inside/gnomAD/v4/v4.1/convertBeds.{exomes,genomes}.sh +# for more information + +# First convert the VCFs to beds: +# Need a jobList for parallel to run: +for c in {1..22} X Y; do + echo "./runVcfToBed.exomes.sh /hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.chr$c.vcf.bgz $WORKDIR/$db/exomes/gnomad.v4.1.chr$c.bed" +done > vcfToBed.jobList +for c in {1..22} X Y; do + echo "./runVcfToBed.genomes.sh /hive/data/outside/gnomAD.4/v4.1/genomes/gnomad.genomes.v4.1.sites.chr$c.vcf.bgz $WORKDIR/$db/genomes/gnomad.v4.1.chr$c.bed" +done >> vcfToBed.jobList +parallel --joblog vcfToBed.run.log --jobs 15 :::: vcfToBed.jobList + +# Average of 105 minutes per job: +tail -n +2 /hive/data/inside/gnomAD/v4/v4.1/vcfToBed.run.log | cut -f4 | awk '{sum += $1} END {print (sum/NR) / 60}' +54.0564 + +# convert the beds into a smaller bed file + a tab sep file with everything else +for c in {1..22} X Y; do + echo "./convertBeds.exomes.sh $WORKDIR/$db/exomes/gnomad.v4.1.chr$c.bed $WORKDIR/$db/exomes/chr$c.bed" +done > convertBed.jobList +for c in {1..22} X Y; do + echo "./convertBeds.genomes.sh $WORKDIR/$db/genomes/gnomad.v4.1.chr$c.bed $WORKDIR/$db/genomes/chr$c.bed" +done >> convertBed.jobList +time (parallel --joblog convertBed.run.log --jobs 25 :::: convertBed.jobList) &> convertBed.log & + +# Now we have one bed file and one external file per chromosome, use bgzip and +# bedJoinTabOffset to combine them together: +# Merge all the tab files together: +mkdir -p joined/{exomes,genomes} +# First exomes: +time sort --merge ${WORKDIR}/hg38/exomes/gnomad.*.tab > gnomad.v4.1.exomes.details.pre.tab +# real 14m24.203s +# user 2m17.384s +# sys 10m16.118s + +time bgzip -iI gnomad.v4.1.exomes.details.tab.gz.gzi -c gnomad.v4.1.exomes.details.pre.tab > gnomad.v4.1.exomes.details.tab.gz +# real 62m28.397s +# user 58m29.676s +# sys 3m31.012s + +# The parallel command changes {} into the input file (example: hg38/exomes/chr1.bed), and +# {/} to the basename (example: chr1.bed) +time (ls -1S hg38/exomes/chr*.bed | + parallel --joblog join.run.log --max-procs 12 \\ + bedJoinTabOffset -verbose=2 -bedKey=3 gnomad.v4.1.exomes.details.pre.tab {} joined/exomes/{/}) &> bedJoinTabOffset.log +# real 155m42.266s +# user 381m15.009s +# sys 50m29.407s +time sort -S 40G --merge joined/exomes/*.bed | grep -v "^#" > gnomad.v4.1.exomes.joined.bed +# real 6m38.608s +# user 3m49.595s +# sys 2m43.676s + +# Then genomes: +time sort --merge ${WORKDIR}/hg38/genomes/gnomad.*.tab > gnomad.v4.1.genomes.details.pre.tab +# real 48m36.383s +# user 7m40.632s +# sys 34m7.595s + +# FIELD_FIX_RESTART_MARK +time bgzip -iI gnomad.v4.1.genomes.details.tab.gz.gzi -c gnomad.v4.1.genomes.details.pre.tab > gnomad.v4.1.genomes.details.tab.gz +# real 230m19.470s +# user 216m16.280s +# sys 11m28.438s + +# The parallel command changes {} into the input file (example: hg38/genomes/chr1.bed), and +# {/} to the basename (example: chr1.bed) +time (ls -1S hg38/genomes/chr*.bed | + parallel --joblog join.run.log --max-procs 12 \\ + bedJoinTabOffset -verbose=2 -bedKey=3 gnomad.v4.1.genomes.details.pre.tab {} joined/genomes/{/}) &> bedJoinTabOffset.log +# real 430m5.112s +# user 1505m15.341s +# sys 286m35.990s + +time sort -S 40G --merge joined/genomes/*.bed | grep -v "^#" > gnomad.v4.1.genomes.joined.bed +# real 26m31.856s +# user 14m24.696s +# sys 9m57.459s + +# and lastly turn the merged beds into bigBeds: +time (bedToBigBed -type=bed9+16 -tab -as=exomes.as -extraIndex=name,rsId,_displayName gnomad.v4.1.exomes.joined.bed /hive/data/genomes/hg38/chrom.sizes exomes.bb) &> bigBed.exomeslog & +# pass1 - making usageList (24 chroms): 239308 millis +# pass2 - checking and writing primary data (183717261 records, 25 fields): 1695174 millis +# Sorting and writing extra index 0: 258529 millis +# Sorting and writing extra index 1: 172276 millis +# Sorting and writing extra index 2: 73698 millis + +# real 43m44.216s +# user 40m4.371s +# sys 1m54.306s + +time (bedToBigBed -maxAlloc=250000000000 -type=bed9+16 -tab -as=genomes.as -extraIndex=name,rsId,_displayName gnomad.v4.1.genomes.joined.bed /hive/data/genomes/hg38/chrom.sizes genomes.bb) &> bigBed.genomeslog & +# pass1 - making usageList (24 chroms): 720685 millis +# pass2 - checking and writing primary data (759336320 records, 25 fields): 8127679 millis +# Sorting and writing extra index 0: 1711243 millis +# Sorting and writing extra index 1: 1627587 millis +# Sorting and writing extra index 2: 436891 millis +# +# real 227m12.682s +# user 204m24.791s +# sys 20m16.421s + +# Make symlinks to /gbdb: +cd /gbdb/hg38/gnomAD/ +mkdir v4.1 +cd v4.1 +mkdir exomes +mkdir genomes +cd $WORKDIR +ln -s `pwd`/genomes.bb /gbdb/hg38/gnomAD/v4.1/genomes/ +ln -s `pwd`/gnomad.v4.1.genomes.details.tab.gz /gbdb/hg38/gnomAD/v4.1/genomes/ +ln -s `pwd`/gnomad.v4.1.genomes.details.tab.gz.gzi /gbdb/hg38/gnomAD/v4.1/genomes/ +ln -s `pwd`/exomes.bb /gbdb/hg38/gnomAD/v4.1/exomes/ +ln -s `pwd`/gnomad.v4.1.exomes.details.tab.gz /gbdb/hg38/gnomAD/v4.1/exomes/ +ln -s `pwd`/gnomad.v4.1.exomes.details.tab.gz.gzi /gbdb/hg38/gnomAD/v4.1/exomes/ + +# Woops the conversion script had the wrong header listing, change it up: +sed -i -e '1,24s/segdup/segdup\tvariant_type\tallele_type/' gnomad.v4.1.genomes.details.pre.tab +# and then re-do the genomes steps from the FIELD_FIX_RESTART_MARK mark