材料和方法

本研究中所有的统计分析都是使用R 的3.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\)。