2020-05-01 多基因风险评分PRS分析

参考学习资料：

https://mp.weixin.qq.com/s/ou9MfogMJJV00S6lZpNGqQ

多基因风险评分PRS分析

目前的主流的 GWAS 分析使用单一变异关联分析，但这一般需要非常大的样本量才可以进行。

相比之下，PRS（多基因风险评分）分析并非仅识别单个 SNP，而是将整个基因组中的遗传风险汇总到一个多基因评分中（ 简化示例请参见图1 ）。

在这种方法中，需要大量的检测样本才能确定每个 SNP 的权重。随后，在样本量较小的独立目标样品中，基于遗传 DNA 图谱和这些权重计算多基因得分。根据经验，需要大约 2,000 名受试者的目标样本才有足够的 power 来识别所解释的很大一部分方差。此外，观测数据集和目标样本应具有相同数量的样本。如果有更多样本可用，这些样本应该包括在观测样本中，以最大程度地估计效应量。

虽然 PRS 不足以在个体水平上预测疾病风险，但它已成功用于显性性状内和性状间的重要关联。

例如，对精神分裂症的 PRS 分析首次发现，根据常见 SNP（来自观测样本）的影响估算出的发展为精神分裂症的遗传风险的总体表现为与疾病发生风险显著相关，同时在目标样本中也表现与精神分裂症风险相关。尽管可用的样本量太小无法检测出全基因组范围内的重要SNP，但通过这种分析仍发现了显著的相关性。

此外，用于精神分裂症的 GWAS（观测样本）分析已被用于预测具有多种表型的目标样本中的遗传风险，例如双相情感障碍等等。

image-20200429151402950.png

图1 将三个单核苷酸多态性（SNP）汇总为一个多基因风险评分（PRS）的例子。权重可以是 beta 或比值比的对数，具体取决于分析的是数量性状还是二元性状。

为了进行 PRS 分析，从发现集 GWAS 获得特定表型的权重（数量性状的 beta 值和二元性状的比值比的对数）。在目标样本中，将根据他或她携带的风险等位基因数量的加权总和乘以特定性状权重，为每个个体计算 PRS。对于许多复杂的性状，SNP 效应大小已经公开（例如，参阅 https://www.med.unc.edu/pgc/downloads）。

尽管原则上所有常见的 SNP 都可以用于 PRS 分析中，但习惯上先将 GWAS 结果 clump（ 参见 clumping ），然后再计算风险评分。p 值的阈值通常用于删除几乎没有或没有关联证据的 SNP（ 例如，仅保留 p 值<0.5 或 <0.1 的SNP 。通常，我们会执行多个 PRS 分析，选取多个阈值。

进行多基因风险预测分析

一旦为目标样品中的所有个体计算了 PRS，就可以将这些分数用于（逻辑）回归分析中，以预测任何预期与目标性状表现出遗传相似的性状。预测精度可以通过回归分析的 (pseudo‐)R^2 表示。重要的是，在回归分析中至少包含一些 MDS 成分作为协变量，以校正群体分层。为了估计有多少变化可以由 PRS 所解释，可以把包括协变量（例如，MDS组分）模型的R^2 和包括协变量+PRS 的模型 R^2进行比较。若由于 PRS 而导致 R^2 增加就说明遗传风险因素增加了预测的精度。

PRS 的预测准确性主要取决于所分析特征的（共）遗传性，SNP 数量和discovery 样本的大小。目标样本的大小仅影响 R^2 的可靠性，如果感兴趣特征的（共）遗传性以及使用的 discovery 样本的大小足够大，即使是样本量为几千的目标样本也可以获得一个显著的 R^2。PRS 计算 power 的代码参考：https://sites.google.com/site/fdudbridge/software 。

PRSice[1] 是一个进行 PRS 分析的程序。它会进行 clumping，设定 p 值阈值，MDS 成分，并绘制图表。下面的教程将介绍用 PRSice 进行 PRS 分析，参考 https://github.com/MareesAT/GWA_tutorial/ （4_PRS.doc）。其他可进行 PRS 分析的软件，例如 PLINK（‐‐score）和 LDpred 。

PRSice-2 安装

#wget https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_linux.nightly.zip
#unzip PRSice_linux.nightly.zip
OSX 64-bit: 
wget https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_mac.nightly.zip
# Windows 32-bit: https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_win32.nightly.zip
# Windows 64-bit: https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_win64.nightly.zip

