不确定基因型的可加模型及变量选择

Additive Model and Variable Selection for Uncertain Genotypes

摘要: 全基因组关联分析(GWAS)是研究复杂疾病相关位点的有效方法.在基因不确定情形下，传统方法利用基因填补方式估计基因概率，继而展开后续基因关联分析。我们对大样本基因考虑一个非参数可加模型对可加分量维数大而非零加性分量数目小的基因数据进行建模，其中加性分量利用B样条基函数的线性组合工具来近似拟合基因概率对性状表征的效应关系；选择非零分量是利用组Lasso惩罚来获得初始估计量。最后我们利用蒙特卡洛模拟证明，可加模型的组lasso方法在基因表达样本中的效果良好。

Abstract: Genome-wide association analysis (GWAS) is an effective method to study the associated loci of complex diseases. In the case of genetic uncertainty, the traditional method uses the gene filling method to estimate the gene probability, and then carries out the subsequent gene association analysis. We used a nonparametric additive model to model the data of large samples of genes with large additive component dimensions but small non-zero-additive component numbers. The additive component was used as a linear combination tool of B-spline basis function to approximate the effect relationship of gene probability on trait characterization. The group Lasso penalty was used to obtain the initial estimator for selecting the non-zero component. Finally, Monte Carlo simulation was used to demonstrate that the group Lasso method of the additive model performed well in gene expression samples.

1. 引言

近年来，全基因组关联研究(GWAS) [1] 被广泛关注，成功地应用于识别与复杂疾病相关的位点及农业畜牧业等重要经济性状相关的遗传变异基因。复杂疾病是环境和遗传变异因子的共同结果，一般涉及多个基因，而大部分疾病的相关基因位点未知，这就要求科学家开发合理的方法识别相关位点。GWAS就是寻找致病基因的方法之一，主要用于数量性状的分析，应用基因组中数以百万计的单核苷酸多态(SNP)，在全基因组水平上进行大样本分析找出影响复杂性状的基因变异位点。也就是说GWAS的主要目标是筛选出与特定疾病相关的基因、确定SNP与疾病表型的关联 [2]。Klein et al. (2005)最初利用该方法找出老年黄斑变性的相关致病基因 [3]，Sladek et al. (2007)确定4个与2型糖尿病有关的致病基因 [4]。相关实例应用陆续开发，在风湿性关节炎 [5]、食管癌 [6]、心血管疾病 [7] 和前列腺癌 [8] 等疾病以及农业(芝麻等)畜牧业等重要经济作物中有了相应研究成果，GWAS快速发展。

随着测序数据的发展，多国协作的HapMap计划的完成，海量的基因型数据被用于GWAS [9]，适用于大数据关联分析的算法也在不断更新。在大型样本的基因测序研究中，有时出现基因型不确定的情况 [10]，在这种情况下，利用“填补基因型方法”确定每个样本对应的基因型 [11] [12] [13]，继而展开后续的基因关联分析 [14] [15] [16] [17]。对于特定的SNP基因分型，通常采用的“最大概率法”和“剂量法”。前者指定概率最大的基因型作为样本基因进行分析估计，后者是用线性组合的方式计算基因型概率，计算方法分别为： ${\tilde{X}}_{ij}=\left\{g:p_{ijg}=\max \left\{p_{ij0},p_{ij1},p_{ij2}\right\}\right\}$ ； ${\bar{X}}_{ij}=p_{ij1}+2p_{ij2}$。

传统的关联分析主要考虑单个基因型对疾病表征关系的关联分析 [18]，但多个基因位点的基因型与特定疾病的关联分析及变量选择就很少 [19]。Li Q. [20]、Loley [21] 分别对多个基因位点对性状表征的传统线性回归模型、广义线性模型。在高维基因数据方面，成青 [22] 和刘璐 [19] 考虑高维数据中n个样本的p个基因位点的基因数据建立线性回归模型并结合Lasso罚函数进行变量选择。而这些方法没有考虑基因型的不确定，且认为基因对疾病表征关系是线性的。在遗传基因表达中，基因位点对疾病的效应关系可以是非线性的。对高维数据处理，一般传统的高维数据模型加罚函数可以初步对变量降维，在这方面已有大量的研究工作 [23] [24] [25] [26]。这里我们将在考虑高维基因数据下，部分基因位点不确定，探索一种利用B样条拟合基因位点对基本性状响应的光滑效应关系，进而建立可加模型 [27]，结合grouplasso罚函数对高维变量降维及非参估计以确定基因型与性状表征的关系。

