GWAS中的各种模型

GWAS 为什么总是在做回归

如果把 GWAS 想成一句话，它其实就是：

对每个变异位点，问一次 这个位点的基因型变化，是否会系统性地伴随表型变化。

最自然的数学语言就是回归。 设有 $n$ 个样本、$m$ 个位点。对第 $j$ 个位点，记

• $y \in \mathbb{R}^n$ 为表型向量
• $g_j \in \mathbb{R}^n$ 为第 $j$ 个位点的基因型向量
• $C \in \mathbb{R}^{n \times q}$ 为协变量矩阵，比如年龄、性别、批次、前几个 ancestry PCs
• $\beta_j$ 为我们真正关心的 SNP 效应

那么最经典的 GWAS 单位点模型可以统一写成

也就是

或者在二分类表型时写成

所以从统计本质上看，GWAS 更像是：

把回归模型重复跑到全基因组的每一个位点上，然后整理出哪一些位点的 $\beta_j$ 明显不为 0。

GWAS 的数据对象

在真正进入模型之前，最好先把 GWAS 里最常见的几个量想清楚。

1）基因型矩阵

通常记作 $G \in \mathbb{R}^{n \times m}$，第 $i$ 行是第 $i$ 个个体，第 $j$ 列是第 $j$ 个位点。

最常见的加性编码是

• 0 表示参考等位基因纯合
• 1 表示杂合
• 2 表示效应等位基因纯合

如果是 imputation dosage，则 $g_{ij}$ 不一定是整数，而是 $[0,2]$ 之间的期望等位基因剂量。

2）表型向量

• 定量性状，比如身高、血压、表达量，常用线性回归
• 二分类性状，比如病例/对照，常用逻辑回归
• 极端不平衡二分类，往往更适合 Firth 或 logistic mixed model
• 稀有变异常常不做单个位点，而做 burden / SKAT 一类的区域或基因水平回归

3）协变量矩阵

最常见的协变量包括

• age
• sex
• batch
• recruitment center
• genotype array
• ancestry PCs
• 其他已知混杂因素

GWAS 的一个关键思想是：

我们希望估计的是在控制这些混杂因素以后，基因型本身对表型的额外贡献（即：遗传对表型的贡献）。

GWAS 的回归

下面这张图先给一个全景图。

接下来我们一个个展开。

基因型编码方式决定在检验什么

很多人一上来就看回归公式，但真正最容易忽略的是：

同一个位点，换一种编码方式，你检验的问题就已经变了。

1）加性模型 additive

最常见，也是现代 GWAS 默认模型。

解释是每增加 1 个效应等位基因，表型平均改变 $\beta$。

这是最常见的原因有两个。

• 它参数少，检验功效通常较高。
• 很多真实效应即便不是完美加性，也常常能被加性模型近似抓到。

2）显性模型 dominant

把 1 和 2 合并。

这时 $\beta$ 问的是：只要带有至少一个效应等位基因，是否就会改变表型。

3）隐性模型 recessive

把只有 2 拿出来。

这时 $\beta$ 问的是：只有效应等位基因纯合时，表型是否改变。

4）基因型模型 genotypic 2 df

用两个 dummy 变量而不是一个加性剂量。

这时你不再强行假设 0→1→2 是线性的，而是让杂合和纯合有各自效应，最后做一个 2 df 联合检验。

5）dosage 模型

如果是 imputed 数据，$g$ 常取期望 alt allele count，比如 0.13、1.72 这种连续值。这时模型仍然是加性回归，但变量不是硬调用的 0/1/2，而是更平滑的剂量值。

线性回归

对于定量性状，最常见的单位点模型是

这里

• $\alpha$ 是截距
• $C\gamma$ 是协变量效应
• $g\beta$ 是 SNP 的固定效应
• $\varepsilon$ 是残差

我们最关心的是零假设

如果拒绝 $H_0$，就说明在控制协变量以后，这个位点与表型仍然显著相关。

矩阵写法

设设计矩阵

那么 OLS 估计量为

其中 $\hat{\theta}$ 最后一项就是 $\hat{\beta}$。

残差方差估计为

其中 $p$ 是参数个数。

于是 $\hat{\beta}$ 的标准误为

接着就有

在自由度 $n-p$ 下计算双侧 $p$ 值。

最小可读 Python 实现

import numpy as np
from scipy import stats

