基因组预测模型：洞悉遗传潜力的统计利器

在基因组选择育种中，核心任务是利用全基因组分子标记数据来准确预测个体的遗传潜力，即基因组估计育种值（GEBV）。为了实现这一目标，科学家们开发了多种统计模型。本节将深入探讨其中一些关键模型，包括岭回归模型、贝叶斯模型，并详细阐述 GBLUP 及其与 AlphaBeta 方法的内在联系。

---

2.2.1 岭回归基因组最佳线性无偏预测（RR-BLUP）

概念解析：
岭回归基因组最佳线性无偏预测 (RR-BLUP) 是一种直接估计每个 SNP 标记（Single Nucleotide Polymorphism，单核苷酸多态性）效应的基因组预测方法。它基于这样的假设：所有 SNP 都对性状有贡献，并且它们的效应服从同一个正态分布。RR-BLUP 通过引入岭回归（Ridge Regression）的惩罚项来处理基因组数据中 SNP 数量远大于个体数量时可能出现的共线性问题，并对 SNP 效应进行收缩（shrinkage），从而防止模型过拟合，提高预测的稳健性。

通俗理解：
想象一下，一个性状的形成是由无数个基因（SNP）共同决定的，每个基因都有自己的“贡献值”。RR-BLUP 就像一个“侦探”，它试图找出每个基因的具体贡献值。但由于基因太多，很多基因的贡献可能微乎其微，甚至有些基因之间“串通一气”（共线性），会干扰侦探的判断。RR-BLUP 的“岭回归”就像给侦探戴上了一副“约束眼镜”，强制它不要给任何一个基因过高的“贡献值”，让那些看似作用很大的基因效应也适度“收缩”回平均水平，从而确保整体判断更准确、更稳定。

数学模型：
RR-BLUP 的核心模型可以表示为：
$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{Ms}+\mathbf{e}$
其中：

• $\mathbf{y}$ 是表型观测值向量。
• $\mathbf{X}$ 是固定效应的设计矩阵，$\mathbf{\beta }$ 是固定效应的参数向量。
• $\mathbf{M}$ 是基因型矩阵，其行代表个体，列代表 SNP 标记（通常编码为 0, 1, 2 分别代表纯合子、杂合子、另一纯合子）。
• $\mathbf{s}$ 是 SNP 效应向量，假设 $\mathbf{s}\sim N(0,\mathbf{I}\sigma \_s^2)$，即所有 SNP 效应独立同分布，且服从均值为 $0$、方差为 $\sigma \_s^2$ 的正态分布。
• $\mathbf{e}$ 是残差向量，$\mathbf{e}\sim N(0,\mathbf{I}\sigma \_e^2)$。

为了估计 $\mathbf{s}$，RR-BLUP 通常通过求解以下岭回归混合模型方程获得：
$\hat{\mathbf{s}}=({\mathbf{M}}^{\prime }\mathbf{M}+\lambda \mathbf{I}{)}^{-1}{\mathbf{M}}^{\prime }(\mathbf{y}-\mathbf{X}\mathbf{\beta })$
其中，$\lambda =\frac{\sigma \_e^2}{\sigma \_s^2}$ 是岭回归惩罚参数，它控制着 SNP 效应收缩的程度。

计算 GEBV：
一旦估算出每个 SNP 的效应 $\hat{\mathbf{s}}$，任何个体的 GEBV 就可以通过将其基因型乘以对应的 SNP 效应并求和得到：
${\text{GEBV}}_i={\sum }_{j=1}^m{M}_{ij}{\hat{s}}_j$
其中，$M\_ij$ 是个体 $i$ 在第 $j$ 个 SNP 上的基因型编码，$\hat{s}\_j$ 是第 $j$ 个 SNP 的估计效应值，$m$ 是 SNP 标记的总数。

R 代码示例： (见上文，rrBLUP 包的 mixed.solve 函数通过 Z = geno_data 参数实现 RR-BLUP)

---

2.2.2 贝叶斯模型（Bayesian Models）

概念解析：
贝叶斯模型 是一类更灵活的基因组预测方法，它们通过贝叶斯统计框架来估计 SNP 效应。与 RR-BLUP 假设所有 SNP 效应服从单一正态分布不同，贝叶斯模型允许 SNP 效应具有不同的先验分布。这种灵活性使得贝叶斯模型能够更好地处理性状受少数大效应基因和大量微效基因共同控制的情况（即“稀疏效应”）。

