多基因风险分数 PRS( Polygenic risk score)系列之十: PRS-CSx 跨祖先PRS的构建

本文内容:

  1. PRS-CSx简介
  2. PRS-CSx使用方法
  3. PRScsx实例应用
  4. 参考

回顾

  1. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之一:概念入门
  2. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)
  3. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之三:使用PRSice计算PRS(C+T方法)
  4. ldpred
  5. GWASLab:多基因风险分数 PRS(Polygenic risk score)系列之五:使用PRS-CS计算PRS(beta-shrinkage方法)
  6. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之六:metaGRS介绍
  7. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之七:Pathway-based PRS 通路PRS
  8. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之八:PGS Catalog 数据库
  9. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之九: 使用PLINK2分染色体计算PRS并加和

PRS-CSx简介

先前的文章中介绍了PRS研究中的一大问题便是在A群体中的构建的PRS难以直接转移应用到B群体中。为了解决这一问题,Yunfeng Ruan等人开发了PRS-CSx。

PRS-CSx是一个贝叶斯多基因模型构建与预测的框架,通过整合多个族裔的GWAS概括性统计数据来提升跨群体PRS的预测能力。该方法为PRS-CS的扩展 (参考:GWASLab:多基因风险分数 PRS(Polygenic risk score)系列之五:使用PRS-CS计算PRS(beta-shrinkage方法))。

原理上,PRS-CSx利用了一个共同的连续收缩先验分布来整合各个群体中SNP的效应,该方法通过在GWAS概括性统计数据之间共享先验分布,利用不同群体间的LD信息,来达到更准确的效应估计量。这个共享的先验分布考虑到了效应估计量在不同群体中相互关联但又存在差异的特点,保持了模型框架的灵活性。

PRS-CSx使用的先验分布 其中全局与局部收缩系数不随群体k变化

给定GWAS概括性统计数据,以及相应群体的LD参考面板,PRS-CSx可以对每个群体计算分别的PRS,并通过最优线性组合来得出最终的PRS.

PRS-CSx使用方法

https://github.com/getian107/PRScsx

PRScsx是一个基于Python的命令行工具,需要安装scipy与h5py这两个依赖包。从github上下载PRS-CSx:

git clone https://github.com/getian107/PRScsx.git

LD 参考面板与 PRS-CS 所使用文件相同 (参考:GWASLab:多基因风险分数 PRS(Polygenic risk score)系列之五:使用PRS-CS计算PRS(beta-shrinkage方法))。

下载链接(国内可用的FTP):https://personal.broadinstitute.org/hhuang//public//PRS-CSx/Reference

记得同时下载对应面板的snp list:snpinfo_mult_1kg_hm3 (1kg),或是 snpinfo_mult_ukbb_hm3(ukbb)

选项

python PRScsx.py \
--ref_dir=PATH_TO_REFERENCE \
--bim_prefix=VALIDATION_BIM_PREFIX \
--sst_file=SUM_STATS_FILE \
--n_gwas=GWAS_SAMPLE_SIZE \
--pop=POPULATION \
--out_dir=OUTPUT_DIR \
--out_name=OUTPUT_FILE_PREFIX \
--a=PARAM_A \
--b=PARAM_B \
--phi=PARAM_PHI \
--n_iter=MCMC_ITERATIONS \
--n_burnin=MCMC_BURNIN \
--thin=MCMC_THINNING_FACTOR \
--chrom=CHROM \
--meta=META_FLAG \
--seed=SEED

必须的参数:

  • PATH_TO_REFERENCE:LD参考面板的路径,路径下应包含相应群体的参考面板以及snp list. 例如,纳入群体为EUR以及EAS,指定路径为:./ldref ,那么该路径下应该有 ldblk_1kg_eas,ldblk_1kg_eur 这两个文件夹, 以及snpinfo_mult_1kg_hm3这个文件。
  • VALIDATION_BIM_PREFIX:目标数据集的bim文件。
  • SUM_STATS_FILE:sumstats的完整路径,由逗号分隔。
  • GWAS_SAMPLE_SIZE:sumstats的样本量大小,由逗号分隔,顺序与SUM_STATS_FILE一致。
  • POPULATION:对应的群体,可以为 AFR, AMR, EAS, EUR, SAS,由逗号分隔,顺序与SUM_STATS_FILE一致。
  • OUTPUT_DIR: 输出的路径
  • OUTPUT_FILE_PREFIX:输出文件前缀

其余为可选参数:

META_FLAG : 如果为True,则输出inverse-variance-weighted meta-analysis of the population-specific posterior effect size estimates。

