GWAS分析全流程：基因、遗传力、表型与质控解读

手上有一份两三百人的芯片数据，一份 Excel 表型表，老板丢过来一句话：“看看这个性状跟哪些位点有关系。”这大概是每个刚接触统计遗传的人都会遇到的场面。GWAS（全基因组关联分析）这个词听上去唬人，各种教程一上来就是曼哈顿图、QQ 图、P 值、λ 膨胀因子，把人绕得云里雾里。可真正拆开来看，它只干一件事：把测出来的基因型数据和观测到的表型数据对齐，逐个位点做回归，看哪个位点的等位基因频率在表型高低两组之间存在系统性差异。而遗传力这个概念，则是判断这件事值不值得做、能做到什么程度的那把尺子。

这篇东西写给两类人：一类是刚进组、被安排跑第一个 GWAS 的学生，手上数据还没理清楚，不知道从哪儿下手；另一类是做过一点分析、但结果老是不显著或者 QQ 图膨胀得离谱、想搞清楚问题出在哪的人。我会把基因、遗传力、表型这三个词先掰开讲透，再顺着一条完整的实操链路走一遍——从原始文件到质控、到关联分析、到结果解读和下游分析，中间夹着参数怎么定、命令怎么写、坑在哪里。所有参数和经验都是基于常见实践给出的参考值，你自己项目里该调的还得按实际情况调。

1. 先把三个词说透：基因、遗传力、表型到底在讲什么

很多人跑 GWAS 跑得很熟，但被问一句“你这个性状的遗传力是多少”就卡壳了。这三个概念不是背定义用的，它们直接决定了你后面每一步该怎么设计。花十分钟把它们捋清楚，比盲目跑十遍流程有用得多。

1.1 基因型数据长什么样：从碱基到 0/1/2

先说最基础的。人的基因组有三十亿个碱基对，其中大部分在人群里是一模一样的，真正有差异的位置叫多态性位点，最常见的是单碱基替换，也就是SNP。一个 SNP 上，人群里通常只有两种碱基在流传，比如某一位点上是 A 或者 G，那么每个人在这个位置就有三种可能的组合：AA、AG、GG。

GWAS 分析里，我们会把这种组合编码成一个数字。约定俗成是数“次等位基因”（minor allele）的个数：AA 记 0，AG 记 1，GG 记 2。这个 0/1/2 的编码就是所谓的加性编码，它背后假设的是“每多带一个次等位基因，表型平均变化 β”，这也是 GWAS 里最主流、最省事的模型。除此之外还有显性编码、隐性编码，但除非你有明确的生物学理由，一般不用。

实际拿到手的基因型数据，来源主要有三种：SNP 芯片、全外显子测序、全基因组测序。芯片便宜、通量大，一般覆盖几十万到两百多万个位点，是目前大规模队列最常用的；测序数据位点全，但贵，样本量上不去。不管哪种，通常都会做一步基因型填充（imputation），用参考面板把没测到的位点推测出来，把位点数量补到几百万甚至上千万。填充后的数据一般以 dosage（0 到 2 之间的小数）形式存储，因为你没法确定一个没测到的位点到底是 0 还是 1，只能给个概率期望。

提醒：填充质量有专门的指标，叫 INFO 分数或者 Rsq，一般建议保留 INFO > 0.8 的位点。低质量的填充位点会引入大量假阳性，尤其是 rare variant（低频位点），这部分数据宁可丢掉，别舍不得。

这里还有个绕不开的概念叫LD（连锁不平衡）。相邻的 SNP 不是独立遗传的，它们往往成块地一起传下去，所以一个位点显著，它周围一片位点往往也显著。这就是曼哈顿图上那些“峰”形成的原因，也是后面 fine-mapping 要处理的问题。LD 的强度用 r² 衡量，r² = 1 表示完全共变，r² = 0 表示独立。做位点剪枝（pruning）时常用的阈值是 r² < 0.2 或者 0.1，目的就是把这种冗余去掉。

1.2 遗传力：为什么身高能遗传，酒量却很难

遗传力回答的问题是：一个性状在人群中的变异，有多大比例是基因差异造成的。用公式写出来就是

