Sillanpää2012 功能QTL定位

作者: 董八七 | 来源:发表于2019-03-14 17:57 被阅读5次

Sillanpää MJ, Pikkuhookana P, Abrahamsson S, et al (2012) Simultaneous estimation of multiple quantitative trait loci and growth curve parameters through hierarchical Bayesian modeling. Heredity 108:134–146. doi: 10.1038/hdy.2011.56

通过分层贝叶斯模型同时估计多个数量性状基因座和生长曲线参数

提出了一种新的分层数量性状基因座（QTL）定位方法，该方法使用多项式生长函数和多重QTL模型（在时间上没有依赖性）。该方法考虑基于群体的样品，其中个体已经针对某些动态性状进行表型分析（随着时间的推移）并且在给定的一组基因座处进行基因分型。所提出的方法的特定特征在于，代替平均功能曲线，每个个体具有其自己的功能曲线。此外，每个QTL可以通过其对一个或多个生长曲线参数的影响来修改个体的性状值的动态特征。该方法的明显优点包括：（1）假设时间无关的QTL和环境影响，（2）减少残差的自回归协方差结构的必要性和（3）使用变量选择方法的灵活性。作为该方法的副产品，还可以估计个体生长曲线参数的遗传力和遗传相关性，其被认为是潜在性状。为了在模型中选择与性状相关的基因座，我们使用了众所周知的贝叶斯自适应收缩技术的修改版本。我们通过分析来自模拟QTLMAS 2009数据集的500个个体的子样本，以及模拟重复和真实的苏格兰松（Pinus sylvestris）数据集，使用高度的时间测量作为感兴趣的动态特征来说明我们的方法。
关键词：功能定位;苏格兰松树; QTL; multitrait;贝叶斯模型; MCMC

介绍

已经提出了几种方法来定位影响动态性状的数量性状基因座（QTL）（即，表达随时间变化的性状）（参见Wu和Lin，2006年的综述）。即使在不同时间点测量的表型可以由不同的QTL组控制，但是时间点上的表型值通常是高度相关的。因此，已经提出重复测量框架用于随时间的性状测量的QTL分析（Lynch和Walsh，1998）。或者，可以将在不同时间点测量的性状视为单独的性状，并在多性状框架中共同分析。这里，就协方差函数而言，多性状框架的有效参数化提供了可行的方法（Macgregor等，2005; Lund等，2008）。然而，最常见的做法是使用一些数学函数来描述动态特征行为，然后定位QTL，这些QTL使用单个或多变量QTL定位来影响这个特殊函数。例如，逻辑增长函数（Ma等，2002; Wu等，2002,2003,2004），以及多项式函数（即多元回归模型）（Gee等，2003）和勒让德多项式（ Yang等人，2006; Yang和Xu，2007）已被提出用于此目的。逻辑斯蒂生长函数在生物学上是合理的（West等，2001）。作为一种批评，使用逻辑函数的回归仅适合于S形的增长轨迹，即单调递增时间函数（Yang和Xu，2007）。 Legendre和其他正交多项式拟合（用于协方差函数）也受到Pletcher和Geyer（1999）的批评。通常，功能的选择应基于特征轨迹的复杂性。即使已经提出了几种方法，但大多数方法仅限于单QTL模型或双QTL模型。例外情况包括Yang和Xu（2007），Min等人，2011年以及Heuven和Janss（2010）的方法。
对于生长的QTL进行了单独的年龄 - 年龄分析，以解决木本树木中QTL稳定性问题（Verhaegen等，1997; Conner等，1998; Kaya等，1999; Lerceteau等，2001）。。然而，据我们所知，Ma等（2004）是唯一一项将功能性QTL作图用于研究森林树种生长轨迹的研究。在他们的工作中，Ma等人（2004）指出，与替代QTL时间点分析相比，基于功能定位方法的QTL检测的统计功效增加。
功能性QTL作图方法通常用时间特异性QTL和环境效应模拟平均曲线行为（Yang和Xu，2007和Min等，2011）。这些方法中的个体特定变化被描述为与平均曲线行为的偏差，并且这些偏差取决于相邻时间点。 Gee等人（2003）和Heuven和Janss（2010）提供了这一共同主题的例外情况，其中所有与时间相关的行为都是通过个体特定的曲线参数来描述的，这些参数允许对QTL效应进行分层建模。这两个世界（等级和非等级）在概念上彼此非常不同。在Gee等人（2003）的参数化中，QTL效应不依赖于时间，它们影响曲线的形状而不是在特定时间点具有特定效应。为了描述随时间变化的功能曲线，我们在此考虑Gee等人（2003）的方法。作为对其方法的改进（以及Heuven和Janss，2010的方法），我们将整个问题表述为单一的层次模型。在我们的公式中，我们同时使用多重多QTL模型和模型选择，同时在贝叶斯框架中估计功能曲线和其他模型参数。

模型

