1- Overview
This is a guide for an introductory analysis to 1) construct a polygenic risk score (PRS) using the base data (GWAS summary statistics, particularly effect-sizes and P-values, generally public available) via a clumping and thresholding method (C + T); and 2) test the constructed PRS for prediction using the target data (PLINK binary data format). In general, it is the ‘user’ data).
2- Learning Objectives
- Apply quality control measures to base/target sample prior to PRS analysis;
- Perform PRS analysis (hands-on);
- Understand the graphs and outputs (hands-on)
3- Material
Ideally for PRS analyses, you would be using the genome-wide genotype data. Here we use base data containing summary statistics and target data containing genotypes for chromosome 16 as an example to demonstrate the workflow for prediction of simulated body mass index (BMI) data. The procedure described below will be the same for the genome-wide dataset. All the materials required for this workshop are attached here. Relevant materials for this workshop are as follows:
BMI_1kgph3_chr16_snps_summarystat.txt
3.1- Base Data
We will use as base data part of GWAS Anthropometric 2015 BMI summary statistics ( https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4382211/), made available by the GIANT consortium and were extracted from their online portal
3.1.1- QC for Base data:
1)- Check SNPs Heritability: h2SNP>0.05
LD Score Regression: https://github.com/bulik/ldsc (add ref)
SumHer: http://dougspeed.com/snp-heritability/ (add ref)
2)- The effect allele must be known. Both Base and Target datasets should have the same effect allele.
3)- Filter out SNPs with MAF < 0.01 and INFO < 0.8.
# BMI_1kgph3_chr16_snps_summarystat.txt:
# With the following columns: “SNP” = marker ID, “A1” = minor allele or effect allele, “A2” = major allele or reference allele, “Freq1.Hapmap” = allele frequency according to Hapmap, “b” = effect estimate for the minor allele, “se” = standard error of the effect estimate, “p” = p-value, “N” = sample size
SNP A1 A2 Freq1.Hapmap b se p N
1 rs1000014 G A 0.8000 -0.0055 0.0049 0.2617 233462
2 rs1000047 C T 0.6917 0.0019 0.0040 0.6348 233959
3 rs1000077 C G 0.4833 0.0006 0.0038 0.8745 220681
4 rs1000078 A G 0.7000 -0.0031 0.0041 0.4496 232588
5 rs1000100 A T 0.4167 0.0006 0.0040 0.8808 233753
6 rs1000174 A G 0.3417 -0.0025 0.0052 0.6307 150418
Tip: The base data may come in different formats. For example, if marker IDs (rs IDs) are not available, we may have to derive them based on available chromosome and map positions. The effect estimate may come in the form of an odds ratio (OR), in which case we will calculate beta using the formula beta=ln(OR). It is also ideal to have the allele calls based on the positive strand. Make sure to identify the effect allele.
3.2- Target Data
Our target data for this tutorial consist of genotype data on the subset of European samples for chromosome 16 from the 1000 Genomes datasets in PLINK format. Phenotype consists of simulated BMI data. For demonstration purposes, the BMI in the target dataset was simulated randomly from a normal distribution with no reference to the genotype data.
Tip: Base and target data should share the same ancestry background.
3.2.1- QC for Target data:
1)- Check if your Target data is in the same genome build as Base data.
LiftOver: https://genome.ucsc.edu/cgi-bin/hgLiftOver
2)- Basic QC: e.g., geno >0.99, mind <0.02, HWE P>1x10-6, 3SD HET, MAF >0.01, INFO >0.8. Also, remove indels and multi-allelic SNPs.
3)- Avoid sample overlap, as well as high degree of relatedness between individuals of Base and Target data.
4)- Check strand and allele calls
# 1kgph3_chr16.bed:
# binary plink pedigree genotype file; do not try to open this file
# 1kgph3_chr16.bim (as appearing in R):
# with 277,663 biallelic markers (SNPs or Ins/Del markers) – 6 columns without column names: “V1”=chromosome, “V2”=marker ID, “V3”=genetic distance, “V4”=chromosomal position in bp, “V5”=Allele1, “V6”=Allele2
V1 V2 V3 V4 V5 V6
1 16 rs185537431 0 60778 A G
2 16 rs377548396 0 62569 A G
3 16 rs368745239 0 66640 T G
4 16 rs187053456 0 70765 A C
5 16 rs193118147 0 70767 A G
6 16 rs201639477 0 75246 C CTTTTTTT
# 1kgph3_chr16.fam (as appearing in R):
# with 476 individuals – 6 columns without column names: “V1”=Family ID, “V2”=Individual ID, “V3”=Father ID, “V4”=Mother ID, “V5”=Sex (Male=1, Female=2), “V6”=Affected status (Yes=2, No=1, Unknown=-9)
V1 V2 V3 V4 V5 V6
1 HG00096 HG00096 0 0 1 -9
2 HG00097 HG00097 0 0 2 -9
3 HG00099 HG00099 0 0 2 -9
4 HG00101 HG00101 0 0 1 -9
5 HG00102 HG00102 0 0 2 -9
6 HG00103 HG00103 0 0 1 -9
# 1kgph3_dummybmi20200804.csv (as appearing in R):
# with the columns with column names: “V1”=Family ID, “V2”=Individual ID, “V3”=Father ID, “V4”=Mother ID, “V5”=Sex, “V6”=Affective status, “dummybmi”=phenotype of interest, “Sex” (male=1, female=0)
V1 V2 V3 V4 V5 V6 dummybmi sex
HG00096 HG00096 0 0 1 -9 13.42457034 1
HG00097 HG00097 0 0 2 -9 11.29648997 0
HG00099 HG00099 0 0 2 -9 10.98411446 0
HG00101 HG00101 0 0 1 -9 8.268308741 1
HG00102 HG00102 0 0 2 -9 11.1241505 0
3.3- Required Softwares
PLINK 1.9 (http://www.cog-genomics.org/plink2/). Use the stable version (19 Oct 2020 (b6.21)).
R (version 3.2.3+) (https://cran.r-project.org/)