通俗理解：
如果说 RR-BLUP 认为所有基因都“或多或少”有点用，只是贡献大小有区别，那么贝叶斯模型就更“聪明”一些。它认为，有些基因可能真的对性状有很大影响，而绝大多数基因可能根本没影响，或者影响可以忽略不计。贝叶斯模型会根据数据，自动调整对不同基因效应的“信任度”，甚至能把那些没用的基因效应直接“归零”，从而更精准地找出那些真正的“关键基因”。

常用贝叶斯模型：

• BayesA： 假设 SNP 效应服从独立的 $t$ 分布（或正态分布，但方差各异），允许少数 SNP 具有较大效应。
• BayesB： 引入了一个混合分布先验，假设一部分 SNP 效应为零（不影响性状），而另一部分 SNP 效应服从 $t$ 分布，这更符合大效应 QTL（数量性状基因座）存在的生物学现实。
• BayesC / BayesC$\pi$： 类似于 BayesB，但假定非零 SNP 效应服从一个共同的正态分布（BayesC）或允许非零效应的比例 $\pi$ 被估计（BayesC$\pi$）。

优势：

• 处理稀疏效应： 能够识别少数具有大效应的基因，适用于复杂性状的遗传结构。
• 灵活性高： 允许更复杂的先验分布假设，更好地拟合生物学现实。
• 预测精度： 在某些情况下，特别是性状受大效应 QTL 影响时，贝叶斯模型的预测精度可能高于 RR-BLUP。

R 代码示例：使用 BGLR 包实现贝叶斯模型

BGLR 包提供了多种贝叶斯基因组预测模型。

# 首先，安装并加载 BGLR 包
# install.packages("BGLR")
library(BGLR)

# 沿用之前的 geno_data 和 pheno_data$trait1

# 定义模型参数
ETA <- list(
  list(M = geno_data, model = "BL") # "BL" 对应 BayesB，其他如 "BRR" 对应RR-BLUP
)

# 运行贝叶斯模型 (这里以BayesB为例，BGLR中"BL"默认接近BayesB)
# nIter: MCMC迭代次数, burnIn: 燃烧期
# verbose = FALSE: 不显示详细的迭代过程
bayesian_results <- BGLR(
  y = pheno_data$trait1,
  ETA = ETA,
  nIter = 1000,
  burnIn = 500,
  verbose = FALSE
)

# 提取 GEBV
# BGLR 会返回一个包含预测结果的列表，bayesian_results$yHat 是预测的表型值
# 对于简单的模型，这通常接近GEBV
cat("\nBayesian Model (BL) 预测的 GEBV（部分）：\n")
print(head(bayesian_results$yHat))

# 实际的GEBV需要从模型的SNP效应累加获得，或从特定输出变量提取。
# 例如，如果 model="BL" (BayesB), SNP效应在ETA[[1]]$s中，可以这样计算GEBV:

# s\_hat\_bayesian <- bayesian\_results$ETA\[\[1]]$s

# GEBVs\_bayesian <- geno\_data %\*% s\_hat\_bayesian

# print(head(GEBVs\_bayesian))

---

2.2.3 基因组最佳线性无偏预测（GBLUP）

概念解析：
GBLUP (Genomic Best Linear Unbiased Prediction) 是一种基于**基因组亲缘关系矩阵（Genomic Relationship Matrix, G）**的 BLUP 模型。它不直接估计每个 SNP 标记的效应，而是通过构建个体间的基因组相似度矩阵（$G$ 矩阵）来捕捉所有 SNP 的累加效应，从而预测个体的基因组育种值（GEBV）。

通俗理解：
GBLUP 就像一个“族谱学家”。它不关注每个基因的具体功能，而是通过分析你的全部基因组，来判断你和哪些亲缘关系近的“族人”在基因组层面最相似。如果你和家族里那些高产的“族人”基因组非常相似，那么你的遗传潜力也很可能高。GBLUP 利用的是整体的基因组相似性（而非单个基因的效应），来推断个体的遗传潜力。

数学模型：
GBLUP 的核心模型是标准的混合线性模型，可以表示为：
$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{Zu}+\mathbf{e}$$
其中：