让我们考虑基于群体的个体样本，其中数据样本已经针对某些动态特征进行了表型分类，并且在给定的一组标记基因座处进行了基因分型。虽然这代表了基于群体的单核苷酸多态性关联研究中的典型设计，但所提出的方法可直接应用于近交系中的回交和双单倍体，以及由近交系杂交产生的后代种群。处理缺失值时，我们完全忽略父母（连锁）信息，以便独立处理标记。这也意味着只有标记位置被认为是推定的QTL位置。有关替代方案，请参阅有关缺失基因型数据的小节。

表型模型随着时间的推移

对于每个个体 $i$ ，让我们假设 $y_{i, t}$ 是在时间点 $t$ ， $t=1, \ldots, T$ 测量的表型值。我们使用以下回归模型来描述一段时间内的表型行为：
$y_{i, t}=\beta_{0, i}+\beta_{1, i} a_{t, i}+\beta_{2, i} a_{t, i}^{2}+e_{i, t}$
式中， $\beta_{i}=\{\beta_{0, i}, \beta_{1, i}, \beta_{2, i}\}$ 是个体 $i$ 的曲线参数，假设误差项 $e_{i, t}$ 独立并且正态分布，在所有时间点都具有平均零和方差 $\sigma_{e}^{2}$ 。因为个体的曲线参数不同，我们预先指定 $\sigma_{e}^{2}$ 以改善我们的分层模型中的参数可识别性（如下所述）。请注意， $\sigma_{e}^{2}$ 描述了允许每个时间点的测量值偏离个体特定曲线（即数据和生长函数之间的一致性水平）。 $\sigma_{e}^{2}$ 的合适值取决于数据的类型。例如，对于增长数据，在我们的小模拟示例中使用 $\sigma_{e}^{2}=0.1$ ，对于实际数据分析和QTLMAS 2009数据分析 $\sigma_{e}^{2}=0.01$ （ $\sigma_{e}^{2}$ 的选定值不应太大，因为这可能导致残余误差错误地解释了所有QTL变化。数量 $a_{i, t}$ 是个体 $i$ 在时间点 $t$ 的年龄（在日历时间中，可以表示为与平均年龄的偏差;参见Gee等人，2003）。为简单起见，我们考虑所有个体的共同时间点和相同的年龄，对于所有时间点 $t=1, \ldots, T$ 和全部个体 $i$ ， $a_{i, t}=t$ 。

多性状QTL模型

我们将 $\beta_{i}$ 中的三个曲线参数视为三个潜在性状，并假设以遗传效应为条件，曲线参数是先验，且彼此相关。通过做出这样的假设，我们可以对曲线参数 $\beta_{i}$ 分层次地拟合多重QTL模型。对于每个个体 $i$ ，假设存在 $p$ 个加性标记基因座，其具有基因型值 $x_{i, j}$ ， $j=1, \ldots, p$ ，两个纯合子编码为0或1，杂合子编码为0.5。给定标记效应 $\left(B_{(k)}=\left\{B_{1(k)}, \ldots, B_{p(k)}\right\}, k=0,1,2\right)$ ，每个曲线参数 $\left(\beta_{k, i}, k=0,1,2\right)$ 建模为不同基因座的基因型 $x_{i, j}$ 的影响的线性组合（加权和）。
$\beta_{0, i}=\mu_{0}+\sum_{j=1}^{p} I_{j(0)} B_{j(0)} x_{i, j}+\epsilon_{i(0)}$
$\beta_{1, i}=\mu_{1}+\sum_{j=1}^{p} I_{j(1)} B_{j(1)} x_{i, j}+\rho_{10} \beta_{0, i}+\epsilon_{i(1)}$
$\beta_{2, i}=\mu_{2}+\sum_{j=1}^{p} I_{j(2)} B_{j(2)} x_{i, j}+\rho_{20} \beta_{0, i}+\rho_{21} \beta_{1, i}+\epsilon_{i(2)}$
式中， $\mu=\left\{\mu_{0}, \mu_{1}, \mu_{2}\right\}$ 是基线参数， $\epsilon_{i(k)}$ 是残差，服从均值为0，方差为 $\sigma^{2}_{\epsilon(k)}$ 的独立同正态分布。不同的残差方差表示在向量 $\sigma^{2}_{\epsilon}=\{\sigma_{\epsilon(0),}^{2} \sigma_{\epsilon(1)}^{2}, \sigma_{\epsilon(2)}^{2}\}$ 中。自回归项 $\rho=\{\rho_{10}, \rho_{20}, \rho_{21}\}$ 包含在模型中以考虑特征之间的残差依赖性，以便可以假设实际残差是独立的。自回归模型通常用于对时间序列数据中不同时间点之间的协方差进行建模。在这里使用相同的原理来模拟性状协方差（参见Bonney的D类模型，1986）。请注意，即使在模型中可见单向依赖性，也会自动引入双向依赖关系，因为 $\beta_{0, i}$ 和 $\beta_{1, i}$ 是模型参数而不是模型中的观测量。虽然假设具有非结构协方差矩阵的多变量正态分布残差的模型将是模拟这种现象的常用方法，但我们决定在计算基础上使用这种自回归模型。
在上述多性状QTL模型中，我们对每个基因座和每个性状 $I_{j(k)}$ 使用自己的指示变量，其中k=0,1或2。尽管这些指标提供了一种监测QTL后验占有的自然方法，如Pikkuhookana和Sillanpa（2009）所示，在模型中使用它们的真正原因是提高遗传力估计。有关详细信息，请参阅处理遗传性和遗传协方差/相关性的小节。
虽然这里没有明确显示，但是像区组效应这样的环境因素可以很容易地作为协变量包含在每个性状的QTL模型中（2-4）。在这些模型中，环境因素可能对不同的曲线特征产生不同的影响。或者，在QTL分析之前，可以尝试首先调整（恒定的或时间依赖的）环境因子的表型数据。这意味着初步分析的残差被作为连续QTL分析的表型。但是，这种调整很可能只能预先校正对截距（ $\beta_{0, i}$ ）有影响。一般来说，这种预矫正实践可能存在许多问题（参见Martinez等，2005），这些问题对于时间依赖的协变量可能更严重。

