多基因风险评分（PRS）分析教程

翻译自：https://choishingwan.github.io/PRS-Tutorial/

需要安装：

1.R[1] (version 3.2.3+)2.PLINK 1.9[2]

多基因风险评分（Polygenic Risk Score）分析过程概览。PRS 分析需要两个输入数据集：i）base data（GWAS）：全基因组范围内遗传变异的基因型-表型关联的摘要统计信息（例如 beta，P值）；ii）target data：目标样本中个体的基因型和表型。基于 base data 得到的 SNP 效应值计算 target data 中样本的 PRS。

教程分为四个部分：

1.Base Data 的质控2.Target Data 的质控3.计算和分析 PRS4.可视化 PRS 结果

示例数据下载地址（需翻墙）：

# base datahttps://drive.google.com/file/d/1RWjk49QNZj9zvJHc9X_wyZ51fdy6xQjv/view?usp=sharing# target datahttps://drive.google.com/file/d/1uhJR_3sn7RA8U5iYQbcmTp6vFdQiF4F2/view?usp=sharing

我也把两个数据集传到了坚果云上，方便大家下载：https://www.jianguoyun.com/p/DQ53W3sQ757XCBik57AD

Base Data 的质控

PRS 分析的第一步就是得到 base data 的 GWAS 结果数据。查看 Height.gwas.txt.gz示例文件：

gunzip -c Height.gwas.txt.gz | head

CHR    BP    SNP    A1    A2    N    SE    P    OR    INFO    MAF1    756604    rs3131962    A    G    388028    0.00301666    0.4831710.997886915712657    0.890557941364774    0.3693895927649211    768448    rs12562034    A    G    388028    0.00329472    0.8348081.00068731609353    0.895893511351165    0.3368457540962891    779322    rs4040617    G    A    388028    0.00303344    0.428970.997603556067569    0.897508290615237    0.3773680109408141    801536    rs79373928    G    T    388028    0.00841324    0.8089991.00203569922793    0.908962856432993    0.4832122453740951    808631    rs11240779    G    A    388028    0.00242821    0.5902651.00130832511154    0.893212523690488    0.4504095589995871    809876    rs57181708    G    A    388028    0.00336785    0.714751.00123165786833    0.923557624081969    0.4997439326567591    835499    rs4422948    G    A    388028    0.0023758    0.7108840.999119752645202    0.906437735120596    0.4810160058161681    838555    rs4970383    A    C    388028    0.00235773    0.1509930.996619945289758    0.907716506801574    0.3271640296727541    840753    rs4970382    C    T    388028    0.00207377    0.1999670.99734567895614    0.914602590137255    0.498936220426316

包含以下几列：

•CHR：SNP 所在的染色体•BP：SNP 的座标•SNP：SNP ID，通常采用 rs-ID 的形式•A1：SNP 的 effect allele•A2：SNP 的 non-effect allele•N：用于计算效应值的样本数•SE：效应量的标准差•P：SNP 基因型与表型关联的 p 值•OR：SNP 的效应量（case-control）。若为连续的表型，那这一列通常为 beta•INFO：imputation information score•MAF：minor allele frequency

检查遗传度

我们建议 base data 计算得到的 h2 > 0.05（遗传度可用 LDSC 计算）。

Effect allele

一些 GWAS 结果文件没有明确哪个等位基因是效应等位基因，哪个是非效应等位基因。如果在计算 PRS 时进行了错误的假设，那么 PRS 对 target data 的效应估算将是错误的。

检查 GWAS 结果文件的完整性

另一个常见的问题是，下载的 base data 文件可能在下载过程中损坏，这可能导致 PRS 软件崩溃或在产生错误的结果。我们可用 md5sum检查文件的完整性：

md5sum Height.gwas.txt.gz

参考基因组

我们还需要检查 base data 和 target data 是否使用了相同的参考基因组。若使用了不同的参考基因组，则需要用 LiftOver[3] 进行坐标转换，使他们保持一致。

GWAS 结果 QC

Base data 和 target data 都应执行严格的质量控制步骤。低 MAF 和 INFO 的 SNPs 在执行下游分析之前通常会被去除。我们建议移除 MAF < 1% 和 INFO < 0.8 的 SNP。这一步可通过以下代码实现:

gunzip -c Height.gwas.txt.gz |awk 'NR==1 || ($11 > 0.01) && ($10 > 0.8) {print}' |gzip  > Height.gz