h² = V_G / V_P

其中 V_G 是遗传因素贡献的方差，V_P 是表型总方差（等于遗传方差加环境方差）。这个数在 0 到 1 之间，越接近 1 说明这个性状越“由基因决定”。

但这里有个非常容易被混淆的地方，得说清楚。经典的遗传力来自双胞胎研究和家系研究，叫狭义遗传力，它估的是所有遗传效应（包括那些我们根本没测到的位点）的总贡献。而 GWAS 里算出来的是SNP 遗传力，只用芯片或测序实际捕获到的那些常见位点来估，数值通常明显低于经典遗传力。举个例子，身高的经典遗传力在 0.8 左右，但基于常见 SNP 估出来的 SNP 遗传力大概只有 0.4 到 0.5；BMI 的经典遗传力约 0.4 到 0.7，SNP 遗传力大致在 0.2 上下。这两者之间的差距，学名叫“缺失的遗传力”，一部分是低频位点和结构变异没被芯片捕获，一部分是基因间交互和非加性效应没被线性模型捕捉到。

计算 SNP 遗传力目前主流是两条路。一条是GREML，也就是用全基因组关系矩阵（GRM）配合限制性最大似然估计，把个体间的遗传相似度和表型相似度做回归，软件实现是 GCTA 那一套。另一条是LD Score 回归，只需要汇总统计数据（summary statistics），不用个体数据，特别适合拿公开的 GWAS 结果来算，软件是 LDSC。两条路的假设不同，结果也会有出入，GREML 在样本量小的时候更稳，LDSC 对样本量的要求宽松但对 LD 参考面板敏感。

实操心得：如果你手上只有汇总统计，别硬套 GREML，直接用 LDSC；如果有个体数据、样本量超过两千，GREML 的结果更可信。另外算 SNP 遗传力时，一定要控制好前几个主成分，不然群体分层会把遗传力估得虚高，我见过把 h² 算到 0.9 的结果，一查是分层没处理干净。

1.3 表型：定义错了，后面全白干

三个词里最容易被轻视、但实际上最要命的是表型。基因型数据是客观的，测出来是什么就是什么；表型却充满主观选择，你怎么定义，直接决定你能找到什么。

先分类。连续型表型，比如身高、体重、血压、某个代谢物的浓度，用线性回归。二分类表型，比如病例/对照、有无某种症状，用逻辑回归。此外还有计数型、生存时间型、有序分类，各有对应的模型，别一把梭全用线性回归。

连续型表型有几个处理细节，值得单独说。第一是离群值，一两个极端值能把回归系数带偏，常规做法是把超过均值 ± 5 个标准差的值设为缺失，或者做 winsorize 缩尾。第二是分布形态，很多生化指标严重右偏（比如某些酶浓度、抗体滴度），直接做线性回归会违反正态假设，通用的处理是rank-based inverse normal transformation（INT），也就是先排秩再转成正态分位数，这一步几乎是标准动作。第三是协变量，年龄、性别、BMI、测序批次、前若干个遗传主成分，这些该放进去的一定要放，漏掉一个混杂因素，出来的峰可能全是假的。

二分类表型更麻烦。病例的定义标准要统一，如果一部分病例是按临床诊断、一部分是按自报，那这个表型的噪声会非常大，效力直接打对折。病例和对照的比例也别太失衡，1:1 到 1:4 之间比较理想，1:10 以上效力会明显下降。还有一点，很多表型存在用药影响，比如服降压药的人血压被压下来了，直接拿测出来的血压做表型，会把真实的遗传效应稀释掉，正规做法是对服药人群的血压做校正（比如加一个固定增量），这个细节在论文方法学里经常被忽略，但对结果影响不小。

表型文件整理成什么样？最省事的格式是三列：样本 ID、表型值、协变量（可以多列）。缺失值用一个统一标记（比如 NA 或 -9），别用 0，因为 0 是合法数值，混在一起会让软件把它当成真实观测。样本 ID 的命名规则务必和基因型文件完全一致，这听起来像废话，但我见过太多人卡在这一步，两边 ID 一个是数字、一个是带前缀的字符串，对不上，白折腾半天。

