计算亲缘关系矩阵（Genetic Relationship Matrix，GRM）

定义

亲缘关系矩阵（GRM）是一种用于量化个体间遗传相似性的统计矩阵。在基因组学和遗传学研究中，GRM 用于评估个体之间的遗传相关性，并帮助控制由于遗传结构导致的统计偏差。

一些背景知识

  • **IBD (Identical by descent)**：血缘同源
  • **IBS (Identical by state)**：状态同源

通常我们会讲，如果一段DNA在两个或者多个个体中均有一致的核苷酸序列，那么可以定义该DNA片段属性为状态同源，即 IBS 。

在两个或者多个个体中的一个IBS片段，如果遗传自一个共同的祖先个体，且没有发生重组，那么该片段也是血缘同源的，即 IBD 。

一段DNA，如果是 IBD，那么肯定也是 IBS；如果片段不是 IBD，也有可能是 IBS，因为不同个体中发生了相同的突变，或者重组没有改变片段等原因所致。

个体间的亲缘关系(relationship)的估计

  • IBD 在公式中体现为基于系谱信息的亲缘关系估计，强调的是共享的祖先DNA片段的概率。
  • IBS 在公式中则表现为通过分子标记计算的相似性，基于基因型是否在特定位点上相同。

1. 基于IBD的亲缘关系估计：

  • 系谱法：通过系谱分析，估计两个个体间的亲缘关系，通常通过计算两个个体共享的祖先DNA片段的比例来实现。这里，亲缘关系矩阵 $\mathbf{R}$ 的元素 $R_{ij}$ 表示个体 $i$ 和 $j$ 之间的亲缘关系，可以通过如下公式估计：

$$
R_{ij} = \text{IBD}(i, j) = \sum_{l=1}^{L} p_l
$$

  • $L$ 是所有可能的祖先路径数。
  • $p_l$ 是路径 $l$ 的概率，即在这条路径下，个体 $i$ 和 $j$ 共享同一段祖先DNA的概率。

2. 基于IBS的相似性估计：

  • 分子标记法：通过分子标记数据，估计两个个体间的基因相似性，常用的方式是计算两个个体间处于IBS（Identical by State）状态的等位基因的比例。在公式上，个体 $i$ 和 $j$ 之间的相似性可以表示为：

$$
S_{ij}=\frac{1}{M}\sum_{k=1}^{M} \delta(x_{ik},x_{ik})
$$

  • $S_{ij}$ 表示个体 $i$ 和 $j$ 之间的相似性。
  • $M$ 是标记位点的总数。
  • $x_{ik}$ 和 $x_{jk}$ 分别是个体 $i$ 和 $j$ 在第 $k$ 个位点上的基因型值。
  • $\delta(x_{ik}, x_{jk})$ 是一个指示函数，当 $x_{ik}$ 和 $x_{jk}$ 处于相同状态（即IBS）时，值为1，否则为0。

用途

  • 遗传关联分析: 在全基因组关联研究（GWAS）中，GRM 有助于识别与特定表型相关的遗传变异。
  • 混合模型分析: 用于控制遗传背景和家系关系，减少由于遗传结构导致的假阳性。
  • 遗传多样性研究: 评估群体内部和群体之间的遗传多样性。

计算过程

步骤一、获取基因型数据

GRM的构建始于基因型数据，通常以矩阵的形式表示。假设基因型矩阵为$X$，其中每一行代表一个个体，每一列代表一个SNP（单核苷酸多态性位点）。基因型数据通常编码为 0, 1, 2，分别表示该位点的等位基因个数。这些数据通常以 SNP（单核苷酸多态性）为基础。数据格式可以是 PLINK 格式、VCF 格式等。

步骤二、标准化基因型数据

这里假设 $Z$ 矩阵是对 $M$ 矩阵进行了中心化，对 $M$ 矩阵中心化的方法：
$$ Z = M - P$$

