# gnomaAD for hg38 # RM 21782 # Genome Aggregation Database at the Broad # https://gnomad.broadinstitute.org # (2019-10-07 kate) # March 6, 2019: gnomAD 2.1.1 mkdir vcf; cd vcf # hg38 mkdir hg38; cd hg38 mkdir -p /gbdb/hg38/gnomAD/vcf # Download from GRCh38 liftover section of downloads page: https://gnomad.broadinstitute.org/downloads#variants-grch38-liftover wget https://storage.googleapis.com/gnomad-public/release/2.1.1/README.txt # Exome variants # vcf # 85.31 GiB, MD5: cff8d0cfed50adc9211d1feaed2d4ca7 wget https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/exomes/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz >&! wget.exomes.log & # from log: 20 minutes download time # --2019-10-07 16:40:48-- https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/exomes/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz # 2019-10-07 17:00:52 (72.6 MB/s) - 'gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz' saved [91604348731/91604348731] du -sh *.bgz # 171G gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz md5sum gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz # cff8d0cfed50adc9211d1feaed2d4ca7 gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz ln -s `pwd`/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz /gbdb/hg38/gnomAD/vcf/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.gz # index file wget https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/exomes/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz.tbi ln -s `pwd`/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.bgz.tbi /gbdb/hg38/gnomAD/vcf/gnomad.exomes.r2.1.1.sites.liftover_grch38.vcf.gz.tbi # Genomes Variants # 743.06 GiB, MD5: 83de3d5b52669f714e810d4fcf047c18 wget https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/genomes/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz >&! wget.genomes.log & # from log: 3 hrs download time #--2019-10-07 17:03:11-- https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/genomes/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz # 2019-10-07 20:02:56 (70.6 MB/s) - 'gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz' saved [797851733355/797851733355] ln -s `pwd`/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz /gbdb/hg38/gnomAD/vcf/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.gz wget https://storage.googleapis.com/gnomad-public/release/2.1.1/liftover_grch38/vcf/genomes/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz.tbi ln -s `pwd`/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.bgz.tbi /gbdb/hg38/gnomAD/vcf/gnomad.genomes.r2.1.1.sites.liftover_grch38.vcf.gz.tbi ################# Use/adapt trackDb and HTML from hg19 track ################# # Add V3. Native hg38 analysis (not lifted) # Released 10/16/2019 by MacArthur lab, announced here: # https://macarthurlab.org/2019/10/16/gnomad-v3-0/ # (2019-10-21 kate) mkdir /hive/data/outside/gnomAD.3 cd /hive/data/outside/gnomAD.3 wget https://storage.googleapis.com/gnomad-public/release/3.0/vcf/genomes/gnomad.genomes.r3.0.sites.vcf.bgz.tbi ln -s `pwd`/gnomad.genomes.r3.0.sites.vcf.bgz.tbi /gbdb/hg38/gnomAD/vcf/gnomad.genomes.r3.0.sites.vcf.gz.tbi (date; wget https://storage.googleapis.com/gnomad-public/release/3.0/vcf/genomes/gnomad.genomes.r3.0.sites.vcf.bgz; date) >&! wget.genomes.log & #Mon Oct 21 11:26:58 PDT 2019 #Mon Oct 21 12:27:14 PDT 2019 # ~1 hr md5sum gnomad.genomes.r3.0.sites.vcf.bgz ln -s `pwd`/gnomad.genomes.r3.0.sites.vcf.bgz /gbdb/hg38/gnomAD/vcf/gnomad.genomes.r3.0.sites.vcf.gz ############################################################################## # gnomAD v3.1 Update - Nov 24, 2020 - ChrisL - DONE # See /hive/data/outside/gnomAD.3/v3.1/make.txt for how to download # the new version WORKDIR=/hive/data/inside/gnomAD/v3.1 cd $WORKDIR db="hg38" mkdir -p $db/genomes # get the headers: bcftools view -h /hive/data/outside/gnomAD.3/v3.1/genomes/gnomad.genomes.v3.1.sites.chr1.vcf.bgz > gnomad.v3.1.header # The VEP field has the important information, format: bcftools view -h /hive/data/outside/gnomAD.3/v3.1/genomes/gnomad.genomes.v3.1.sites.chr1.vcf.bgz | grep vep | grep -o "Format: .*$" | cut -d' ' -f2- | tr '|' '\t' | tl > gnomad.v3.1.vep.fields # the specific fields (same as v2.1.1 except '_' has been swapped to '-' EXCEPT in popmax # fields -___-), also add some new fields like CADD scores: fields="AC,AN,AF,faf95,nhomalt,vep,popmax,AC_popmax,AN_popmax,AF_popmax,AC-afr,AN-afr,AF-afr,nhomalt-afr,AC-ami,AN-ami,AF-ami,nhomalt-ami,AC-amr,AN-amr,AF-amr,nhomalt-amr,AC-asj,AN-asj,AF-asj,nhomalt-asj,AC-eas,AN-eas,AF-eas,nhomalt-eas,AC-fin,AN-fin,AF-fin,nhomalt-fin,AC-mid,AN-mid,AF-mid,nhomalt-mid,AC-nfe,AN-nfe,AF-nfe,nhomalt-nfe,AC-sas,AN-sas,AF-sas,nhomalt-sas,AC-oth,AN-oth,AF-oth,nhomalt-oth,cadd_phred,revel_score,splice_ai_max_ds,splice_ai_consequence,primate_ai_score" # Don't use parallel anymore, use parasol cause this still takes far too long and bogs down hgwdev: # make script to run each job: cat <<'_EOF' > run.sh #!/bin/bash # don't do the following, as vcfToBed will correctly exit with an error #set -beEu -o pipefail fields="AC,AN,AF,faf95,nhomalt,vep,popmax,AC_popmax,AN_popmax,AF_popmax,AC-afr,AN-afr,AF-afr,nhomalt-afr,AC-ami,AN-ami,AF-ami,nhomalt-ami,AC-amr,AN-amr,AF-amr,nhomalt-amr,AC-asj,AN-asj,AF-asj,nhomalt-asj,AC-eas,AN-eas,AF-eas,nhomalt-eas,AC-fin,AN-fin,AF-fin,nhomalt-fin,AC-mid,AN-mid,AF-mid,nhomalt-mid,AC-nfe,AN-nfe,AF-nfe,nhomalt-nfe,AC-sas,AN-sas,AF-sas,nhomalt-sas,AC-oth,AN-oth,AF-oth,nhomalt-oth,cadd_phred,revel_score,splice_ai_max_ds,splice_ai_consequence,primate_ai_score" infile=$1 finalBed=$2 # get the suffix name after the final '/': suffix=${infile#/hive/data/outside/gnomAD.3/v3.1/genomes/} suffix=${suffix%.vcf.bgz} outBed=/hive/data/inside/gnomAD/v3.1/hg38/genomes/${suffix}.bed vcfToBed -fields="${fields}" ${infile} ${outBed} 2>/dev/null /hive/data/inside/gnomAD/v3.1/gnomadVcfBedToBigBed ${outBed} stdout | sort -k1,1 -k2,2n > ${finalBed} _EOF # make a jobList file, here's a sample line: ./run.sh /hive/data/outside/gnomAD.3/v3.1/genomes/gnomad.genomes.v3.1.sites.chr1.vcf.bgz {check out exists+ /hive/data/inside/gnomAD/v3.1/hg38/genomes/chr1.bed} ssh ku "cd ${WORKDIR}; para create jobList; para try; exit" # woops errors in the gnomadVcfBedToBigBed script caused the first 10 jobs fail, fix up the script # and run those jobs separately: # the regular expression is for extracting chr*.bed out of the input file and making just that name # the output file # first start the rest of the jobs: ssh ku "cd ${WORKDIR}; para push; exit" # then finish off: parallel -j10 --joblog secondhalf.run.log --plus "./gnomadVcfBedToBigBed {} hg38/genomes/{= s/.*(chr[0-9]+\.bed$)/\$1/ =}" ::: hg38/genomes/gnomad.genomes.v3.1.sites.chr{1,10,11,12,13,14,15,16,17,18}.bed # See para.time and secondhalf.run.log for timing info time sort -k1,1 -k2,2n hg38/genomes/chr*.bed > gnomad.v3.1.genomes.bed # real 26m7.992s # user 31m40.053s # sys 14m31.118s nohup time bedToBigBed -type=bed9+64 -tab -as=genomes.as gnomad.v3.1.genomes.bed /hive/data/genomes/hg38/chrom.sizes genomes.bb &> bigBed.log & mkdir -p /gbdb/hg38/gnomAD/v3.1/variants/ ln -s `pwd`/genomes.bb /gbdb/hg38/gnomAD/v3.1/variants/ ############################################################################## # gnomAD v3.1.1 Update - Jun 29, 2021 - ChrisL - DONE # See /hive/data/outside/gnomAD.3/v3.1.1/make.txt for how to download # the new version ############################################################################## WORKDIR=/hive/data/inside/gnomAD/v3.1.1 cd $WORKDIR db="hg38" mkdir -p $db/genomes # get the headers: bcftools view -h /hive/data/outside/gnomAD.3/v3.1.1/genomes/gnomad.genomes.v3.1.1.sites.chr1.vcf.bgz > gnomad.v3.1.1.header # The VEP field has most of the important information, format: bcftools view -h /hive/data/outside/gnomAD.3/v3.1.1/genomes/gnomad.genomes.v3.1.1.sites.chr1.vcf.bgz | grep vep | grep -o "Format: .*$" | cut -d' ' -f2- | tr '|' '\t' | tl > gnomad.v3.1.1.vep.fields # Now make two control scripts for specifying fields, see # /hive/data/inside/gnomAD/v3.1.1/runVcfToBed.sh and # /hive/data/inside/gnomAD/v3.1.1/convertBeds.sh # for more information # First convert the VCFs to beds: # Need a jobList for parallel to run: today=`date +%F` for c in {1..22} X Y M; do echo "./runVcfToBed.sh /hive/data/outside/gnomAD.3/v3.1.1/genomes/gnomad.genomes.v3.1.1.sites.chr$c.vcf.bgz $WORKDIR/$db/genomes/gnomad.v3.1.1.chr$c.bed" done > vcfToBed.jobList parallel --joblog vcfToBed.run.log --jobs 10 :::: vcfToBed.jobList # WOOPS wrong path to chrM, run separately: ./runVcfToBed.sh /hive/data/outside/gnomAD.3/v3.1.1/genomes/gnomad.genomes.v3.1.sites.chrM.vcf.bgz /hive/data/inside/gnomAD/v3.1.1/hg38/genomes/gnomad.v3.1.chrM.bed # uh oh chrM went wrong again. The chrM variant calling pipeline is different than # the standard variant calling pipeline gnomAD used for the rest of the genome. Thus # we need to get specific fields out of the chrM file, and append empty fields # into the other tab files (the chrM specific data will be missing for the other # chromosomes while the regular genome data will be mostly missing for the chrM # data). This also brings up a bug in the original where I was missing Total counts # for the population data, so just do a full re-run which unifies the tab files # for both chrM and the rest of the genome: ./runVcfToBed.chrM.sh /hive/data/outside/gnomAD.3/v3.1.1/genomes/gnomad.genomes.v3.1.sites.chrM.vcf.bgz /hive/data/inside/gnomAD/v3.1.1/hg38/genomes/gnomad.v3.1.chrM.bed for c in {1..22} X Y; do echo "./convertBeds.sh $WORKDIR/$db/genomes/gnomad.v3.1.1.chr$c.bed $WORKDIR/$db/genomes/chr$c.bed" done > convertBed.jobList echo "./convertBeds.chrM.sh $WORKDIR/$db/genomes/gnomad.v3.1.chrM.bed $WORKDIR/$db/genomes/chrM.bed" >> 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: time bgzip -iI gnomad.v3.1.1.details.tab.gz.gzi -c gnomad.v3.1.1.details.pre.tab > gnomad.v3.1.1.details.tab.gz & # real 255m53.801s # user 238m3.982s # sys 16m7.130s time (ls -1S hg38/genomes/chr*.bed | parallel --joblog join.run.log --max-procs 12 \\ bedJoinTabOffset -verbose=2 -bedKey=3 gnomad.v3.1.1.details.pre.tab {} joined/{/}) &> bedJoinTabOffset.log & time sort --merge joined/*.bed | grep -v "^#" > gnomad.v3.1.1.genomes.joined.bed # real 16m4.876s # user 12m2.699s # sys 6m30.244s # and lastly turn the merged bed into a bigBed: time (bedToBigBed -type=bed9+15 -tab -as=genomes.as -extraIndex=name,rsId,_displayName gnomad.v3.1.1.genomes.joined.bed /hive/data/genomes/hg38/chrom.sizes genomes.bb) &> bigBed.log & # Max Tue Apr 5 04:14:40 PDT 2022 - adding mutation constraint subtrack cd /hive/data/genomes/hg38/bed/gnomad/constraint # Downloaded BED file from https://www.biorxiv.org/content/10.1101/2022.03.20.485034v1.supplementary-material # got chrX scores from konrad.j.karczewski@gmail.com because they feel unsure about them cat Supplementary_Data_2.bed constraint_z_genome_1kb_chrx.bed > mutConstraint.bed bedSort mutConstraint.bed mutConstraint.bed bedGraphToBigWig *.bed ../../../chrom.sizes mutConstraint.bw bedGraphToBigWig *.bed ../../../chrom.sizes mutConstraint.bw ############################################################################## # gnomAD v4 - Nov 1, 2023 - ChrisL - DONE ############################################################################## # See /hive/data/outside/gnomAD.4/make.txt for how to download # the new version # Make the work dir: mkdir -p /hive/data/inside/gnomAD/v4 cd /hive/data/inside/gnomAD/v4 # for each vcf in the download dirs, install files into /gbdb and load up a table with the pointers gbdbPath="/gbdb/hg38/gnomAD/v4/" dataPath="/hive/data/outside/gnomAD.4/" for tbl in exomes genomes do mkdir -p ${gbdbPath}${tbl} ln -s ${dataPath}${tbl}/*.vcf.bgz* ${gbdbPath}${tbl} cp /dev/null ${tbl}.txt if [ ${tbl} == "exomes" ]; then for c in {1..22} X Y do f=${gbdbPath}${tbl}/gnomad.exomes.v4.0.sites.chr${c}.vcf.bgz echo -e "${f}\tchr${c}" >> ${tbl}.txt done else for c in {1..22} X Y do f=${gbdbPath}${tbl}/gnomad.genomes.v4.0.sites.chr${c}.vcf.bgz echo -e "${f}\tchr${c}" >> ${tbl}.txt done fi tName="${tbl^}" hgLoadSqlTab hg38 gnomad${tName}V4 ~/kent/src/hg/lib/bbiChroms.sql genomes.txt done ############################################################################## # Add cancer/non-cancer filter to gnomAD v3.1.1 - Feb 16, 2024 - ChrisL - DONE ############################################################################## # see /hive/data/inside/gnomAD/v3.1.1/2024-02-12/README # for the steps ############################################################################## ############################################################################## # Update transcript contraint scores for gnomAD v4 - Mar 1, 2024 - ChrisL - DONE ############################################################################## # The data is not the same as the hg19 version, many less fields, but we mostly threw a lot of them # away before anyways. Also this file had both ENST and NM/XM transcripts # Here are the fields used in the previous version: zcat ../../gnomAD.2/constraint/gnomad.v2.1.1.lof_metrics.by_transcript.txt.bgz | head -1 | tawk '{print $76,$77,$78,$65,$66,$1,$2,$4,$5,$6,$34,$13,$14,$15,$33,$18,$21,$22,$25,$26,$27,$28,$29,$30,$31}' chromosome start_position end_position gene_id transcript_level gene transcript obs_mis exp_mis oe_mis mis_z obs_syn exp_syn oe_syn syn_z obs_lof exp_lof pLI oe_lof oe_syn_lower oe_syn_upper oe_mis_lower oe_mis_upper oe_lof_lower oe_lof_upper # Right off the bat we will be missing the chrom,start and end # Get the transcripts to get the coordinates and exon-intron boundaries hgsql -Ne "select * from wgEncodeGencodeCompV39" hg38 \ | cut -f2- | genePredToBed -tab stdin stdout | sed -Ee 's/\.[0-9]+//' \ | sort -k4 > hg38.gencodeCompV39.bed12 hgsql -Ne "select * from ncbiRefSeq" hg38 \ | cut -f2- | genePredToBed -tab stdin stdout \ | sort -k4 > hg38.refSeq.bed12 cat hg38.gencodeCompV39.bed12 hg38.refSeq.bed12 | sort -k4 > transcripts.coords f=gnomad.v4.0.constraint_metrics.tsv # I don't think this command will work just copying and pasting like the above will 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