Evaluation of predictive accuracy of gene expression in pigs using machine learning models
-
摘要:目的
通过对比不同机器学习模型利用基因顺式单核苷酸多态性(Single-nucleotide polymorphisms, SNP)预测猪的基因表达量的效果，探究基因顺式遗传力(Cis-heritability, cis-h2)和顺式SNP(Cis-SNP)数量与不同模型预测准确性间的关系。
方法基于PigGTEx项目猪的肌肉组织样本的蛋白编码基因，使用18种不同机器学习模型，将基因转录起始位点±1 Mb范围内的cis-SNP用于训练，评估每种模型的预测准确性。
结果机器学习模型的预测准确性与基因cis-h2间存在正相关，弹性网络回归模型和Lasso回归模型整体预测准确性最高，R2平均值分别为
0.0362 和0.0358 ；一定范围内，模型预测准确性与基因cis-SNP数量间存在正相关。结论使用机器学习模型预测猪基因表达的准确性受基因cis-h2和基因cis-SNP数量影响较大，根据不同基因的cis-h2和cis-SNP数量选择合适的机器学习模型预测猪的基因表达量有利于提高预测准确性。
Abstract:ObjectiveBy comparing the performance of various machine learning models in predicting gene expression in pigs utilizing single nucleotide polymorphisms (SNP), we investigated the relationship among cis-heritability (cis-h2), the number of cis-SNPs and the prediction accuracy.
MethodBased on the protein encoding genes of pigs derived from muscle tissue of the PigGTEx project, we trained 18 distinct machine learning models by employing cis-SNPs located within a ±1 Mb window from the transcription start sites of genes. Subsequently, we evaluated the prediction accuracy of each model.
ResultThere wasa positive correlation between the prediction accuracy of machine learning models and the cis-h2 values of genes. Notably, the Elastic Net regression model and the Lasso regression model exhibited the highest overall prediction accuracy, with mean R2 values of
0.0362 and0.0358 , respectively. Furthermore, there was a positive correlation between the prediction accuracy of these machine learning models and the number of cis-SNPs around the genes within certain range.ConclusionThe accuracy of utilizing machine learning models to predict gene expression in pigs is largely influenced by both cis-h2 and the number of cis-SNPs. Therefore, selecting an appropriate machine learning model tailored to the specific cis-h2 values and the number of cis-SNPs of different genes of pigs can be advantageous in enhancing prediction accuracy.
-
猪是重要的农业动物，其重要经济性状的遗传解析对于提高猪生产效率至关重要。目前，全基因组关联分析 (Genome-wide association studies, GWAS)已成为挖掘畜禽性状重要分子标记和功能基因的主要手段之一[1]，人们使用该方法陆续进行了猪乳头数[2]、背标厚度[3-4]、内脏重量[5]、产仔数[6]等重要生长、繁殖性状的关键位点分析。随着多组学技术普及以及猪基因型−组织表达(Pig genotype-tissue expression, PigGTEx)项目的开展[7]，我们能够将基因组和转录组等多组学数据进行联合分析，找到与分子表型相关的数量性状位点[8]，更有助于解析猪复杂性状的分子遗传机制。
作为GWAS后分析的重要方法，全转录组关联研究(Transcriptome-wide association study, TWAS)被提出后获得了广泛的应用。TWAS能将基因表达、GWAS检测出的显著变异以及性状关联起来。TWAS分析的其中关键一步是在参考群体中构建预测基因表达的预测模型，以用于预测GWAS群体中个体的基因表达水平[9]。提高TWAS预测模型的准确性有助于提升TWAS分析的统计效力。在现有的TWAS软件中，PrediXcan软件使用弹性网络的线性回归模型拟合基因顺式单核苷酸多态性(Single nucleotide polymorphism, SNP)以预测组织中的基因表达水平[10]，Fusion软件使用贝叶斯稀疏线性混合模型(Bayesian sparse linear mixed model, BSLMM)预测组织中的基因表达水平[11]。PigGTEx项目使用S-PrediXcan和S-MultiXcan中的嵌套交叉验证弹性网络模型，利用基因的顺式SNP(cis-SNP)预测表达量并用于后续TWAS分析[7]。然而，目前在猪等畜禽群体中仍缺乏对基因表达量预测的模型效果评估。
在本研究中，我们以猪群体数据为对象，比较了多种机器学习回归模型(包括线性回归模型和非线性回归模型)在基因表达量预测上的效果，并且探究了基因cis-SNP数量和基因顺式遗传力(Cis-heritability, cis-h2)大小对预测准确性的影响，为后续TWAS模型选择和模型优化提供参考。
1. 材料与方法
1.1 试验数据
本研究选取的
1321 头猪的肌肉组织基因型数据和基因表达量数据均来自PigGTEx项目(https://piggtex.farmgtex.org)。基因型数据只保留小等位基因频率(Minor allele frequency, MAF)> 0.05的SNP位点，质控后保留了2928814 个SNP。猪肌肉组织基因表达数据中共包含13454 个基因的TPM(Transcripts per million)表达量矩阵。使用R包edgeR[12]中的TMM(Trimmed mean of M-value)方法对TPM标准化，并进行逆正态转换(Inverse normal transformation, INT)。后续分析均使用INT的TMM基因表达矩阵。为了校正数据中的群体分层，我们使用R包SNPRelate(v 1.26.0)[13]计算基因型主成分，选取前10个基因型主成分作为协变量。为了消除基因表达数据的潜在混杂效应，本研究首先使用PEER(Probabilistic estimation of expression residual)方法[14]估计基因表达量数据中潜在的混杂因素，最终选择前10个PEER因子用于作为协变量。本研究利用一般线性回归模型校正表达量数据中的群体分层和潜在混杂因素，模型为：
$$ \mathit{y}={\mathit{X}}_{\mathit{c}\mathit{o}\mathit{v}}\mathit{\beta }+\mathit{\varepsilon } \text{，} $$ (1) 式中，
$ \mathit{y} $ 是基因表达量向量；$ {\mathit{X}}_{\mathit{c}\mathit{o}\mathit{v}} $ 是协变量矩阵；$ \mathit{\beta } $ 是回归系数向量；$ \mathit{\varepsilon } $ 是误差项向量。我们使用上述计算得到的10个基因型主成分和10个PEER因子构建协变量矩阵，对于第$ i $ 个基因，使用最小二乘法估计回归系数$ {\hat {\mathit{\beta }}}_{\mathit{i}} $ ，再计算校正后的基因表达量($ {\mathit{y}}_{\mathit{c}} $ )：$$ {\mathit{y}}_{\mathit{c}}=\mathit{y}-{\mathit{X}}_{\mathit{c}\mathit{o}\mathit{v}}\hat {\mathit{\beta }} 。 $$ (2) 本研究依据Ensembl数据库中的猪基因注释文件(Sscrofa 11.1 v100)，选取猪1号染色体上的
1091 个蛋白编码基因用于后续分析。1.2 基因表达量的机器学习预测模型
本研究所使用的机器学习回归模型均来自Python(v 3.12.3)中的scikit-learn包，一共18种监督学习模型，可以分为最小二乘法回归模型、正则化回归模型、岭回归模型、支持向量机回归模型、广义线性模型、K最近邻回归模型和树回归模型七类。
1.2.1 最小二乘法回归模型
包括普通最小二乘法回归(Ordinary least squares, OLS)和偏最小二乘法回归(Partial least squares, PLS)，是线性回归模型，其基本模型如下：
$$ \mathit{y}=\mathit{X}\mathit{\beta }+\mathit{\varepsilon } \text{，} $$ (3) 式中，
$ \mathit{X} $ 是个体基因型矩阵；$ \mathit{\beta } $ 是cis-SNP效应值向量。假设$ {y}_{1},{y}_{2},\cdots ,{y}_{n} $ 为$ n $ 个样本在特定组织中的基因表达量，则效应值向量求解公式(优化函数)如下：$$ \hat {\mathit{\beta }}=arg\underset{\beta }{\mathrm{min}}\left(1/2\right)\displaystyle\sum _{i=1}^{n}{\left({y}_{i}-{{x}_{i}}^{T}\mathit{\beta }\right)}^{2} 。 $$ (4) PLS同时考虑因变量和自变量之间的协方差求解回归系数，虽然计算更复杂，但在处理高维和共线性数据时能够提供比OLS更稳健的模型[15]。
1.2.2 正则化回归模型
正则化回归模型包括岭回归(Ridge)、Lasso回归、LassoLars回归和弹性网络回归(Elastic Net)。这些模型通过在OLS基础上引入罚函数(即正则化项)来防止过拟合，提高模型的泛化能力，是线性模型，其优化函数如下：
$$ \hat {\mathit{\beta }}=arg\underset{\beta }{\mathrm{min}}\left(1/2\right)\displaystyle\sum _{i=1}^{n}{\left({y}_{i} - {{x}_{i}}^{T}\mathit{\beta }\right)}^{2} + \lambda \left(\left(1 - \alpha \right){\|\mathit{\beta }\|}_{2}^{2} + \alpha {\|\mathit{\beta }\|}_{1}\right) \text{，} $$ (5) 式中，
$ {\|\mathit{\beta }\|}_{1} $ 和$ {\|\mathit{\beta }\|}_{2} $ 分别是$ \mathit{\beta } $ 的L1范数(即$ \mathit{\beta } $ 系数的绝对值和)和L2范数(即$ \mathit{\beta } $ 系数的平方和)，$ \lambda $ 是正则化强度参数。Ridge回归只包含L2范数[16]，即上式中$ \alpha =0 $ 。Lasso回归和LassoLars回归只包含L1范数，即上式中$ \alpha =1 $ ，两者优化函数相同，由于Lars回归是一种逐步回归方法，每一步选择一个变量进入模型，LassoLars回归会在每一步中使用L1正则化调整回归系数，与Lasso回归相比能更高效生成稀疏解，适用于高维数据[17]。Elastic Net回归包含L1和L2范数，可以调整L1和L2正则化的权重比例系数$ \alpha $ 使模型同时具备稀疏性和稳定性[18]，在本研究中取$ \alpha =0.5 $ ，即L1和L2正则化权重相同。Lasso回归、LassoLars回归和Elastic Net回归受$ \lambda $ 影响大，且存在收敛问题，即无法从已有数据中预测出结果。经过前期验证，针对本研究的数据，我们在使用这3种回归模型时，先对每个基因表达量预测进行5折交叉验证(将数据集分为5个子集，循环使用每个子集作为验证集，其余子集作为训练集)，从$ \lambda =0.1 $ 、$ \lambda =0.05 $ 、$ \lambda =0.01 $ 、$ \lambda =0.005 $ 、$ \lambda =0.001 $ 中选择每个基因预测效果最好的$ \lambda $ ，用于后续分析。1.2.3 岭回归模型
岭回归模型包括核岭回归模型(Kernel Ridge)和贝叶斯岭回归(Bayesian Ridge)。Kernel Ridge在OLS基础上引入了岭正则化和核函数，其优化函数如下：
$$ \hat {\mathit{\beta }}=arg\underset{\beta }{\mathrm{min}}\left(1/2\right)\displaystyle\sum _{i=1}^{n}{\left({y}_{i}-\mathit{K}\mathit{\beta }\right)}^{2}+\lambda {\|\mathit{\beta }\|}^{2} \text{，} $$ (6) 式中，
$ \mathit{K} $ 是$ n\times n $ 的核矩阵。最开始核函数在支持向量机(Support vector machine, SVM)中使用，可以大大降低其非线性变换的计算量，之后进一步将核函数引入到传统的机器学习算法中，可以方便地把线性算法转换为非线性算法[19]。其中核函数有4种：线性核(Linear kernel)，形式为$ K\left({x}_{i},{x}_{j}\right)={x}_{i}\cdot {x}_{j} $ ，适用于线性数据或特征数远大于样本数的数据；多项式核(Polynomial kernel，后文简称为poly)，形式为$ K\left({x}_{i},{x}_{j}\right)={\left({x}_{i}\cdot {x}_{j}+r\right)}^{d} $ ，其中$ r $ 是常数，$ d $ 是多项式的度数，适用于存在非线性关系的数据；径向基函数核 (RBF kernel)，形式为$ K\left({x}_{i},{x}_{j}\right)=exp\left(-\dfrac{{\|{x}_{i}-{x}_{j}\|}^{2}}{2{\sigma }^{2}}\right) $ ，其中$ \sigma $ 是核函数参数，同样适用于存在非线性关系的数据；Sigmoid核 (Sigmoid kernel)，形式为$ K\left({x}_{i},{x}_{j}\right)=\mathrm{tanh}\left(\alpha \cdot \left({x}_{i}\cdot {x}_{j}\right)+c\right) $ ，其中$ \alpha $ 和$ c $ 是可调参数，通常用于神经网络的训练。Bayesian Ridge与传统岭回归相比通过考虑后验分布来量化模型系数的不确定性。1.2.4 支持向量机回归模型
支持向量机回归模型(Support vector regression, SVR)的目标是找到一个函数
$ f\left(x\right)=\left\langle{\mathit{\beta },x}\right\rangle+b $ ，使对于所有$ \left({x}_{i},{y}_{i}\right) $ ，误差$ \left|{y}_{i}-f\left({x}_{i}\right)\right| $ 大部分位于给定的$ \varepsilon $ 内，目标优化函数如下：$$ \hat {\mathit{\beta }}=arg\underset{\beta ,b}{\mathrm{min}}\left(1/2\right){\|\mathit{\beta }\|}^{2}+loss \text{，} $$ (7) 式中，
$ loss $ 为损失函数。与核岭回归类似，SVR也可引入核函数，在本研究中，经过前期验证，我们在SVR中只使用linear、poly和RBF这3种核函数构建预测模型。SVR适用于具有复杂分布的数据或者高维数据。1.2.5 广义线性模型
广义线性模型(Generalize linear model, GLM) 主要由3个部分组成：随机部分
$ {Y}_{1},{Y}_{2},\cdots ,{Y}_{n} $ 服从的分布来自指数分布族；系统部分由线性预测器组成；存在一个单调可微的连接函数，用于连接随机部分的期望和系统部分的线性预测值，即：$$ g\left(\mu \right)=\eta \text{，} $$ (8) 式中，
$ \mu $ 是$ \mathit{Y} $ 的期望，$ \eta $ 是线性预测器。GLM存在多种连接函数，受$ Y $ 的分布决定。在本研究中，我们假定$ Y $ 的分布属于正态分布，则连接函数为$ \mu =\eta $ 。GLM适用于处理误差向量存在非恒等方差的情况，也适用于存在非线性关系的数据。1.2.6 K最近邻回归模型
K最近邻回归(K-nearest neighbors regression, KNN)是一种简单且直观的非参数回归方法，不需要对数据进行假设，直接利用数据中的实例进行预测。对于一个待预测基因表达量的样本，计算与训练集中每个样本间的欧式距离，找出距离最小的K个样本，计算基因表达量均值用于预测待预测样本的基因表达量。KNN无需训练，受异常值影响较小，且适用于存在非线性关系的数据，但不适用于处理高维数据。
1.2.7 树回归模型
树回归包括决策树回归(Decision tree regression)和随机森林回归(Random forest regression)，是非参数回归方法。决策树递归地将基因型矩阵
$ X $ 分割成更小的子集来构建树结构，通过最小化基因表达量$ y $ 均方误差(Mean squared error, MSE)选择分割点$ s $ ，对于第$ j $ 个基因型和选择的分割点$ s $ ，可以将基因型数据分成如下两个部分：$$ {K}_{1}\left(j,s\right)=\left\{\left(X,y\right)|{X}_{j}\in s\right\} $$ $$ {K}_{2}\left(j,s\right)=\left\{\left(X,y\right)|{X}_{j}\notin s\right\} $$ (9) 式中，
$ K $ 是叶节点包含的样本集合。对每个子集递归使用类似的分割过程，直至达到数的最大深度。对于待预测基因表达量的个体基因型向量$ {x}_{i} $ ，沿决策树从根节点基于分割点决策，最终到达一个叶节点，计算基因表达量均值用于预测待预测样本的基因表达量，即：$$ \hat {y}=\dfrac{1}{\left|K\right|}\displaystyle\sum _{i\in K}{y}_{i} 。 $$ (10) 随机森林回归通过构建多棵决策树进行回归。从基因型数据中有放回地抽取多个子集，每个子集用于构建一棵决策树，对于每棵决策树，递归随机选择一个特征子集用于分割，直至达到树的最大深度。随机森林模型虽然结构复杂，但所需的假设条件少，大多数情况下模型参数的缺省设置可以给出接近最优的结果[20]。随机森林回归集成多个决策树模型，提高了模型的稳定性和预测精度，能有效减少单棵决策树所导致的过拟合问题。这两种模型均能够捕捉基因型特征与基因表达量之间的复杂非线性关系。
1.3 模型的预测准确性评估
对于每个基因，选取转录起始位点(Transcription start site, TSS)上下游1 Mb范围内的cis-SNP用于训练，在
1321 个样本中进行10折交叉验证并重复5次。在一次10折交叉验证中，将1321 个样本平均分成10份，循环选取其中的1份(132个样本)作为验证集，剩余9份(1189 个样本)作为训练集。我们通过计算每个验证集中真实基因表达量和模型预测基因表达量之间的皮尔逊相关系数的平方(R2)来衡量模型预测基因表达量的准确性。1.4 基因顺式遗传力计算
为探究基因的遗传力与模型预测准确性的关系，我们选取每个基因TSS上下游1 Mb范围内的cis-SNP，建立混合线性模型估算基因的cis-h2，统计模型为：
$$ \mathit{y}={\mathit{X}}_{\mathit{c}\mathit{o}\mathit{v}}\mathit{\beta }+\mathit{Z}\mathit{g}+\mathit{\varepsilon } \text{，} $$ (11) 式中，
$ \mathit{Z} $ 是基因型矩阵；$ \mathit{g}\sim N\left(0,{\sigma }_{g}^{2}\mathit{G}\right) $ ,是与cis-SNP相关的随机效应，服从以基因组遗传关系矩阵(Genetic relationship matrix, GRM)为协方差的正态分布。该分析基于PLINK(v 1.90)提取每个基因的cis-SNP，利用GCTA(v 1.94.1)软件[21]中计算GRM矩阵，将其作为其随机效应，并将基因表达量作为表型，使用AIREML(Average information restricted maximum likelihood)估算基因的cis-h2。1.5 基因分组与筛选
为了研究基因cis-h2与模型预测准确性的关系，我们将基因按照cis-h2值大小分成cis-h2≤ 0.01、0.01 <cis-h2≤ 0.1、0.1 <cis-h2≤ 0.2、cis-h2> 0.2四组，计算每组基因中不同模型的预测准确性R2平均值。为了研究基因的cis-SNP数量对模型预测准确性的影响，我们计算所有基因cis-SNP数量的上下四分位数和中位数，分别选取cis-SNP数量在[0,下四分位数]、(下四分位数,中位数]、(中位数,上四分位数]、大于上四分位数这四个范围内的基因，计算它们在不同模型预测准确性R2与cis-h2的皮尔逊相关系数(r)。
2. 结果与分析
2.1 不同机器学习模型预测基因表达量总体准确性
不同机器学习模型预测的总体预测效果如图1所示。ElasticNet和Lasso两个模型的准确性R2平均值最高，分别为
0.0362 和0.0358 ；使用RBF非线性核的SVR和KernelRidge的R2平均值其次，分别为0.0354 和0.0347 ；GLM_Gaussian和OLS模型的R2平均值最低，分别为0.0089 和0.0072 。所有模型中98%以上的基因R2小于0.2，2%以内的基因R2大于0.2。![]()
图 1 不同机器学习模型的总体预测准确性(R2)柱子中的红色圆点表示R2平均值，右上方数字表示平均值具体数值。Figure 1. Overall prediction accuracy (R2) of different machine learning modelsRed points in the columns indicate means of R2, the numbers in the upper right corner indicate the values of means.2.2 不同模型预测准确性与基因顺式遗传力关系
18种机器学习模型预测准确性与基因cis-h2的关系如图2所示。在所有模型中，机器学习模型的预测准确性与基因cis-h2间存在正相关，两者间的皮尔逊相关系数大于0.60；对于整体预测准确性R2平均值较高的模型，皮尔逊相关系数大于0.75。
![]()
图 2 不同机器学习模型预测表达量准确性(R2)与基因顺式遗传力关系Figure 2. Relationship between the accuracy (R2) of predicting expression levels using different machine learning models and the cis-heritability of genes我们计算出所有基因的cis-h2中位数为
0.0267 ，平均值为0.0762 。基因按照cis-h2值大小分组情况如表1所示。其中，cis-h2≤ 0.01的基因占全部基因的39.54%，cis-h2> 0.2的基因占全部基因的12.88%。表 1 四组基因的顺式遗传力数据Table 1. Cis-heritability data of four groups of genes顺式遗传力
cis-h2基因个数
Number of genes中位数
Median value平均值
Mean value标准误
Standard error≤0.01 433 1.00e-6 1.00e-3 1.02e-4 (0.01,0.1] 372 0.0433 0.0475 1.35e-3 (0.1,0.2] 139 0.139 0.144 2.49e-3 > 0.2 141 0.281 0.315 9.21e-3 我们选取每组基因预测准确性R2平均值排名前7的机器学习模型，结果如图3所示。对于cis-h2≤ 0.01的基因，Lasso和RandomForest的预测准确性最高，R2平均值分别为
0.0152 和0.0146 ；对于cis-h2范围在(0.01, 0.1]的基因，ElasticNet和Lasso这两个正则化回归模型的预测准确性最高，R2平均值分别为0.0275 和0.0272 ；对于cis-h2范围在(0.1, 0.2]的基因，使用RBF非线性核的KernelRidge和SVR的预测准确性最高，R2平均值分别为0.0530 和0.0526 ；对于cis-h2> 0.2的基因，ElasticNet和Lasso的预测准确性最高，R2平均值分别为0.115和0.107。![]()
图 3 四组不同顺式遗传力基因排名前7的模型的预测准确性(R2)Figure 3. The prediction accuracy (R2) of the top 7 models with four different cis-heritability gene rankings2.3 顺式SNP数量对不同模型预测准确性的影响
不同机器学习回归模型由于算法存在差异，不同的特征数量(即cis-SNP数量)也会对预测准确性有影响，一般情况模型预测准确性与cis-SNP数量间存在正相关。我们计算出cis-SNP数量与cis-h2的皮尔逊相关系数为0.236(P = 3.05×10−15)。基因cis-SNP数量与cis-h2存在相关性，同样表明cis-SNP数量对模型预测准确性有一定影响。
我们计算得到基因cis-SNP数量的下四分位数是
1564 ，中位数是3097 ，上四分位数是5384 。我们选取cis-SNP数量在[0,1564 ]、(1564 ,3097 ]、(3097 -5384 ]、>5384 范围的基因，计算它们在不同模型预测准确性与cis-h2的皮尔逊相关系数，得到结果如表2所示。除了OLS、Ridge、KernelRidge_sigmoid、GLM_Gaussian四个模型，其余14种模型在基因cis-SNP数量在(3097 -5384 ]范围时预测准确性最高，说明cis-SNP数量过少或过多会使这些模型训练过程产生欠拟合和过拟合现象。OLS和GLM_Gaussian随基因cis-SNP数量增多，预测准确性提升，与基因cis-SNP数量在[0,1564 ]范围内相比，基因cis-SNP数量大于5384 时两个模型预测准确性分别提升8.78%和10.87%。表 2 不同顺式SNP数量的基因在不同模型预测准确性与顺式遗传力的相关系数(r)Table 2. 2The correlation coefficient (r) between the prediction accuracy and cis-heritability of genes with different numbers of cis-SNPs in different models模型 0≤顺式SNP个数≤ 1564 1564 <顺式SNP个数≤3097 3097 <顺式SNP个数≤5384 顺式SNP个数> 5384 OLS 0.638 0.647 0.673 0.694 PLS 0.656 0.668 0.698 0.647 Ridge 0.649 0.621 0.639 0.626 Lasso 0.782 0.78 0.795 0.774 LassoLars 0.796 0.792 0.801 0.767 ElasticNet 0.806 0.802 0.811 0.792 KernelRidge_linear 0.651 0.613 0.654 0.625 KernelRidge_poly 0.787 0.803 0.814 0.808 KernelRidge_RBF 0.794 0.805 0.81 0.768 KernelRidge_sigmoid 0.657 0.65 0.653 0.615 BayesianRidge 0.777 0.778 0.786 0.75 SVR_linear 0.617 0.617 0.634 0.632 SVR_poly 0.742 0.77 0.788 0.774 SVR_RBF 0.795 0.801 0.814 0.772 GLM_Gaussian 0.64 0.671 0.704 0.709 KNN 0.71 0.723 0.754 0.732 DecisionTree 0.618 0.644 0.68 0.631 RandomForest 0.756 0.766 0.783 0.773 3. 讨论与结论
个体的基因表达水平可以被分解为由遗传决定的成分(Genetically regulated expression，GReX)、一个由性状本身改变的成分和归因于环境和其他因素的剩余成分。对于大多数基因，其表达量被认为具有稀疏结构，即由几个主要SNP调控，每个SNP都对基因表达量有较大影响[22]，而剩余的SNP对基因表达量都有微小的影响。利用SNP预测基因表达量是一个多元回归任务，包含多个SNP作为特征，预测准确性会直接影响TWAS后续进行基因表达量与表型之间的关联分析。经过前期研究发现，若使用所有SNP而非使用基因的cis-SNP进行基因表达量预测会产生较多噪音同时运行速度缓慢[10]，所以本研究使用了18种不同机器学习回归模型利用cis-SNP预测猪基因表达量的预测准确性。
结果表明，不同模型预测准确性(R2)与基因cis-h2间存在正相关，但R2一般不高于cis-h2。使用机器学习回归模型预测的基因表达量实际上是预测基因表达水平的GReX，无法预测环境及其他因素对基因表达水平的影响。因此，基因的遗传力越大，即基因表达受遗传因素影响较大，模型的预测准确性越高。在本研究中所有基因cis-h2的平均值为
0.0762 ，所有模型整体预测准确性(R2)的平均值范围为0.0072 -0.0362 ，符合预期结果。整体预测准确性最高的模型是Elastic Net回归模型和Lasso回归模型，这与Gamazon等人[10]的研究结果一致。Elastic Net回归模型结合了Lasso回归的L1正则化项和Ridge回归的L2正则化项，通过灵活调整L1和L2正则化的权重比例系数
$ \alpha $ ，可以防止模型过拟合，且有利于生成稀疏模型，相比Lasso回归对于输入 SNP之间的微小差异具有鲁棒性，且也在许多研究中被验证在预测人的某些复杂表型表现较好[23, 24]。Elastic Net回归模型、Lasso回归模型、使用非线性核poly和RBF的支持向量机回归模型和核岭回归模型、LassoLars回归模型、随机森林回归模型这八个机器学习模型对所有cis-h2区间的基因预测准确性较高。一些非线性模型也被验证在某些时候具有比线性模型高的基因表达预测准确性，如K最近邻回归、随机森林回归[25]，但在本研究中未发现这些非线性模型总体上具有高于Elastic Net回归的预测能力，但是在个别基因上，它们具有较好的预测能力，且在回归模型中使用非线性核RBF和poly相比于线性核能提高模型的预测能力。可能原因是部分基因表达量的SNP效应存在一定的非线性关系。对于cis-h2极低(cis-h2< 0.01)的基因，Lasso回归模型和随机森林回归模型预测能力较强；对于cis-h2中高(cis-h2> 0.2)的基因，Elastic Net回归模型和Lasso回归模型预测能力较强。尽管如此，预测效果最好的Elastic Net回归模型R2均值小于0.1，预测准确性较低。Fryett等[26]在多个人类数据集上同样使用多种机器学习模型预测基因表达量，将皮尔逊相关系数(R)用于评估预测准确性，发现预测效果最好的BSLMM的R均值为
0.0743 ，预测准确性同样较低。我们认为预测准确性低的原因主要有：虽然我们将基因表达量数据进行校正，但仍包含环境因素等非遗传成分，直接将其与预测得到的GReX计算相关性会造成偏差；同时猪群体具有复杂的群体和遗传结构，我们使用的机器学习模型并未考虑个体间的亲缘关系。Wainberg等[27]发现预测的基因表达量与实际表达量相关性较低的情况下，TWAS分析的z分数仍然较高，表明预测的表达量在检测关联性上比真实表达量具有优势，能够捕捉到遗传变异、基因表达和性状之间的调控关系。这表明即使预测的基因表达量准确性较低，在TWAS分析中具有实质意义。在未来，我们考虑构建混合线性模型，将个体亲缘关系矩阵纳入随机效应，以提升预测准确性。在回归问题中，不同回归模型所适用的最佳特征数会因其内部算法不同而不同，在本研究中，我们选取不同cis-SNP数量的基因进行分析，发现一般情况下，cis-SNP数量在一定范围内，随cis-SNP数量即特征数量增加，回归模型拟合效果也会越好，预测准确性提升，但是SNP数量过少或过多，也会使大部分模型在训练时产生欠拟合和过拟合现象，这与Fan等人[28]的研究结果相似。所有模型中，线性模型预测准确性受cis-SNP数量影响大于非线性模型。在Guyon等[29]的研究中提到在数据是高维的情况下，线性模型往往更依赖于特征选择，意味着线性模型的预测准确性会更大程度上受到特征数量的影响，故结果符合预期。
综上所述，相同条件下不同机器学习模型对基因表达量的预测能力存在差异，使用同种机器学习模型的预测准确性受基因cis-h2和cis-SNP数量影响较大。所以根据基因cis-h2和cis-SNP数量选择合适的机器学习模型预测基因表达量，并在未来根据已有数据结构开发新的回归预测模型，能够有效提高预测准确性，进而提高TWAS分析的统计效力。
-
![]()
图 1 不同机器学习模型的总体预测准确性(R2)
柱子中的红色圆点表示R2平均值，右上方数字表示平均值具体数值。
Figure 1. Overall prediction accuracy (R2) of different machine learning models
Red points in the columns indicate means of R2, the numbers in the upper right corner indicate the values of means.
![]()
图 2 不同机器学习模型预测表达量准确性(R2)与基因顺式遗传力关系
Figure 2. Relationship between the accuracy (R2) of predicting expression levels using different machine learning models and the cis-heritability of genes
![]()
图 3 四组不同顺式遗传力基因排名前7的模型的预测准确性(R2)
Figure 3. The prediction accuracy (R2) of the top 7 models with four different cis-heritability gene rankings
表 1 四组基因的顺式遗传力数据
Table 1 Cis-heritability data of four groups of genes
顺式遗传力
cis-h2基因个数
Number of genes中位数
Median value平均值
Mean value标准误
Standard error≤0.01 433 1.00e-6 1.00e-3 1.02e-4 (0.01,0.1] 372 0.0433 0.0475 1.35e-3 (0.1,0.2] 139 0.139 0.144 2.49e-3 > 0.2 141 0.281 0.315 9.21e-3 
下载: 导出CSV
表 2 不同顺式SNP数量的基因在不同模型预测准确性与顺式遗传力的相关系数(r)
Table 2 2The correlation coefficient (r) between the prediction accuracy and cis-heritability of genes with different numbers of cis-SNPs in different models
模型 0≤顺式SNP个数≤ 1564 1564 <顺式SNP个数≤3097 3097 <顺式SNP个数≤5384 顺式SNP个数> 5384 OLS 0.638 0.647 0.673 0.694 PLS 0.656 0.668 0.698 0.647 Ridge 0.649 0.621 0.639 0.626 Lasso 0.782 0.78 0.795 0.774 LassoLars 0.796 0.792 0.801 0.767 ElasticNet 0.806 0.802 0.811 0.792 KernelRidge_linear 0.651 0.613 0.654 0.625 KernelRidge_poly 0.787 0.803 0.814 0.808 KernelRidge_RBF 0.794 0.805 0.81 0.768 KernelRidge_sigmoid 0.657 0.65 0.653 0.615 BayesianRidge 0.777 0.778 0.786 0.75 SVR_linear 0.617 0.617 0.634 0.632 SVR_poly 0.742 0.77 0.788 0.774 SVR_RBF 0.795 0.801 0.814 0.772 GLM_Gaussian 0.64 0.671 0.704 0.709 KNN 0.71 0.723 0.754 0.732 DecisionTree 0.618 0.644 0.68 0.631 RandomForest 0.756 0.766 0.783 0.773
下载: 导出CSV
-
[1] 牛安然, 张兴, 杨雨婷, 等. 全基因组关联分析在猪育种中的研究进展[J]. 畜牧与兽医, 2023(55): 139-147. [2] LI T, WAN P, LIN Q, et al. Genome-wide association study meta-analysis elucidates genetic structure and identifies candidate genes of teat number traits in pigs[J]. International Journal of Molecular Sciences, 2023(25): 451.
[3] ZENG H, ZHONG Z, XU Z, et al. Meta-analysis of genome-wide association studies uncovers shared candidate genes across breeds for pig fatness trait[J]. BMC Genomics, 2022(23): 1-786.
[4] 窦腾飞, 吴姿仪, 白利瑶, 等. 全基因组关联分析鉴定大白猪生长性状遗传变异及候选基因[J]. 中国畜牧杂志, 2023(59): 264-272. [5] LI X, WU J, ZHUANG Z, et al. Integrated single-trait and multi-trait GWASs reveal the genetic architecture of internal organ weight in pigs[J]. Animals, 2023(13).
[6] 张宇, 周佳伟, 吴俊静, 等. 大白猪繁殖性状全基因组关联分析[J]. 中国畜牧杂志, 2022(58): 94-99. [7] TENG J, GAO Y, YIN H, et al. A compendium of genetic regulatory effects across pig tissues[J]. Nature Genetics, 2024(56): 112-123.
[8] 郑韵頔, 冉雪琴, 牛熙, 等. 全基因组eQTL揭示猪11号染色体肉质性状新候选位点[J]. 农业生物技术学报, 2024(32): 807-819. doi: 10.3969/j.issn.1674-7968.2024.04.007 [9] MAI J, LU M, GAO Q, et al. Transcriptome-wide association studies: recent advances in methods, applications and available databases[J]. Communications Biology, 2023(6): 899.
[10] GAMAZON E R, WHEELER H E, SHAH K P, et al. A gene-based association method for mapping traits using reference transcriptome data[J]. Nature Genetics, 2015(47): 1091-1098.
[11] GUSEV A, KO A, SHI H, et al. Integrative approaches for large-scale transcriptome-wide association studies[J]. Nature Genetics, 2016(48): 245-252.
[12] ROBINSON M D, OSHLACK A. A scaling normalization method for differential expression analysis of RNA-seq data[J]. Genome Biology, 2010(11): R25.
[13] ZHENG X, LEVINE D, SHEN J, et al. A high-performance computing toolset for relatedness and principal component analysis of SNP data[J]. Bioinformatics, 2012(28): 3326-3328.
[14] STEGLE O, PARTS L, PIIPARI M, et al. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses[J]. Nature Protocols, 2012(7): 500-507.
[15] MEHMOOD T, LILAND K H, SNIPEN L, et al. A review of variable selection methods in Partial Least Squares Regression[J]. Chemometrics and Intelligent Laboratory Systems, 2012(118): 62-69.
[16] HOERL A E, KENNARD R W. Ridge regression: biased estimation for nonorthogonal problems[J]. Technometrics, 1970(12): 55-67.
[17] TIBSHIRANI R. Regression shrinkage and selection via the lasso[J]. Journal of the Royal Statistical Society Series B: Statistical Methodology, 1996(58): 267-288.
[18] ZOU H, HASTIE T. Regularization and variable selection via the elastic net[J]. Journal of the Royal Statistical Society Series B: Statistical Methodology, 2005(67): 301-320.
[19] 汪廷华, 陈峻婷. 核函数的选择研究综述[J]. 计算机工程与设计, 2012(33): 1181-1186. doi: 10.3969/j.issn.1000-7024.2012.03.068 [20] 李欣海. 随机森林模型在分类与回归分析中的应用[J]. 应用昆虫学报, 2013(50): 1190-1197. doi: 10.7679/j.issn.2095-1353.2013.163 [21] YANG J, LEE S H, GODDARD M E, et al. GCTA: a tool for genome-wide complex trait analysis[J]. The American Journal of Human Genetics, 2011(88): 76-82.
[22] WHEELER H E, SHAH K P, BRENNER J, et al. Survey of the heritability and sparse architecture of gene expression traits across human tissues[J]. PLoS Genetics, 2016(12): e1006423.
[23] BAE S, CHOI S, KIM S M, et al. Prediction of quantitative traits using common genetic variants: application to body mass index[J]. Genomics Inform, 2016(14): 149-159.
[24] SPILIOPOULOU A, NAGY R, BERMINGHAM M L, et al. Genomic prediction of complex human traits: relatedness, trait architecture and predictive meta-models[J]. Human Molecular Genetics, 2015(24): 4167-4182.
[25] WANG J, GAMAZON E R, PIERCE B L, et al. Imputing gene expression in uncollected tissues within and beyond GTEx[J]. The American Journal of Human Genetics, 2016(98): 697-708.
[26] FRYETT J J, MORRIS A P, CORDELL H J. Investigation of prediction accuracy and the impact of sample size, ancestry, and tissue in transcriptome-wide association studies[J]. Genetic Epidemiology, 2020(44): 425-441.
[27] WAINBERG M, SINNOTT-ARMSTRONG N, MANCUSO N, et al. Opportunities and challenges for transcriptome-wide association studies[J]. Nature Genetics, 2019(51): 592-599.
[28] FAN J, LV J. A selective overview of variable selection in high dimensional feature space[J]. Statistica Sinica, 2010(20): 101.
[29] GUYON I, ELISSEEFF A. An introduction to variable and feature selection[J]. Journal of Machine Learning Research, 2003(3): 1157-1182.
计量
- 文章访问数: 118
- HTML全文浏览量: 9
- PDF下载量: 29