def ols_snp(y, g, C=None):
    """
    y: (n,) quantitative phenotype
    g: (n,) genotype dosage / 0,1,2
    C: (n,q) covariates, can be None
    """
    y = np.asarray(y, dtype=float)
    g = np.asarray(g, dtype=float)

    if C is None:
        X = np.column_stack([np.ones(len(y)), g])
    else:
        C = np.asarray(C, dtype=float)
        X = np.column_stack([np.ones(len(y)), C, g])

    XtX_inv = np.linalg.inv(X.T @ X)
    beta_hat_all = XtX_inv @ X.T @ y
    y_hat = X @ beta_hat_all
    resid = y - y_hat

    n, p = X.shape
    sigma2 = (resid @ resid) / (n - p)
    cov_beta = sigma2 * XtX_inv

    beta = beta_hat_all[-1]
    se = np.sqrt(cov_beta[-1, -1])
    t_stat = beta / se
    pval = 2 * stats.t.sf(np.abs(t_stat), df=n - p)

    return {
        "beta": beta,
        "se": se,
        "t": t_stat,
        "p": pval
    }

解释方式

如果 $g$ 用的是 0/1/2 编码，那么 $\hat{\beta}$ 表示：每增加 1 个效应等位基因，表型均值平均改变多少单位。

比如

• 身高单位是 cm，那么 $\beta=0.35$ 表示每个 effect allele 对应平均身高增加 0.35 cm
• 表达量单位是 log2(TPM+1)，那么它就是在这个刻度下的效应

逻辑回归

标准逻辑回归

如果表型是病例/对照，线性回归就不合适了。因为线性回归可能给出小于 0 或大于 1 的概率预测，而且误差分布假设也不匹配。

最常见的模型是逻辑回归：

也可以写成

这里的 $\beta$ 不再是表型均值变化，而是 log-odds 的变化。

因此

表示每增加 1 个 effect allele，患病 odds 乘上多少倍。

比如

• $\beta = 0.182$，则 $OR \approx 1.20$
• $\beta = -0.223$，则 $OR \approx 0.80$

也就是风险升高 20% 或下降到 80%。

最常见的检验

通常输出的是 Wald 统计量

再用标准正态分布计算 $p$ 值。

最小可读 Python 实现

import numpy as np
import statsmodels.api as sm
from scipy import stats

def logistic_snp(y, g, C=None):
    y = np.asarray(y, dtype=float)
    g = np.asarray(g, dtype=float)

    if C is None:
        X = np.column_stack([g])
    else:
        C = np.asarray(C, dtype=float)
        X = np.column_stack([C, g])

    X = sm.add_constant(X, has_constant="add")
    fit = sm.Logit(y, X).fit(disp=0)

    beta = fit.params[-1]
    se = fit.bse[-1]
    z = beta / se
    pval = 2 * stats.norm.sf(np.abs(z))
    odds_ratio = np.exp(beta)
    ci_low = np.exp(beta - 1.96 * se)
    ci_high = np.exp(beta + 1.96 * se)

    return {
        "beta": beta,
        "se": se,
        "z": z,
        "p": pval,
        "or": odds_ratio,
        "or_ci95": (ci_low, ci_high)
    }

Firth 逻辑回归

Firth 逻辑回归本质上是给逻辑回归的似然加了一个偏差修正型惩罚项，用来缓解小样本偏差和 separation 问题。

惩罚后的目标函数可以写成

其中

• $\ell(\theta)$ 是普通逻辑回归的 log-likelihood
• $I(\theta)$ 是 Fisher information matrix

它的直观理解是： 普通逻辑回归在样本小、MAC 很低、病例对照极不平衡时，估计会往无穷大方向跑。Firth 会给它加一点“刹车”，让估计回到有限值。

所以在 GWAS 中，Firth 特别适合这些情形：

• 低 MAC 变异
• 病例对照严重不平衡
• 某个位点出现 complete / quasi separation
• 普通 logistic 不收敛

下面这段代码将展示 Firth 的核心思想。

import numpy as np
from scipy.special import expit