unzip PRSice_mac.nightly.zip

# 下载需要安装的 R 包
Rscript PRSice.R --dir .

使用 PRSice 进行 PRS 分析

这里用 PRSice 自带的示例数据为例。

  • 二元性状
Rscript PRSice.R --dir . \
    --prsice ./PRSice \
    --base TOY_BASE_GWAS.assoc \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --stat OR \
    --binary-target T

按照教程进行探索：

(base) bogon:PRSice_mac chelsea$ Rscript PRSice.R --dir . --prsice ./PRSice_mac   --base TOY_BASE_GWAS.assoc     --target TOY_TARGET_DATA     --thread 1     --stat OR     --binary-target T
PRSice 2.2.13 (2020-03-10)
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)
2020-04-29 16:01:27
./PRSice_mac \
    --bar-levels 0.001,0.05,0.1,0.2,0.3,0.4,0.5,1 \
    --base TOY_BASE_GWAS.assoc \
    --binary-target T \
    --clump-kb 250kb \
    --clump-p 1.000000 \
    --clump-r2 0.100000 \
    --interval 5e-05 \
    --lower 5e-08 \
    --num-auto 22 \
    --or  \
    --out PRSice \
    --seed 2741587729 \
    --stat OR \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --upper 0.5

Initializing Genotype file: TOY_TARGET_DATA (bed)

Start processing TOY_BASE_GWAS
==================================================

Reading 100.00%
Base file: TOY_BASE_GWAS.assoc
91062 variant(s) observed in base file, with:
2226 variant(s) located on haploid chromosome
88836 total variant(s) included from base file

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

2000 people (1024 male(s), 976 female(s)) observed
2000 founder(s) included

2226 variant(s) not found in previous data
88836 variant(s) included

There are a total of 1 phenotype to process

Start performing clumping

Clumping Progress: 100.00%

Number of variant(s) after clumping : 88836

Processing the 1 th phenotype
Phenotype is a binary phenotype
1000 control(s)
1000 case(s)

Preparing Output Files

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.2.12
Plotting Bar Plot
Plotting the high resolution plot
(base) bogon:PRSice_mac chelsea$ ls
PRSice.R                            PRSice_mac
PRSice.best                         PRSice_mac.Rproj
PRSice.log                          TOY_BASE_GWAS.assoc
PRSice.prsice                       TOY_TARGET_DATA.bed
PRSice.summary                      TOY_TARGET_DATA.bim
PRSice_BARPLOT_2020-04-29.png       TOY_TARGET_DATA.fam
PRSice_HIGH-RES_PLOT_2020-04-29.png TOY_TARGET_DATA.pheno

这个运行结果非常赞，写论文的时候可以直接使用这里的描述了！！！

  • 数量性状
Rscript PRSice.R --dir . \
    --prsice ./PRSice \
    --base TOY_BASE_GWAS.assoc \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --stat BETA \
    --beta \
    --binary-target F

运行结果如下：

(base) bogon:PRSice_mac chelsea$ Rscript PRSice.R --dir . --prsice ./PRSice_mac --base TOY_BASE_GWAS.assoc --target TOY_TARGET_DATA --thread 1 --stat BETA --beta
PRSice 2.2.13 (2020-03-10)
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)
2020-04-29 16:04:46
./PRSice_mac \
    --bar-levels 0.001,0.05,0.1,0.2,0.3,0.4,0.5,1 \
    --base TOY_BASE_GWAS.assoc \
    --beta  \
    --binary-target F \
    --clump-kb 250kb \
    --clump-p 1.000000 \
    --clump-r2 0.100000 \
    --interval 5e-05 \
    --lower 5e-08 \
    --num-auto 22 \
    --out PRSice \
    --seed 2105946812 \
    --stat BETA \
    --target TOY_TARGET_DATA \
    --thread 1 \
    --upper 0.5

Error:
Execution halted

不知出了什么问题，程序运行过程中断，也没有读懂提示，可能是缺了什么？

base 参数是指具有基本 base 样本（也称为发现样本或训练样本）的统计信息文件。 这些汇总统计信息包含每个遗传变异的效应值和 p 值。target 参数是指 target 样本二进制 plink 格式的基因型数据前缀，也支持 BGEN 格式。 base 样本和 target 样本也称为验证样本和测试样本，target 样本应独立于 base 样本。