2. GWAS 到底在算什么：模型、选型与效力

明白了三个词，接下来看 GWAS 的核心机制。它本质上就是一个大规模并行的回归分析，几百万个位点，每个都跑一遍，看谁显著。听起来暴力，但里面有几处关键选择会决定结果的可靠性。

2.1 一个位点一次回归，检验的到底是什么

以连续型表型为例，对每一个 SNP，我们拟合这样一个模型：

y = μ + β·G + γ₁C₁ + γ₂C₂ + … + e

y 是表型，G 是该位点的加性编码（0/1/2），C 是协变量（年龄、性别、主成分等），e 是残差。我们要检验的原假设是β = 0，也就是“这个位点的等位基因剂量和表型无关”。检验统计量是 β 除以它的标准误，得到一个 z 值，再转成 P 值。

把上面这个过程重复几百万次，就得到几百万个 P 值。多重检验校正是必须的，因为按 α = 0.05 去判定，几百万次检验里光靠随机就能冒出十几万个“显著”位点。全基因组显著性阈值通行的做法是 Bonferroni 校正：0.05 除以有效独立检验次数。人类基因组里独立位点大概一百万量级，所以阈值取5×10⁻⁸。这个数字现在基本成了行业默认，看到 P 值小于它就认为是全基因组显著。另外还有个“提示性显著”的说法，通常取 1×10⁻⁵，只当作线索，不能当结论用。

对于二分类表型，把线性回归换成逻辑回归：

logit(P(y=1)) = μ + β·G + γC

解释方式变成比值比（OR），β 取指数就是每多一个次等位基因对应的患病风险倍数。OR 一般都很小，1.05 到 1.3 是常见区间，别指望看到 2 以上的数，那基本说明你数据有问题。

2.2 模型怎么选：从最简单的 OLS 到混合模型

这一步是新手最容易迷路的地方。各种软件、各种模型，名字一堆，其实按复杂度排队就清楚了。

最朴素的是普通最小二乘 / 普通逻辑回归，不考虑任何样本间的相关性。它的优点是快、结果好解释，缺点是遇到群体分层或者隐性亲缘关系时会产生大量假阳性。如果你做的是同质人群、样本间没有亲缘关系、并且已经用主成分校正过分层，那 OLS 交出的结果其实是可以用的。

往上一档是线性混合模型（LMM），把样本间的遗传相似性作为随机效应塞进模型：

y = μ + β·G + γC + u + e，其中 u ~ N(0, σ²_g·K)

K 就是前面提到的 GRM。这一项的作用是把“这两个人长得像，是因为他们基因像”这层信息吸收掉，从而把分层和亲缘造成的虚假关联压下去。LMM 是目前中等样本量研究的主流选择，GCTA、GEMMA、BOLT-LMM 都属于这一类。BOLT-LMM 用了贝叶斯方法加速，几万样本也能跑得动；GEMMA 在多性状、多模型上更灵活。

最近几年更常用的是SAIGE这一类基于 saddle-point approximation 的方法，它专门解决二分类表型在病例对照极度不平衡时的 P 值失准问题。传统逻辑回归在病例很少（比如几千人里只有几十个病例）的时候，P 值会偏小，假阳性一大堆，SAIGE 把这部分偏差校正掉了，代价是计算稍微慢一点。

怎么选？我的判断顺序是这样：先看表型类型，二分类且病例少，直接上 SAIGE 或 REGENIE；连续型、样本量几万以内，LMM 类工具足够；样本量不到一千、人群又比较同质，PLINK 加主成分校正完全能应付，别为了用高级模型把流程搞得过于复杂。模型不是越复杂越好，够用就行，因为越复杂的模型，出错的地方也越多。

2.3 样本量和统计效力：为什么你的 P 值总是不显著

这是被问得最多的问题：“我跑完了，一个显著位点都没有，是不是代码写错了？”大概率不是代码的问题，是效力不够。

单位点关联的检验效力，近似取决于一个叫非中心参数的量：