等位基因不匹配

在 base data 和 target data 中不匹配的 SNP 可通过“链翻转”进行匹配，例如某个 SNP 在 base data 中为 A/C，target data 中为 G/T，亦或者是一些不可解析的 SNP，例如在 base data 中为 C/G，target data 中为 C/T。大多数 PRS 软件会对可解析的 SNPs 进行自动链翻转，并删除不可解析的 SNPs。

因为我们需要 target data 来知道哪些 SNPs 具有不匹配的等位基因，所以我们将在 target data 中执行这种链翻转。

重复的 SNPs

大多数 PRS 软件不允许 base data 中存在重复的 SNPs，我们可以用下面的命令删除它们：

# 输出重复的 SNP IDgunzip -c Height.gz |awk '{ print $3}' |sort |uniq -d > duplicated.snp# 删除重复的 SNPgunzip -c Height.gz  |grep -vf duplicated.snp |gzip - > Height.nodup.gz

去除一些可能产生歧义的 SNPs

可以使用以下方法保留无歧义的 SNPs:

gunzip -c Height.nodup.gz |awk '!( ($4=="A" && $5=="T") ||         ($4=="T" && $5=="A") ||         ($4=="G" && $5=="C") ||         ($4=="C" && $5=="G")) {print}' |    gzip > Height.QC.gz

检查样本性别

有时样品标签会出现错误，最典型的问题就是通过 plink 计算出的性别与记录的性别有差。这部分请参阅 target data 的质控部分。

样本重复

由于示例的 target data 是模拟产生的，所以这里的 base data 和 target data 之间没有重叠的样本。

亲缘关系

base data 和 target data 中有亲缘关系的样本可能导致过拟合的结果。此部分可参考 target data 的质控部分。

最终，Height.QC.gz 文件可用于下游分析。

Target Data 的质控

# 解压数据unzip EUR.zip

样本量

我们建议 target data 的样本量大于 100。示例数据的样本量是 503 人。

基因型数据质控

plink     --bfile EUR     --maf 0.01     --hwe 1e-6     --geno 0.01     --mind 0.01     --write-snplist     --make-just-fam     --out EUR.QC

主要参数的含义如下：

•--bfile：输入基因型文件•--maf：删除 MAF 小于 0.01 的 SNPs•--hwe：删除 HWE p 值低于 1e-6 的 SNPs•--geno：排除大部分样本中缺失的 SNPs•--mind：删除基因型缺失率高的样本

执行 prunning 来删除高度相关的 SNPs:

plink     --bfile EUR     --keep EUR.QC.fam     --extract EUR.QC.snplist     --indep-pairwise 200 50 0.25     --out EUR.QC

主要参数的含义如下：

•--keep：只保留 EUR.QC.fam 中存在的样本•--extract：只保留 EUR.QC.snplist 中的 SNPs•--indep-pairwise：执行 prunning 的窗口大小为 200kb，每次步长为 50 个 SNP，并过滤掉 LD r2 大于 0.25 的 SNPs

这行命令将生成两个文件1) EUR.QC.prune.in 和2) EUR.QC.prune.out。

下一步我们需要排除样本杂合率均值偏离 ±3 SD 的个体。用 plink 计算杂合率：

plink     --bfile EUR     --extract EUR.QC.prune.in     --keep EUR.QC.fam     --het     --out EUR.QC

这行代码会生成 EUR.QC.het 文件。用 R 代码进行过滤:

library(data.table)# Read in filedat <- fread("EUR.QC.het")# Get samples with F coefficient within 3 SD of the population meanvalid <- dat[F<=mean(F)+3*sd(F) & F>=mean(F)-3*sd(F)] # print FID and IID for valid samplesfwrite(valid[,c("FID","IID")], "EUR.valid.sample", sep="t") q() # exit R

等位基因不匹配的 SNPs

在 base data 和 target data 中等位基因不匹配的 SNPs，可通过将等位基因翻转到它们的互补等位基因来解决。R 代码如下。

1.首先读入数据：

