多基因风险分数 PRS( Polygenic risk score)系列之五：使用PRS-CS计算PRS（beta-shrinkage方法）

本文内容

PRS系列回顾
PRS-CS简介
概念框架
使用方法
参考

系列回顾：

PRS-CS简介

本文将介绍一种使用GWAS数据与外部LD参考面板来计算SNP后验效应量的方法， PRS-CS。

该方法利用了高维贝叶斯回归的框架，与之前类似研究最大的不同点便是对于beta使用了连续收缩（continuous shrinkage， CS）先验分布，适应更多样的遗传结构，并大幅提高了计算效率，使局部LD特征的多变量建模成为可能。

概念框架

考虑一个贝叶斯高维回归框架
常用的beta的先验分布通常可以表示为多个正态分布的比例混合
为了构建更为灵活的模型以更好反映遗传结构，一些方法采用了两个或多个质点或分布的离散混合，相比单纯的正太先验分布，这类方法通常能够构建更广的效应量分布。例如LDpred。
但此类方法后验推断需要极大的计算量，可能不能准确的建模出局部的LD结构。
为了解决这些问题，PRS-CS的作者采用了一种与先前离散混合概念不同的先验分布，也就是连续收缩先验分布，也就是正态分布的全局-局部比例混合
PRS-CS对于全局比例系数采取两种不同算法进行估计。一是搜索固定的系数，而是通过指定全局系数的先验分布来在给定数据中进行估计。

详细推导请参考原论文：T Ge, CY Chen, Y Ni, YCA Feng, JW Smoller. Polygenic Prediction via Bayesian Regression and Continuous Shrinkage Priors. Nature Communications, 10:1776, 2019.

使用方法

PRS-CS由Python语言写成，是python2与3同时兼容的。需要安装scipy与h5py这两个包。

PRS-CS 本体与外部LD参考面板的下载：

https://github.com/getian107/PRScs

git clone <https://github.com/getian107/PRScs.git>

下载完成后，检验是否安装成功：

./PRScs.py --help

下载对应人群的LD参考面板，作者提供了千人基因组以及UKB数据的参考面板：

https://github.com/getian107/PRScs#getting-started

输入准备

上面下载的LD参考面板 plink格式
目标群体的bim文件：该文件只是用来提供target数据集中的SNP列表
GWAS sumstats

sumstats的格式如下

SNP          A1   A2   BETA      P
rs4970383    C    A    -0.0064   4.7780e-01
rs4475691    C    T    -0.0145   1.2450e-01
rs13302982   A    G    -0.0232   2.4290e-01

或是

SNP          A1   A2   OR        P
rs4970383    A    C    0.9825    0.5737                 
rs4475691    T    C    0.9436    0.0691
rs13302982   A    G    1.1337    0.0209
...

注意A1是effect allele， A2是non-effect allele

PRS模型计算

使用方法如下

python PRScs.py \
  --ref_dir=PATH_TO_REFERENCE \ # LD参考面板的文件夹路径
  --bim_prefix=VALIDATION_BIM_PREFIX \ #目标群体plink的bim格式文件的路径
  --sst_file=SUM_STATS_FILE \  # GWAS sumstats的路径
  --n_gwas=GWAS_SAMPLE_SIZE \ # GWAS的样本量大小
  --out_dir=OUTPUT_DIR  #输出文件夹

#以下为可选项
  #--a=PARAM_A  \ gamma-gamma prior中的参数a，默认为1
  #--b=PARAM_B  \ gamma-gamma prior中的参数b，默认为0.5
  #--phi=PARAM_PHI \ 全局比例系数，不指定时自动估计（PRS-CS-auto），也可以小规模网格搜索（phi=1e-6, 1e-4, 1e-2, 1）
  #--n_iter=MCMC_ITERATIONS \ MCMC迭代次数，默认1000
  #--n_burnin=MCMC_BURNIN \ MCMC中burnin的次数，默认500
  #--thin=MCMC_THINNING_FACTOR \ 马尔科夫链的thinning factor，默认为5
  #--chrom=CHROM \ 可以只单独计算一条染色体
  #--beta_std=BETA_STD \ 若为True，则输出标准化的后验SNP效应量 
  #--seed=SEED 随机数种子

PRS-CS使用scipy，会自动占用所有可用的cpu核心，使用服务器时可能会出现干扰，可以通过以下shell脚本指定使用核心数：

export MKL_NUM_THREADS=$N_THREADS
export NUMEXPR_NUM_THREADS=$N_THREADS
export OMP_NUM_THREADS=$N_THREADS

例如，只用一个核时，N_THREADS=1

输出结果

输出文件包含五列，rsID，碱基位置，A1，A2 与后验效应量估计值

然后使用该文件和plink可以对目标人群计算PRS： https://www.cog-genomics.org/plink/1.9/score

分染色体计算时，最后加和各个染色体的PRS即可

测试例子

使用eur的参考面板与test_data文件夹里的gwas数据，以及一个有22号染色体上1000个SNP的bim文件：

python PRScs.py \
--ref_dir=path_to_ref/ldblk_1kg_eur \
--bim_prefix=path_to_bim/test \
--sst_file=path_to_sumstats/sumstats.txt \
--n_gwas=200000 \
--chrom=22 \
--phi=1e-2 \
--out_dir=path_to_output/eur

参考

T Ge, CY Chen, Y Ni, YCA Feng, JW Smoller. Polygenic Prediction via Bayesian Regression and Continuous Shrinkage Priors. Nature Communications, 10:1776, 2019.