本文的其余部分安排如下：在第2节中，我们介绍在基因不确定情形下的加性模型，我们将在第3节中详细介绍选择一致性证明结果；在第4节将我们的方法应用于模拟GWAS数据的分析。最后，我们作一些结束语。

2. 符号与模型

2.1. B样条曲线

在介绍我们的模型前，首先回顾B样条曲线概念。1972年Riesenfeld等人 [28] 提出B样条曲线，它可以通过一组标准B样条基函数的线性组合来有效逼近上述非线性函数。假定 $a={\eta }_0<{\eta }_1<\cdots <{\eta }_{K-1}<{\eta }_K=b$ 将区间 $\left[a,b\right]$ 分割为K个子区间，记 $I_{Kt}=\left[{\eta }_t,{\eta }_{t+1}\right),t=0,\cdots ,K-2$ 为前K − 1个子区间，对于第k个子区间，其p次样条基函数 ${\varphi }_{k,p}$ 定义如下

${\varphi }_{k,0}=\{\begin{array}{l}1,x\in \left[{\eta }_k,{\eta }_{k+1}\right]\\ 0,x\notin \left[{\eta }_k,{\eta }_{k+1}\right]\end{array}$

其递归公式为

${\varphi }_{k,p}\left(x\right)=\frac{x-{\eta }_k}{{\eta }_{k+p}-{\eta }_k}{\varphi }_{k,p-1}\left(x\right)+\frac{{\eta }_{k+p+1}-x}{{\eta }_{k+p+1}-{\eta }_{k+p}}{\varphi }_{k+1,p-1}\left(x\right),k>0$

假设 $S_n$ 为 $p\ge 1$ 次多项式样条空间， $\left\{{\varphi }_k,k\le d_n\right\}$ 是定义在 $S_n$ 上的一组标准基向量，则有

$f_{nj}\left(x\right)={\displaystyle \underset{k=1}{\overset{d_n}{\sum }}{\varphi }_k\left(x\right){\beta }_{jk}},1\le j\le p$

其中 $d_n=K+p$，在适当的条件下， $f_{nj}\left(x\right)$ 可以和好的逼近 $f_j\left(x\right)$ [29]。记 $\beta ={\left({\beta }_1^{\text{T}},\cdots ,{\beta }_p^{\text{T}}\right)}^{\text{T}}$， ${\beta }_j={\left({\beta }_{j1},\cdots ,{\beta }_{j{d}_n}\right)}^{\text{T}}$。

2.2. 基于基因不确定的可加模型

考虑等位基因为a和A的单核苷酸多态性。假设A是导致疾病的高风险等位基因。三种基因型分别表示为aa、Aa、AA，对应编码为0、1、2。定义 $X_{ij}$ 为第i个个体第j个点位的基因型，则对应的基因型概率定义为 $p_{ij}=P\left(X_{ij}=g\right),g=0,1,2$， $Y_i$ 为第i个个体的性状观测值。这里基因不确定数据，我们的 $X_{ij}$ 可采取传统填补基因型方法： ${\tilde{X}}_{ij}=\left\{g:p_{ijg}=\max \left\{p_{ij0},p_{ij1},p_{ij2}\right\}\right\}$ ； ${\bar{X}}_{ij}=p_{ij1}+2p_{ij2}$。

假设我们的数据向量 $\left(X_i,Y_i\right)$ 独立同分布服从 $\left(X,Y\right)$，Y是响应变量， $X={\left(X_1,\cdots ,X_p\right)}^{\text{T}}$ 是p维基因数据，考虑如下非线性可加模型：

$Y_i=\mu +{\displaystyle {\sum }_{j=1}^p{f}_j\left(X_{ij}\right)+{\epsilon }_i},i=1,\cdots ,n;j=1,\cdots ,p$

其中 $\mu$ 为截距项， $f_j(\ast )$ 为未知光滑函数， ${\epsilon }_i$ 为服从均值为0方差为 ${\sigma }^2$ 的随机误差。在高维稀疏性的假定下，这里的大部分 $f_j(\ast )\equiv 0,j=1,\cdots ,p$，我们的目的是找出不恒等于0的部分。结合B样条拟合光滑函数的工具，表达式等价于：

$Y_i\triangleq \mu +{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^{d_n}{\varphi }_k\left(X_{ij}\right){\beta }_{jk}}}+{\epsilon }_i\text{,}i=1,\cdots ,n;j=1,\cdots ,p;k=1,\cdots ,d_n$