NCP ≈ 2·N·p(1−p)·β² / σ²_y

N 是样本量，p 是次等位基因频率，β 是每个等位基因的效应大小，σ²_y 是表型方差（β 用标准化单位时 σ²_y = 1）。要在全基因组显著水平（5×10⁻⁸）下达到 80% 的检出率，NCP 大致需要 30 左右。把这个式子反过来解，就能估出你需要多少样本。

举个具体的：假设某个位点 MAF = 0.2，每等位基因效应 β = 0.02 个标准差（这在复杂性状里已经算中等偏上了），代入得到 2 × 0.2 × 0.8 × 0.0004 = 1.28×10⁻⁴，那么 N ≈ 30 / 1.28×10⁻⁴ ≈ 23 万。也就是说，二十多万人才有八成把握抓到这个位点。如果你手上只有两千人，那这个位点对你的研究来说基本等于不存在。

这解释了为什么顶刊上的 GWAS 动辄几十万、上百万人。也解释了为什么你复现不出来别人的结果——不是方法错了，是量级差了两三个数量级。

提示：低频位点（MAF < 0.05）的效力下降得特别快，因为 p(1−p) 这一项变小了。这也是为什么罕见变异关联分析通常要换一套方法（比如 burden test、SKAT），而不是硬套 SNP 级别的回归。

理解了效力这件事，你对结果的预期就会现实很多。一个两千人的小样本 GWAS，能挖出一两个提示性显著（P < 1×10⁻⁵）的位点，再用独立样本验证一下，就已经是不错的结果了。非要追 5×10⁻⁸，容易走上“调参数调到显著”的歪路。

3. 从原始文件到结果表：一套可复现的实操流程

理论讲完，进入动手环节。这一节我按实际执行顺序走一遍，命令都是可以直接改改就用的，参数给的是常规起点，你按自己数据的情况微调。

3.1 数据准备：格式、ID 与对齐

基因型数据最常见的两种格式：PLINK 的 bed/bim/fam三件套，和VCF。两者可以互相转：

# VCF 转 PLINK plink --vcf raw.vcf.gz --make-bed --out raw --double-id --vcf-half-call m # PLINK 转 VCF（需要 --recode） plink --bfile raw --recode vcf-iid --out raw_vcf

--double-id这个参数容易踩坑，它把样本 ID 和家系 ID 设成一样，如果你的 VCF 里只有样本名没有家系信息，加这个参数能省掉一堆报错。--vcf-half-call处理半缺失基因型（比如只有一条链上测到的位点），设成 m 表示当作缺失处理。

接下来最关键的一步：对齐表型和基因型的样本 ID。别跳过，这步出问题的概率比你想象的高得多。常见的情况是基因型文件用数字 ID，表型表用带前缀的字符串，或者两边顺序不同但你以为一样。稳妥的做法是导出两边的 ID 列表，用comm或sort求交集：

awk '{print $2}' raw.fam | sort > geno_id.txt awk 'NR>1{print $1}' pheno.txt | sort > pheno_id.txt comm -12 geno_id.txt pheno_id.txt | wc -l # 看交集有多少

交集数量应该等于你预期的样本数。如果差得多，先停下来查 ID 规则，别急着往下走。另外，表型的行顺序和 fam 文件的行顺序必须一致，PLINK 默认按 fam 文件顺序读表型，如果你直接拿 Excel 排过序的表进去，结果会整体错位，而且不会报错，这是最阴险的一类 bug。

3.2 质控：这一步偷懒，后面全是假的

质控是整个流程里最枯燥、也最不能省的一步。我见过太多人拿到数据直接跑关联，出来一堆漂亮结果，最后一查质控发现一半的位点 call rate 不到 90%，全部作废。

质控分两个层面：位点层面和样本层面。

位点层面主要看三个指标。call rate（也叫 missingness），即这个位点有多大比例的样本没测出来，常规阈值是剔除缺失率高于 2% 到 5% 的位点。MAF，次等位基因频率，低于 1% 的位点一般剔除，因为频率太低的位点统计上不稳，P 值容易失真。HWE（哈迪-温伯格平衡），检验位点基因型分布是否偏离期望，偏离太厉害说明可能存在基因分型错误，常见阈值是 P < 1×10⁻⁶ 剔除。注意 HWE 过滤在病例对照研究里需要分群体做，因为疾病本身就会让病例组偏离平衡，混在一起做会把真实信号筛掉。