def firth_logistic_core(y, X, max_iter=100, tol=1e-8):
    """
    Firth logistic 核心迭代
    y: (n,)
    X: (n,p) design matrix, should include intercept
    """
    y = np.asarray(y, dtype=float)
    X = np.asarray(X, dtype=float)
    beta = np.zeros(X.shape[1])

    for _ in range(max_iter):
        eta = X @ beta
        p = expit(eta)
        W = p * (1 - p)

        # Fisher information
        XtWX = X.T @ (W[:, None] * X)
        XtWX_inv = np.linalg.inv(XtWX)

        # leverage h_i
        H_diag = np.sum((X @ XtWX_inv) * X, axis=1) * W

        # adjusted score
        a = (0.5 - p) * H_diag
        U_star = X.T @ (y - p + a)

        step = XtWX_inv @ U_star
        beta_new = beta + step

        if np.max(np.abs(beta_new - beta)) < tol:
            beta = beta_new
            break
        beta = beta_new

    return beta

混合模型

为什么普通回归不总是够用

如果样本彼此独立、群体结构很弱、协变量控制得也不错，那么普通 OLS / logistic regression 很好用。但真实数据经常不是这样。常见问题有两个。

1）population stratification

不同 ancestry 子群体的等位基因频率不同，而表型均值也不同，这会制造虚假的 SNP 关联。

2）sample relatedness

样本之间有亲缘关系，比如兄弟姐妹、堂表亲，或者来自同一家系。此时样本不再独立，普通回归的标准误会出问题。

这时就进入了混合模型的世界。

线性混合模型 LMM

线性混合模型常写成

其中

这里

• $K$ 是 genetic relationship matrix，简称 GRM 或 kinship matrix
• $u$ 是由全基因组背景和样本相关性带来的随机效应
• $\tau$ 和 $\sigma^2$ 是方差分量

于是表型协方差变成

它的直觉是：样本之间并不是相互独立，而是“越像亲戚、越像同 ancestry”的个体，残差越相关。

这样做的好处是可以同时更稳地控制

• relatedness
• cryptic relatedness
• population structure
• 一部分 polygenic background

这也是为什么现代大样本定量性状 GWAS 里，LMM 非常常见。

一个最简的“已知协方差矩阵”教学版实现

为了把核心想法讲清楚，我们暂时假设 $V=\tau K+\sigma^2 I$ 已经给定。此时可以对白化后的数据做 GLS：

import numpy as np
from scipy import stats

def gls_snp_with_known_V(y, g, C, V):
    """
    GLS / LMM 核心思想
    假设 V 已知，先白化，再做 OLS
    """
    y = np.asarray(y, float)
    g = np.asarray(g, float)
    C = np.asarray(C, float)

    X = np.column_stack([np.ones(len(y)), C, g])

    # Cholesky whitening: V = L L^T
    L = np.linalg.cholesky(V)
    Linv = np.linalg.inv(L)

    y_t = Linv @ y
    X_t = Linv @ X

    XtX_inv = np.linalg.inv(X_t.T @ X_t)
    beta_all = XtX_inv @ X_t.T @ y_t
    resid = y_t - X_t @ beta_all
    n, p = X_t.shape
    sigma2 = (resid @ resid) / (n - p)

    beta = beta_all[-1]
    se = np.sqrt(sigma2 * XtX_inv[-1, -1])
    t_stat = beta / se
    pval = 2 * stats.t.sf(np.abs(t_stat), df=n - p)

    return {"beta": beta, "se": se, "t": t_stat, "p": pval}

上面这段不是完整的 LMM，因为完整 LMM 还要估计 $\tau$ 和 $\sigma^2$。

工业级工具

• GEMMA
• GCTA / fastGWA
• BOLT-LMM
• FaST-LMM
• EMMAX

其中 BOLT-LMM 特别适合大样本定量性状，GEMMA 非常经典，fastGWA 则在资源效率上很有代表性。

logistic mixed model 与 SAIGE

对于二分类表型，如果你既有

• case-control imbalance
• sample relatedness
• 大样本 biobank

那么普通 logistic regression 就常常不够稳，普通 LMM 又不是专门为 binary trait 设计的。

这时常见做法是 logistic mixed model：

问题在于，它比线性模型更难算，尤其是在几十万样本、上千万变异时。

SAIGE 的思想很重要，因为它解决了一个现实痛点：二分类表型特别不平衡时，普通渐近正态近似的 p 值会失真，尤其在低频变异处更明显。

因此 SAIGE 的典型思路是

