多基因风险分数 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

需要以下输入文件：

base样本的GWAS概括性数据文件 Height.QC.Transformed
一个包含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)系列之二：使用PLINK计算PRS（C+T方法）》有8个想法

Pingback：多基因风险分数 PRS( Polygenic risk score)系列之三：PRSice – GWAS实验室 – GWASLab – GWASラボ
Pingback：多基因风险分数 PRS( Polygenic risk score)系列之三：PRSice – GWAS实验室 – GWASLab – GWASラボ
Pingback：多基因风险分数 PRS( Polygenic risk score)系列之七：Pathway-based PRS 通路PRS – GWAS实验室 – GWASLab
Liyan说道：

2022年4月1日下午9:19

请问，你有遇到过这样的报错嘛，谢谢：Error: –clump-p1 must be used with –clump.

赞赞

回复
1. Yunye说道：
  
  2022年4月2日下午9:20
  
  你好，看样子是没指定–clump，这里–clump-p1是指定p值阈值的，要用–clump 指定输入sumstats的文件名
  
  赞赞
  
  回复
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)系列之二：使用PLINK计算PRS（C+T方法）》有8个想法

发表评论 取消回复

发表评论取消回复