This page provides four example analyses; three using GigaPRS and one using MegaPRS. Please note that the examples provide only a brief description of what each command does. If you wish to fully understand the process, you should follow the links or view the help pages for GigaPRS and MegaPRS. Note also that in the MegaPRS example, we use a reference panel containing 404 individuals, but this ideally should be larger (e.g., at least 2000 individuals).
When run using a single CPU, Examples 1 and 2 should complete in less than five minutes, while Examples 3 and 4 take about 30 and 15 minutes, respectively (all times can be reduced by using multiple CPU). All examples require less than 8Gb memory.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Three examples of running GigaPRS (requires only summary statistics)
Before we can perform the analyses, we must do some preparation; namely we obtain the data used in the examples, make the correlations details file and format the summary statistics.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Download and extract data:
The GigaPRS examples require us to download nine (compressed) files. The total size of these files is about 7Gb (so it helps to have a decent internet connection).
We will analyze summary statistics produced by European, African and East Asian GWAS, and therefore require European and East Asian SNP-SNP correlations. We can download these from a terminal using the following commands
wget https://genetics.ghpc.au.dk/doug/Correlations/HAPMAP.UK.zip
wget https://genetics.ghpc.au.dk/doug/Correlations/HAPMAP.CARAFR.zip
wget https://genetics.ghpc.au.dk/doug/Correlations/HAPMAP.EAS.zip
then extract them using the commands
unzip HAPMAP.UK.zip
unzip HAPMAP.CARAFR.zip
unzip HAPMAP.EAS.zip
The first two GigaPRS examples will use summary statistics from an East Asian GWAS of Type 2 Diabetes performed by Sakaue et al. (using data from Biobank Japan). We can download and extract these using the commands
wget https://humandbs.dbcls.jp/files/hum0197/hum0197.v3.BBJ.T2D.v1.zip
unzip hum0197.v3.BBJ.As.v1.zip
The second example also uses the European bulk summary statistics available in Resources. We can download and extract them using the commands
wget https://genetics.ghpc.au.dk/doug/HAPMAP.bulk.summaries.gz
gunzip hum0197.v3.BBJ.T2D.v1.zip
The third GigaPRS example will use summary statistics from European, African and East Asian GWAS of height performed by Yengo et al. (using data from the GIANT Consortium). In particular, we will use results from the GWAS that excluded the UK Biobank, which can be downloaded using the commands
wget https://giant-consortium.web.broadinstitute.org/images/8/8e/GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_EUR_excluding_UKB.gz
wget https://giant-consortium.web.broadinstitute.org/images/c/c5/GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_AFR_excluding_UKB.gz
wget https://giant-consortium.web.broadinstitute.org/images/a/ad/GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_EAS_excluding_UKB.gz
The following commands will download and extract genetic data for nine individuals from the 1000 Genome Project, which will be used as target individuals for all three GigaPRS analyses
wget https://genetics.ghpc.au.dk/doug/target.zip
unzip target.zip
Finally, the following command downloads a file containing details of HapMap3 SNPs (while not required, it can save time to reduce the summary statistics to HapMAP3 SNPs, and will not have any impact on accuracy)
wget https://genetics.ghpc.au.dk/doug/hapmap3.snps
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Make correlations details file and format summary statistics:
The following command makes a file containing the prefixes of the three available SNP-SNP correlations
echo "HAPMAP.UK
HAPMAP.CARAFR
HAPMAP.EAS" > cors.txt
Next we convert the summary statistics to the format required by LDAK (see Summary Statistics for details). I generally do this using the tool AWK, which is very efficient at processing large files and is usually installed by default on any UNIX operating system (you can read more about AWK here).
We can convert the Type 2 Diabetes summary statistics using the command
gunzip -c hum0197.v3.BBJ.T2D.v1/GWASsummary_T2D_Japanese_SakaueKanai2020.auto.txt.gz | awk '(NR==FNR){arr[$1];next}(FNR==1){print "Predictor", "A1", "A2", "Z", "n", "A1Freq"}($4 in arr && $ 12>0){print $4, $6, $5, $11/$12, $10, $8}' hapmap3.snps - > T2D.EAS.summaries
Note that the above command only prints out details for the SNPs present in the file hapmap3.snps (this will not impact the results, because the SNP-SNP correlations were constructed using only HapMap3 SNPs, so only these will contribute towards the PRS).
We can convert the height summary statistics using three commands
gunzip -c GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_EUR_excluding_UKB.gz | awk '(NR==1){print "Predictor", "A1", "A2", "Z", "n", "A1Freq"}(NR>1 && $8!="NA" && $9>0){print $2, $5, $6, $8/$9, $11, $7}' > height.EUR.summaries
gunzip -c GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_AFR_excluding_UKB.gz | awk '(NR==1){print "Predictor", "A1", "A2", "Z", "n", "A1Freq"}(NR>1 && $8!="NA" && $9>0){print $2, $5, $6, $8/$9, $11, $7}' > height.AFR.summaries
gunzip -c GIANT_HEIGHT_YENGO_2022_GWAS_SUMMARY_STATS_EAS_excluding_UKB.gz | awk '(NR==1){print "Predictor", "A1", "A2", "Z", "n", "A1Freq"}(NR>1 && $8!="NA" && $9>0){print $2, $5, $6, $8/$9, $11, $7}' > height.EAS.summaries
Note that there is no need to reduce to HapMap3 SNPs (because the GIANT GWAS results have already been filtered).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Summary of files:
If the downloads and extractions were successful, we should have three sets of correlations, with prefixes HAPMAP.UK, HAPMAP.CARAFR and HAPMAP.EAS. In addition, we should have target genotypes stored in target.bed, target.bim and target.fam, and the file cors.txt, which has three lines providing details of the SNP-SNP correlations.
If the formatting commands were successful, we should have four sets of formatted summary statistics, the tops of which look as follows
head -n 3 T2D.EAS.summaries
Predictor A1 A2 Z n A1Freq
rs3094315 A G -1.67198 177415 0.842874803896686
rs1048488 T C -1.65707 177415 0.844555807724612
head -n 3 height.AFR.summaries
Predictor A1 A2 Z n A1Freq
rs6678176 T C 0.944917 100691 5.328e-01
rs12069019 A G 0.223503 100692 1.359e-01
head -n 3 height.EAS.summaries
Predictor A1 A2 Z n A1Freq
rs6678176 T C -0.233064 260452 3.973e-01
rs12069019 A G -0.589327 262239 3.591e-01
head -n 3 height.EUR.summaries
Predictor A1 A2 Z n A1Freq
rs6678176 T C -1.59522 1175334 3.080e-01
rs12069019 A G -0.931926 1175334 1.432e-01
If everything looks correct, you are ready to run the three GigaPRS examples
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
GigaPRS Example 1:
This example produces East Asian estimates of SNP effect sizes for Type 2 Diabetes, then uses these to compute PRS for the target individuals. We run GigaPRS with one set of summary statistics (those providing results from the East Asian GWAS of Type 2 Diabetes).
We make the file sums1.txt, which contains details of the summary statistics
echo "T2D.EAS.summaries F 3" > sums1.txt
This file specifies that there is one set of summary statistics, which is stored in the file T2D.EAS.summaries, and that these summary statistics are focal and match to the third set of SNP-SNP correlations (those with prefix HAPMAP.EAS).
We can now run GigaPRS Step 1 using the command
./ldak.out --giga-prs fit1 --corslist cors.txt --sumslist sums1.txt
Note that if you have multiple CPUs, you can speed up the analysis by adding (say) --max-threads 4. The estimated SNP effect sizes are stored in Column 5 of fit1.effects. Note that because we only use one set of summary statistics and one ancestry, it would be equivalent to run the command
./ldak.out --giga-prs fit1 --cors HAPMAP.EAS --summary T2D.EAS.summaries
We can run GigaPRS Step 2 using the command
./ldak.out --calc-scores prs1 --bfile target --scorefile fit1.effects
The estimated PRS for the target individuals are stored in prs1.profile, which contains the following lines
ID1 ID2 Phenotype Covariates Profile_1 BLANK
HG00246 EUR.GBR NA 0.000000 0.416744 NA
HG00473 EAS.CHS NA 0.000000 0.023841 NA
HG01107 AMR.PUR NA 0.000000 1.375774 NA
HG02816 AFR.GWD NA 0.000000 3.203354 NA
HG03012 SAS.BEB NA 0.000000 -0.139241 NA
HG03926 SAS.BEB NA 0.000000 1.027284 NA
HG04225 SAS.ITU NA 0.000000 2.169381 NA
NA18563 EAS.CHB NA 0.000000 0.936203 NA
NA19782 AMR.MXL NA 0.000000 3.903008 NA
Note that because Step 1 only estimated East Asian effect sizes, the PRS will likely have poor accuracy for individuals of different ancestries (e.g., European, Hispanic, African and South Asian).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
GigaPRS Example 2:
This example also produces East Asian estimates of SNP effect sizes for Type 2 Diabetes, then uses these to compute PRS for the target individuals. We again run GigaPRS with one set of focal summary statistics (those providing results from the East Asian GWAS of Type 2 Diabetes), but now also provide 52 sets of secondary summary statistics (those in the bulk summary statistics file).
We make the file sums2.txt, which contains details of the summary statistics
echo "T2D.EAS.summaries F 3" > sums2.txt
Note that this matches the file used for the first example (i.e., provides details for one set of summary statistics, which are focal and matched to the third set of SNP-SNP correlations).
We can now run GigaPRS Step 1 using the command
./ldak.out --giga-prs fit2 --corslist cors.txt --sumslist sums2.txt --bulk-summaries HAPMAP.bulk.summaries
This command differs from that used in the first GigaPRS example, because we now provide bulk summary statistics (results from 52 European GWAS of heritable quantitative traits that are treated as secondary summary statistics). As part of its analysis, GigaPRS estimates the informativeness of each set of secondary summary statistics (which depends mainly on the genetic correlation with the Type 2 Diabetes summary statistics). From the screen output, we see that GigaPRS determines that three sets of secondary summary statistics are sufficiently informative to include when estimating effect sizes, namely those for body mass index (Trait 10), SHBG (Trait 20) and HbA1c (Trait 31), whose inclusions are estimated to increase the effective sample size by 54%, 37% and 123%, respectively. The estimated SNP effect sizes are stored in Column 5 of fit2.effects.
We can run GigaPRS Step 2 using the command
./ldak.out --calc-scores prs2 --bfile target --scorefile fit2.effects
The estimated PRS for the target individuals are stored in prs2.profile, which contains the following lines
ID1 ID2 Phenotype Covariates Profile_1 BLANK
HG00246 EUR.GBR NA 0.000000 -2.985358 NA
HG00473 EAS.CHS NA 0.000000 -0.677189 NA
HG01107 AMR.PUR NA 0.000000 -2.297894 NA
HG02816 AFR.GWD NA 0.000000 -0.313217 NA
HG03012 SAS.BEB NA 0.000000 -1.165671 NA
HG03926 SAS.BEB NA 0.000000 0.538093 NA
HG04225 SAS.ITU NA 0.000000 0.290979 NA
NA18563 EAS.CHB NA 0.000000 1.115676 NA
NA19782 AMR.MXL NA 0.000000 1.087831 NA
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
GigaPRS Example 3:
This example produces European, African and East Asian estimates of SNP effect sizes for height, then uses these to compute PRS for the target individuals. We now run GigaPRS with three sets of focal summary statistics (those providing results from the European, African and East Asian GWAS of height).
We make the summary statistics details file
echo "height.EUR.summaries F 1
height.AFR.summaries F 2
height.EAS.summaries F 3" > sums3.txt
This file specifies that there are three set of focal summary statistics, stored in the files height.EUR.summaries, height.AFR.summaries and height.EAS.summaries, which match to the first, second and third SNP-SNP correlations (those with prefix HAPMAP.EUR, HAPMAP.AFR and HAPMAP.EAS, respectively).
We can now run GigaPRS Step 1 using the command
./ldak.out --giga-prs fit3 --corslist cors.txt --sumslist sums3.txt
The European, African and East Asian estimates of SNP effect sizes are stored in Columns 5, 7 and 9 of fit3.combined.effects (meanwhile, Columns 4, 6 and 8 contain the average A1 allele counts for each ancestry).
We can run GigaPRS Step 2 using the command
./ldak.out --calc-scores prs3 --bfile target --multi-scorefile fit3.combined.effects
The estimated PRS for the target individuals are stored in prs3.profile, which contains the following lines
ID1 ID2 Phenotype Covariates Combined_Profile BLANK
HG00246 EUR.GBR NA 0.000000 -1.085577 NA
HG00473 EAS.CHS NA 0.000000 0.853548 NA
HG01107 AMR.PUR NA 0.000000 -3.927584 NA
HG02816 AFR.GWD NA 0.000000 0.317691 NA
HG03012 SAS.BEB NA 0.000000 -3.768587 NA
HG03926 SAS.BEB NA 0.000000 -3.088692 NA
HG04225 SAS.ITU NA 0.000000 -1.729275 NA
NA18563 EAS.CHB NA 0.000000 -0.299295 NA
NA19782 AMR.MXL NA 0.000000 -2.522715 NA
Meanwhile, the estimated ancestry proportions are stored in prs3.ancestry, which contains the following lines
ID1 ID2 Raw1 Raw2 Raw3 Processed1 Processed2 Processed3
HG00246 EUR.GBR 0.9851 0.0042 0.0109 1.0000 0.0000 0.0000
HG00473 EAS.CHS -0.0143 -0.0010 1.0212 0.0000 0.0000 1.0000
HG01107 AMR.PUR 0.6646 0.2617 0.0729 0.6651 0.2619 0.0730
HG02816 AFR.GWD -0.0512 1.0493 0.0134 0.0000 1.0000 0.0000
HG03012 SAS.BEB 0.4856 0.1007 0.4014 0.4916 0.1020 0.4064
HG03926 SAS.BEB 0.4764 0.0994 0.4182 0.4793 0.1000 0.4207
HG04225 SAS.ITU 0.5742 0.0935 0.3163 0.5836 0.0950 0.3214
NA18563 EAS.CHB 0.0107 0.0002 0.9917 0.0000 0.0000 1.0000
NA19782 AMR.MXL 0.6034 0.0763 0.3105 0.6094 0.0770 0.3136
Reassuringly, we see that the European individual (ID HG00246 EUR.GBR) is estimated to be 100% from Ancestry 1 (European), the African individual (ID HG02816 AFR.GGWD) is estimated to be 100% from Ancestry 2 (African), while the East Asian individuals (IDs HG00473 EAS.CHS and NA18563 EAS.CHB) are both estimated to be 100% from Ancestry 3 (East Asian). Meanwhile the Hispanic and South Asian individuals (those with AMR or SAS in their IDs) are estimated to be a mixture of the three ancestries.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Example of running MegaPRS (requires summary statistics and a reference panel):
Here we use summary statistics from the 2014 meta-analysis of height from the GIANT Consortium, while the reference panel uses genotype data for non-Finnish Europeans from the 1000 Genomes Project. The scripts use the tool AWK, which is very efficient at processing large files and is usually installed by default on any UNIX operating system. You can read more about AWK here.
When running LDAK, always read the screen output, which suggests arguments and estimates memory usage.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
1 - Preparation
Download and extract a non-Finnish European reference panel constructed from 1000 Genomes Project data (the scripts we used to make this data match those in the example for Reference Panel, except we reduced to HapMap3 SNPs). Note that we normally recommend using a larger panel (at least 2000 samples).
wget https://genetics.ghpc.au.dk/doug/RefPanels/reduced.1000G.404EUR.tar.gz
tar -xzvf reduced.1000G.404EUR.tar.gz
Download and format the summary statistics. Note that the second command retains only SNPs with non-matching, single-character alleles, MAF above 0.01 and sample size above 200k, then excludes rows that contain NAs (for more details, see Summary Statistics).
wget https://portals.broadinstitute.org/collaboration/giant/images/0/01/GIANT_HEIGHT_Wood_et_al_2014_publicrelease_HapMapCeuFreq.txt.gz
gunzip -c GIANT_HEIGHT_Wood_et_al_2014_publicrelease_HapMapCeuFreq.txt.gz | awk '(NR==1){print "Predictor A1 A2 Z n A1Freq"}(NR>1){snp=$1;a1=$2;a2=$3;freq=$4;maf=freq;if(freq>0.5){maf=1-freq};effect=$5;se=$6;z=effect/se;n=$8; if(a1!=a2 && (a1=="A"||a1=="C"||a1=="G"||a1=="T") && (a2=="A"||a2=="C"||a2=="G"||a2=="T") && maf>0.01 && n>200000){print snp, a1, a2, z, n, freq}}' - | grep -v NA > height.txt
Download the breakpoints file from Resources.
wget https://genetics.ghpc.au.dk/doug/berisa.txt
Download and unzip the file BaselineLD.zip from Resources.
wget https://genetics.ghpc.au.dk/doug/BaselineLD.zip
unzip BaselineLD.zip
At this point, we should have the following files:
1000G.404EUR.hapmap3.bed
1000G.404EUR.hapmap3.bim
1000G.404EUR.hapmap3.fam
height.txt
berisa.txt
96 files with prefix BaselineLD
We can check we have all the required files by running these two commands
ls 1000G.404EUR.hapmap3.{bed,bim,fam} height.txt berisa.txt | grep XXX
ls BaselineLD{1..96} | grep XXX
The screen output from these commands will report if any files are missing (so if there is no output, all files are present).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2 - Construct the prediction model.
Calculate predictor-predictor correlations (see MegaPRS). Note that if you have multiple CPUs, you can speed up the analysis by adding (say) --max-threads 4.
./ldak.out --calc-cors cors --bfile 1000G.404EUR.hapmap3 --break-points berisa.txt
Estimate the effect sizes (see MegaPRS). Again, you can speed up the analysis by using multiple CPUs.
./ldak.out --mega-prs mega --cors cors --summary height.txt --annotation-number 86 --annotation-prefix BaselineLD
The estimated effect sizes are saved in mega.effects.