为简单起见，我们在此分析中未包括主成分以及其他协变量，但在进行您自己的分析时，我们强烈建议添加协变量进行分析。

解释结果

默认情况下，PRSice 会输出两个图和几个文本文件。第一个图是PRSice_BARPLOT_.png （图2）。此图显示了不同阈值得到的关联结果对应的 R2 分布。 如图2 所示，使用 p 值 0.4463 的模型在目标样品中的 p 值最高，为 4.7*10-18 。 但通常样本量较小的多基因风险评分分析的预测值相对较低（ R2 约为5％）。

PRSice_BARPLOT_2020-04-29.png

图2 PRSice 条形图

第二个图是 PRSice_HIGH-RES_PLOT_ .png（图3）显示了许多不同的 p值阈值，以及对应 R2 的 p 值（黑点），绿线为趋势线。y 轴最高点代表最优解。

PRSice_HIGH-RES_PLOT_2020-04-29.png

图3 PRSice 高分辨率图

这两个图均表明，许多影响 base 样本性状的 SNP 可用于预测 target 样品中的性状。 注意，这两个性状可以相同也可以不同。 如果使用相同的性状，则预测值与该性状的遗传性（以及 base 样本的样本量）有关。 如果分析不同的性状，则预测值与两个性状之间的遗传重叠有关。无论哪种方式，多基因风险评分分析通常都表明，宽松 p 值阈值的模型通常比更严格阈值的模型有更好的预测能力，这表明许多统计学上无关紧要的 SNP 在多基因性状上具有预测价值。

文本文件则详细记录了每个 p 值的 R2。

# 不同阈值对应的 SNP 数量，SNP 解释度，回归系数
(base) bogon:PRSice_mac chelsea$ head PRSice.prsice
Set Threshold   R2  P   Coefficient Standard.Error  Num_SNP
Base    0.00025005  0.0133696   8.43169e-06 -0.197266   0.0442903   2
Base    0.00030005  0.00824473  0.000456434 -0.225204   0.0642503   3
Base    0.00040005  0.0089725   0.000256089 -0.350267   0.0958035   5
Base    0.00045005  0.0101339   0.000102845 -0.445497   0.114707    6
Base    0.00065005  0.00532975  0.004775    -0.402003   0.142462    8
Base    0.00070005  0.00876654  0.00030122  -0.549246   0.151967    9
Base    0.00080005  0.00233607  0.061455    -0.369219   0.197422    13
Base    0.00085005  0.00153157  0.129826    -0.342923   0.226384    15
Base    0.00095005  0.000124324 0.665873    -0.100725   0.233258    16

# 最佳阈值的分析结果
(base) bogon:PRSice_mac chelsea$ head PRSice.summary
Phenotype   Set Threshold   PRS.R2  Full.R2 Null.R2 Prevalence  Coefficient Standard.Error  P   Num_SNP
-   Base    0.4463  0.0520082   0.0520082   0   -   86.288  9.96331 4.69368e-18 36759

# 每个个体的分值
(base) bogon:PRSice_mac chelsea$ head PRSice.best
FID IID In_Regression PRS
CAS_1 CAS_1 Yes -0.00599501328
CAS_2 CAS_2 Yes -0.00631017938
CAS_3 CAS_3 Yes -0.00227495325
CAS_4 CAS_4 Yes -0.00204360007
CAS_5 CAS_5 Yes -0.000830676955
CAS_6 CAS_6 Yes -0.00224943517
CAS_7 CAS_7 Yes -0.000687589983
CAS_8 CAS_8 Yes -0.00413102565
CAS_9 CAS_9 Yes 0.00256661049