样本层面看四个指标。样本 call rate，低于 95% 到 98% 的样本剔除。杂合度，个体的杂合位点比例如果偏离群体均值太远（一般超过 ±3 个标准差就值得警惕），可能是样本污染、重复或者近亲。性别核查，用 X 染色体上的杂合率和 Y 染色体上的位点密度判断性别，和登记性别不符的样本要处理。重复样本，用 IBD 或者 KING 的亲缘系数找出重复个体，同一个人测了两次只留一份。

# 第一轮：位点层面 plink --bfile raw \ --geno 0.02 \ --maf 0.01 \ --hwe 1e-6 \ --make-bed --out qc_snp # 第二轮：样本层面 plink --bfile qc_snp \ --mind 0.02 \ --het \ --make-bed --out qc_ind # 查看杂合度分布，找出偏离的个体 awk 'NR>1{print $6}' qc_ind.het | sort -n | awk '{a[NR]=$1} END{print a[int(NR*0.5)]}'

参数建议表：

注意：质控顺序会影响最终保留的位点数。一般先做位点层面的过滤，再做样本层面的过滤，最后再补一轮位点过滤。因为剔除样本后，某些位点的 call rate 可能重新变合格。这个迭代不是必须的，但做一遍心里更有底。

还有一步容易被忽略：只保留常染色体（1 到 22 号）。性染色体上的位点在男女之间剂量不同，处理方式特殊，除非你有明确的研究目的，否则第一遍分析把性染色体排除掉，省心得多。PLINK 里加--autosome就行。

3.3 群体分层与亲缘关系：假阳性的头号来源

群体分层是这么回事：假设你要研究身高，样本里混了一部分北方人（平均身高偏高）和南方人（平均身高偏低），同时这两个群体在等位基因频率上也有系统性差异。那么任何一个频率在北高南低的位点，都会和身高“显著相关”——但这个相关完全是人口结构造成的，跟身高的生物学机制毫无关系。这就是分层导致的假阳性，而且它能在 QQ 图上留下非常典型的上翘尾巴。

处理方式主要是三步。

第一步，PCA 校正。先用剪枝后的独立位点做主成分分析，取前 10 个（或者前 20 个）主成分作为协变量放进回归模型。这一步能吸收掉大部分由群体结构引起的变异。

# 先生成剪枝列表（去掉强 LD 的冗余位点） plink --bfile qc_ind \ --indep-pairwise 50 5 0.2 \ --out prune # 用剪枝后的位点做 PCA plink --bfile qc_ind \ --extract prune.prune.in \ --pca 10 \ --out pca

第二步，看 PCA 图。把前两个主成分画出来，按自报人群标签着色，看看有没有明显分簇。如果有一簇人孤零零地待在外面，那可能是另一个群体的混入，或者干脆是样本标签错了。这一步建议一定要做，图上能看出很多数据质量问题的苗头。

第三步，处理亲缘关系。如果样本里有大量家系、或者有明显的一二级亲属，普通 LMM 也压不住的时候，一般两种做法：一是剔除其中一个（保留表型更完整的），二是用 LMM 的 GRM 把亲缘信息整体吸收掉。前者简单粗暴但会损失样本量，后者更优雅但计算量大。

# 计算 GRM（GCTA 示例） gcta64 --bfile qc_ind \ --autosome \ --maf 0.01 \ --make-grm \ --out grm # 查看亲缘系数大于 0.0884 的配对（对应二级亲属及以上） gcta64 --grm grm --grm-cutoff 0.0884 --make-grm --out grm_cut

判断分层的另一个直观指标是λGC（基因组控制膨胀因子）。它等于观测到的检验统计量中位数除以理论中位数。理想情况 λGC 接近 1.0（比如 1.02 到 1.05 之间很正常，大样本下会略微偏大）。如果 λGC 跑到 1.2 以上，八成是分层没控制好，别急着往下跑，回去查 PCA。反过来，λGC 明显小于 1（比如 0.95）也不正常，可能是数据重复、P 值被过度校正，或者表型被过度整理过。

