【翻译】RAINBOW:采用新型SNP-set方法的基于单倍型的全基因组关联分析【第二部分:材料和方法】

材料和方法

本研究中所有的统计分析都是使用R3.6.0版[31]进行的,图像是使用R包ggplot 3.2.1版[32]制作的。我们的R包RAINBOWR是使用R包Rcpp 1.0.2版[33-35]和RcppEigen 0.3.3.5.0版[36]实现的,以减少下面描述的多核混合效应模型所需的计算时间。本研究的整体模拟框架在S1图中以流程图的形式展示。

RAINBOW的方法

在这个部分中,我们描述了RAINBOW的基础思想。

RAINBOW模型

RAINBOW使用的模型(model)如下:

\[{y}= {X}{\beta}+ {Z}_c{u}_c+ {Z}_{r_i}{u}_{r_i}+ \epsilon \tag{1} \]

其中,\({y}\)是一个\(n\times1\)的表型值向量,\({X}{\beta}\)是一个\(n\times1\)的固定效应(fixed effect)向量,包括一个截距(intercept)、一个修正种群结构(correct the population structure)的项和其他协变量(covariates),\({Z}_c{u}_c\)\({Z}_{r_i}{u}_{r_i}\)\(n\times1\)的随机效应(random effect)向量,\(\epsilon_,\)是一个\(n\times1\)的残差向量。这里\({\beta}\)是一个\(p\times1\)的固定效应向量,其中\(p\)是固定效应的数量。\({u}_c\)\({u}_{r_i}\)分别是基因型值的\(m_c\times1\)\(m_{r_i}\times1\)向量,其中\(m_c\)是加性多遗传效应的基因型数量,\(m_{r_i}\)是感兴趣的第\(i\)个SNP-set的基因型数量。\({X}\)\({Z}_c\)\({Z}_{r_i}\)分别对应于\({\beta}\)\({u}_c\)\({u}_{r_i}\)的,大小尺寸分别为\(n\times p\)\(n\times m_c\)\(n\times m_{r_i}\)的矩阵。如以下方程2所示,我们假设多遗传效应\(u_c\)遵循多元正态分布,其方差-协方差矩阵与加性分子关系矩阵(multivariate normal distribution)\({K}_c\)成正比。

\[{u}_c \sim MVN({0},{K}_c\sigma^2_c) \tag{2} \]

其中\(\sigma^2_c\)是在“方差分量的估计”(Estimation of variance components)部分中要估计的加性遗传方差,这里的\(m_c\times m_c\)矩阵\({K}_c={A}\),其中\({A}\)是根据标记基因型数据估计的已知加性遗传关系矩阵\(W_c\)[37]。

我们还假设,来自第\(i\)个感兴趣的SNP-set的随机效应\({u}_{r_i}\)遵循多变量正态分布,其方差-协方差矩阵与Gram矩阵\({K}_{r_i}\)

\[{u}_{r_i} \sim MVN({0},{K}_{r_i}\sigma^2_{r_i}) \tag{3} \]

其中\(\sigma^2_{r_i}\)是在“方差分量的估计”一节中要估计的第\(i\)个SNP-set的遗传方差,\({K}_{r_i}\)是由属于第\(i\)个SNP-set的标记基因型数据\({W}_{r_i}\)估计的已知的\(m_{r_i}\times m_{r_i}\)Gram矩阵。我们为Gram矩阵提供了一个线性核、一个指数核和一个高斯核,在使用线性核的情况下,我们可以实现更快的计算(S1附录的补充说明)[24]。

最后,假设残差项完全独立地遵循正态分布,如以下公式所示:

\[\epsilon \sim MVN({0},{I}_n\sigma^2_e) \tag{4} \]

其中\({I}_n\)是一个\(n\times n\)的识别矩阵(identity matrix),\(\sigma_e^2\)将在“方差分量的估计”一节中进行估计。

方差分量的估计

方差分量(variance components)通过最大似然(maximum-likelihood, ML)[26, 38]和受限最大似然(resitricted maximum-likelihood)[39]来估计。在这里,我们解释了如何为一般\({K}_{r_i}\)获得方程1的ML和REML估计值。