• 第一步拟合 null logistic mixed model
• 第二步对每个变异做 score test
• 再用 saddlepoint approximation 修正 p 值

从统计直觉上说，SAIGE 不是把 $\beta$ 的估计“变得更漂亮”，而是把极端不平衡情况下的 p 值校准 做得更靠谱。

稀有变异

如果一个基因里有很多 rare variants，那么逐个位点单独检验经常功效很差。因为每个位点太稀有，单独拿出来信号很弱。 于是就有了 set-based regression。

1）burden test

先把一组变异压缩成一个 burden score

然后做

或者二分类版本的逻辑回归。它适合的前提是：同一组变异大致朝同一个方向起作用。

2）SKAT

SKAT 不强迫所有变异同方向、同大小，而是把变异效应当作随机效应，最后检验一个方差分量是否为 0。它更像在问：这个基因里的变异整体上有没有贡献，不管它们的方向是否一致。

3）SKAT-O

折中 burden 和 SKAT，常常是实战里很常见的选择。

最简 burden 实现

import numpy as np
from scipy import stats

def burden_test_linear(y, G_set, weights=None, C=None):
    """
    y: (n,)
    G_set: (n,p_set)
    weights: (p_set,)
    """
    y = np.asarray(y, float)
    G_set = np.asarray(G_set, float)

    if weights is None:
        weights = np.ones(G_set.shape[1])

    burden = G_set @ weights

    return ols_snp(y=y, g=burden, C=C)

最简 SKAT 思想版

严格的 SKAT p 值通常要用混合卡方分布近似或 Davies / Liu 方法。可以先抓住它的核心统计量：

其中

• $r$ 是在零模型 $y \sim C$ 下的残差
• $K_G = G W W G^\top$ 是由该基因内变异构造的核矩阵

你甚至可以先用 permutation 去近似 $p$ 值。

import numpy as np

def skat_toy_permutation(y, G_set, C, weights=None, n_perm=2000, seed=42):
    rng = np.random.default_rng(seed)
    y = np.asarray(y, float)
    G_set = np.asarray(G_set, float)
    C = np.asarray(C, float)

    if weights is None:
        weights = np.ones(G_set.shape[1])

    # null model residuals
    X0 = np.column_stack([np.ones(len(y)), C])
    beta0 = np.linalg.inv(X0.T @ X0) @ X0.T @ y
    r = y - X0 @ beta0

    GW = G_set * weights
    K = GW @ GW.T
    Q_obs = r.T @ K @ r

    count = 0
    for _ in range(n_perm):
        r_perm = rng.permutation(r)
        Q_perm = r_perm.T @ K @ r_perm
        if Q_perm >= Q_obs:
            count += 1

    pval = (count + 1) / (n_perm + 1)
    return {"Q": Q_obs, "p_perm": pval}

数据的标准化

标准化之前

GWAS 里最容易混淆的不是“该不该标准化”，而是这三种操作经常被混在一起：

1）中心化 centering

只减去均值。

2）标准化 standardization / z-score

既减均值，又除以标准差。

3）残差化 residualization

先把 $x$ 对协变量 $C$ 回归，再取残差

这三者的数学效果并不一样。

表型标准化

设原始模型是

现在把表型做 z-score：

则新模型中的 SNP 效应会变成

直觉上就是：原来 $\beta$ 的单位是“表型原始单位每 allele”，标准化以后就变成“表型标准差单位每 allele”。

更重要的是：

• $\beta$ 会缩放
• $SE(\beta)$ 也会按同样比例缩放
• $t=\beta/SE$ 不变
• $p$ 值不变
• 截距会变化

所以对线性回归来说，只标准化 $y$ 通常不会改变显著性，只改变效应量单位。 这也是为什么很多 quantitative trait GWAS 会把表型标准化后再回归，这样不同表型之间效应大小更可比较。

基因型标准化

如果把基因型也做 z-score：

其中对加性编码 $g \in \{0,1,2\}$ 来说，在 Hardy-Weinberg 平衡近似下

这里的 $p$ 是 effect allele frequency。

这时模型变成

有

也就是说，原来是 每增加 1 个 allele 的效应，现在变成 每增加 1 个 genotype SD 的效应。

因此：

• $\beta$ 的解释发生了变化
• $SE$ 也按同样比例变化
• Wald / t / z 统计量通常不变
• $p$ 值通常不变
• 但 OR 或 BETA 的生物学解释已经不是“per allele”了