# magrittr allow us to do piping, which help to reduce the # amount of intermediate data typeslibrary(data.table)library(magrittr)# Read in bim file bim <- fread("EUR.bim") %>%    # Note: . represents the output from previous step    # The syntax here means, setnames of the data read from    # the bim file, and replace the original column names by     # the new names    setnames(., colnames(.), c("CHR", "SNP", "CM", "BP", "B.A1", "B.A2")) %>%    # And immediately change the alleles to upper cases    .[,c("B.A1","B.A2"):=list(toupper(B.A1), toupper(B.A2))]# Read in summary statistic data (require data.table v1.12.0+)height <- fread("Height.QC.gz") %>%    # And immediately change the alleles to upper cases    .[,c("A1","A2"):=list(toupper(A1), toupper(A2))]# Read in QCed SNPsqc <- fread("EUR.QC.snplist", header=F)

2.识别需要翻转 DNA 链的 SNPs

# Merge summary statistic with targetinfo <- merge(bim, height, by=c("SNP", "CHR", "BP")) %>%    # And filter out QCed SNPs    .[SNP %in% qc[,V1]]# Function for calculating the complementary allelecomplement <- function(x){    switch (x,        "A" = "T",        "C" = "G",        "T" = "A",        "G" = "C",        return(NA)    )} # Get SNPs that have the same alleles across base and targetinfo.match <- info[A1 == B.A1 & A2 == B.A2, SNP]# Identify SNPs that are complementary between base and targetcom.snps <- info[sapply(B.A1, complement) == A1 &                     sapply(B.A2, complement) == A2, SNP]# Now update the bim filebim[SNP %in% com.snps, c("B.A1", "B.A2") :=        list(sapply(B.A1, complement),             sapply(B.A2, complement))]

3.识别需要在 target data 中重新编码的 SNPs

# identify SNPs that need recodingrecode.snps <- info[B.A1==A2 & B.A2==A1, SNP]# Update the bim filebim[SNP %in% recode.snps, c("B.A1", "B.A2") :=        list(B.A2, B.A1)]# identify SNPs that need recoding & complementcom.recode <- info[sapply(B.A1, complement) == A2 &                     sapply(B.A2, complement) == A1, SNP]# Now update the bim filebim[SNP %in% com.recode, c("B.A1", "B.A2") :=        list(sapply(B.A2, complement),             sapply(B.A1, complement))]# Write the updated bim filefwrite(bim, "EUR.QC.adj.bim", col.names=F, sep="t")

4.确定在 base 和 target 上有不同等位基因的 SNPs (通常是由于使用不同参考基因组或 Indel造成的)

mismatch <- bim[!(SNP %in% info.match |                    SNP %in% com.snps |                    SNP %in% recode.snps |                    SNP %in% com.recode), SNP]write.table(mismatch, "EUR.mismatch", quote=F, row.names=F, col.names=F)q() # exit R

5.将 EUR.bim 替换为 EUR.QC.adj.bim：

# Make a back upmv EUR.bim EUR.bim.bkln -s EUR.QC.adj.bim EUR.bim

注意：大多数 PRS 软件将自动执行这种转换，所以一般不需要这一步。

重复的 SNPs

确保删除 target data 中重复的 SNPs (示例 target data 是模拟产生的，因此不包括重复的 SNPs)。

检查性别

我们可根据 X 染色体杂合/纯合率检查数据集中记录的性别与真实性别之间的差异。男性 X 染色体纯合度估计值 > 0.8，女性 <0.2 。

在执行性别检查之前，应先进行 pruning：

plink     --bfile EUR     --extract EUR.QC.prune.in     --keep EUR.valid.sample     --check-sex     --out EUR.QC

这将生成 EUR.QC.sexcheck 的文件，根据这个文件进行过滤，R 代码如下：

library(data.table)# Read in filevalid <- fread("EUR.valid.sample")dat <- fread("EUR.QC.sexcheck")[FID%in%valid$FID]fwrite(dat[STATUS=="OK",c("FID","IID")], "EUR.QC.valid", sep="t") q() # exit R

样本重叠

由于 target data 是模拟的，所以这里的 base data 和 target data 之间没有重叠的样本。

亲缘关系

Target data 中密切相关的个体可能导致过拟合，需要进行过滤。

plink     --bfile EUR     --extract EUR.QC.prune.in     --keep EUR.QC.valid     --rel-cutoff 0.125     --out EUR.QC

生成最终 QC 后的 target data 文件

plink     --bfile EUR     --make-bed     --keep EUR.QC.rel.id     --out EUR.QC     --extract EUR.QC.snplist     --exclude EUR.mismatch