实操心得：有人喜欢用 genomic control 直接把 λ 校正到 1，这招能压住膨胀，但也可能把真实的弱信号一起压没，不建议作为常规操作。它适合用在 meta 分析阶段处理单个队列的波动，不适合用来掩盖自己流程里的分层问题。

3.4 关联分析：跑起来，并把结果存干净

质控和分层处理完，进入正式关联。这里给两套常见命令，按前面选定的模型来。

用 PLINK 做基础线性/逻辑回归：

# 连续型表型，前 10 个主成分当协变量 plink --bfile qc_ind \ --pheno pheno.txt --pheno-name height \ --covar pca.eigenvec --covar-name PC1-PC10 \ --covar pheno.txt --covar-name age,sex \ --linear \ --out gwas_height # 二分类表型 plink --bfile qc_ind \ --pheno pheno.txt --pheno-name case \ --covar pca.eigenvec --covar-name PC1-PC10 \ --logistic \ --out gwas_case

用 SAIGE 处理不平衡的病例对照：

library(SAIGE) # 第一步：拟合空模型 fitNULLGLMM( plinkFile = "qc_ind", phenoFile = "pheno.txt", phenoCol = "case", covarColList = c("age", "sex", paste0("PC", 1:10)), sampleIDColinfile = "ID", traitType = "binary", outputPrefix = "null_model" ) # 第二步：逐位点关联 SPAGMMATtest( vcfFile = "qc_ind.vcf.gz", GMMATmodelFile = "null_model.rda", varianceRatioFile = "null_model.varianceRatio.txt", SAIGEOutputFile = "assoc_case.txt", minMAF = 0.01, minMAC = 5 )

结果文件一般包含这些列：CHR、BP、SNP（rsID）、A1、A2、BETA 或 OR、SE、P、N。拿到结果先别急着画图，做两件事：一是按 MAF 和 P 值筛掉明显的异常值（比如 P = 0 或者 P = 1 的可疑记录），二是检查 β 的符号方向是不是和等位基因的定义对得上。A1 是次等位基因还是主等位基因，不同软件的定义不同，PLINK 默认 A1 是次等位基因，SAIGE 的规则也不完全一样，符号搞反了，后面所有解释都是错的。

命令里有个细节单独说一下：PLINK 的--covar-name可以写多次，但每次都要跟一个--covar指明文件来源，如果主成分和协变量在同一个文件里，写一次就够了。上面分开写是为了演示，实际用的时候建议把表型、协变量、主成分合并成一张宽表，减少出错概率。

4. 结果解读：从 P 值到生物学含义

跑出结果只是开始，接下来要判断这些结果可不可信，以及它们意味着什么。这一节讲三个层面的工作：可视化、精细定位、下游应用。

4.1 曼哈顿图与 QQ 图怎么读

曼哈顿图是最直观的展示方式，横轴是染色体上的位置，纵轴是 −log₁₀(P)。图上一簇簇的高峰就是候选信号。看这个图的时候，重点看三件事：峰的绝对高度（有没有越过 5×10⁻⁸ 这条线）、峰的形态（是单个尖峰还是宽宽的一大片）、以及峰出现在哪条染色体上。

峰的形态很有意思。一个窄而尖的峰，说明信号集中，可能就是一个强效位点；一个很宽、覆盖几十万到上百万碱基的峰，说明这一整段区域都处在一个强 LD 块里，真正起作用的位点可能只是其中一小个，这种就得靠精细定位去分辨。

library(qqman) gwas <- read.table("assoc_case.txt", header = TRUE) manhattan(gwas, chr = "CHR", bp = "BP", snp = "SNP", p = "P", col = c("grey30", "grey60"), suggestiveline = -log10(1e-5), genomewideline = -log10(5e-8), annotatePval = 5e-8)