$M$ 矩阵包括了每个样本、每个位点次等位基因的含量。为了计算方便，其中每个元素都减掉了1，即 $M=M^*-1$，由0、1、2变为-1、0和1形式。这样做的目的是为了将每列的均值标准化为0。

$$ P = 2(p_i-0.5) $$

$P$ 矩阵包括了每个位点次等位基因的频率。减去0.5的原因在于，$M$ 矩阵中每个 SNP 位点的基因型值被减去了1，因此原来2倍的次等位基因频率需要减去1（即 $2\times0.5$）。这是因为每个位点的等位基因含量是由两个等位基因组成，所以总和应该是 $2 \times p_i-1$。通过这样的方式，$P$ 矩阵实际上反映了每个 SNP 位点的次等位基因频率偏离0.5的程度。

以一个简单的例子来说明：

1. 构建M矩阵

首先，我们有一个基因型矩阵，它记录了每个个体在每个SNP位点上的基因型。

个体 SNP位点
Ind1 AA
Ind2 AA
Ind3 TT

在这个例子中，小写字母 ‘t’ 代表出现频率较低的等位基因(MAF，Minor Allele Frequency)。转换成数字表示法，我们可以将基因型转换为次等位基因的数量：

个体 SNP位点
Ind1 0
Ind2 0
Ind3 2

这样，我们就得到了M矩阵的一列：

2. 中心化 $M$ 矩阵

为了便于计算，通常会将M矩阵中的每个元素减去1，使得基因型数据变为-1、0、1的形式。这样做之后，M矩阵变为：

$$ M = \begin{bmatrix} 0 & 0 & 2 \end{bmatrix}^\top $$

3. 构建 $P$ 矩阵

$P$ 矩阵用于记录每个位点次等位基因的频率。在这个例子中，次等位基因 t （一般用小写字母表示次等位基因）的频率为 $ p = \frac{2}{6} = \frac{1}{3} $。根据公式$ P_i = 2(p_i - 0.5) $，我们计算得到：

$$ P = 2\left(\frac{1}{3} - 0.5\right) = 2\left(\frac{1}{3} - \frac{3}{6}\right) = 2\left(-\frac{1}{6}\right) = -\frac{1}{3} $$

因此，$P$ 矩阵为：

$$ P = \begin{bmatrix} -\frac{1}{3} & -\frac{1}{3} & -\frac{1}{3} \end{bmatrix}^\top $$

4. 计算Z矩阵

Z矩阵是M矩阵减去P矩阵的结果，即：

$$ Z = M - P $$

因此，

$$ Z = \begin{bmatrix} -1 \ -1 \ 1 \end{bmatrix} - \begin{bmatrix} -\frac{1}{3} \ -\frac{1}{3} \ -\frac{1}{3} \end{bmatrix} = \begin{bmatrix} -1 - (-\frac{1}{3}) \ -1 - (-\frac{1}{3}) \ 1 - (-\frac{1}{3}) \end{bmatrix} = \begin{bmatrix} -\frac{2}{3} \ -\frac{2}{3} \ \frac{4}{3} \end{bmatrix} $$

5. 中心化的效果

通过减去P矩阵，我们得到了 $Z$ 矩阵。此时， $Z$ 矩阵每一列的基因型效应和为0，这确保了每列的均值为0，实现了中心化的效果。

具体而言，对于 $Z$ 矩阵中的每个值，我们计算其平均值：

$$ \text{Mean of } Z = \frac{-\frac{2}{3} - \frac{2}{3} + \frac{4}{3}}{3} = \frac{0}{3} = 0 $$

这就意味着，通过中心化处理，我们确保了 $Z$ 矩阵的每一列的均值为0，达到了中心化的目标。

步骤三、计算 GRM 矩阵

这里根据VanRaden (2008) 提出的三种计算GRM的方法：

方法1：基于等位基因频率的标准化方法

通过对基因型进行标准化处理，并以次要等位基因频率来调整计算的权重，从而得到个体间的基因关系。