这里的待估参数 $\beta =\left(\mu ,{\beta }_{11},\cdots ,{\beta }_{1d_n},\cdots ,{\beta }_{p1},\cdots ,{\beta }_{p{d}_n}\right)$。这里利用Bspline基函数线性组合工具，可视为其将非线性模型转化线性关系模型。

由于参数 $\beta$ 的维数是 $p{d}_n+1$，利用Bspline线性组合近似非线性函数 $f_j(\ast )$ 时， $d_n$ 很大则系数的稀疏性是我们考虑的范围，因而这里用稀疏组lasso选择重要分量和非零系数，这一部分内容我们在下一小节详细介绍。

2.3. grouplasso变量选择及参数估计

特定疾病对应与之相关的基因位点往往只有少部分，假定真实模型是稀疏的，也就是重要的协变量就是少数，只有小部分加性分量 $f_j$ 是非零的，记真实的非零分量指标集为 $M^{\ast }=\left\{j:f_j\ne 0\right\}$，这里指标集中的数目大小为 $s=\Vert M^{\ast }\Vert$。我们的目标是利用罚函数对高维基因数据降维，选择分量 $\hat{M}=\left\{j:\hat{\beta }\ne 0\right\}$，这里 $\hat{s}=\Vert M^{\ast }\Vert$，且 $\hat{s}\ge s$。即当样本量增大时，重要变量选入模型的概率趋于1，我们的模型具有选择一致性。

在高维基因数据的情况下，在模型上加groupLasso惩罚项有选择一致的特性，这里我们考虑惩罚最小二乘法：

$S_n\left(\mu ,{\beta }_n\right)={\displaystyle {\sum }_{i=1}^n{\left[y_i-\mu -{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^{d_n}{\varphi }_k}\left(x_{ij}\right){\beta }_{jk}}\right]}^{\text{2}}+{\lambda }_n}{\displaystyle {\sum }_{j=1}^p{\Vert {\beta }_{nj}\Vert }_2}$

记 ${\beta }_n={\left({\beta }_{n1}^{\text{T}},\cdots ,{\beta }_{np}^{\text{T}}\right)}^{\text{T}}={\left({\beta }_{11},\cdots ,{\beta }_{1d_n},\cdots ,{\beta }_{p1},\cdots ,{\beta }_{p{d}_n}\right)}^{\text{T}}$ 、 $1\le j\le p$， ${\lambda }_n$ 是惩罚参数。

为了模型的可识别性，加入约束条件 ${\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{k=1}^{d_n}{\varphi }_k\left(x_{ij}\right){\beta }_{jk}}}=0,j=1,\cdots ,p$。对 $S_n\left(\mu ,{\beta }_n\right)$ 极小化求参数：

$\hat{\beta }=\arg \min \left\{{\displaystyle {\sum }_{i=1}^n{\left(y_i-\mu -{\displaystyle {\sum }_{j=1}^p{\displaystyle {\sum }_{k=1}^{d_n}{\varphi }_k}\left(x_{ij}\right){\beta }_{jk}}\right)}^2}+{\lambda }_n{\displaystyle {\sum }_{j=1}^p{\Vert {\beta }_{nj}\Vert }_2}\right\}$

通过GroupLasso估计出每组参数的估计，根据稀疏性假定，其中只有少部分 ${\Vert {\hat{\beta }}_{nj}\Vert }_2\ne 0,j=1,\cdots ,p$，

等价于 $f_j(\ast )$ 不恒等于，也就是说对应的第j个基因位点是与疾病存在显著关系的。

3. 模拟结果

