多基因风险评分分析（Polygenic risk score）
多基因风险评分是一个用来评估个体患某种疾病风险的方法，它是通过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
- GitHub - choishingwan/PRSice: A software package for calculating, applying, evaluating and plotting the results of polygenic risk scores https://github.com/choishingwan/PRSice/
- https://choishingwan.github.io/PRSice/step_by_step/