$$
\mathbf{G} = \frac{\mathbf{Z}\mathbf{Z}^T}{2 \sum p_j(1 - p_j)}
$$

其中：

  • $\mathbf{Z}$ 是标准化后的基因型矩阵， $\mathbf{Z}{ij} = x{ij} - 2p_j$，其中 $x_{ij}$ 是第 $i$ 个个体在第 $j$ 个SNP的基因型值，$p_j$ 是第 $j$ 个SNP的等位基因频率。
  • $p_j$ 是次要等位基因频率（MAF）。

方法2：简单加和法

该方法直接使用基因型矩阵来计算个体间的关系，不考虑等位基因频率的调整，是最简单的计算GRM的方法。但是未考虑SNP的等位基因频率差异，可能导致计算结果偏差较大，特别是在MAF不均匀时。

$$
\mathbf{G} = \frac{\mathbf{M}\mathbf{M}^T}{N}
$$

其中：

  • $\mathbf{M}$ 是未经标准化的基因型矩阵。
  • $N$ 是 SNP 的总数。

方法3：基于等位基因频率的方差化方法

第二种方法是第一种方法的一个变体，利用每个位点期望方差的倒数，对标记效应值进行加权。

$$ \mathbf{G} = \mathbf{ZDZ’} $$

其中，$\mathbf{D}$ 是一个对角矩阵，其中 $D_{ii}$ 等于：

$$ D_{ii} = \frac{1}{m[2p_i(1 - p_i)]} $$

$\mathbf{D}$ 是一个对角矩阵，其中每个对角元素 $D{ii}$ 是根据每个 SNP 位点的等位基因频率 $p$ 计算的。具体来说，每个对角元素 $D_{ii}$ 是该位点期望方差的倒数，即 $\frac{1}{m[2p_i(1 - p_i)]}$ ，其中 $m$ 是SNP的总数。

这种方法通过将每个SNP位点的效应值乘以该位点期望方差的倒数，对标记效应值进行加权。这样可以确保每个SNP位点对遗传关系矩阵 $\mathbf{G}$ 的贡献与其方差成反比，从而实现对不同SNP位点的加权处理。

要使用 Plink 计算遗传关系矩阵 (GRM)，并指定上述参数，可以运行以下命令：

1
2
3
4
5
plink --bfile example_data \
      --make-grm-gz \
      --numRandomMarkerforSparseKin 2000 \
      --relatednessCutoff 0.125 \
      --out example_grm

--numRandomMarkerforSparseKin=2000

在处理大规模数据集的时候，考虑到内存占用等问题，从所有的标记中随机选择一部分来构建一个简化版本的稀疏矩阵（SparseGRM）。

通常用 --numRandomMarkerforSparseKin=2000 指定了用于计算稀疏亲缘矩阵的随机标记数量，这样做的目的是减少计算量。

--relatednessCutoff=0.125

以--relatednessCutoff=0.125为例，设置亲缘关系的截断值，用于定义个体间亲缘关系的最小可接受水平，任何个体对之间的亲缘关系小于等于 0.125 的都会被视为无关或相关性不足。

亲缘关系系数（或亲缘矩阵的元素） $\phi_{ij}$ 用来衡量个体 $i$ 和 $j$ 之间的遗传相似度。如果 $\phi_{ij}$ 小于或等于 0.125，则认为个体 $i$ 和 $j$ 之间没有足够的遗传相关性。这意味着如果个体 $i$ 和 $j$ 之间的亲缘关系系数小于或等于 0.125，则它们被认为是没有显著亲缘关系的个体。

在遗传学中，0.125 通常对应于近亲关系，比如第一代堂兄弟姐妹或祖孙之间的亲缘关系。因此，这个截断值可以帮助过滤掉那些遗传上关系较远的个体对，使得进一步的分析更加集中于那些遗传上更接近的个体

Reference

VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008 Nov;91(11):4414-23. doi: 10.3168/jds.2007-0980. PMID: 18946147.