计算和分析 PRS

用 plink 计算 PRS

根据前面的质控步骤，我们得到了以下六个文件：

•Height.QC.gz•EUR.QC.bed•EUR.QC.bim•EUR.QC.fam•EUR.height•EUR.covariate

更新效应量

这里为了简化计算，我们使用 OR 的自然对数，这样可以用求和来计算 PRS。

library(data.table)dat <- fread("Height.QC.gz")fwrite(dat[,BETA:=log(OR)], "Height.QC.Transformed", sep="t")q() # exit R

由于数值的四舍五入，使用 awk 来进行转换可能会导致不太精确的结果。因此，我们建议在 R 中执行转换，或者用 PRS 软件直接执行转换。

Clumping

尽管原则上所有常见的 SNP 都可以用于 PRS 分析中，但习惯上先将 GWAS 结果 clump，然后再计算风险评分。所谓 clumping 就是识别并选择每个 LD block 中最显著的 SNP（即 p 值最低）以进行进一步分析。这样可以减少 SNP 之间的相关性，同时保留具有最强统计证据的 SNP。

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

主要参数的含义如下：

•--clump-p1：index SNP 的 P 值，取 1 即是纳入所有的 SNPs•--clump-r2：r2 > 0.1 的 SNP 将被删除•--clump-kb：范围取 250kb•--clump：base data 的 summary statistic 结果文件•--clump-snp-field：指定包含 SNP ID 的列名•--clump-field：指定 P 值的列名

上面的命令将得到 EUR.clumped。我们可以通过执行以下命令提取 index SNP ID:

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

如果你的 target data 很小(比如 n < 500) ，那么你可以用千人基因组计划的数据来计算 LD。确保使用最接近反映 base 样本的人群数据。

计算 PRS

在 plink 中我们可使用 --score 和 --q-score-range 计算 PRS。

这里我们需要三个文件：

1.base data 文件：Height.QC.Transformed2.一个包含 SNP id 及其对应的 p 值的文件：

awk '{print $3,$8}' Height.QC.Transformed > SNP.pvalue

1.一个包含不同 p 值阈值的文件：

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

三列分别表示：

Threshold 名称	下限	上限

例如对于 0.05 的阈值，我们包括 p 值从 0 到 0.05 的所有 SNPs，包括 p 值等于 0.05 的 SNPs。

然后我们就可以用 plink 计算 PRS 了：

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

参数含义如下:

•--score Height.QC.Transformed 3 4 12 header：输入 Height.QC.Transformed 文件，第三列为 SNP ID，第四列为 effective allele，第十二列为效应量，且该文件包含 header。•--q-score-range range_list SNP.pvalue：定义 PRS 阈值文件，和 SNP P 值文件。

上面的命令和 range_list 将生成7个文件:

1.EUR.0.5.profile2.EUR.0.4.profile3.EUR.0.3.profile4.EUR.0.2.profile5.EUR.0.1.profile6.EUR.0.05.profile7.EUR.0.001.profile

考虑人群分层因素

人口结构是 GWAS 中混杂因素的主要来源，我们通常将主成分(PC)作为协变量进行校正。

用 plink 计算主成分：

# First, we need to perform prunningplink     --bfile EUR.QC     --indep-pairwise 200 50 0.25     --out EUR# Then we calculate the first 6 PCsplink     --bfile EUR.QC     --extract EUR.prune.in     --pca 6     --out EUR

我们可以用 LDSC 来挑选合适的主成分数量。使用不同数量的主成分进行校正，LDSC 截距最接近 1 的即是最佳的主成分数量。如果 base 样本和 target 样本是从世界各地不同的人群中收集的，PRS 分析的结果可能存在偏差。

寻找最佳的 PRS 阈值

为了接近“最佳拟合”的 PRS，我们可在一定 P 值范围内计算 PRS，然后选择能够解释最高表型差异的 PRS。这可以通过 R来实现：