所以在 GWAS summary statistics 里，最常见仍然是使用 0/1/2 或 dosage 的原始加性编码，而不是把每个 SNP 单独 z-score 以后再输出结果。因为后者不利于跨研究比较 per-allele effect。

什么时候会标准化基因型

• 构建 GRM 时，通常会按等位基因频率中心化/缩放
• ridge / BLUP / whole-genome prediction 里常常会内部标准化
• PCA 里也常做中心化与按位点方差缩放

所以你会发现：

“GWAS 最终单位点回归结果常不用标准化 genotype 输出，但很多中间矩阵计算会对 genotype 做标准化。”

表型和基因型都标准化

若 $y$ 和 $g$ 都 z-score，那么线性模型中的 $\beta$ 可以看成一种标准化效应大小。

在没有协变量时，它和相关系数的关系尤其直接。

在有协变量时，它更接近 partial standardized effect，也就是在控制协变量后，基因型每增加 1 个标准差，表型改变多少个标准差。

这对于比较不同表型、不同位点的相对效应大小很有用，但它会牺牲 per-allele 的直观解释。

对逻辑回归来说，标准化影响更要小心解释

对于逻辑回归

如果你只把 $g$ 缩放成 $g^*=g/s_g$，那么会有

于是

这已经不再是“每增加 1 个 allele 的 OR”，而是“每增加 1 个 genotype SD 的 OR”。所以二分类 GWAS summary statistics 中，最常见做法仍然是输出 per allele OR，也就是保留 0/1/2 或 dosage 的原始尺度。

协变量该怎么处理

这一块非常关键，因为很多“显著位点”其实死在协变量没处理好上。 最常见的完整模型是

1）不要轻易漏掉 ancestry PCs

因为 population stratification 是 GWAS 假阳性的大户。

2）对于定量性状，残差化是可以讲清楚的

在线性回归中，如果你把 $y$ 和 $g$ 都先对 $C$ 残差化，再用残差做一元回归，得到的 SNP 系数和直接在全模型里拟合是一致的。这是 Frisch–Waugh–Lovell 定理。

也就是说：

• 先回归掉协变量，再做 SNP 回归
• 直接把协变量和 SNP 一起放进模型

在线性模型里本质等价。

3）但二分类表型不要把 logistic 问题偷换成“先残差化 y 再做 OLS”

这在教学里经常见到，但并不是标准的病例对照 GWAS 处理方式。二分类表型还是应该在 logistic / mixed logistic 的框架里分析。

4）协变量也可以标准化，但含义要清楚

• age 标准化后，系数变成每 1 个 age SD 的效应
• sex 通常不要 z-score，二值变量保留 0/1 更直观
• PCs 可以原样放，通常不需要额外解释其单位

GWAS summary 数据的各列

GWAS 最终输出表通常会有这些列：

• SNP / ID
• CHR
• POS
• A1 / effect allele
• A2 / other allele
• BETA 或 OR
• SE
• STAT
• P
• N
• A1FREQ 或 EAF
• INFO
• TEST

下面把它们逐个对应到回归里。

1）BETA

线性回归里就是 $\hat{\beta}$。 逻辑回归里通常也是 $\hat{\beta}$，只是很多软件会额外给出

2）SE

标准误，反映 $\hat{\beta}$ 的不确定性。

3）STAT

• 线性模型通常是 $t$ 或近似 $z$
• 逻辑模型通常是 $z$
• Score test 则可能输出 score statistic
• LRT 则可能输出 deviance 或 $\chi^2$

4）P

根据统计量在零假设分布下计算出来。

5）95% CI

对线性回归

对逻辑回归更常报告 OR 区间

6）N 为什么有时每个 SNP 不一样

因为不同 SNP 可能缺失模式不同，某些个体在这个位点没有可靠基因型，所以参与该位点回归的样本数会变化。

OLS 的 Python 示例

下面给一个简单的单 SNP scan，对定量性状使用 OLS。

import numpy as np
import pandas as pd

