多基因风险评分是一个用来评估个体患某种疾病风险的方法，它是通过GWAS统计数据的基因型效应值来计算的。计算公式如下：

image.png

n: 纳入的SNP的数量；
Xi：SNP I的基因型（0，1，2）
Wi：GWAS发现的SNP的权重，其中：

• 连续型变量的表型用Effect size表示，beta值
• 分类变量的表型用OR值表示

分析之前要对GWAS数据做质控，e.g. removing SNPs according to low genotyping rate, minor allele frequency>1%, removing SNPs with linkage disequilibrium with an R2 value greater than 0.1 in 500-kb
windows or imputation ‘info score>0.6’ and individuals with low genotyping rate.

plink --bfile ($target) \
    --maf 0.05 \
    --mind 0.1 \
    --geno 0.1 \
    --hwe 1e-6 \
    --make-just-bim \
    --make-just-fam \
    --out ($target).qc

做完以上步骤，将--keep (target).qc.bim 命令添加到PRSice 软件命令行中过滤SAMPLE和SNP。

安装PRSice软件：

git clone https://github.com/choishingwan/PRSice.git
cd PRSice
git checkout beta_testing
g++ --std=c++11 -I inc/ -isystem lib/ -lz -DNDEBUG -O2 -pthread src/*.cpp -o PRSice

# 安装完成后测试下是否成功
./PRSice
# 自动下载需要的包
Rscript PRSice.R --dir .

出现以下界面，说明安装成功：

image.png

需要输入文件：

• Base Dataset：以空格分割的SNP关联分析文件（association analysis results for SNPs）
  必须要有的几列是：SNP, A1, OR or BETA, P，其他几列不是必须要的，如下图：
  对于分类型表型提供OR值

  
  
  image.png
  
  

  对于连续型表型提供beta值

  
  
  
  如果表头没有以上述名称命名，必须用参数--chr, --A1, --A2, --stat, --snp, --bp, --se, --pvalue指明表头。
  比如下表：
  
  image.png
  
  参数可以这样设置：--snp SNP --chr CHR --bp BP --A1 A1 --A2 A2 --stat OR --se SE --pvalue P
  

  或者：--snp 0 --chr 1 --bp 2 --A1 3 --A2 4 --stat 5 --se 6 --pvalue 7 --index

如果文件中含有imputation info score 或者MAF，可以用--info-base <Info Name>,<Info Threshold> 和 --maf-base <MAF Name>,<MAF Threshold>参数对SNP进行过滤。
也可以对case/control分别进行过滤 --maf-base <Name for Control>,<Threshold for Control>:<Name for Case>,<Threshold for Case>

• Target Dataset 是plink产生的二进制文件.bed, .bim, 和 .fam file

• 如果fam文件没有表型数据，可以另外提供Phenotype files，加上参数--pheno-file就可以了。
  表型文件格式为FID IID pheno，以空格或tab键隔开，缺失值用NA或-9表示。也可以用--pheno-col A,B,C,D指定多个表型。

• LD reference
  当目标样本数量较少(e.g. < 500 samples)时，可以用reference panel来提高LD评估水平，参数为--ld <LD refernce>

• Clumping
  PRSice会自动评估SNP的距离并移除高LD的SNP，参数：--clump-kb, --clump-r2， --clump-p、--no-clump

• PRS calculation
  PRS计算可以有不同的model如 (e.g. add, dom, het, rec)

• 样本基因型缺失，参数--missing;
  If --missing NO_MEAN_IMPUTATION is set, the SNP for the missing samples will be excluded. Alternatively, if --missing CENTER is set, all PRS calculated will be minused by the MAF of the SNP (therefore, missing samples will have PRS of 0).

输出结果：

• Bar Plot，一般命名为[Name]BARPLOT[date].png，bar plot displaying the model fit of the PRS at P-value threshold as indicated by --bar-levels， The tallest bar indicates the best fit polygenic risk score (PRS) for the disease，only variants exceeding a P-value threshold, PT, are included，本例中最高的bar对应的P的阈值为4.7*10-18，该阈值可以解释约5%的疾病变异。

image.png

• High Resolution Plot, 一般命名为[Name]HIGH-RES_PLOT[date].png，This plot present the model fit of PRS calculated at all P-value thresholds. 最高点代表best-fit PRS
  

  image.png

• Quantile Plots，[Name]QUANTILE_PLOT[date].png，The quantile plot provide an illustration of the effect of increasing PRS on predicted risk of phenotype.
  

  image.png

• 示例：对于分类变量来说：

Rscript PRSice.R --dir . \
    --prsice ./PRSice \
    --base TOY_BASE_GWAS.assoc \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --stat OR \
    --binary-target T

• 对于连续型变量

Rscript PRSice.R --dir . \
    --prsice ./PRSice \
    --base TOY_BASE_GWAS.assoc \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --stat BETA \
    --beta \
    --binary-target F

参考：1. A guide to performing Polygenic Risk Score analyses

2. GitHub - choishingwan/PRSice: A software package for calculating, applying, evaluating and plotting the results of polygenic risk scores https://github.com/choishingwan/PRSice/
3. https://choishingwan.github.io/PRSice/step_by_step/