分层模型

上面给出的所有模型（1-4）同时被认为是更大的分层模型的一部分。我们将表型和标记数据分别表示为Y和X. 我们将模型参数联合表示为
$\theta=\{\beta_{1}, \ldots, \beta_{N}, \mu, B_{(0)}, B_{(1)}, B_{(2)}, I_{(0)}, I_{(1)}, I_{(2)},\rho, \sigma_{\epsilon}^{2}, \tau_{(0)}^{2}, \tau_{(1)}^{2}, \tau_{(2)}^{2}\}$
请注意，此向量包括模型中所需的所有未知参数（1-4）。后验分布 $P(\theta | X, Y)$ 与数据和参数的联合分布 $P(X, Y, \theta)$ 成比例。该联合分布可以描述为似然 $P(Y | \theta)$ 和先前 $P(\theta | X)$ 的乘积，其中似然（具有预选值 $\sigma^{2}_e$ ）是
$P(Y | \theta)=\prod_{i=1}^{N} \prod_{t=1}^{T} \frac{1}{\sqrt{2 \pi \sigma_{e}^{2}}} \exp \left(-\frac{1}{2 \sigma_{e}^{2}}\left(y_{i, t}-\beta_{0, i}-\beta_{1, i} a_{t, i}-\beta_{2, i} a_{t, i}^{2}\right)^{2}\right)$
先验是
$P(\theta | X)=\prod_{k=0}^{2} \prod_{i=1}^{N} P\left(\beta_{k, i} | \beta_{<k, i}, X, \mu_{k}, B_{(k)}, I_{(k)}, \rho, \sigma_{\in(k)}^{2}\right) \times$
$\times \prod_{k=0}^{2}\left[P\left(\mu_{k}\right) P\left(\sigma_{\epsilon(k)}^{2}\right) \prod_{j=1}^{p}\left[P\left(I_{j(k)}\right) P\left(B_{j(k)} | \tau_{j k}^{2}\right) P\left(\tau_{j k}^{2}\right)\right]\right] \times$
$\times P\left(\rho_{10}\right) P\left(\rho_{20}\right) P\left(\rho_{21}\right)$
这里，符号 $\beta_{<k, i}$ 指的是 $\beta_{i}$ 中出现在 $\beta_{k, i}$ 前的所有项。例如 $\beta_{1, i}$ 前面那一项是 $\beta_{0, i}$ 。先验的功能形式 $\mathrm{P}\left(\beta_{k, i} | \beta_{<k, i} X_{t} \mu_{k}, B_{(k)}, I_{(k)}, \rho, \sigma_{\in(k)}^{2}\right)$ 是模型(2–4)的残差的正态密度，均值为0，方差 $\sigma_{\in(k)}^{2}$ 。对于截距，
$\begin{array}{l}{P\left(\beta_{0, i} | \beta_{<0, i}, X, \mu_{0}, B_{(0)}, I_{(0)}, \rho, \sigma_{\epsilon(0)}^{2}\right)} \\ {=\prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi \sigma_{\epsilon(0)}^{2}}} \exp \left(-\frac{1}{2 \sigma_{\epsilon(0)}^{2}}\left(\beta_{0, i}-\mu_{0}-\sum_{j=1}^{p} I_{j(0)} B_{j(0)} x_{i, j}\right)^{2}\right)}\end{array}$
假设 $P(\sigma_{\in(k)}^{2})$ 中的每一个先验是Inverse-Gamma (0.001, 0.001)。每一个 $P\left(\mu_{k}\right)$ ， $P\left(\rho_{10}\right), P\left(\rho_{20}\right)$ 和 $P\left(\rho_{21}\right)$ 是 $N(0,100)$ 。Inverse-Gamma分布支持正范围内的值，并且上述正态分布相当平坦。因此，它们提供了适用于许多数据集的实用先验，而没有标准化。先验 $P\left(I_{j(k)}\right)$ ， $P\left(B_{j(k)} | \tau_{j k}^{2}\right)$ 和 $P\left(\tau_{j k}^{2}\right)$ 在下一节中介绍。