多基因风险分数 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的计算涉及到参数优化，那么必然会有过拟合的问题。

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

在独立的检验样本中评估
交叉验证
计算经验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)系列之三：使用PRSice计算PRS（C+T方法）》有4个想法

Pingback：多基因风险分数 PRS( Polygenic risk score)系列之七：Pathway-based PRS 通路PRS – GWAS实验室 – GWASLab
Pingback：多基因风险分数 PRS( Polygenic risk score)系列之五：使用PRS-CS计算PRS（beta-shrinkage方法） – GWAS实验室 – GWASLab
Pingback：多基因风险分数 PRS( Polygenic risk score)系列之一：概念入门 – GWAS实验室 – GWASLab
Pingback：多基因风险分数 PRS( Polygenic risk score)系列之八：PGS Catalog – GWAS实验室 – GWASLab

发表评论取消回复

请在此输入您的评论... ...

Fill in your details below or click an icon to log in:

电子邮箱地址（必填） (Address never made public)

显示名称（必填）

网站地址

您正在使用您的 WordPress.com 账号评论。 ( 注销 / 更改 )

您正在使用您的 Facebook 账号评论。 ( 注销 / 更改 )

取消

Connecting to %s

共享此文章：

赞过：

相关

《多基因风险分数 PRS( Polygenic risk score)系列之三：使用PRSice计算PRS（C+T方法）》有4个想法

发表评论 取消回复

发表评论取消回复