基因组选择育种基础

简介

基因组选择是利用覆盖全基因组的成千上万个分子标记（主要是单核苷酸多态性，SNPs）来估计个体对目标性状的总遗传贡献，并据此选择优良的下一代育种材料或种畜。

基因组选择育种 (GS) 与传统育种、MAS 的核心区别

连锁不平衡性

连锁指的是位于同一染色体上的两个或多个基因在世代遗传过程中倾向于一起遗传的现象。 两个位点之间的距离越近，它们在减数分裂时通过重组而被分开的概率就越低，它们结合在一起遗传的趋势就越强。而连锁不平衡指的是两个或多个位点上的等位基因非随机的关联远高于预期随机关联的现象。

假设 A、B 两个基因座，各有两个等位基因 a、A 和 b、B：

• 如果是连锁平衡，那么出现单倍型 AB 的预期概率就是 $P_{AB,exp}=P_{A} \times P_{B}$
• 然而实际情况是单倍型 AB 出现的概率 $P_{AB}$与预期的概率偏差很大，连锁不平衡指数 D 为实际概率与期望概率之差：

$$
D = P(A \cap B) - P(A)P(B)
$$

但是 D 值会受等位基因频率影响，这使得我们无法比较不同频率的等位基因对之间连锁不平衡的大小，因此引入标准化连锁不平衡指标 $D'$：

$$
D'=D/\lvert D_{max} \rvert
$$

• $D'$的值域为$(-1,1)$
• $D'=0$ 表示完全连锁平衡，$D'=1$ 或 $D'=-1$ 表示完全连锁不平衡或者某单倍型缺失

在基因组选择中，使用平方相关系数 $R^2$来衡量连锁不平衡指数：

$$
R^2 = \frac{D^2}{P_A P_a P_B P_b}
$$

$R^2$衡量的是 SNP 标记对 QTL 位点的遗传变异信息的捕捉程度或替代程度。

• $R^2=0$ 表示完全连锁平衡，两个位点完全独立
• $R^2=1$ 表示两个位点具有完全相同的遗传信息，可以完美替代

常见的训练方法

• GBLUP(Genomic Best Linear Unbiased Prediction)
• BAYESIAN ALPHABET
• RR-BLUP(Ridge Regression BLUP)
• Machine Learning
• Deep Learning

基因组估计育种值 (GEBV)

GEBV (Genomic Estimated Breeding Value) 是指基于个体全基因组的分子标记信息所预测出的该个体的加性遗传价值。GEBV 是对个体所有微效 QTL 累积效应的综合评估。它比传统的估计育种值（EBV）更准确，因为它能够更全面地捕捉到个体基因组中对性状有贡献的每一个微效基因（QTL）的效应。

• 育种值(BV)：代表一个体将其优良基因传递给其后代的能力。BV 是后代性能的直接预测指标。
• 估计(E)：因为这是通过统计模型预测出来的数值，而不是真实观测值。
• 基因组(G)：强调这种估计是基于全基因组高密度分子标记（SNP）的信息，而非传统的系谱和表型数据。

GEBV 的计算原理是：将全基因组所有分子标记的效应值累加起来，从而得到个体的总遗传价值。

• $\text{GEBV}_j$： 个体 $j$ 的基因组估计育种值。
• $M$：全基因组中使用的 SNP 标记总数。
• $Z_{ij}$：个体 $j$ 在第 $i$ 个 SNP 位点上的基因型编码（通常为 0, 1 或 2，代表某个等位基因的数量）。
• $\hat{\alpha}_i$：第 $i$ 个 SNP 标记的效应值。这是模型训练的关键产物，它代表了该 SNP 及其连锁的 QTL 对性状的贡献大小和方向。

杨树测试