def linear_gwas_scan(y, G, C=None, snp_names=None):
    """
    y: (n,)
    G: (n,m)
    C: (n,q)
    """
    y = np.asarray(y, float)
    G = np.asarray(G, float)

    m = G.shape[1]
    results = []

    if snp_names is None:
        snp_names = [f"SNP_{j}" for j in range(m)]

    for j in range(m):
        g = G[:, j]

        # 简单 missing 过滤
        mask = np.isfinite(y) & np.isfinite(g)
        if C is not None:
            mask &= np.all(np.isfinite(C), axis=1)

        y_j = y[mask]
        g_j = g[mask]
        C_j = None if C is None else C[mask]

        # 跳过方差为 0 的位点
        if np.var(g_j) == 0:
            results.append({
                "SNP": snp_names[j],
                "BETA": np.nan,
                "SE": np.nan,
                "STAT": np.nan,
                "P": np.nan,
                "N": len(y_j)
            })
            continue

        out = ols_snp(y_j, g_j, C_j)
        results.append({
            "SNP": snp_names[j],
            "BETA": out["beta"],
            "SE": out["se"],
            "STAT": out["t"],
            "P": out["p"],
            "N": len(y_j)
        })

    return pd.DataFrame(results)

二分类表型也可以类似循环，只不过把 ols_snp() 换成 logistic_snp()。

GWAS 最小 workflow 示例

从模拟基因型、计算 PCs、做 GWAS，到整理结果。

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.decomposition import PCA

rng = np.random.default_rng(42)

# ---------------------------
# 1. 模拟 genotype matrix
# ---------------------------
n = 500
m = 200

maf = rng.uniform(0.05, 0.5, size=m)
G = np.column_stack([rng.binomial(2, p, size=n) for p in maf]).astype(float)

# ---------------------------
# 2. 模拟 covariates
# ---------------------------
age = rng.normal(50, 10, size=n)
sex = rng.binomial(1, 0.5, size=n)

# 对 genotype 做中心化后提 PCA
G_std = (G - G.mean(axis=0)) / G.std(axis=0, ddof=1)
pcs = PCA(n_components=5).fit_transform(G_std)

C = np.column_stack([age, sex, pcs])

# ---------------------------
# 3. 模拟 quantitative phenotype
#    让第 17 个 SNP 真的有因果效应
# ---------------------------
causal_j = 17
beta_g = 0.8

y = (
    0.03 * age
    - 0.5 * sex
    + beta_g * G[:, causal_j]
    + rng.normal(0, 1.0, size=n)
)

# ---------------------------
# 4. 跑教学版 GWAS
# ---------------------------
res = linear_gwas_scan(y=y, G=G, C=C)

# ---------------------------
# 5. 加一点注释列
# ---------------------------
res["CHR"] = 1
res["POS"] = np.arange(1, m + 1) * 10000
res["A1FREQ"] = G.mean(axis=0) / 2.0

# 看最显著结果
print(res.sort_values("P").head())

这段代码做了几件很接近真实分析逻辑的事：

• 用 0/1/2 模拟 genotype
• 用 genotype 计算 PCs
• 把 age、sex、PCs 都放进协变量
• 对每个 SNP 做线性回归
• 输出 BETA、SE、P、N

当然，真实研究中你不会用这段代码直接跑 UK Biobank，但它非常适合理解“GWAS 从回归角度到底在干什么”。

大型 GWAS 工具

• PLINK 2 跑单 SNP 线性/逻辑/Firth
• GEMMA / BOLT-LMM / fastGWA 跑 LMM
• SAIGE 跑大样本不平衡二分类
• REGENIE 跑大型 biobank two-step regression
• Python / R 负责整理结果、过滤、注释、作图、下游分析

PLINK2示例

plink2 \
  --pfile mydata \
  --pheno pheno.tsv \
  --covar covar.tsv \
  --glm hide-covar cols=+a1freq,+nobs \
  --out gwas_demo

然后在 Python 里读结果：

import pandas as pd

glm = pd.read_csv("gwas_demo.trait.glm.linear", sep="\t")
glm = glm[glm["TEST"] == "ADD"]  # 只保留加性模型
glm = glm.sort_values("P")
print(glm.head())

回归由专用工具做，解释和展示由 Python/R 做。

QC 和预处理会直接影响回归

很多人把 QC 当成“回归前的行政流程”，其实不是。QC 会直接决定你的 $\beta$、$SE$、$P$ 有没有意义。

1）MAF / MAC 过滤

MAC 太低时，单位点 logistic 非常容易不稳，线性模型的标准误也会很大。

2）HWE 过滤

