多基因风险分数 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)系列之二:使用PLINK计算PRS(C+T方法)》有8个想法

发表评论

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

WordPress.com 徽标

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

Twitter picture

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

Facebook photo

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

Connecting to %s