library(data.table)library(magrittr)p.threshold <- c(0.001,0.05,0.1,0.2,0.3,0.4,0.5)phenotype <- fread("EUR.height")pcs <- fread("EUR.eigenvec", header=F) %>%    setnames(., colnames(.), c("FID", "IID", paste0("PC",1:6)) )covariate <- fread("EUR.cov")pheno <- merge(phenotype, covariate) %>%        merge(., pcs)null.r2 <- summary(lm(Height~., data=pheno[,-c("FID", "IID")]))$r.squaredprs.result <- NULLfor(i in p.threshold){    pheno.prs <- paste0("EUR.", i, ".profile") %>%        fread(.) %>%        .[,c("FID", "IID", "SCORE")] %>%        merge(., pheno, by=c("FID", "IID"))    model <- lm(Height~., data=pheno.prs[,-c("FID","IID")]) %>%            summary    model.r2 <- model$r.squared    prs.r2 <- model.r2-null.r2    prs.coef <- model$coeff["SCORE",]    prs.result %<>% rbind(.,        data.frame(Threshold=i, R2=prs.r2,                     P=as.numeric(prs.coef[4]),                     BETA=as.numeric(prs.coef[1]),                    SE=as.numeric(prs.coef[2])))}print(prs.result[which.max(prs.result$R2),])q() # exit R

用 PRSice-2 计算 PRS

在运行 PRSice-2 之前，我们还需要一个协变量文件，整合协变量信息和主成分信息。R 代码如下:

library(data.table)covariate <- fread("EUR.covariate")pcs <- fread("EUR.eigenvec", header=F)colnames(pcs) <- c("FID","IID", paste0("PC",1:6))cov <- merge(covariate, pcs)fwrite(cov,"EUR.cov", sep="t")q()

生成 EUR.cov 文件。

然后可以运行 PRSice-2，计算 PRS：

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

这行代码会自动寻找最佳的参数并进行可视化。

可视化 PRS 结果

我们可以用 ggplot2 进行可视化：

# ggplot2 is a handy package for plottinglibrary(ggplot2)# generate a pretty format for p-value outputprs.result$print.p <- round(prs.result$P, digits = 3)prs.result$print.p[!is.na(prs.result$print.p) &                       prs.result$print.p == 0] <-    format(prs.result$P[!is.na(prs.result$print.p) &                            prs.result$print.p == 0], digits = 2)prs.result$print.p <- sub("e", "*x*10^", prs.result$print.p)# Initialize ggplot, requiring the threshold as the x-axis (use factor so that it is uniformly distributed)ggplot(data = prs.result, aes(x = factor(Threshold), y = R2)) +    # Specify that we want to print p-value on top of the bars    geom_text(        aes(label = paste(print.p)),        vjust = -1.5,        hjust = 0,        angle = 45,        cex = 4,        parse = T    )  +    # Specify the range of the plot, *1.25 to provide enough space for the p-values    scale_y_continuous(limits = c(0, max(prs.result$R2) * 1.25)) +    # Specify the axis labels    xlab(expression(italic(P) - value ~ threshold ~ (italic(P)[T]))) +    ylab(expression(paste("PRS model fit:  ", R ^ 2))) +    # Draw a bar plot    geom_bar(aes(fill = -log10(P)), stat = "identity") +    # Specify the colors    scale_fill_gradient2(        low = "dodgerblue",        high = "firebrick",        mid = "dodgerblue",        midpoint = 1e-4,        name = bquote(atop(-log[10] ~ model, italic(P) - value),)    ) +    # Some beautification of the plot    theme_classic() + theme(        axis.title = element_text(face = "bold", size = 18),        axis.text = element_text(size = 14),        legend.title = element_text(face = "bold", size =                                        18),        legend.text = element_text(size = 14),        axis.text.x = element_text(angle = 45, hjust =                                       1)    )# save the plotggsave("EUR.height.bar.png", height = 7, width = 7)q() # exit R

我们还可以看看 PRS 与表型之间的关系：

library(ggplot2)# Read in the filesprs <- read.table("EUR.0.3.profile", header=T)height <- read.table("EUR.height", header=T)sex <- read.table("EUR.cov", header=T)# Rename the sexsex$Sex <- as.factor(sex$Sex)levels(sex$Sex) <- c("Male", "Female")# Merge the filesdat <- merge(merge(prs, height), sex)# Start plottingggplot(dat, aes(x=SCORE, y=Height, color=Sex))+    geom_point()+    theme_classic()+    labs(x="Polygenic Score", y="Height")q() # exit R

引用链接

[1] R: https://www.r-project.org/ [2] PLINK 1.9: https://www.cog-genomics.org/plink2 [3] LiftOver: https://genome.ucsc.edu/cgi-bin/hgLiftOver