利用蒙特卡罗模拟数值，比较基因型数据 $X_i$ 采用剂量法与最大可能概率法在线性模型与非线性可加模型的表现。这里 $X_{ij}$ 中的三个概率( $p_{ij0},p_{ij1},p_{ij2}$ )通过将基因不确定率的情况考虑在内的三项分布生成。这里的是三项分布可以通过三种频率对模拟生成特定频率生成，在模拟设定基因确定率和不确定基因频率则为r和1 − r。有不确定率的存在，一个基因位点会有三个基因型(aa、Aa、AA)的概率。基因型的检测频率，如表1所示。

Table 1. Frequency of genotype detection

表1. 基因型的检测频率

下面我们通过两个例子比较说明可加模型的groupLasso罚函数的变量选择的优良性。

例1：从下面的模型中产生数据：

$Y_i=\mu +{\displaystyle {\sum }_{j=1}^p{f}_j\left(X_{ij}\right)+{\epsilon }_i},i=1,\cdots ,n$

这里取 $\mu =2$， $f_1\left(x\right)=-x$ ； $f_2\left(x\right)=2x^2$ ； $f_3\left(x\right)=-3\log \left(x\right)$ ； $f_j\left(x\right)\equiv 0,4\le j\le 100=p$ ；这里非零分量个数 $s=3$， $p=100$，样本量 $n=1000$，随机误差服从标准正态分布。模型超参数 $\lambda$ 的选取采用十折交叉验证进行筛选，数据按照7:3比例随机分为训练集与测试集。统计变量正确选择平均个数(TM)，错误选择平均个数(FM)以及测试集均方误差(TMSE)如表2。

Table 2. Results of variable selection

表2. 变量选择结果

根据表1结果可以知道，在模型一样的情况下，剂量法相较于最大概率法能更多的识别出重要变量以及拥有更小的均方误差，与此同时识别变量错误的情况也偏多；当然线性模型应用范围比较局限，在存在非线性部分时，其变现结果很差。

例2：从下面的模型中产生数据：

$Y_i=\mu +{\displaystyle {\sum }_{j=1}^p{f}_j\left(X_{ij}\right)+{\epsilon }_i,}i=1,\cdots ,n$

这里取 $f_1\left(x\right)=5x$， $f_2\left(x\right)=3{\left(2x-1\right)}^2$ ； $f_3\left(x\right)=4\sin \left(2\pi x\right)/\left(2-\sin \left(2\pi x\right)\right)$ ；

$f_4\left(x\right)=6\left(0.1\sin \left(2\pi x\right)+0.2\cos \left(2\pi x\right)+0.3\sin {\left(2\pi x\right)}^2+0.4\cos {\left(2\pi x\right)}^3+0.5\sin {\left(2\pi x\right)}^3\right)$ ； $f_5=\cdots =f_p=0$。这里

非零分量个数 $s=\text{4}$， $p=100\text{0}$ 和三种不同的样本容量：n=100，200和1000，随机误差服从标准正态分布。这里的数据生成模型和Huang的一样，但是我们这里基于基因不确定情形。针对这个模型比较GroupLasso和GroupSCAD惩罚的优越性，得到下面表3。

Table 3. Results of variable selection

表3. 变量选择结果

通过对表3两种变量选择方法进行比较，采用最大概率法时，GroupLasso方法与GroupSCAD方法变量选择的效果差距不大；从错选变量(FM)个数看，Group Lasso方法倾向于筛选出更少的伪变量，特别在用剂量法时，Group Lasso错选变量个数明显小于Group SCAD方法。

4. 结束语

复杂疾病与基因位点通常存在非线性关系，全基因组测序建立复杂疾病与基因位点之间关系是一种有效途径，但是全基因组的准确测序也是一个耗时耗力的庞大工程。在部分基因位点基因型不确定下，采用“最大概率法”或“剂量法”来填补基因数据，GroupLasso以及GroupSCAD方法对非线性可加模型进行重要基因位点识别是一种有效途径。

本文在基因不确定情形下，基于GWAS创新一种基因表达方法，不同于传统的基因填补估计概率方法，结合了B样条工具拟合加性分量建立非参数可加模型，考虑了大样本基因的环境下进行组Lasso变量选择以获得初始估计；最后利用数值模拟对比证明了我们方法的优良性。在基因测序中，我们非参加性模型更灵活，对模型的假设要求更低，适用范围更广。