要利用杨树的 SNP 和株高数据实施基因组选择，核心过程是：首先，将已测量株高（表型）和已分型 SNP（基因型）的杨树作为参考群体，对数据进行预处理和标准化。然后，选择并训练一个线性混合模型，通过计算基因组关系矩阵或标记效应，建立基因型与株高之间的关联。模型的调参主要通过 k 折交叉验证来评估和优化预测准确性。模型训练完成后，保存关键参数，即可用于对只进行 SNP 分型、未测量株高的幼龄杨树计算基因组估计育种值（GEBV），并根据 GEBV 直接选择优良个体，从而极大地缩短育种周期。

参考群体数据处理

• 目标变量：株高，通过方差分析去除固定效应，得到校正表型值 Y
• 输入变量：SNP，筛选并移除低质量以及偏离哈迪-温伯格平衡的 SNP 标记和个体，将 SNP 转换为数值：AA 标记为 2，Aa 标记为 1，aa 标记为 0

基因组关系矩阵（G 矩阵）

传统育种中，依靠谱系关系来构建亲缘关系矩阵，即 A 矩阵，从理论上可以计算其基因组的相似性，但是实际上由于减数分裂过程中存在的随机组合等过程会导致遗传信息传递时会与理论有较大差异。而在基因组选择育种之中使用全基因组的 SNP 来计算个体之间实际的遗传相似度，由于它捕捉到了由于减数分裂过程中染色体随机分离和重组所带来的变异，因此准确率比理论值更高。

假设有 10 颗杨树（Tree_001 ... Tree_100），测了 10 个 SNP 和每株杨树的株高，利用 SNP 数据来构建 10x10 的 G 矩阵。

# 模拟株高表型
heights = np.array([15.2, 16.8, 14.5, 17.1, 15.9, 16.3, 14.8, 17.5, 15.6, 16.0])

pheno_df = pd.DataFrame({
    'Tree_ID': tree_names,
    'Height': heights
})

print("\n杨树株高表型:")
print(pheno_df)

首先需要将基因型数据数值化，例如对于一对等位基因 A/G，G 为备选育种基因，那么基因型 AA、AG、GG 分别为 0、1、2，那么前三颗杨树在前五个 SNP 位点上的基因型应该如下：

由此，我们得到一个 10 行（个体） x 10 列（SNP） 的数值矩阵 M。将矩阵输入到 Python 中：

import numpy as np
import pandas as pd
np.random.seed(42)

# 10株杨树，10个SNP
n_trees = 10
n_snps = 10

# 基因型矩阵 M
# 0=AA, 1=AG, 2=GG
M = np.array([
    [0, 1, 2, 1, 0, 2, 1, 0, 1, 2],  # Tree_001
    [1, 2, 1, 0, 1, 1, 2, 1, 0, 1],  # Tree_002
    [0, 1, 0, 2, 2, 0, 1, 2, 1, 0],  # Tree_003
    [2, 0, 1, 1, 0, 2, 0, 1, 2, 1],  # Tree_004
    [1, 1, 2, 0, 1, 1, 2, 0, 1, 2],  # Tree_005
    [0, 2, 1, 1, 2, 0, 1, 1, 0, 1],  # Tree_006
    [1, 0, 0, 2, 1, 2, 0, 2, 1, 0],  # Tree_007
    [2, 1, 2, 1, 0, 1, 1, 0, 2, 1],  # Tree_008
    [0, 1, 1, 0, 2, 2, 1, 1, 0, 2],  # Tree_009
    [1, 2, 0, 1, 1, 0, 2, 1, 1, 0]   # Tree_010
], dtype=float)

tree_names = [f"Tree_{i+1:03d}" for i in range(n_trees)]
snp_names = [f"SNP_{i+1:03d}" for i in range(n_snps)]

df_M = pd.DataFrame(M, index=tree_names, columns=snp_names)

然后要计算每个 SNP 位点中等位基因 G 的频率：

$$
p_j = \frac{\sum_{i=1}^{n} M_{ij}}{2n}
$$

• $p_j$ 是第 $j$ 个 SNP 的等位基因频率。
• $M_{ij}$ 是个体 $i$ 在 SNP $j$ 上的基因型编码值。
• $n$ 是个体总数（10）。
• 分母是 $2n$ 因为每个个体在该位点贡献两个等位基因。