名词解释

  • Clumping ：识别并选择每个LD block中最显著的 SNP（即 p 值最低）以进行进一步分析。这样可以减少 SNP 之间的相关性，同时保留具有最强统计证据的 SNP。

  • Co-heritability ：疾病之间遗传关系的量度。 基于 SNP 的共遗传性是 SNP 对疾病效应之间协方差的比例。

  • Heterozygosity杂合性：对于特定 SNP 的两种不同等位基因的携带。个体的杂合率是杂合基因型的比例。个体内高水平的杂合性可能表明样品质量低，而低水平的杂合性可能是近亲繁殖所致。

  • Individual-level missingness 个体级别的缺失：特定个体缺少的 SNP 数量。较高的缺失表明 DNA 可能质量较差或存在技术问题。

  • Linkage disequilibrium（LD）连锁不平衡：给定种群中同一染色体上不同基因座等位基因之间非随机关联的一种度量。当其等位基因的关联频率高于随机分类下的预期频率时，SNP 位于 LD 中。LD 涉及 SNP 之间的模式。

  • Minor allele frequency （MAF）次要等位基因频率：这是在特定位置上出现频率最低的等位基因的频率。大多数研究的 power 不足以检测表型与 MAF 低的 SNP 的关联，因此需要过滤这些 SNP。

  • Population stratification 群体分层：研究中是否存在多个亚人群（例如，具有不同种族背景的个人）。由于等位基因频率在亚群之间可能不同，因此群体分层可能导致假阳性关联和/或掩盖真实关联。筷子基因就是一个很好的例子，由于群体分层的现象而导致得到 SNP 可以用来解释用筷子吃饭的习惯的结论[2]。

  • Pruning：选择近似连锁平衡的标记子集的方法。在 PLINK 中，此方法使用染色体特定窗口（区域）内 SNP 之间的 LD 强度，并根据用户指定的 LD 阈值仅选择不相关的SNP。与clumping 相比，这种方法过滤时不考虑 SNP 的 p 值。

  • Relatedness ：个体之间在遗传上有多强的关联性。常规的 GWAS 研究假定所有受试者都是无关的（即，没有任何一个个体比二级亲属更接近）。若数据集中包括亲属关系，不进行适当校正的话，可能对 SNP 效应大小的标准误的估计导致一定偏差。目前已有用于分析家族数据的特定工具。

  • Sex discrepancy 性别差异：表型信息中的性别与根据基因型确定的性别之间的差异。这个差异可以验证实验室中的样品是否混淆。注意，只有在计算了性染色体（X和Y）上的 SNP 时，才能进行此测试。

  • Single nucleotide polymorphism （SNP）单核苷酸多态性：在基因组中特定位置发生的单核苷酸（即A，C，G或T）变异。SNP 通常以两种不同的形式存在（例如，A与T）。这些不同的形式称为等位基因。包含两个等位基因的 SNP 有三种不同的基因型（例如，AA，AT和TT）。

  • SNP‐heritability：这是分析中一定集合内 SNP 解释的性状的表型变异分数。

  • SNP‐level missingness：这是样本中缺少特定 SNP信息的个体数量。具有高度缺失的SNP 可能导致误差。

  • 摘要统计数据：这些是进行 GWAS 后得到的结果，包括染色体位置， SNP 的位置，SNP（rs）标识符，MAF，效应大小（比值/β），标准误 和 p 值的信息。GWAS的摘要统计信息通常可以在研究人员之间免费访问或共享。

  • 哈迪-温伯格定律：指在理想状态下，各等位基因的频率在遗传中是稳定不变的，即保持着基因平衡。该定律运用在生物学、生态学、遗传学。条件：①种群足够大；②种群个体间随机交配；③没有突变；④没有选择；⑤没有迁移；⑥没有遗传漂变[3]。违反 HWE 法则表明基因型频率与预期值显着不同（例如，如果等位基因A 的频率= 0.20，等位基因T的频率= 0.80；基因型AT的预期频率为2 * 0.2 * 0.8 = 0.32），并且观察到的频率不应有显着差异。在 GWAS 中，通常假设与 HWE 的差异是基因分型错误的结果。病例中的 HWE 阈值通常不如对照组中的阈值严格，因为在病例中违反HWE法则可表明。

引用参考：

[1] http://www.prsice.info/step_by_step/
[2] https://www.nature.com/articles/4000662
[3] https://baike.baidu.com/item/%E5%93%88%E8%BF%AA-%E6%B8%A9%E4%BC%AF%E6%A0%BC%E5%AE%9A%E5%BE%8B/6265394?fr=aladdin

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 202,980评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,178评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 149,868评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,498评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,492评论 5赞 364
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,521评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,910评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,569评论 0赞 256
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,793评论 1赞 296
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,559评论 2赞 319
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,639评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,342评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,931评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,904评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,144评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,833评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,350评论 2赞 342