• $\mathbf{y}$ 是表型观测值向量。
• $\mathbf{X}$ 是固定效应的设计矩阵，$\mathbf{\beta }$ 是固定效应的参数向量。
• $\mathbf{Z}$ 是连接个体表型和其育种值的关联矩阵（通常是对角矩阵，如果每个个体只有一个表型记录）。
• $\mathbf{u}$ 是基因组育种值（GEBV）的向量，且其协方差结构为 $Var(\mathbf{u})=\mathbf{G}{\sigma }_u^2$。这里的 $\mathbf{G}$ 是基因组亲缘关系矩阵，${\sigma }_u^2$ 是加性遗传方差。
• $\mathbf{e}$ 是残差向量，其协方差结构为 $Var(\mathbf{e})=\mathbf{I}{\sigma }_e^2$。

基因组亲缘关系矩阵 (G) 的构建：
$\mathbf{G}$ 矩阵通常由个体间的 SNP 基因型数据构建，最常用的是 VanRaden 方法：
$$\mathbf{G}=\frac{\mathbf{M}{\mathbf{M}}^{\prime }}{2{\sum }_{j=1}^m{p}_j(1-p_j)}$$
其中，$\mathbf{M}$ 是中心化后的基因型矩阵（每个 SNP 位点的基因型值减去 $2p_j$，其中 $p_j$ 为该位点的等位基因频率），${\mathbf{M}}^{\prime }$ 是 $\mathbf{M}$ 的转置，$m$ 是 SNP 标记的总数。

优势：

• 计算效率高： 对于大规模群体，GBLUP 的计算相对高效，因为它处理的是个体间的关系矩阵而非每个 SNP 的效应。
• 易于理解和实现： 模型形式与传统的系谱 BLUP 相似，易于理解和推广。
• 整合所有 SNP 信息： 通过 $G$ 矩阵，隐含地整合了所有 SNP 的累加效应。

R 代码示例： (见上文，rrBLUP 包的 mixed.solve 函数通过 K = G_matrix 参数实现 GBLUP)

GBLUP 与 AlphaBeta 方法的内在联系

GBLUP 模型实际上可以从一个基于 SNP 效应的模型（如 RR-BLUP）推导出来。这种联系通常被称为 AlphaBeta 方法，它揭示了从 SNP 效应到基因组育种值的转换机制。

Alpha 方法：
指的是直接估算每个 SNP 的效应值（$\hat{\mathbf{s}}$），例如 RR-BLUP 和贝叶斯模型。
$$\hat{\mathbf{s}}=\text{function}(\mathbf{y},\mathbf{X},\mathbf{M},\text{parameters})$$
然后通过 ${\text{GEBV}}_i={\sum }_{j=1}^m{M}_{ij}{\hat{s}}_j$ 来计算 GEBV。

Beta 方法：
指的是直接估算个体的基因组育种值（GEBV，即 GBLUP 中的 $\hat{\mathbf{u}}$），而不需要显式地估计每个 SNP 的效应。
$$\hat{\mathbf{u}}=\text{function}(\mathbf{y},\mathbf{X},\mathbf{G},\text{parameters})$$
例如，在 GBLUP 中，$\hat{\mathbf{u}}$ 是通过求解 BLUP 方程得到的。

等价性：
在假设所有 SNP 效应服从相同方差的正态分布时（RR-BLUP 的核心假设），GBLUP 和 RR-BLUP 在数学上是等价的。这意味着：

• 从 Alpha 到 Beta： 通过 RR-BLUP 估算出 SNP 效应 $\hat{\mathbf{s}}$，然后用 $\mathbf{M}\hat{\mathbf{s}}$ 计算出的 GEBV，就等同于 GBLUP 直接预测的 GEBV。
• 从 Beta 到 Alpha： 理论上，如果从 GBLUP 获得了 GEBV，也可以反推出隐含的 SNP 效应，但这不是标准操作，且可能需要额外的假设。

这种等价性表明，无论我们是直接关注每个基因的贡献（Alpha 方法），还是关注个体基因组整体的相似性（Beta 方法），只要基础统计假设一致，我们都能殊途同归地获得相同的最佳基因组育种值预测。这为基因组预测提供了强大的理论支撑和灵活的实现途径。