b85941107d295df78e38d8020dffa07896fe7cbb kate Tue Dec 22 13:07:58 2020 -0800 Adding make doc for R4 track. refs #26616 diff --git src/hg/makeDb/doc/covid/covidHgiGwas.txt src/hg/makeDb/doc/covid/covidHgiGwas.txt index a9da78a..9e010aa 100644 --- src/hg/makeDb/doc/covid/covidHgiGwas.txt +++ src/hg/makeDb/doc/covid/covidHgiGwas.txt @@ -1,173 +1,369 @@ # GWAS from the COVID-19 Host Genetics Initiative (HGI) # (2020-07-02 kate) # From: covid19hg.org/results # Contacts: Rachel Liao, Juha Karjalainen (Broad) juha.karjalainen@helsinki.fi # Create build dir cd /hive/data/outside/covid/covidHostGenetics # GWAS meta-analyses file format 1 #CHR chromosome 2 POS chromosome position in build 37 3 REF non-effect allele 4 ALT effect allele (beta is for this allele) 5 SNP #CHR:POS:REF:ALT {STUDY}_AF_Allele2 allele frequency in {STUDY} -{STUDY}_AF_fc allele frequency in {STUDY} / allele frequency in gnomAD v3 (1000000 if frequency in gnomAD is 0). Calculated based on each study's ancestry in gnomAD +{STUDY}_AF_fc allele frequency in {SeUDY} / allele frequency in gnomAD v3 (1000000 if frequency in gnomAD is 0). Calculated based on each study'ss ancestry in gnomAD {STUDY}_N 6 + (X = #studies * 3) all_meta_N number of studies that had the variant after AF and INFO filtering and as such were used for the meta 7 + X all_inv_var_meta_beta effect size on log(OR) scale 8 + X all_inv_var_meta_sebeta standard error of effect size 9 + X all_inv_var_meta_p p-value 10 + X all_inv_var_het_p p-value from Cochran's Q heterogeneity test # additional columns: 11 + X "all_meta_sample_N" 12 + X "all_meta_AF" 13 + X "rsid" # additional for hg19 liftover. Values in hg38. 14 + X "anew_chr" 15 + X "anew_pos" 16 + X "REF.1" 17 + X "ALT.1" # Studies 1 Genetic determinants of COVID-19 complications in the Brazilian population BRACOVID 2 Genetic modifiers for COVID-19 related illness BelCovid 3 deCODE DECODE 4 FinnGen FinnGen 5 GEN-COVID, reCOVID GENCOVID 6 UK 100,000 Genomes Project genomicsengland100kgp_EUR 7 Genes & Health GNH 8 Generation Scotland GS 9 COVID19-Host(a)ge HOSTAGE 10 Helix Exome+ COVID-19 Phenotypes Helix 11 UK Blood Donors Cohort INTERVAL 12 LifeLines CytoSNP LifelinesCyto 13 LifeLines Global Screening Array LifelinesGsa 14 Netherlands Twin Register NTR 15 Partners Healthcare Biobank PHBB 16 Qatar Genome Program QGP 17 UK Biobank UKBB ##################### # Consult with Ana on Aug 12 1. Restrict to analyses with enough power (based on plots): B2: hospitalized covid vs. population (3199 cases, 8 studies) C2: covid vs. population (6696 cases, 18 studies) 2. Label options: rsID and/or ALT/REF (use ALT/REF if no label) 3. Mouseover: pValue, effect size, #studies 4. Filters: pValue (default=5), #studies, ###################### # Make with new .as (input from Juha), and hg38 tracks # (2020-09-04 kate) # Analysis B2: hospitalized covid vs. population # cases: 3199 # studies: 8 # Analysis C2: covid vs. population # cases: 6696 # studies: 18 # download hg19 files wget https://storage.googleapis.com/covid19-hg-public/20200619/results/build_37/COVID19_HGI_ANA_B2_V2_20200701.b37.txt.gz wget https://storage.googleapis.com/covid19-hg-public/20200619/results/build_37/COVID19_HGI_ANA_C2_V2_20200701.b37.txt.gz # download hg38 files wget https://storage.googleapis.com/covid19-hg-public/20200619/results/COVID19_HGI_ANA_B2_V2_20200701.txt.gz wget https://storage.googleapis.com/covid19-hg-public/20200619/results/COVID19_HGI_ANA_C2_V2_20200701.txt.gz gunzip *.z wc -l covidHgiGwas*.hg38.txt 15392647 covidHgiGwas.B2.hg38.txt 24600933 covidHgiGwas.C2.hg38.txt # rename ln -s COVID19_HGI_ANA_B2_V2_20200701.txt covidHgiGwas.B2.hg38.txt ln -s COVID19_HGI_ANA_C2_V2_20200701.txt covidHgiGwas.C2.hg38.txt ln -s COVID19_HGI_ANA_B2_V2_20200701.b37.txt covidHgiGwas.B2.hg19.txt ln -s COVID19_HGI_ANA_C2_V2_20200701.b37.txt covidHgiGwas.C2.hg19.txt wc -l covidHgiGwas* 15392647 covidHgiGwas.B2.hg38.txt 24600933 covidHgiGwas.C2.hg38.txt mkdir -p /gbdb/hg19/covidHgiGwas /gbdb/hg38/covidHgiGwas cat > makeHgi.csh << 'EOF' set bin = ~/kent/src/hg/makeDb/outside/covid foreach d (B2.hg19 B2.hg38 C2.hg19 C2.hg38) set db = $d:e set a = $d:r set sizes = /hive/data/genomes/$db/chrom.sizes set b = covidHgiGwas set f = $b.$a.$db set bb = $b$a.$db echo $b.$a.$db if ($a == "B2") then set studies = 8 else set studies = 18 endif perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt > $f.bed bedSort $f.bed $f.sorted.bed bedToBigBed -type=bed9+12 -as=$bin/covidHgiGwas.as -tab $f.sorted.bed $sizes $bb.bb ln -s `pwd`/$bb.bb /gbdb/hg19/covidHgiGwas end 'EOF' make makeHgi.csh >&! makeHgi.out & ################## # Get effect size distribution to inform color gradient selection % textHistogram -real -minVal=-25 -binSize=1 -maxBinCount=50 effectSizes.txt -13.000000 1 -12.000000 1 -11.000000 0 -10.000000 2 -9.000000 3 -8.000000 6 -7.000000 6 -6.000000 * 19 -5.000000 *** 53 -4.000000 ********************************************************* 944 -3.000000 ************************************************************ 1000 -2.000000 ************************************************************ 1000 -1.000000 ************************************************************ 1000 0.000000 ************************************************************ 1001 1.000000 ************************************************************ 1000 2.000000 ************************************************************ 1000 3.000000 ************************************************************ 1000 4.000000 ************************************************************ 994 5.000000 ********************************************************** 971 6.000000 **************************************************** 866 7.000000 ******************************************* 717 8.000000 ****************************** 504 9.000000 ********************** 371 10.000000 *********** 183 11.000000 ***** 78 12.000000 *** 43 13.000000 ** 37 14.000000 * 13 15.000000 4 16.000000 6 17.000000 2 18.000000 0 19.000000 0 20.000000 1 + +################################## +# Release 4 (October 2020) +# Dec. 2020 Kate + +#23 contributors world-wide, including 23 and me +#14 European contributors in 23andme leave-out analysis + +# download, all leaving out 23andme (to allow full download) +# using 4 analyses this time, as per Ana recommendation (enough power now in A2 and C1) + +# Note the analysis methods have changed slightly. Effect size ranges are much smaller. + +B2: Hospitalized Covid vs population +Cases: 8638 - 23andme (613+140) = 7885 +Controls: 1736547 - 23andme (680416+94327) = 961804 + +C2: Covid vs population +Cases: 30937 - 23andme (506+9913+2553) = 17965 +Controls: 1471815 - 23andme (3110+85072+13086) = + +A2: Very severe respratory confirmed covid vs population +Cases: 4933 - 23andme (495 + 102) = 4336 +Controls: 1398672 (680440 + 94330) = 623902 + +C1: Covid vs. lab/self-reported negative +Cases: 24057 - 23andme (506 + 9913 + 2553) = 11085 +Controls: 218062 - 23andme (85072 + 13086 + 3110) = 116794 +Studies: 23 - 23andme (3) = 20 + +# hg38 +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_B2_ALL_leave_23andme_20201020.txt.gz +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_C2_ALL_leave_23andme_20201020.txt.gz + +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_A2_ALL_leave_23andme_20201020.txt.gz +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_C1_ALL_leave_23andme_20201020.txt.gz + +# hg19 +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_B2_ALL_leave_23andme_20201020.b37.txt.gz +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_C2_ALL_leave_23andme_20201020.b37.txt.gz + +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_A2_ALL_leave_23andme_20201020.b37.txt.gz +wget https://storage.googleapis.com/covid19-hg-public/20200915/results/20201020/COVID19_HGI_C1_ALL_leave_23andme_20201020.b37.txt.gz + +gunzip *.z + +# rename +ln -s COVID19_HGI_B2_ALL_leave_23andme_20201020.b37.txt covidHgiGwasR4.B2.hg19.txt +ln -s COVID19_HGI_B2_ALL_leave_23andme_20201020.txt covidHgiGwasR4.B2.hg38.txt +ln -s COVID19_HGI_C2_ALL_leave_23andme_20201020.b37.txt covidHgiGwasR4.C2.hg19.txt +ln -s COVID19_HGI_C2_ALL_leave_23andme_20201020.txt covidHgiGwasR4.C2.hg38.txt + +ln -s COVID19_HGI_A2_ALL_leave_23andme_20201020.txt covidHgiGwasR4.A2.hg38.txt +ln -s COVID19_HGI_A2_ALL_leave_23andme_20201020.b37.txt covidHgiGwasR4.A2.hg19.txt +ln -s COVID19_HGI_C1_ALL_leave_23andme_20201020.txt covidHgiGwasR4.C1.hg38.txt +ln -s COVID19_HGI_C1_ALL_leave_23andme_20201020.b37.txt covidHgiGwasR4.C1.hg19.txt + +wc -l covidHgiGwas*.txt + +# previous (release 3) track +#15392647 covidHgiGwas.B2.hg38.txt +# 15M +#24600933 covidHgiGwas.C2.hg38.txt + + 11830413 covidHgiGwasR4.A2.hg19.txt + 11846556 covidHgiGwasR4.A2.hg38.txt + 11434012 covidHgiGwasR4.B2.hg19.txt + 11449434 covidHgiGwasR4.B2.hg38.txt +# 11M + 13060711 covidHgiGwasR4.C1.hg19.txt + 13075183 covidHgiGwasR4.C1.hg38.txt + 12706061 covidHgiGwasR4.C2.hg19.txt + 12722892 covidHgiGwasR4.C2.hg38.txt + +# Note more filtered than release 3 data + + +# Columns +# +# head -2 covidHgiGwasR4.B2.hg19.txt + +CHR POS REF ALT SNP +all_meta_N all_inv_var_meta_beta all_inv_var_meta_sebeta all_inv_var_meta_p all_inv_var_het_p +all_meta_sample_N all_meta_AF rsid + +# 13 fields + +# These are same as Release 3, except: +# 1. hg19 lacks trailing columns from lift +# 2. there are no study-specific values (between SNP and all_meta_N) + +# determine range of effect size for filter settings, excluding low significance (log pvalue<3 items) + +# parse and generate bigBed files + +cat > makeHgi.csh << 'EOF' + set bin = ~/kent/src/hg/makeDb/outside/covid +foreach d (B2.hg19 B2.hg38 C2.hg19 C2.hg38 A2.hg19 A2.hg38 C1.hg19 C1.hg38) + set db = $d:e + set a = $d:r + set sizes = /hive/data/genomes/$db/chrom.sizes + set b = covidHgiGwasR4 + set f = $b.$a.$db + set bb = $b$a.$db + echo $b.$a.$db + if ($a == "B2") then + set studies = 21 + else if ($a == "C2") then + set studies = 33 + else if ($a == "A2") then + set studies = 12 + else if ($a == "C1") then + set studies = 20 + endif + + # colored by pvalue (effect size determines lolly height) + # NOTE: this track will not be on RR (perhaps will host on genome-test for comparison with + # earlier data + perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt R4 > $f.bed + bedSort $f.bed $f.sorted.bed + bedToBigBed -type=bed9+12 -as=$bin/covidHgiGwas.as -tab $f.sorted.bed $sizes $bb.bb + ln -s `pwd`/$bb.bb /gbdb/$db/covidHgiGwas + + # colored by effect size (pvalue determines lolly height) + #perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt R4 .025 > $f.EC.bed +# Use medians: positive theshold .178, negative threshold .119 + perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt R4 .066 .115 > $f.EC.bed + bedSort $f.EC.bed $f.sorted.EC.bed + bedToBigBed -type=bed9+12 -as=$bin/covidHgiGwas.as -tab $f.sorted.EC.bed $sizes $bb.EC.bb + ln -s `pwd`/$bb.EC.bb /gbdb/$db/covidHgiGwas + end +'EOF' + +# At Ana's suggestion, remake track with per-track effect size coloring thresholds + +cat >getEffectSizes.csh << 'EOF' +foreach a (A2 B2 C1 C2) + echo -n "$a\t" + awk '$13 >= 3.0 {print}' covidHgiGwasR4.$a.hg38.EC.bed > $a.highPvalues.bed + awk '$10 < 0 {print $21}' $a.highPvalues.bed | sort -n > negEffects.$a.txt + awk '$10 > 0 {print $21}' $a.highPvalues.bed | sort -n > posEffects.$a.txt + set neg = `ave negEffects.$a.txt | grep median | awk '{print $2}'` + set pos = `ave posEffects.$a.txt | grep median | awk '{print $2}'` + echo "$neg\t$pos" +end +'EOF' + +cat getEffectsizes.out +A2 0.130000 0.276000 +B2 0.118000 0.206000 +C1 0.094000 0.169000 +C2 0.066000 0.115000 + + +cat > reColor.csh << 'EOF' + set bin = ~/kent/src/hg/makeDb/outside/covid +foreach d (A2.hg38 B2.hg38 C1.hg38 C2.hg38 C2.hg19 A2.hg19 C1.hg19 B2.hg19 ) + set db = $d:e + set a = $d:r + set sizes = /hive/data/genomes/$db/chrom.sizes + set b = covidHgiGwasR4 + set f = $b.$a.$db + set bb = $b.$a.$db + echo $b.$a.$db + if ($a == "A2") then + set studies = 12 + # median negative and positive effect sizes + set neg = 0.130000; set pos = 0.276000 + else if ($a == "B2") then + set studies = 21 + set neg = 0.118000 ; set pos = 0.206000 + else if ($a == "C1") then + set studies = 20 + set neg = 0.094000; set pos = 0.169000 + else if ($a == "C2") then + set studies = 33 + set neg = 0.066000; set pos = 0.115000 + endif + + # colored by pvalue (effect size determines lolly height) + #perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt R4 > $f.beta.bed + #bedSort $f.beta.bed $f.beta.sorted.bed + #bedToBigBed -type=bed9+12 -as=$bin/covidHgiGwas.as -tab $f.beta.sorted.bed $sizes $bb.beta.bb + #ln -s `pwd`/$bb.beta.bb /gbdb/$db/covidHgiGwas + + # colored by effect size (pvalue determines lolly height) + perl $bin/makeCovidHgiGwas.pl $db $studies $f.txt R4 $neg $pos > $f.bed + bedSort $f.bed $f.sorted.bed + bedToBigBed -type=bed9+12 -as=$bin/covidHgiGwas.as -tab $f.sorted.bed $sizes $bb.bb + ln -s `pwd`/$bb.bb /gbdb/$db/covidHgiGwas + end +'EOF' + +csh reColor.csh >&! reColor.out & +tail -100f reColor.out +