PARAM_A, PARAM_B, PARAM_PHI,MCMC_ITERATIONS,MCMC_BURNIN,MCMC_BURNIN,SEED与CHROM 使用方法与PRScs一致。(参考:GWASLab:多基因风险分数 PRS(Polygenic risk score)系列之五:使用PRS-CS计算PRS(beta-shrinkage方法)

示例代码

python PRScsx.py \
--ref_dir=path_to_ref \
--bim_prefix=path_to_bim/test \
--sst_file=path_to_sumstats/EUR_sumstats.txt,path_to_sumstats/EAS_sumstats.txt \
--n_gwas=200000,100000 \
--pop=EUR,EAS \
--chrom=22 \
--phi=1e-2 \
--out_dir=path_to_output \
--out_name=test

注意:将路径替换为自己的路径

大约一分钟即可完成计算。

运行log如下:

*** 2 discovery populations detected ***

##### process chromosome 22 #####
... parse reference file: /home/heyunye/tools/prscs/ldref/snpinfo_mult_1kg_hm3 ...
... 18944 SNPs on chromosome 22 read from /home/heyunye/tools/prscs/ldref/snpinfo_mult_1kg_hm3 ...
... parse bim file: /home/heyunye/tools/prscsx/PRScsx/test_data/test.bim ...
... 1000 SNPs on chromosome 22 read from /home/heyunye/tools/prscsx/PRScsx/test_data/test.bim ...
... parse EUR sumstats file: /home/heyunye/tools/prscsx/PRScsx/test_data/EUR_sumstats.txt ...
... 1000 SNPs read from /home/heyunye/tools/prscsx/PRScsx/test_data/EUR_sumstats.txt ...
... 1000 common SNPs in the EUR reference, EUR sumstats, and validation set ...
... parse EAS sumstats file: /home/heyunye/tools/prscsx/PRScsx/test_data/EAS_sumstats.txt ...
... 1000 SNPs read from /home/heyunye/tools/prscsx/PRScsx/test_data/EAS_sumstats.txt ...
... 901 common SNPs in the EAS reference, EAS sumstats, and validation set ...
... parse EUR reference LD on chromosome 22 ...
... parse EAS reference LD on chromosome 22 ...
... align reference LD on chromosome 22 across populations ...
... 1000 valid SNPs across populations ...
... MCMC ...
--- iter-100 ---
--- iter-200 ---
--- iter-300 ---
--- iter-400 ---
--- iter-500 ---
--- iter-600 ---
--- iter-700 ---
--- iter-800 ---
--- iter-900 ---
--- iter-1000 ---
--- iter-1100 ---
--- iter-1200 ---
--- iter-1300 ---
--- iter-1400 ---
--- iter-1500 ---
--- iter-1600 ---
--- iter-1700 ---
--- iter-1800 ---
--- iter-1900 ---
--- iter-2000 ---
... Done ...

输出为EUR以及EAS的PRS:

test_EAS_pst_eff_a1_b0.5_phi1e-02_chr22.txt

test_EUR_pst_eff_a1_b0.5_phi1e-02_chr22.txt

head test_EAS_pst_eff_a1_b0.5_phi1e-02_chr22.txt
22      rs9605903       17054720        C       T       8.694291e-04
22      rs5746647       17057138        G       T       -1.005430e-03
22      rs5747999       17075353        C       A       -2.499230e-04
22      rs2845380       17203103        A       G       6.037999e-04
22      rs2247281       17211075        G       A       4.780305e-04
22      rs2845346       17214252        C       T       7.767527e-04
22      rs2845347       17214669        C       T       1.671207e-03
22      rs1807512       17221495        C       T       -1.778397e-03
22      rs5748593       17227461        T       C       9.849030e-04
22      rs9606468       17273728        C       T       1.442600e-04

使用该文件便可以利用plink进行PRS计算:

GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之九: 使用PLINK2分染色体计算PRS并加和

PRScsx实例应用

PRScsx的通讯作者以第一作者的身份,将PRScsx应用于二型糖尿病的跨族裔PRS研究中, 文中使用PRScsx和European, African American,以及East Asian的GWAS数据,构建了二型糖尿病的跨族裔PRS。

https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-022-01074-2

参考

Ruan, Y., Lin, Y. F., Feng, Y. C. A., Chen, C. Y., Lam, M., Guo, Z., … & Ge, T. (2022). Improving polygenic prediction in ancestrally diverse populations. Nature Genetics54(5), 573-580.

Ge, T., Irvin, M. R., Patki, A., Srinivasasainagendra, V., Lin, Y. F., Tiwari, H. K., … & Karlson, E. W. (2022). Development and validation of a trans-ancestry polygenic risk score for type 2 diabetes in diverse populations. Genome medicine14(1), 1-16.

多基因风险分数 PRS( Polygenic risk score)系列之九: 使用PLINK2分染色体计算PRS并加和

本文内容:

  1. 下载PGS分数文件
  2. 第一步 计算每条染色体的原始分数
  3. 第二步 将各个染色体的结果提取
  4. 第三步 将分数进行简单加和

本文主要演示分染色体计算PRS并求和的过程。

下载PGS分数文件 Download PGS score

演示PGS分数文件下载自PGScatalog

下载harmonise后的文件:

wget https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS000012/ScoringFiles/Harmonized/PGS000012_hmPOS_GRCh37.txt.gz
gunzip PGS000012_hmPOS_GRCh37.txt.gz
head PGS000012_hmPOS_GRCh37.txt
###PGS CATALOG SCORING FILE - see <https://www.pgscatalog.org/downloads/#dl_ftp_scoring> for additional information
#format_version=2.0
##POLYGENIC SCORE (PGS) INFORMATION
#pgs_id=PGS000012
#pgs_name=GRS49K
#trait_reported=Coronary artery disease
#trait_mapped=coronary artery disease
#trait_efo=EFO_0001645
#genome_build=hg19
#variants_number=49310

#下载后去除文件头部的comment,方便Plink2使用
awk '$1!~/#/{print $0}' PGS000012_hmPOS_GRCh37.txt > PGS000012_hmPOS_GRCh37_plink2.txt

head PGS000012_hmPOS_GRCh37_plink2.txt
rsID    effect_allele   effect_weight   hm_source       hm_rsID hm_chr  hm_pos  hm_inferOtherAllele
rs1333045       C       0.187251        ENSEMBL rs1333045       9       22119195        T
rs1537370       T       0.17296 ENSEMBL rs1537370       9       22084310        C
rs9863247       T       0.136328        ENSEMBL rs9863247       3       161125373       C
rs11203077      T       0.131725        ENSEMBL rs11203077      10      91097085        G
rs10268558      C       0.12966 ENSEMBL rs10268558      7       18924927        T
rs7148203       C       0.128366        ENSEMBL rs7148203       14      45062275        T
rs12747328      T       0.127016        ENSEMBL rs12747328      1       159429702       C
rs17810947      G       0.126626        ENSEMBL rs17810947      18      70754886        A/T
rs2412710       A       0.125677        ENSEMBL rs2412710       15      42683787        G

基因组文件使用 https://choishingwan.github.io/PRS-Tutorial/plink/ 所提供的的EUR文件

EUR.QC.bed
EUR.QC.bim
EUR.QC.fam

注:本教程计算结果无实际意义,仅为演示计算过程使用。

第一步 计算每条染色体的原始分数 step1 calculate the raw scores for each chromosome

score=PGS000012_hmPOS_GRCh37_plink2.txt

for chr in $(seq 1 22)
do
plink2 \
    --bfile ./EUR.QC \
    --score ${score} 1 2 3 header list-variants cols=+scoresums \
    --chr ${chr} \
    --out EUR_CAD.chr${chr} 
done

—score 后跟 分数文件路径, 并指定 文件的 ID, effect allele 以及 effect_weight的列数。

header :分数文件有header,计算时会自动跳过第一行

list-variants : 将计算所用的所有variant ID写进plink2.sscore.vars

—chr 指定所计算的染色体编号

注1: 这里要加上cols=+scoresums 来获取原始结果,默认输出的仅为平均后的结果

--score applies one or more linear scoring systems to each sample, and reports results to plink2.sscore. More precisely, if G is the full genotype/dosage matrix (rows = alleles, columns = samples) and a is a scoring-system vector with one coefficient per allele, --score computes the vector-matrix product aTG, and then divides by the number of variants when reporting score-averages.

注2:有时候分数文件中不包括某些染色体的variants,所以有些染色会计算失败,求和时应注意

所得文件

EUR_CAD.chr10.log          EUR_CAD.chr15.log          EUR_CAD.chr1.log           EUR_CAD.chr3.log          EUR_CAD.chr8.log
EUR_CAD.chr10.sscore       EUR_CAD.chr15.sscore       EUR_CAD.chr1.sscore        EUR_CAD.chr3.sscore       EUR_CAD.chr8.sscore
EUR_CAD.chr10.sscore.vars  EUR_CAD.chr15.sscore.vars  EUR_CAD.chr1.sscore.vars   EUR_CAD.chr3.sscore.vars  EUR_CAD.chr8.sscore.vars
EUR_CAD.chr11.log          EUR_CAD.chr16.log          EUR_CAD.chr20.log          EUR_CAD.chr4.log          EUR_CAD.chr9.log
EUR_CAD.chr11.sscore       EUR_CAD.chr16.sscore       EUR_CAD.chr20.sscore       EUR_CAD.chr4.sscore       EUR_CAD.chr9.sscore
EUR_CAD.chr11.sscore.vars  EUR_CAD.chr16.sscore.vars  EUR_CAD.chr20.sscore.vars  EUR_CAD.chr4.sscore.vars  EUR_CAD.chr9.sscore.vars
EUR_CAD.chr12.log          EUR_CAD.chr17.log          EUR_CAD.chr21.log          EUR_CAD.chr5.log          
EUR_CAD.chr12.sscore       EUR_CAD.chr17.sscore       EUR_CAD.chr21.sscore       EUR_CAD.chr5.sscore       
EUR_CAD.chr12.sscore.vars  EUR_CAD.chr17.sscore.vars  EUR_CAD.chr21.sscore.vars  EUR_CAD.chr5.sscore.vars  
EUR_CAD.chr13.log          EUR_CAD.chr18.log          EUR_CAD.chr22.log          EUR_CAD.chr6.log          
EUR_CAD.chr13.sscore       EUR_CAD.chr18.sscore       EUR_CAD.chr22.sscore       EUR_CAD.chr6.sscore       
EUR_CAD.chr13.sscore.vars  EUR_CAD.chr18.sscore.vars  EUR_CAD.chr22.sscore.vars  EUR_CAD.chr6.sscore.vars  
EUR_CAD.chr14.log          EUR_CAD.chr19.log          EUR_CAD.chr2.log           EUR_CAD.chr7.log
EUR_CAD.chr14.sscore       EUR_CAD.chr19.sscore       EUR_CAD.chr2.sscore        EUR_CAD.chr7.sscore
EUR_CAD.chr14.sscore.vars  EUR_CAD.chr19.sscore.vars  EUR_CAD.chr2.sscore.vars   EUR_CAD.chr7.sscore.vars

分数文件内容

head EUR_CAD.chr1.sscore
#FID    IID     ALLELE_CT       NAMED_ALLELE_DOSAGE_SUM SCORE1_AVG      SCORE1_SUM
HG00096 HG00096 1046    486     0.00757595      7.92444
HG00097 HG00097 1046    484     0.00783352      8.19386
HG00099 HG00099 1046    505     0.00798495      8.35226
HG00101 HG00101 1046    491     0.00791584      8.27997
HG00102 HG00102 1046    474     0.00769178      8.04561
HG00103 HG00103 1046    502     0.00799328      8.36097
HG00105 HG00105 1046    500     0.00793847      8.30364
HG00107 HG00107 1046    468     0.00770316      8.0575
HG00108 HG00108 1046    469     0.00732639      7.66341

分染色体的计算的目的是,实际操作中,通常会使用巨大的pgen或bgen进行计算,合并文件并不现实, 这是一般需要分染色体计算并手动求和:

for chr in $(seq 1 22)
do
plink2 \
    --pfile chr${chr} \
    --score ${score} 1 2 3 header list-variants cols=+scoresums \
    --out chr${chr}
done

第二步 将各个染色体的结果提取 step2 extract raw scores into a single file

scoreFile=EUR_CAD

count=1
for chr in $(seq 1 22)
do
    file=${scoreFile}.chr${chr}.sscore
    if [ -s "$file" ] ; then 
        if [ ${count} -eq 1 ] ;then
        echo "ID        chr${chr}" > ./${scoreFile}.chrALL.score
        awk 'NR>1{print $1,$6}' ./${file} >> ./${scoreFile}.chrALL.score
        else
        paste ./${scoreFile}.chrALL.score <(cut -f6 ./${file} | awk -v chr=chr${chr} 'NR==1 {$0=chr} 1') >./${scoreFile}.temp&&mv ./${scoreFile}.temp ./${scoreFile}.chrALL.score
        fi
        count=`expr $count + 1`
     fi
done


所得文件

head EUR_CAD.chrALL.score
ID      chr1    chr2    chr3    chr4    chr5    chr6    chr7    chr8    chr9    chr10   chr11   chr12   chr13   chr14   chr15   chr16   chr17   chr18   chr19   chr20       chr21   chr22
HG00096 7.92444 6.41628 4.9392  3.9937  5.08024 7.50241 4.61922 3.59464 4.24131 4.56312 5.03334 5.55227 2.95436 2.32917 3.39488 3.61625 3.40913 2.02032 1.95575 2.74394     0.79493 1.33408
HG00097 8.19386 6.33598 5.0713  4.0514  4.92887 7.30242 4.19193 3.66247 4.514   4.42716 5.26378 5.11748 2.7589  2.47905 3.82365 3.83448 3.44519 2.09788 1.80506 2.7079      0.811121        1.27481
HG00099 8.35226 6.61542 5.199   3.94708 4.88824 7.72775 4.31098 3.58159 4.30785 4.81232 5.2327  6.05063 2.8803  2.3651  3.36827 3.54378 3.29725 2.08449 1.88188 2.52587     0.767397        1.21001
HG00101 8.27997 6.13397 5.36465 4.04592 4.88809 7.86323 4.49118 3.60872 4.27056 4.6258  5.42121 5.26046 2.70375 2.52093 3.56581 3.45794 3.30534 2.1327  1.70916 2.66579     0.818805        1.3808
HG00102 8.04561 6.32975 5.20324 3.86051 4.67724 7.78206 4.69662 3.38095 4.40059 4.67174 5.2312  5.04913 2.69084 2.48935 3.62454 3.61557 3.3044  2.12043 1.96028 2.82496     0.808453        1.46265
HG00103 8.36097 6.1202  5.31402 4.24136 4.91736 7.34715 4.55809 3.60716 4.15705 4.66235 5.04477 5.53387 2.70537 2.36537 3.76842 3.64108 3.33273 2.11081 2.00393 2.68089     0.763437        1.32735
HG00105 8.30364 6.13937 5.15716 4.19656 5.06392 7.8889  4.33995 3.78158 4.24975 4.58436 5.32105 5.72797 2.84577 2.43968 3.39822 3.60722 3.44641 1.9986  1.775   2.7736      0.793725        1.41736
HG00107 8.0575  6.51259 5.19973 4.06067 4.78711 7.14954 4.66686 3.81154 4.42637 4.45144 5.26257 5.46242 2.82289 2.28362 3.33811 3.58244 3.25041 2.11758 2.00922 2.70618     0.917279        1.43502
HG00108 7.66341 6.18669 4.98761 4.05088 4.7294  7.39041 4.08688 3.53173 4.43728 4.80822 5.02485 5.52436 3.12156 2.35199 3.54311 3.49076 2.97043 2.12883 1.91326 2.51056     0.817221        1.47096

第三步 将分数进行简单加和 step3 combine the scores using awk

scoreFile=EUR_CAD
#sum the scores from each chr
awk 'BEGIN{t=0}
        {if(NR==1){
                print "ID","score"
                  }
                else{
                        for(i=2;i<=NF;i++) t+=$i 
                        print $1,t
                    }
        }' ./${scoreFile}.chrALL.score > ./${scoreFile}.summary.score

所得文件

head EUR_CAD.summary.score
ID score
HG00096 88.013
HG00097 176.112
HG00099 265.062
HG00101 353.577
HG00102 441.807
HG00103 530.37
HG00105 619.62
HG00107 707.931
HG00108 794.672

求和完成,获得了原始分数(没有进行任何操作,没有标准化,没有平均,只是所有variant的效应量与dosage乘积的和),可以用于后续PRS分析

参考:

https://www.cog-genomics.org/plink/2.0/score

https://www.pgscatalog.org/browse/scores/

https://choishingwan.github.io/PRS-Tutorial/plink/

多基因风险分数 PRS( Polygenic risk score)系列之八:PGS Catalog

本文内容

  1. PGS Catalog 简介
  2. PGS Catalog的纳入标准
  3. 从PGS Catalog寻找PGS
  4. PGS分数文件格式与下载
  5. 下载PGS后使用PLINK计算PGS
  6. 参考

回顾

PGS Catalog 简介

本文简要介绍PGS Catalog的基本信息与使用方法。该数据库对于未来的PGS相关研究可以说是必不可少的,或多或少都需要通过此数据库查询,下载已有的PGS或上传自己的PGS模型。

与 GWAS catalog 类似, PGS Catalog 是一个已发表多基因风险分数 (polygenic scores)的公开数据库。在PGS Catalog中的每个PGS都被统一地标注了相关的元信息:包括分数文件(variants, effect alleles/weights),PGS如何构建与应用的注释,以及其预测表现的评价等。PGS对应的表型会被连接到相应的EFO(Experimental Factor Ontology,https://www.ebi.ac.uk/efo/)以保持研究间的统一(GWAS catalog 也使用EFO)。

PGS Catalog旨在对PGS构建索引,并以标准化的形式分发每个PGS的关键信息(variants,结果,实验设计等),以促进对PGS分析有效性的评价。

该数据库由剑桥大学Michael Inouye(在推上很活跃的大佬,建议关注)组的Samuel Lambert与HDR UK及NHGRI-EBI (GWAS Catalog)合作开发。

PGS Catalog 的主页

PGS Catalog的纳入标准

纳入PGSCatalog的标准主要有两大块:

  1. 新近开发的PGS,包含其分数与预测能力的必要基础信息 (需要在独立样本中评估
  2. 对已开发的PGS在新的群体中进行评估。

纳入后每一个PGS都被赋予了识别编号, 例如 PGS000001

从PGS Catalog寻找PGS

查询PGS时, 可以通过搜索框直接搜索关键词查询PGS,或是通过表型,发表的文献等方式浏览数据库中的PGS.

以breast cancer 为例,查询后可以看到数据库中目前有112个乳腺癌相关的PGS:

点击后,可以查看这些PGS的汇总信息:

可以通过ancestry对PGS进行过滤,列表中的ancestry distribution表示的是,所用样本群体中各个族裔的构成。

选取感兴趣的PGS后,可以点击进入查看详细信息,或是直接下载PGS模型文件。

每个PGS的页面包括了PGS的详细信息,构建方法与参数,原始GWAS数据,评价指标,评价时所用样本信息等等。

PGS分数文件格式与下载

PGS Catalog数据库中的文件格式说明可以参考:https://www.pgscatalog.org/downloads/

如下所示,Scoring File Format由两部分组成,header和数据。

Header部分主要为该文件版本信息,PGS的基础信息,以及原始研究的信息,数据部分则包括了variant和计算PGS的allele与权重,大多可以来直接使用。

###PGS CATALOG SCORING FILE - see <https://www.pgscatalog.org/downloads/#dl_ftp_scoring> for additional information
#format_version=2.0
##POLYGENIC SCORE (PGS) INFORMATION
#pgs_id=PGS000001
#pgs_name=PRS77_BC
#trait_reported=Breast Cancer
#trait_mapped=breast carcinoma
#trait_efo=EFO_0000305
#weight_type=NR
#genome_build=NR
#variants_number=77
##SOURCE INFORMATION
#pgp_id=PGP000001
#citation=Mavaddat N et al. J Natl Cancer Inst (2015). doi:10.1093/jnci/djv036
rsID	chr_name	effect_allele	other_allele	effect_weight	locus_name	OR
rs78540526	11	T	C	0.16220387987485377	CCND1	1.1761
...

下载PGS后使用PLINK计算PGS

确认基因组版本等信息无误后,结合手头的基因型数据,可以通过PLINK来计算PGS。(与手头基因型文件的variant ID不一致时需要重新匹配)

https://www.cog-genomics.org/plink/2.0/scoreGWASLab:多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)的后半部分介绍了计算方法。

plink2 --score <filename> [i] [j] [k] [{header | header-read}]
                   [{center | variance-standardize | dominant | recessive}]
                   ['no-mean-imputation'] ['se'] ['zs'] ['ignore-dup-ids']
                   [{list-variants | list-variants-zs}]
                   ['cols='<column set descriptor>]

需要注意的是,有时我们手里的插补后的dosage文件是分染色体的,而PGS模型文件通常包括所有染色体上的variants,这种情况下一般需要分染色体进行计算单纯的分数加和 (使用—score里的cols=+scoresums 选项),然后再把22条染色体的分数再次加和算得总分数。

PLINK2 score的输出文件表头

Header	Column set	Contents
FID	maybefid, fid	Family ID
IID	(required)	Individual ID
SID	maybesid, sid	Source ID
PHENO1	pheno1	All-missing phenotype column, if none loaded
<Pheno name>, ...	pheno1, phenos	Phenotype value(s) (only first if just 'pheno1')
ALLELE_CT	nallele	Number of alleles across scored variants
DENOM	denom	Denominator used for score average
NAMED_ALLELE_DOSAGE_SUM	dosagesum	Sum of named allele dosages
<Score name>_AVG, ...	scoreavgs	Score averages
<Score name>_SUM, ...	scoresums	Score sums #分染色体时使用这一列求和

参考

https://www.pgscatalog.org/

Samuel A. Lambert, Laurent Gil, Simon Jupp, Scott C. Ritchie, Yu Xu, Annalisa Buniello, Aoife McMahon, Gad Abraham, Michael Chapman, Helen Parkinson, John Danesh, Jacqueline A. L. MacArthur and Michael Inouye.

The Polygenic Score Catalog as an open database for reproducibility and systematic evaluation

Nature Geneticsdoi: 10.1038/s41588-021-00783-5 (2021).

https://www.ebi.ac.uk/efo/

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

本文内容

  1. PRS系列回顾
  2. PRS-CS简介
  3. 概念框架
  4. 使用方法
  5. 参考

系列回顾:

PRS-CS简介

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

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

概念框架

  1. 考虑一个贝叶斯高维回归框架
  2. 常用的beta的先验分布通常可以表示为多个正态分布的比例混合
  3. 为了构建更为灵活的模型以更好反映遗传结构,一些方法采用了两个或多个质点或分布的离散混合,相比单纯的正太先验分布,这类方法通常能够构建更广的效应量分布。例如LDpred。
  4. 但此类方法后验推断需要极大的计算量,可能不能准确的建模出局部的LD结构。
  5. 为了解决这些问题,PRS-CS的作者采用了一种与先前离散混合概念不同的先验分布,也就是连续收缩先验分布,也就是正态分布的全局-局部比例混合
  6. 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

输入准备

  1. 上面下载的LD参考面板 plink格式
  2. 目标群体的bim文件:该文件只是用来提供target数据集中的SNP列表
  3. 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.

多基因风险分数 PRS( Polygenic risk score)系列之七:Pathway-based PRS 通路PRS

本文是多基因风险分数PRS系列的第七篇文章,主要介绍基于通路的PRS(Pathway-based PRS)的概念,以及简要的计算方法(仅供参考)。

回顾:

基于通路的PRS (Pathway-based PRS )简介

目前主流PRS计算模型的局限

当前主流的PRS计算方法都是基于疾病的经典多基因模型,也就是假设个体处在一个从低到高风险的线性分布上,其PRS作为易感性(liability)的估计值,能够总结这个个体的遗传概况。

易感性:易感性-阈值模型 Genetic liability, Threshold model

但这种方法也存在显著的缺陷,那就是这种概括性的估计值不能完全捕捉个体遗传特征,存在明显的信息丢失,一个最显著的例子就是无法判别不同生物学过程与通路之间遗传风险大小的变化,但这类信息通常在PRS的应用中(例如病人分层,或是治疗反应预测)起到重要的作用。

为了解决上述问题,近年来基于通路的PRS方法逐渐获得关注,本文将主要介绍Paul O’Reilly课题组开发的Pathway-based PRS方法PRSet。PRSet方法可以解释基因组中的亚结构(genomic sub-structure),该模型在经典多基因模型中进行扩展,使其能更好反映疾病的异质性。

基于通路的PRS : Pathway-based PRS

经典PRS模型计算整个基因组上risk allele与其效应估计乘积的加和,而pathway-based PRS只分别纳入k个相关通路上的SNP。也就是说,使用pathway-based PRS方法时,一个个体会有k个不同通路的通路PRS,每个通路PRS都对应基因组上一个具体的通路。(如图a)

在这里,通路可以由多种方式定义,可以是一些经典的通路数据库(如:KEGG,REACTOME等),也可以是通过基因共表达,蛋白质互作,实验扰动的功能产出等分析得到的模组。

理想状态下,通路能够反映不同生物学功能的编码。这些功能是可分割的,就像流行病学研究中,不同环境风险因素(抽烟,饮食,等等)也是分割的一样。基于这种观点,我们可以吧GWAS的结果也视为不同通路编码的生物学功能的信号的综合。(如图b)

PRSet的通路PRS方法简要原理

PRSet采用传统的C+T方法,分别对k条通路分别计算PRS,mk为每条通路上SNP数量,对于每个个体i,计算后会得到对应k个通路的k个PRS:

Clumping

与传统C+T对基因组上所有SNP进行clumping的方法不同,PRSet独立地对每条通路上的SNP进行clumping,确保能最大程度的保留信号。PRSet采用如下所示的bit-flag方法,对不同通路独立clumping。简单来说就是只有indexSNP包含在该gene set的时候才会去除掉需要被clump的SNP,否则会保留。

PRSet软件链接:https://www.prsice.info/prset_detail/

注:该方法的优势是简单快速,而且已经在软件中实现,作为探索可以尝试,但具体效果还有待检验。

(20211211:目前该文章还没有正式发表,文章里的检验方法不能令人信服,所以仅供参考)

参考:https://www.prsice.info/prset_detail/

https://europepmc.org/article/ppr/ppr362752

多基因风险分数 PRS( Polygenic risk score)系列之六:metaGRS介绍

目录

  1. PRS回顾
  2. 为什么要构建metaGRS
  3. metaGRS构建方法
  4. 总结
  5. 参考

回顾

  1. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之一:概念入门
  2. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)
  3. GWASLab:多基因风险分数 PRS( Polygenic risk score)系列之三:使用PRSice计算PRS(C+T方法)
  4. ldpred2 (预留链接)
  5. prs-cs(预留链接)

为什么要构建metaGRS

以往的PRS研究中,PRS模型通常以单一的GWAS概括性数据为基础进行计算,但此方法往往不能达到较高的预测能力,主要是因为

  1. 有限的样本量
  2. 目标表型较为显著的异质性
  3. 对基因组覆盖的不完全
  4. 基因定型与插补存在较大不确定性

等原因。

在存在上述问题的情况下使用单一PRS进行预测会降低准确度,为了解决以上描述的问题,Michael Inouye等人最早在冠状动脉疾病(CAD)的预测中,引入了metaGRS的概念,类似于meta分析,meGRS将多个GRS进行合并得到metaGRS,以提升其风险预测能力。

metaGRS构建方法

本文以 Abraham,G等人在2019年发表的论文为例,简要介绍metaGRS的构建方法:

第一步:常规的单一PRS的计算

第二步:使用弹性网络来构建metaGRS,具体步骤如下所示

例如,Abraham,G等人构建的对IS缺血性中风的metaGRS中,纳入了如下的GRS,包括了AS(所有中风),以及中风的各个亚型(IS, SVS, LAS, CES),以及各个中风的风险因子(例如血压,BMI,2型糖尿病,冠状动脉疾病等)的GRS

注:弹性网络 elastic net:

在损失函数中同时加入L1 (Lasso回顾)与L2(Ridge回归)的正则项:

该计算可以通过R包glmnet实现

https://cran.r-project.org/web/packages/glmnet/index.html

第三步:最后和单一PRS时一样,需要使用独立样本对metaGRS进行验证


总结

目前使用metaGRS的研究还比较少,且集中于特定的表型(CAD和IS),但可以预计随着PRS研究数量不断增加,在未来metaGRS的应用也会随之增加。

参考

Inouye, Michael, et al. “Genomic risk prediction of coronary artery disease in 480,000 adults: implications for primary prevention.” Journal of the American College of Cardiology 72.16 (2018): 1883-1893.

Abraham, G., Malik, R., Yonova-Doing, E. et al. Genomic risk score offers predictive performance comparable to clinical risk factors for ischaemic stroke. Nat Commun 10, 5819 (2019). https://doi.org/10.1038/s41467-019-13848-1

Friedman, Jerome, Trevor Hastie, and Rob Tibshirani. “glmnet: Lasso and elastic-net regularized generalized linear models.” R package version 1.4 (2009): 1-24.

多基因风险分数 PRS( Polygenic risk score)系列之三:使用PRSice计算PRS(C+T方法)

前情回顾:

本文内容:

1 PRSice-2简介与下载
2 PRSice-2教程
3 输入文件
4 C+T 参数
5 PRS计算
6 经验p值计算
7 示例演示
8 参考

PRSice-2简介与下载

PRSice-2是一款方便快捷的PRS计算软件,其主要功能是可以自动化执行一系列PLINK中用于PRS计算的功能,本文将简要讲解PRSice2并演示C+T方法的PRS计算。

PRSice下载链接:

https://www.prsice.info/

演示数据下载 (由教程原作者提供,与系列之二PLINK教程中使用的示例数据相同):

https://drive.google.com/file/d/1x_G0Gxk9jFMY-PMqwtg6-vdEyUPp5p5u/view


PRSice-2教程:

首先,需要的输入文件包括:

1 .base datatset GWAS概括性数据
2 .target dataset 目标数据集
3 .Phenotype files 表型文件
4. covariate files 协变量文件
5 .LD reference (optional: <500 samples) LD参考面板

1.base datatset

通常为GWAS的概括性数据,示例输入文件如下

如果base dataset的格式与示例不同,则可以通过以下选项自行指定

通过列名指定
--snp SNP --chr CHR --bp BP --A1 A1 --A2 A2 --stat OR --pvalue P
或者通过列数指定
--snp 0 --chr 1 --bp 2 --A1 3 --A2 4 --stat 5 --pvalue 7 --index

2.Target Dataset 目标数据集

目标数据集即为目标群体的基因型或填补后的基因型,目前支持PLINK的二进制格式,与BGEN格式

plink文件通过 -target 来指定,输入bed,bim,fam文件的前缀即可

bgen文件还需额外加上 -type bgen 选项

3 & 4 Phenotype & Covariates files

通过–pheno来指定表型文件,该文件是一个tab或space间隔的文本文件,前两列为FID与IID,第三列开始为表型或协变量,例如

FID	IID	Height
HG00096	HG00096	169.132168767547
HG00097	HG00097	171.256258630279
HG00099	HG00099	171.534379938588
HG00101	HG00101	169.850176470551
HG00102	HG00102	172.788360878389
HG00103	HG00103	169.862973824923
HG00105	HG00105	168.939248611414
HG00107	HG00107	168.972346393861
HG00108	HG00108	171.311736719186

如果表性文件中有多列,可以通过–pheno-col来指定表型的列名。

协变量文件格式与表型的一致。

5 LD reference

当样本量小于500时,建议额外指定LD参考面板, 输入格式与上述Target Dataset一致

--ld <LD refernce>

C+T 参数

Clumping

这里的clumping方法与plink的大体一致,可以参考:

SNP的LD剪枝与聚集 LD pruning & clumping

多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)

Clumping 参数可以通过 –clump-kb, –clump-r2 与–clump-p 选项来改变。

PRS计算

prsice2提供了多种计算PRS的模型与公式,默认的为计算均值,如下图所述

其他参数

可以指定运行所占用的内存,线程数,以及随机种子

--memory
--thread
--seed

经验p值计算

PRS的计算涉及到参数优化,那么必然会有过拟合的问题。

为了解决这一问题,通常有以下三种方法:

  1. 在独立的检验样本中评估
  2. 交叉验证
  3. 计算经验p值

PRSice2采用了第三种方法,通过permutation 来计算经验p值


PRS计算示例演示

使用之前下载的示例文件,执行以下命令:

Rscript PRSice.R \
    --prsice PRSice_linux \
    --base Height.QC.gz \
    --target EUR.QC \
    --binary-target F \
    --pheno EUR.height \
    --cov EUR.covariate \
    --base-maf MAF:0.01 \
    --base-info INFO:0.8 \
    --stat OR \
    --or \
    --out EUR

其中

prsice	PRSice_xxx	指定PRSice二进制文件路径 
base	        Height.QC.gz	 输入的base的GWAS summary statistic target	
EUR.QC	指定target的plink文件名的前缀 
binary-target	F	指定表型是否为二进制表型,F为否 
pheno	EUR.height	提供表型文件 
cov	EUR.covariate	        提供协变量文件 
base-maf	MAF:0.01	使用base文件中MAF列来过滤掉maf<0.01的变异 
base-info	INFO:0.8	使用base文件中INFO列来过滤掉info<0.8的变异 stat	
OR	base文件中包含效应量的列名 
or	-	指定效应量的形式为OR 
out	EUR	指定输出前缀为EUR

执行后log文件如下所示:

PRSice 2.3.3 (2020-08-05) 
<https://github.com/choishingwan/PRSice>
(C) 2016-2020 Shing Wan (Sam) Choi and Paul F. O'Reilly
GNU General Public License v3
If you use PRSice in any published work, please cite:
Choi SW, O'Reilly PF.
PRSice-2: Polygenic Risk Score Software for Biobank-Scale Data.
GigaScience 8, no. 7 (July 1, 2019)
2021-09-06 12:47:56
PRSice_linux \
    --a1 A1 \
    --a2 A2 \
    --bar-levels 0.001,0.05,0.1,0.2,0.3,0.4,0.5,1 \
    --base Height.QC.gz \
    --base-info INFO:0.8 \
    --base-maf MAF:0.01 \
    --binary-target F \
    --bp BP \
    --chr CHR \
    --clump-kb 250kb \
    --clump-p 1.000000 \
    --clump-r2 0.100000 \
    --cov EUR.covariate \
    --interval 5e-05 \
    --lower 5e-08 \
    --num-auto 22 \
    --or  \
    --out EUR \
    --pheno EUR.height \
    --pvalue P \
    --seed 3490962723 \
    --snp SNP \
    --stat OR \
    --target EUR.QC \
    --thread 1 \
    --upper 0.5

Initializing Genotype file: EUR.QC (bed) 

Start processing Height.QC 
================================================== 

Base file: Height.QC.gz 
GZ file detected. Header of file is: 
CHR	BP	SNP	A1	A2	N	SE	P	OR	INFO	MAF 

Reading 100.00%
499617 variant(s) observed in base file, with: 
499617 total variant(s) included from base file 

Loading Genotype info from target 
================================================== 

483 people (232 male(s), 251 female(s)) observed 
483 founder(s) included 

489805 variant(s) included 

Phenotype file: EUR.height 
Column Name of Sample ID: FID+IID 
Note: If the phenotype file does not contain a header, the 
column name will be displayed as the Sample ID which is 
expected. 

There are a total of 1 phenotype to process 

Start performing clumping 

Clumping Progress: 100.00%
Number of variant(s) after clumping : 193758 

Processing the 1 th phenotype 

Height is a continuous phenotype 
11 sample(s) without phenotype 
472 sample(s) with valid phenotype 

Processing the covariate file: EUR.covariate 
============================== 

Include Covariates: 
Name	Missing	Number of levels 
Sex	0	- 
PC1	0	- 
PC2	0	- 
PC3	0	- 
PC4	0	- 
PC5	0	- 
PC6	0	- 

After reading the covariate file, 472 sample(s) included in 
the analysis 

Start Processing
Processing 100.00%
There are 1 region(s) with p-value less than 1e-5. Please 
note that these results are inflated due to the overfitting 
inherent in finding the best-fit PRS (but it's still best 
to find the best-fit PRS!). 
You can use the --perm option (see manual) to calculate an 
empirical P-value. 

Begin plotting
Current Rscript version = 2.3.3
Plotting Bar Plot
Plotting the high resolution plot

默认输出结果中有两张图,以及两个文本文件:

第一张为柱状图可以看到不同P值阈值下PRS模型R2的大小:

第二种则是高精度图,表示所有p值阈值的模型拟合:

第一个后缀为prsice的文本文件各个C+T参数的模型的结果:

head EUR.prsice
Pheno	Set	Threshold	R2	P	Coefficient	Standard.Error	Num_SNP
-	Base	5e-08	0.0192512	0.000644758	563.907	164.163	1999
-	Base	5.005e-05	0.0619554	5.24437e-10	3080.12	485.501	7615
-	Base	0.00010005	0.0713244	2.27942e-11	3816.99	557.044	9047
-	Base	0.00015005	0.0785205	2.00391e-12	4362.01	603.601	10014
-	Base	0.00020005	0.0799241	1.24411e-12	4666.13	639.344	10756
-	Base	0.00025005	0.0820088	6.11958e-13	4947.48	668.217	11365
-	Base	0.00030005	0.0818422	6.47686e-13	5146.63	695.905	11884
-	Base	0.00035005	0.0836689	3.47353e-13	5392.8	720.237	12373
-	Base	0.00040005	0.0849879	2.21299e-13	5589.41	739.972	12817

第二个后缀为best的文本文件则为使用最优模型计算得到的target样本的PRS值

head EUR.best
FID IID In_Regression PRS
HG00096 HG00096 Yes -1.90234673e-05
HG00097 HG00097 Yes 1.27505253e-05
HG00099 HG00099 Yes -3.46628033e-06
HG00101 HG00101 Yes 9.02490969e-06
HG00102 HG00102 Yes 1.69546146e-05
HG00103 HG00103 Yes 1.91180157e-05
HG00105 HG00105 Yes 4.62832238e-06
HG00107 HG00107 Yes 7.59715813e-07
HG00108 HG00108 Yes 1.25453063e-05

参考:

https://www.prsice.info/

多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)