使用 Python 计算：

def calculate_allele_frequencies(M):
    n = M.shape[0]
    p = np.sum(M, axis=0) / (2 * n)
    return p

p = calculate_allele_frequencies(M)
for i, freq in enumerate(p):
    print(f"  {snp_names[i]}: p = {freq:.3f}")

SNP_001: p = 0.400
SNP_002: p = 0.550
SNP_003: p = 0.500
SNP_004: p = 0.450
SNP_005: p = 0.500
SNP_006: p = 0.550
SNP_007: p = 0.550
SNP_008: p = 0.450
SNP_009: p = 0.450
SNP_010: p = 0.500

接下来就是构建 G 矩阵的核心步骤，构建中心化矩阵 $\mathbf{Z}$，其中每个元素是原始基因型值减去该 SNP 的期望值（$2p_j$）。计算公式如下：

$$
Z_{ij} = M_{ij} - 2p_j
$$

• $2p_j$ 是在哈迪-温伯格平衡下，一个随机个体在 SNP $j$ 上基因型值的群体平均。
• $Z_{ij}$ 表示个体 $i$ 在 SNP $j$ 上偏离群体平均的程度，它代表了该个体的"遗传独特性"。

def center_genotype_matrix(M, p):
    Z = M - 2 * p  # 广播机制，每列减去对应的 2p_j
    return Z

Z = center_genotype_matrix(M, p)
df_Z = pd.DataFrame(Z, index=tree_names, columns=snp_names)
print(df_Z.round(3))

          SNP_001  SNP_002  SNP_003  SNP_004  SNP_005  SNP_006  SNP_007  SNP_008  SNP_009  SNP_010
Tree_001     -0.8     -0.1      1.0      0.1     -1.0      0.9     -0.1     -0.9      0.1      1.0
Tree_002      0.2      0.9      0.0     -0.9      0.0     -0.1      0.9      0.1     -0.9      0.0
Tree_003     -0.8     -0.1     -1.0      1.1      1.0     -1.1     -0.1      1.1      0.1     -1.0
Tree_004      1.2     -1.1      0.0      0.1     -1.0      0.9     -1.1      0.1      1.1      0.0
Tree_005      0.2     -0.1      1.0     -0.9      0.0     -0.1      0.9     -0.9      0.1      1.0
Tree_006     -0.8      0.9      0.0      0.1      1.0     -1.1     -0.1      0.1     -0.9      0.0
Tree_007      0.2     -1.1     -1.0      1.1      0.0      0.9     -1.1      1.1      0.1     -1.0
Tree_008      1.2     -0.1      1.0      0.1     -1.0     -0.1     -0.1     -0.9      1.1      0.0
Tree_009     -0.8     -0.1      0.0     -0.9      1.0      0.9     -0.1      0.1     -0.9      1.0
Tree_010      0.2      0.9     -1.0      0.1      0.0     -1.1      0.9      0.1      0.1     -1.0

不同的 SNP 位点包含的遗传信息量是不同的。一个等位基因频率为 0.5 的 SNP，其多态性最高，包含的信息量最大；而一个频率为 0.99 的 SNP，几乎大家都有相同的等位基因，信息量很小。因此需要一个因子来对信息量进行标准化，这个因子就是所有 SNP 的预期方差之和。

$$
\text{标准化因子} = 2 \sum_{j=1}^{m} p_j (1 - p_j)
$$

其中 $p_j (1 - p_j)$ 正比于该 SNP 在随机交配群体中的期望方差。

def calculate_scale_factor(p):
    scale_factor = 2 * np.sum(p * (1 - p))
    return scale_factor

scale_factor = calculate_scale_factor(p)
print(f"标准化因子 = 2 * Σ p_j(1-p_j) = {scale_factor:.4f}")
print("\n各SNP对标准化因子的贡献 p(1-p):")
for i, snp in enumerate(snp_names):
    contrib = p[i] * (1 - p[i])
    print(f"  {snp}: {p[i]:.2f} * {1-p[i]:.2f} = {contrib:.3f}")