首先我们通过以下算法估计遗传方差(\(\sigma^2_c\)\(\sigma^2_{r_i}\))之间的权重(我们将其定义为\(w_c\)\(w_{r_i}\))。

  1. \(w_c\)\(w_{r_i}\)设置初始参数:

    \[w_c = w_{r_i} = \frac{1}{2} \tag{5} \]

  2. 计算下面的\(n\times n\)矩阵\({K}_{s}\)

    \[{K}_{s} = {Z}_c{K}_c{Z}_cw_C \space+\space {Z}_{r_i}{K}_{r_i}{Z}_{r_i}^Tw_{r_i} \tag{6} \]

  3. 通过使用EMMA(efficient mixed model association)或者GEMMA(genome-wide efficient mixed association)[40, 41]来解下面的单核线性混合模型(linear mixed model, LMM)。

    \[{y} = {X}{\beta} + {u}_s + \epsilon \tag{7} \]

    其中

    \[{u}_s \sim MVN({0},{K}_s\sigma_s^2) \tag{8} \]

  4. 使用估计参数\(\hat{\beta}\)\(\hat{\sigma}_s^2\)\(\hat{\sigma}^2_e\)计算方程7中的完整对数似然\(l_F\))或者受限对数似然\(l_R\)):

    \[l_F({y};\hat{\beta},\hat{\sigma}_s,\hat{\delta}) = \frac{1}{2} \left[ - n\log \left( 2\pi\hat{\sigma}_s^2 \right) - \log \left| \hat{{H}} \right| - \frac{1}{\hat{\sigma}_s^2} \left( {y}-{X}\hat{\beta} \right)^T \hat{{H}}^{-1} \left( {y}-{X}\hat{\beta} \right) \right] \tag{9} \]

    \[l_R = l_F({y};\hat{\beta},\hat{\sigma}_s,\hat{\delta}) + \frac{1}{2} \left[ p\log \left( 2\pi\hat{\sigma^2_s} \right) + \log \left| {X}^T{X} \right| - \log \left| {X}^T\hat{{H}}^{-1}{X} \right| \right] \tag{10} \]

    这里的\(\hat{{H}}\)是:

    \[\hat{{H}} = \frac{\hat{{V}}}{\hat{\sigma}^2_s} = {K}_s + \hat{\sigma}{I}_n \tag{11} \]

  5. 通过重复步骤2-4[42],使用L-BFGS优化方法优化\(w_c\)\(w_{r_i}\)来最大化完整/受限对数似然。

在估计出了权重\(w_c\)\(w_{r_i}\)之后,我们通过EMMA/GEMMA方法,使用\(\hat{w_c}\)\(\hat{w}_{r_i}\)估计了模型方程7和8的方差分量(\(\sigma_s^2\)\(\sigma^2_e\))。然后我们得到了\(\hat{\sigma}^2_c=\hat{w_c}\hat{\sigma}_s^2\)\(\hat{\sigma}^2_{r_i}=\hat{w_{r_i}}\hat{\sigma}_s^2\)

如上所述,我们的拟合方法是一种两步法(two-step approach),首先估计遗传方差的权重,然后用估计得到的权重,通过EMMA/GEMMA来估计方程7和8中所示模型的方差分量。在此之外,一些通过AIREML(average information REML,平均信息受限最大似然)[43]直接估计方程1的方差分量的拟合方法也已经被提出,并在一些包/软件[44, 45]中中实现。与通过AIREML的直接估计方法相比,我们的两步法优势在于权重的搜索空间限于区间\([0,1]\),并且即使在遗传性(heritability)太低/太高的情况下,结果的收敛性(convergence)也相对有保证。

GWAS的似然比检验

为了检验每个SNP-set的显著性,我们对\(\sigma^2_{r_i}\)是否为0进行了似然比检验(likelihood ratio test, LR test)。作为零假设(null hypothesis),假设以下模型不包括SNP-set效应项。

\[{y} = {X}{\beta} + {Z}_c{u}_c + \epsilon \tag{12} \]

相反的,作为替代假设模型(alternative hypothesis model),假设方程1的多核线性混合效应模型(multi-kernel linear mixed model, MKLMM)。因此,我们在估计每个SNP-set的方差分量后计算了以下偏差(deviance)。

\[D = 2\times \left( \hat{l}_{R,model} - \hat{l}_{R,null} \right) \tag{13} \]

其中,\(\hat{l}_{R,model}\)是方程1模型的首先对数似然的最大值,\(\hat{l}_{R,null}\)是方程12模型的受限对数似然的最大值。

最后,我们测试了\(\sigma_{r_i}^2\)的显著性,并通过假设方程13中的偏差遵循具有不同*度的两个卡方分布的混合 来计算p值[47,48]。

\[D \sim \pi_0\chi^2_0 + \left( 1-\pi_0 \right) \chi_0^2 \tag{14} \]

其中\(\pi_0\)是混合参数,这里我们使用的值为\(\pi_0=1/2\)

【翻译】RAINBOW:采用新型SNP-set方法的基于单倍型的全基因组关联分析【第二部分:材料和方法】

上一篇:SAP-PP T-Code集合


下一篇:dremio 升级的一个实践