本文内容

简介
第1步 事先准备
第2步 clumping
第3步 使用PLINK计算PRS
第4步 寻找最佳的P值阈值
参考

简介:

PRS基础回顾 :多基因风险分数 PRS( Polygenic risk score)系列之一:概念入门

PLINK是群体遗传学研究中一款非常强大的软件,尽管PLINK并不是专门为计算PRS而开发,但其内置的功能足以使我们完成C+T (clumping + p value thresholding,也称P + T) 计算PRS方法必要的步骤,在初学阶段对于理解PRS计算过程有很大帮助。本文将简要介绍使用PLINK计算PRS的过程 (本教程改编自https://choishingwan.github.io/PRS-Tutorial/plink/ ):

本文所用数据(已经完成QC):

下载链接,由教程原作者提供:

https://drive.google.com/file/d/1x_G0Gxk9jFMY-PMqwtg6-vdEyUPp5p5u/view

也可以使用自己准备的数据,必要的数据包括

1 target样本的plink文件 (bed,bim,fam)
2 base样本的GWAS概括性数据 
3 target样本的表型数据

第1步 事前准备:更新效应量

当表型为二变量表型时,SNP的效应量通常表示为odds ratio(OR),这种情况下计算PRS时要计算乘积,这里为了简化我们将OR转化为beta,就可以只计算加和。

head Height.QC
CHR	BP	SNP	A1	A2	N	SE	P	OR	INFO	MAF
1	756604	rs3131962	A	G	388028	0.00301666	0.483171	0.997886915712650.890557941364774	0.369389592764921
1	768448	rs12562034	A	G	388028	0.00329472	0.834808	1.000687316093530.895893511351165	0.336845754096289
1	779322	rs4040617	G	A	388028	0.00303344	0.42897	0.997603556067569	0.897508290615237	0.377368010940814
1	801536	rs79373928	G	T	388028	0.00841324	0.808999	1.002035699227930.908962856432993	0.483212245374095
1	808631	rs11240779	G	A	388028	0.00242821	0.590265	1.001308325111540.893212523690488	0.450409558999587
1	809876	rs57181708	G	A	388028	0.00336785	0.71475	1.00123165786833	0.923557624081969	0.499743932656759
1	835499	rs4422948	G	A	388028	0.0023758	0.710884	0.999119752645200.906437735120596	0.481016005816168
1	838555	rs4970383	A	C	388028	0.00235773	0.150993	0.996619945289750.907716506801574	0.327164029672754
1	840753	rs4970382	C	T	388028	0.00207377	0.199967	0.997345678956140.914602590137255	0.498936220426316

可以利用简单的R脚本进行转换,在这里我们在文件的最后一列后添加转换后的beta

dat <- read.table(gzfile("Height.QC.gz"), header=T)
dat$BETA <- log(dat$OR)
write.table(dat, "Height.QC.Transformed", quote=F, row.names=F)
q() # exit R

将OR转换成beta后

head Height.QC.Transformed
CHR BP SNP A1 A2 N SE P OR INFO MAF BETA
1 756604 rs3131962 A G 388028 0.00301666 0.483171 0.997886915712657 0.890557941364774 0.369389592764921 -0.00211532000000048
1 768448 rs12562034 A G 388028 0.00329472 0.834808 1.00068731609353 0.895893511351165 0.336845754096289 0.000687079999998082
1 779322 rs4040617 G A 388028 0.00303344 0.42897 0.997603556067569 0.897508290615237 0.377368010940814 -0.00239932000000021
1 801536 rs79373928 G T 388028 0.00841324 0.808999 1.00203569922793 0.908962856432993 0.483212245374095 0.00203362999999797
1 808631 rs11240779 G A 388028 0.00242821 0.590265 1.00130832511154 0.893212523690488 0.450409558999587 0.00130747000000259
1 809876 rs57181708 G A 388028 0.00336785 0.71475 1.00123165786833 0.923557624081969 0.499743932656759 0.00123090000000354
1 835499 rs4422948 G A 388028 0.0023758 0.710884 0.999119752645202 0.906437735120596 0.481016005816168 -0.000880634999999921
1 838555 rs4970383 A C 388028 0.00235773 0.150993 0.996619945289758 0.907716506801574 0.327164029672754 -0.0033857799999997
1 840753 rs4970382 C T 388028 0.00207377 0.199967 0.99734567895614 0.914602590137255 0.498936220426316 -0.00265785000000019

第2步 Clumping (基于P值的LD-pruning)

回顾:SNP的LD剪枝与聚集 LD pruning & clumping

plink \
    --bfile EUR.QC \
    --clump-p1 1 \
    --clump-r2 0.1 \
    --clump-kb 250 \
    --clump Height.QC.Transformed \
    --clump-snp-field SNP \
    --clump-field P \
    --out EUR

其中

–bfile 是输入的PLINK格式文件(target样本)

–clump 则为GWAS的包含P值的概括性数据(base样本)

--clump-p1 1 \
	该选项依据P值,对索引SNP(每个clump中心的SNP)施加过滤,默认值为0.0001
--clump-r2 0.1 \
	与索引SNP的r2大于0.1的位点被归于与索引SNP为同一clump,在计算时去除
--clump-kb 250 \
	进行clumping时只考虑与索引SNP的距离小于250 kb的位点
--clump-snp-field SNP \
	概括性数据文件中变异ID的列名
--clump-field P \
	概括性数据文件中表示P值的列名

详细文档可参考:https://www.cog-genomics.org/plink/1.9/postproc#clump

执行以上命令后我们会得到文件 EUR.clumped, 该文件包含clumping后的索引SNP。

head EUR.clumped
CHR    F         SNP         BP        P    TOTAL   NSIG    S05    S01   S001  S0001    SP2
   6    1   rs3134762   31210866  4.52e-165      180     20      0      0      0    160 rs2523898(1),rs13210132(1),rs28732080(1),rs28732081(1),rs28732082(1),rs2517552(1),rs2517546(1),rs2844647(1),rs2517537(1),rs2844645(1),rs2523857(1),rs2517527(1),rs3131788(1),rs3132564(1),rs3130955(1),rs3130544(1),rs2535296(1),rs2517448(1),rs6457327(1),rs28732096(1),rs2233956(1),rs1265048(1),rs3130985(1),rs13201757(1),rs6909321(1),rs3130557(1),rs17196961(1),rs9501055(1),rs13200022(1),rs3130573(1),rs28732101(1),rs1265095(1),rs2233940(1),rs130079(1),rs746647(1),rs2240064(1),rs3131012(1)...

我们可以通过以下代码提取这些SNP的id, 用于下一步计算

awk 'NR!=1{print $3}' EUR.clumped > EUR.valid.snp

注意:

这里的r2计算是基于 最大似然单倍型频率估计 的,如果你的target样本比较少(<500),那么计算LD的r2时最好使用千人基因组项目中相应群体的数据。


第3步 计算 PRS

plink提供了—score选项,是我们可以方便的计算PRS

需要以下输入文件:

  1. base样本的GWAS概括性数据文件 Height.QC.Transformed
  2. 一个包含SNP ID与相应P值的文件
    使用awk命令从Height.QC.Transformed中提取SNPID与P值即可
    awk ‘{print $3,$8}’ Height.QC.Transformed > SNP.pvalue

在这里我们计算多个P值范围内SNP的PRS

echo "0.001 0 0.001" > range_list 
echo "0.05 0 0.05" >> range_list
echo "0.1 0 0.1" >> range_list
echo "0.2 0 0.2" >> range_list
echo "0.3 0 0.3" >> range_list
echo "0.4 0 0.4" >> range_list
echo "0.5 0 0.5" >> range_list

上述范围的三个数字分别为,名称,范围的下限,范围的上限

例如0.001 0 0.001 ,名称为0.001的这个范围,下限为0,上限为0.0001

下面开始计算,PLINK代码如下

plink \
    --bfile EUR.QC \
    --score Height.QC.Transformed 3 4 12 header \
    --q-score-range range_list SNP.pvalue \
    --extract EUR.valid.snp \
    --out EUR

其中 –q-score-range 即为包含P值范围的文件 与 一个包含SNP ID与相应P值的文件

–score 为计算分数的选项,3, 4,12 分别表示SNP ID,allele code与beta 分别在文件的第 3, 4, 12列, header表示输入文件有表头。

当然在计算时也可以加入协变量,例如PC等,这里就略过

执行后输出对应上述7个范围的PRS的文件:

EUR.0.5.profile
EUR.0.4.profile
EUR.0.3.profile
EUR.0.2.profile
EUR.0.1.profile
EUR.0.05.profile
EUR.0.001.profile

head EUR.0.001.profile
FID       IID  PHENO    CNT   CNT2    SCORE
  HG00096   HG00096     -9  32678   6894 -4.41148e-05
  HG00097   HG00097     -9  32678   6921 5.71558e-05
  HG00099   HG00099     -9  32678   6759 7.23017e-07
  HG00101   HG00101     -9  32678   6875 3.48096e-05
  HG00102   HG00102     -9  32678   6966 4.68491e-05
  HG00103   HG00103     -9  32678   6763 8.78172e-05
  HG00105   HG00105     -9  32678   6960 4.39891e-05
  HG00107   HG00107     -9  32678   6835 2.45655e-05
  HG00108   HG00108     -9  32678   6884 7.68976e-05

第4步 寻找 best-fit PRS

通常情况下,事先我们并不知道最优的P值阈值,所以在计算完成多组PRS后,为了找到最适阈值,需要对PRS进行回归分析,然后选取能够解释最多表型方差的P值阈值。

示例R代码如下 bestfit.R

p.threshold <- c(0.001,0.05,0.1,0.2,0.3,0.4,0.5)
# Read in the phenotype file 
phenotype <- read.table("EUR.height", header=T)
# Read in the PCs
pcs <- read.table("EUR.eigenvec", header=F)
# The default output from plink does not include a header
# To make things simple, we will add the appropriate headers
# (1:6 because there are 6 PCs)
colnames(pcs) <- c("FID", "IID", paste0("PC",1:6)) 
# Read in the covariates (here, it is sex)
covariate <- read.table("EUR.cov", header=T)
# Now merge the files
pheno <- merge(merge(phenotype, covariate, by=c("FID", "IID")), pcs, by=c("FID","IID"))
# We can then calculate the null model (model with PRS) using a linear regression 
# (as height is quantitative)
null.model <- lm(Height~., data=pheno[,!colnames(pheno)%in%c("FID","IID")])
# And the R2 of the null model is 
null.r2 <- summary(null.model)$r.squared
prs.result <- NULL
for(i in p.threshold){
    # Go through each p-value threshold
    prs <- read.table(paste0("EUR.",i,".profile"), header=T)
    # Merge the prs with the phenotype matrix
    # We only want the FID, IID and PRS from the PRS file, therefore we only select the 
    # relevant columns
    pheno.prs <- merge(pheno, prs[,c("FID","IID", "SCORE")], by=c("FID", "IID"))
    # Now perform a linear regression on Height with PRS and the covariates
    # ignoring the FID and IID from our model
    model <- lm(Height~., data=pheno.prs[,!colnames(pheno.prs)%in%c("FID","IID")])
    # model R2 is obtained as 
    model.r2 <- summary(model)$r.squared
    # R2 of PRS is simply calculated as the model R2 minus the null R2
    prs.r2 <- model.r2-null.r2
    # We can also obtain the coeffcient and p-value of association of PRS as follow
    prs.coef <- summary(model)$coeff["SCORE",]
    prs.beta <- as.numeric(prs.coef[1])
    prs.se <- as.numeric(prs.coef[2])
    prs.p <- as.numeric(prs.coef[4])
    # We can then store the results
    prs.result <- rbind(prs.result, data.frame(Threshold=i, R2=prs.r2, P=prs.p, BETA=prs.beta,SE=prs.se))
}
# Best result is:
prs.result[which.max(prs.result$R2),]
q() # exit R

执行后就可以看到本次PRS计算中,最优的阈值为0.3

Rscript bestfit.R 
  Threshold        R2            P     BETA       SE
5       0.3 0.1638468 1.025827e-25 47343.32 4248.152

参考:

Basic Tutorial for Polygenic Risk Score Analyses

多基因风险分数 PRS( Polygenic risk score)系列之一:概念入门

本文将讲解多基因风险分数PRS(Polygenic risk score,或称PGS) 的相关基础概念,

目录

  1. PRS的背景
  2. PRS的概念与定义
  3. PRS的一般性质
  4. 构建PRS的注意事项
  5. PRS的验证与预测
  6. 相关软件
  7. 参考

1. PRS的背景与概念

(首先复习: 遗传结构 Genetic architecture

一般情况下,对于单基因疾病(孟德尔遗传病)来说,只有单个或少数基因对表型有很大的影响,与之相对,对于复杂疾病来,通常有大量的遗传位点对表型有较小的影响,目前GWAS研究多基于此类无限小的假设(详见:解释复杂疾病的四种主流模型 CDCV/RAME/infinitesimal/Broad-sense-heritability ),这种情况下单个变异不足以用来评估个体对某一复杂疾病的风险,所以为了找到一个能够评估个体疾病风险的值,PRS (多基因风险评分)就应运而生,PRS的概念简单说就是,总和多个遗传变异与表型关系的数值。

2. PRS的一般定义

PRS (polygenic risk score, 多基因风险分数) , 对于非疾病的表型也称为 PGS (polygenic score)

实际研究中,PGS 的数学上的定义一般如下所示:

PGS是基因组上与表型相关等位基因的加权线性组合 (权重通常为GWAS中估计的效应量)。

其中: i 表示第i个个体, j 为第j个SNP, wj为该SNP的权重,a则为第i个个体第j个SNP的等位基因拷贝数

这里要注意:

通常PGS假设潜在的模型是加性模型 (additive model)

上述式子是一个概念性的定义式,实际操作中还需要进行额外操作。

3 PRS的一般性质

PRS可以被认为是多个独立的遗传信号的总和,那么根据中心极限定理,PRS也近似服从正态分布。

4.构建PRS的注意事项

4.1 GWAS discovery阶段的样本量要足够大

大的样本量:

好在从GWAS Diversity Monitor (https://gwasdiversitymonitor.com/)上可以看到,GWAS的样本量是在逐年上升的,目前规模最大的GWAS的样本量已经达到了三百万的级别,这将在未来有效促进PRS的构建。

4.2 选择纳入计算的SNP

这包含了两方面因素,1是纳入计算的SNP的数量如何决定,2是对于纳入的SNP如何施加权重。通常情况下这两方面的选择取决于研究的表型,与应用的类型。

目前主要的方法包括 :

GWAS中对SNP的检验通常是逐个进行的,由于LD的存在,这会使得SNP的效应估计值有偏差,继而导致PRS出现偏倚。为了减弱这种偏差目前有两种主流方法:

4.2.1 p值 clumping + thresholding法 (P+T 或 C+T) :

一种常用的方法就是在PRS的计算中只纳入一部分SNP,也就是先进行clumping(基于p值的pruning) (详见:SNP的LD剪枝与聚集 LD pruning & clumping),筛出各个loci里p值最低的SNP,然后再基于p值的某个阈值,选择纳入的SNP。

4.2.2 beta 缩减法

与第一种纳入部分SNP的思路不同,第二种方法是纳入所有的SNP,但在计算时会基于LD信息调整SNP的权重,例如LASSO回归(lassosum),与一些基于贝叶斯方法的算法 (LDpred等)。

5 PRS的验证与预测

5.1 要注意的是,在PRS研究中要使用独立的样本,也就是在GWAS discovery阶段使用的样本与PRS的目标样本之间不应该有重复。这主要是为了避免overfitting 过拟合的问题。只有当样本量足够大时,才可以使用同一样本。

5.2 目标样本应当为同一族裔。

由于不同族裔之间,MAF,局部LD等等的不同,PRS的泛用性较差。例如由BBJ计算而得的二型糖尿病PRS模型,应用到UKBB的人群中时,预测的r2明显降低。

5.3 提升PRS跨族裔泛用性

值得注意的是,目前也有研究致力于提升PRS在不同族裔间的泛用性。例如Amariuta等就基于转录因子介导的细胞特异的调节的位点的功能注释 (Functional annotations marking the precise location of TF-mediated cell-type-specific regulation )来降低群体特异的LD结构(population-specific LD),以提升PRS在不同族裔间的准确度。

6 目前使用较多的PRS软件包括:

PLINK 多基因风险分数 PRS( Polygenic risk score)系列之二:使用PLINK计算PRS(C+T方法)

PRSice 多基因风险分数 PRS( Polygenic risk score)系列之三:使用PRSice计算PRS(C+T方法)

LDpred 1/ 2

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

Lassosum

等等,我会在后续的文章中介绍具体使用方法。

其他PRS相关文章:

7 参考:

Choi, S. W., Mak, T. S. H. & O’Reilly, P. F. Tutorial: a guide to performing polygenic risk score analyses. Nature Protocols 15, 2759–2772 (2020).

McCarthy, M., Abecasis, G., Cardon, L. et al. Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nat Rev Genet 9, 356–369 (2008).

Amariuta, T. et al. Improving the trans-ancestry portability of polygenic risk scores by prioritizing variants in predicted cell-typespecific regulatory elements. Nature Genetics 52, 1346–1354 (2020).