QQ 图是诊断工具，横轴是理论 P 值分布，纵轴是实际观测值。理想情况下点应该贴在对角线上。解读的关键在两端：中间大部分贴合、只在最右上角翘起来，这是好现象，说明模型没系统性问题、且真的抓到了信号；整个图从头到尾系统性偏离对角线，那是膨胀，回去查分层；反过来整条线趴在对角线下方，说明 P 值被过度校正了。

还有一类现象叫“早期翘起”，也就是在 −log₁₀(P) 大概 2 到 3 的地方就开始偏离对角线。这通常不是真实信号，而是混杂没控干净，或者表型做了不该做的变换。看到这种情况，先检查协变量是不是漏了关键项。

提示：画图之前先把 P 值为 0 的位点替换成一个极小值（比如 1×10⁻³⁰⁰），否则取对数会得到无穷大，图直接画崩。这个坑几乎人人都会踩一次。

4.2 精细定位：从一片区域缩到一个位点

如前所述，LD 会让一个真实的因果位点周围一大片都显著。**精细定位（fine-mapping）**的目的就是把这一片区域缩到一个可信的候选集合。

主流方法是SuSiE或者FINEMAP，核心思路是贝叶斯模型：假设这个区域里有若干个因果位点（通常是 1 到 5 个），每个位点的后验概率由它的 Z 值和它与周围位点的 LD 结构共同决定。输入需要三样东西：区域内所有位点的汇总统计（Z 值、样本量、等位基因频率）、一个匹配人群的 LD 参考矩阵、以及位点列表。输出是每个位点的后验包含概率（PIP），以及若干个“可信集”，通常要求可信集里包含真实因果位点的概率达到 95%。

library(susieR) # 提取目标区域的汇总统计和 LD 矩阵 z <- region_stats$Z R <- ld_matrix_local # 从参考面板提取，必须和汇总统计的位点顺序一致 n <- 50000 fit <- susie_rss(z = z, R = R, n = n, L = 5) summary(fit)

这里最常出错的地方是LD 矩阵和汇总统计的位点顺序必须严格一致，少一个位点、多一个位点、顺序打乱一个，结果全废，而且不会报错。建议在跑之前先做一次显式检查：把两边的 rsID 向量做完全相等的比较，不相等就直接停下。

做完精细定位，通常会再叠一层功能注释，看看这些候选位点落在哪些基因上、是不是在调控区域、有没有 eQTL 证据支持。常用工具是 FUMA，它可以自动完成位置映射、染色质状态注释、eQTL 共定位等一整套流程。再往深一步可以做colocalization 分析（比如 coloc），判断“这个位点既和疾病相关、又和某个基因的表达量相关”到底是同一个因果变异在起作用，还是两个独立的信号恰好靠得近。这一步对判断后续实验方向很有价值。

4.3 下游应用：多基因风险评分与因果推断

GWAS 结果最有价值的两个下游方向，一个是多基因风险评分（PRS），一个是孟德尔随机化（MR）。

PRS 的思路很直接：把一个个体的所有风险位点的效应量加起来，得到一个总分，用来预测他/她患某个疾病或具有某个性状的风险。最经典的方法是 C+T（clumping + thresholding），先按 P 值阈值筛选位点、再按 LD 剪枝，剩下的位点按效应量加权求和。现在更常用的是贝叶斯收缩类方法，比如 LDpred2 和 PRS-CS，它们不粗暴地筛位点，而是对所有位点的效应量做收缩估计，效果通常好不少。

library(bigsnpr) obj <- snp_attach("geno.rds") G <- obj$genotypes # 用公开的 GWAS 汇总统计训练 gwas <- readRDS("sumstats.rds") beta <- snp_ldpred2_inf(G, gwas$beta, gwas$N) # 计算每个个体的得分 score <- big_prodVec(G, beta)

PRS 有一个非常关键的限制：跨人群的可迁移性很差。用某一群体训练出来的模型，直接套到另一个群体上，预测效力往往掉一半甚至更多。原因有两个，一是 LD 结构在不同人群里不一样，二是位点频率分布不同。所以如果你要在目标人群里用 PRS，要么用该人群自己的 GWAS 训练，要么做跨人群的迁移学习校正，直接照搬别人的模型是不可靠的。