明显偏离 HWE 的位点可能提示分型错误或批次问题。

3）missingness

样本或位点缺失太高，会让每个位点参与回归的 $N$ 波动很大，结果也不稳。

4）imputation INFO

如果剂量质量很差，拿去回归的其实是高噪声 genotype。

5）表型变换

对高度偏态定量表型，常见做法包括

• log transform
• rank-based inverse normal transform

它们会影响 $\beta$ 的单位解释，但往往能让残差分布更接近模型假设，从而让统计推断更稳。

线性回归、逻辑回归、LMM、REGENIE，到底怎么选

可以粗略按下面这个思路。

情形 1：定量性状，小到中等样本，结构不复杂

• 线性回归就够
• 加上 age、sex、PCs

情形 2：二分类性状，样本比较平衡，常见变异为主

• 逻辑回归通常够用

情形 3：二分类，低频/稀有变异多，或者存在 separation

• Firth 更稳

情形 4：样本有 relatedness，或者 population structure 较强

• 定量性状优先考虑 LMM
• 二分类优先考虑 logistic mixed model / SAIGE 一类

情形 5：几十万样本的大型 biobank

• 定量性状常见 BOLT-LMM、fastGWA、REGENIE
• 二分类常见 SAIGE、REGENIE

情形 6：稀有变异

• 不要只盯着单 SNP
• 基因或区域层面的 burden / SKAT / SKAT-O 更常见

做 PCA 时为什么要标准化 genotype

做 PCA 时，常见做法是先对每个位点中心化，有时还会按位点方差缩放。原因是：

如果不做适当标准化，MAF 高的位点会因为方差更大而主导主成分。

加性编码下，位点方差近似为

因此 MAF 不同的位点天然方差不同。

而 PCA 的目标是刻画群体结构，不是让高方差位点压制一切，所以在构建 PCA 或 GRM 时，genotype 的标准化通常是合理且常见的。

这和“单位点 GWAS 输出是否要保留 per-allele 解释”是两回事，千万不要混为一谈。

结果解释时常犯的错误

错误 1：把 BETA 当成 OR

线性回归里的 BETA 不是 OR。

逻辑回归里的 BETA 也不是 OR，本身是 log-odds ratio，要指数化后才是 OR。

错误 2：没看 effect allele

$\beta>0$ 只是说对 effect allele 而言是正向，不代表对另一条等位基因也是正向。

错误 3：把标准化后的 $\beta$ 当成 per-allele effect

一旦对 genotype 做了 z-score，它就不再是每 allele 解释。

错误 4：以为加了 PCs 就彻底解决了 relatedness

PCs 更像是在固定效应层面校正 ancestry 结构，和 LMM 用 GRM 建模样本相关性不是一回事。

错误 5：binary trait 也用 residualized y 直接做 OLS

这在严格意义上不是标准病例对照 GWAS。

一张速查表

内容小结

如果把 GWAS 回归模型这一整套内容压缩成几句话，我会这样总结。

1) GWAS 的核心是对每个位点，在控制协变量以后，检验 $\beta$ 是否显著偏离 0。

2) 不同模型回答的问题不同

• 线性回归回答的是定量表型均值变化
• 逻辑回归回答的是 log-odds / OR 的变化
• Firth 解决的是小样本和 separation
• LMM 解决的是结构和 relatedness
• SAIGE 解决的是不平衡二分类的大样本问题
• REGENIE 解决的是 biobank 时代的计算可扩展性
• burden / SKAT 解决的是稀有变异单 SNP 功效不足的问题

3) 标准化不会凭空制造显著性，但会直接改变 $\beta$ 的单位解释。

• 标准化 $y$，效应变成表型 SD 单位
• 标准化 $g$，效应变成 genotype SD 单位
• 两者都标准化，变成 standardized effect
• $p$ 值在同一个模型下通常基本不变，但解释会变

参考阅读

下面这些名字你在 GWAS 里会高频遇到，建议形成一个脑图。

• PLINK 2：单 SNP 线性、逻辑、Firth 回归
• GEMMA：经典 LMM / Bayesian mixed model
• GCTA / fastGWA：GRM、GREML、MLM GWAS
• BOLT-LMM：大型定量性状 GWAS
• SAIGE：大样本不平衡二分类 GWAS
• REGENIE：two-step whole-genome regression
• SKAT：rare variant set-based association