标准化因子 = 2 * Σ p_j(1-p_j) = 4.9500

各SNP对标准化因子的贡献 p(1-p):
  SNP_001: 0.40 * 0.60 = 0.240
  SNP_002: 0.55 * 0.45 = 0.248
  SNP_003: 0.50 * 0.50 = 0.250
  SNP_004: 0.45 * 0.55 = 0.248
  SNP_005: 0.50 * 0.50 = 0.250
  SNP_006: 0.55 * 0.45 = 0.248
  SNP_007: 0.55 * 0.45 = 0.248
  SNP_008: 0.45 * 0.55 = 0.248
  SNP_009: 0.45 * 0.55 = 0.248
  SNP_010: 0.50 * 0.50 = 0.250

最后，将中心化矩阵 $\mathbf{Z}$ 与其转置 $\mathbf{Z'}$ 相乘，然后除以标准化因子得到 G 矩阵：

$$
\mathbf{G} = \frac{\mathbf{ZZ'}}{2\sum_{j=1}^{m} p_j(1-p_j)}
$$

• $\mathbf{ZZ'}$ 是一个 $n \times n$ 的矩阵，其元素 $(\mathbf{ZZ'})_{ik}$ 粗略表示个体 $i$ 和个体 $k$ 在所有 SNP 上"遗传独特值"的相似度。
• 除以标准化因子后，G 矩阵的尺度与基于系谱的 A 矩阵相当，对角线平均值约为 1。

def calculate_G_matrix(Z, scale_factor):
    G = (Z @ Z.T) / scale_factor
    return G

G = calculate_G_matrix(Z, scale_factor)
print("\n基因组关系矩阵 G :")
df_G = pd.DataFrame(G, index=tree_names, columns=tree_names)
print(df_G.round(4))

          Tree_001  Tree_002  Tree_003  Tree_004  Tree_005  Tree_006  Tree_007  Tree_008  Tree_009  Tree_010
Tree_001    1.0707   -0.1414   -0.8485    0.2222    0.4848   -0.3232   -0.4040    0.3838    0.2424   -0.6869
Tree_002   -0.1414    0.6667   -0.2424   -0.5859    0.2828    0.2828   -0.6061   -0.2222    0.2424    0.3232
Tree_003   -0.8485   -0.2424    1.4747   -0.4848   -0.8283    0.5859    0.7071   -0.7273   -0.2626    0.6263
Tree_004    0.2222   -0.5859   -0.4848    1.3939   -0.1616   -0.9697    0.7677    0.7475   -0.4040   -0.5253
Tree_005    0.4848    0.2828   -0.8283   -0.1616    0.9091   -0.1010   -0.9899    0.4040    0.2626   -0.2626
Tree_006   -0.3232    0.2828    0.5859   -0.9697   -0.1010    0.9091   -0.3838   -0.6061    0.2626    0.3434
Tree_007   -0.4040   -0.6061    0.7071    0.7677   -0.9899   -0.3838    1.5556   -0.2828   -0.2222   -0.1414
Tree_008    0.3838   -0.2222   -0.7273    0.7475    0.4040   -0.6061   -0.2828    1.1111   -0.6465   -0.1616
Tree_009    0.2424    0.2424   -0.2626   -0.4040    0.2626    0.2626   -0.2222   -0.6465    1.0303   -0.5051
Tree_010   -0.6869    0.3232    0.6263   -0.5253   -0.2626    0.3434   -0.1414   -0.1616   -0.5051    0.9899

GBLUP 模型

GBLUP 是一个线性混合模型，它将利用 G 矩阵计算出的个体间亲缘关系，整合到传统的 BLUP 模型框架中，从而实现对个体遗传值（育种值）的预测。GBLUP 假设所有 SNP 的效应是随机的，并且都来自于同一个正态分布，即每个 SNP 都贡献了一点点效应，没有“主效”基因。这使得模型能够捕捉到由于减数分裂过程中随机的染色体分离和重组所带来的孟德尔抽样变异，从而更精确地估计育种值。

• 目标性状由大量基因共同控制，每个基因的产生的作用都很小。
• 所有 SNP 标记的效应值都来自同一个正态分布，平均值为 $0$，方差为 $\sigma_g^2$。这意味着模型假设没有起决定性作用的主效基因，所有 SNP 都对遗传变异有贡献。

模型的数学公式为：

$$
\mathbf{y} = \mathbf{Xb} + \mathbf{Zg} + \mathbf{e}
$$

• $\mathbf{y}$是观测到的表型值向量。
• $\mathbf{b}$是固定效应向量，包括试验地点、年份、品系或者其他需要校正的环境因素。数据预处理时可以通过计算方差去除。
• $\mathbf{X}$是与固定效应相关联的设计矩阵。
• $\mathbf{g}$是个体产生的随机遗传效应的向量，也就是要求解的 GEBV。
• $\mathbf{Z}$是与随机效应$\mathbf{g}$关联的设计矩阵。
• $\mathbf{e}$是随机残差效应的向量，代表了环境效应和无法被模型解释的其他变异。

这个公式实际上就是一个分解过程，将观测到的表型值分解为三个主要来源：

• $\mathbf{Xb}$:由已知的、系统性的固定效应（如不同地点、年份）造成的影响。
• $\mathbf{Zg}$:由个体的基因组造成的影响。$\mathbf{g}$ 就是我们想求的 GEBV。
• $\mathbf{e}$:由无法解释的随机误差（如微环境差异、测量误差）造成的影响。

使用前面杨树的例子进行计算，已知的值有：

假设在杨树种植地，前五棵树在 A 地点，后五棵树在 B 地点，得到固定效应设计矩阵$\mathbf{X}$：

这里需要使用混合模型方程 (Mixed Model Equations, MME)来进行计算：

$$
\begin{bmatrix}
\mathbf{X'X} & \mathbf{X'Z} \\
\mathbf{Z'X} & \mathbf{Z'Z} + \mathbf{G}^{-1}\lambda
\end{bmatrix}
\begin{bmatrix}
\hat{\mathbf{b}} \\
\hat{\mathbf{g}}
\end{bmatrix} = \begin{bmatrix}
\mathbf{X'y} \\
\mathbf{Z'y}
\end{bmatrix}
$$

这里还需要一个未知数

$$
\lambda = \frac{\sigma_e^2}{\sigma_g^2}
$$

即残差方差与加性遗传方差的比值， $\lambda$ 是一个“惩罚”或“缩放”因子。$\lambda$ 本身也是未知的，需要通过限制性最大似然法 (REML)从数据中迭代估计出来。

• 如果遗传方差很大（基因很重要），$\lambda$ 就很小，模型会更相信基因的作用。
• 如果残差方差很大（环境噪音很大），$\lambda$ 就很大，模型在估计 GEBV 时会更“保守”，使其向 0 收缩。

使用 Python 来进行计算：

import numpy as np
import pandas as pd
import scipy.linalg
from sklearn.model_selection import LeaveOneOut
import matplotlib.pyplot as plt
import seaborn as sns

def calculate_g_matrix(M):
    n, m = M.shape
    p = np.sum(M, axis=0) / (2 * n) # 每个SNP的等位基因频率 p
    Z = M - np.tile(2 * p, (n, 1)) #  中心化的基因型矩阵 Z
    scale_factor = 2 * np.sum(p * (1 - p))
    G = (Z @ Z.T) / scale_factor # 基因组关系矩阵 G
    return G

def solve_gblup(y, X, G, h2):
    n = len(y)
    if h2 <= 0: h2 = 1e-6
    if h2 >= 1: h2 = 1 - 1e-6
    vy = np.var(y)
    vg = vy * h2
    ve = vy * (1 - h2)

    lambda_val = ve / vg # 方差比 λ= (σ_e^2 / σ_g^2)

    # 简化版的MME求解，假设Z是单位矩阵I
    # LHS左侧矩阵
    C11 = X.T @ X
    C12 = X.T
    C21 = X
    C22 = np.eye(n) + np.linalg.inv(G) * lambda_val
    LHS = np.block([[C11, C12], [C21, C22]])

    # RHS右侧向量
    RHS = np.concatenate((X.T @ y, y))

    # 求解
    solution = np.linalg.pinv(LHS) @ RHS

    p = X.shape[1]
    b_hat = solution[:p]
    g_hat = solution[p:]

    return b_hat, g_hat

M = np.array([                                      # 基因型矩阵 M
    [0, 1, 2, 1, 0, 2, 1, 0, 1, 2], [1, 2, 1, 0, 1, 1, 2, 1, 0, 1],
    [0, 1, 0, 2, 2, 0, 1, 2, 1, 0], [2, 0, 1, 1, 0, 2, 0, 1, 2, 1],
    [1, 1, 2, 0, 1, 1, 2, 0, 1, 2], [0, 2, 1, 1, 2, 0, 1, 1, 0, 1],
    [1, 0, 0, 2, 1, 2, 0, 2, 1, 0], [2, 1, 2, 1, 0, 1, 1, 0, 2, 1],
    [0, 1, 1, 0, 2, 2, 1, 1, 0, 2], [1, 2, 0, 1, 1, 0, 2, 1, 1, 0]
], dtype=float)

y = np.array([15.2, 16.8, 14.5, 17.1, 15.9, 16.3, 14.8, 17.5, 15.6, 16.0]) # 表型向量 y

X = np.array([[1, 0]] * 5 + [[1, 1]] * 5, dtype=float) # 固定效应设计矩阵 X

G = calculate_g_matrix(M)
assumed_h2 = 0.6562
b_hat, g_hat = solve_gblup(y, X, G, h2=assumed_h2) # 固定效应解 b, GEBV解 g

results_df = pd.DataFrame({'GEBV': g_hat}, index=[f"Tree_{i+1:03d}" for i in range(10)])
print(results_df.sort_values(by='GEBV', ascending=False))

loo = LeaveOneOut()
y_observed, y_predicted = [], []

for train_index, test_index in loo.split(X):
    y_train, y_test = y[train_index], y[test_index]
    X_train, X_test = X[train_index], X[test_index]
    G_train_train = G[np.ix_(train_index, train_index)]
    G_test_train = G[np.ix_(test_index, train_index)]

    b_hat_train, g_hat_train = solve_gblup(y_train, X_train, G_train_train, h2=assumed_h2)

    g_hat_test = G_test_train @ np.linalg.pinv(G_train_train) @ g_hat_train
    predicted_pheno = (X_test @ b_hat_train) + g_hat_test

    y_observed.append(y_test[0])
    y_predicted.append(predicted_pheno[0])

accuracy = np.corrcoef(y_observed, y_predicted)[0, 1]
print(f"\nPrediction Accuracy (r): {accuracy:.4f}")

              GEBV
Tree_003      0.871591
Tree_005      0.777822
Tree_010      0.251176
Tree_009      0.041215
Tree_004     -0.008354
Tree_001     -0.045589
Tree_002     -0.045901
Tree_007     -0.171184
Tree_008     -0.637850
Tree_006     -1.032927

Prediction Accuracy (r): 0.0157

从结果看，Tree_003 是最值得选择的育种材料，其基因组估计育种值（GEBV）最高（+0.87），表明它拥有最优秀的生长基因潜力，尽管其表观长势不佳。然而，模型的预测准确性极低（相关系数 r ≈ 0.016），这说明虽然该模型能有效评估现有群体的内在遗传价值，但由于样本量过小，不具备预测新未知个体表现的泛化能力。

GBLUP 属于一类简单的机器学习框架，属于监督学习，线性模型，通过一次性求解大型矩阵方程组来得到最优解。