MR 的思路是拿基因型当工具变量，推断某个暴露因素对结局的因果关系。它有三个核心假设：工具变量和暴露强相关、工具变量和结局之间除了通过暴露之外没有其他路径（无多效性）、工具变量和混杂因素无关。实际分析中，第三条最难验证。常用方法有 IVW、MR-Egger、加权中位数，其中 MR-Egger 的截距项可以用来检测方向性多效性，加权中位数对部分无效工具变量更稳健。

实操心得：MR 分析里最容易被审稿人挑的是工具变量的选择。P 值阈值别卡太松（比如 5×10⁻⁸），LD 剪枝一定要做（r² < 0.001），并且要报告 F 统计量（一般要求大于 10，太低的工具变量是弱工具，结果不可信）。这三条做到位，基本能挡住大部分质疑。

5. 常见问题与排查实录

最后这一节是我这些年踩过的坑的汇总。很多问题软件不会给你明确报错，结果看起来还挺正常，但就是不能用，这类问题最耗时间。

5.1 现象与原因速查表

关于最后那条，补充一句。KING 这个工具判断亲缘关系比 PLINK 的 IBD 更稳，尤其是在人群有分层的时候。它输出的是 kinship 系数，和 IBD 的对应关系大致是：kinship 0.177 对应一级亲属（亲子、全同胞），0.0884 对应二级亲属，0.0442 对应三级亲属。如果你看到一堆样本的 kinship 在 0.4 以上，那基本就是同一个人被测了两次，直接去重。

5.2 几条真正来自实战的经验

第一条：先把小流程跑通，再上全量数据。拿到一个几百万位点、几万样本的数据集，第一反应不应该是直接扔进去跑。先抽 500 个样本、5 万个位点，把整条链路跑一遍，确认每一步的输入输出都对得上，再扩到全量。这样出问题的时候，排查范围小得多。

第二条：每个中间文件都留一份。质控后的、剪枝后的、PCA 后的、关联结果，全部保留，别为了省磁盘空间覆盖掉。我吃过一次亏，跑到最后发现结果不对，想回头查是哪一步出的问题，结果中间文件全被覆盖了，只能从头再来。磁盘现在便宜得很，这点空间别省。

第三条：用已知位点做 sanity check。每个性状几乎都有几个研究得很透的位点，比如身高上的某个高频位点、某些代谢指标上的经典位点。跑完自己的数据，先去看这些位点在不在、方向对不对、效应量级合不合理。这一步能在十分钟内告诉你整套流程有没有系统性错误，比慢慢排查高效得多。

第四条：注意等位基因的“链方向”。不同芯片、不同参考面板用的可能是不同的链（正链或负链），同一个位点的 A/T 可能在这份数据里写作 T/A。做 meta 分析或者用公开汇总统计的时候，一定要先做等位基因对齐（harmonization），把 A/C 和 T/G 这类互补的情况单独标记出来处理。这个环节出错，meta 分析的结果会莫名其妙地互相抵消。

第五条：不要过度解读提示性显著的结果。P 值在 1×10⁻⁵ 到 5×10⁻⁸ 之间的位点，只能当作“值得进一步看的线索”，不能写进结论。真想追，就找独立样本做验证，验证通过了再说。这一点在写论文的时候尤其重要，审稿人几乎一定会问验证样本的事。

第六条：记录每一次参数变更。GWAS 流程里参数多，很容易出现“我记得改过但想不起来改成什么了”的情况。用一个简单的日志文件，每次跑都记下时间、输入文件、命令、参数、输出文件，出问题的时候能快速回溯。这个习惯看似琐碎，实际上能省下大量时间。

第七条：表型比方法重要。这句话我想放在最后强调。一个定义清晰、测量一致、噪声低的表型，配上最朴素的线性回归，往往比一个乱七八糟的表型配上最先进的混合模型能出更多真结果。花在表型清洗上的时间，回报率远高于花在模型调参上的时间。可惜的是，大部分人（包括早年的我）都更愿意折腾模型，因为这看起来更“技术”。