基于气相色谱-质谱的尿液代谢组学技术结合化学计量学用于戊二酸血症Ⅰ型早期检测研究

肖 雯，牛芊芊，孙智勇，杨 琴*，吴本清

（1．长江大学 物理与光电工程学院，湖北 荆州 434023；2．深圳爱湾医学检验实验室 深圳罕见病代谢组学精准医学工程研究中心，广东 深圳 518000；3．中国科学院大学 深圳医院，广东 深圳 518000）

戊二酸血症Ⅰ型（GA-Ⅰ）作为一种常染色体隐性遗传病，因编码基因发生突变导致对应的戊二酰辅酶A脱氢酶（GCDH）活性降低或缺失，有毒代谢产物（如戊二酰肉碱、戊二酸、3-羟基戊二酸、戊烯二酸等）在体液与组织中异常蓄积，形成酸中毒而造成严重的不可逆神经系统受损［1-2］。患者临床表现复杂，易与其他神经系统疾病混淆，缺乏特异性［3］。因此，建立高效的GA-Ⅰ早期筛查模型，对缓解新生儿遗传代谢病公共健康难题具有积极的促进作用。随着新生儿筛查项目的推广和普及，基于气相色谱-质谱（GC-MS）的尿液代谢组学技术在GA-Ⅰ早期检测中展示出独特优势，能全面揭示代谢产物的浓度差异［4-6］。化学计量学中各种建模算法的挖掘能力可将数据内在信息有效转化为可理解的知识，已广泛应用于疾病的早期检测和辅助诊断［7-9］。然而，质谱技术的不断发展产生高维小样本的疾病代谢组学图谱，使得应用单个数据分析学习模型进行特征筛选不够稳健，对样本变异（增加或删除一些样本）敏感，不仅影响疾病早期检测模型的建模性能，也会降低各领域专家对分析结果的可信度［10-11］。

自助抽样法（Bootstrap）是集成特征筛选中一种常用的重抽样技术，先通过有放回的抽样方式产生大量样本变异较小且与原始数据结构相同的扰动数据集，再使用相同的数据分析学习算法建立多个特征选择器［12-13］。该法能够在多个扰动数据集中挑选出持续被筛选的变量，其稳健性和解释样本组别间差异的优势已得到实验和理论上的证实［14］。偏最小二乘判别分析（PLS-DA）是处理高维数据中共线性和信息冗余问题的有力统计分析工具，其数据解释能力强，通过计算系列信息向量如载荷（LW）、变量投影重要性（VIP）、显著性多元相关（sMC）等直接刻画各变量对模型的贡献大小［15-16］。本文应用GC-MS技术采集新生儿尿液的代谢组学图谱信息，结合自助抽样法和PLS-DA实现真正解释组别间差异且与疾病机理密切相关的稳健特征变量筛选，建立了高性能的GA-Ⅰ早期检测模型，并进行分析和验证。

1 实验部分

1.1 样本材料

尿液样本采集于2015~2020年，由中国科学院大学深圳医院提供。疾病组由30例GA-Ⅰ患者（＜28天）组成，对照组由32例年龄匹配的健康新生儿组成。本研究经中国科学院大学深圳医院伦理委员会同意，尿液代谢组学图谱检测均获得参与者监护人的同意。尿液样本保存于-20℃，不添加任何防腐剂。

1.2 仪器与试剂

GCMS-QP2010型气相色谱-质谱联用仪（日本岛津公司）。乙酸乙酯（色谱纯）、甲基硅烷化衍生试剂（99∶1，上海安谱实验科技股份有限公司），十七烷酸、尿素酶、盐酸羟胺（美国Sigma-Aldrich公司），盐酸（优级纯，湖南凯信公司），苦味酸（分析纯，成都西亚化学工业有限公司），氢氧化钠（分析纯，美国Sigma-Aldrich公司）。以十七烷酸作为内标，用乙酸乙酯配制质量浓度为0．5 mg/mL的十七烷酸溶液。

1.3 GC-MS代谢组学图谱数据采集

尿液样本用蒸馏水稀释定容至2 mL，其中肌酐质量为0．2 mg。样本中加入20 μL尿素酶分解尿素，并与0．02 mg十七烷酸混合，然后经盐酸羟胺、氢氧化钠和盐酸处理，再加入3 mL乙酸乙酯对有机酸萃取2次。将离心（4 000 r/min，5 min）后的有机酸层移至干净离心管中，用氮气在60℃下吹干。加入100 μL甲基硅烷化衍生试剂（99∶1），在70℃衍生反应30 min后，取1 μL样本采用GC-MS测定尿液中的有机酸。测试条件为：进样口温度280℃，色谱柱温度100~280℃，以4℃/min逐步升温，质谱接口温度280℃，离子源温度200℃。电子束能量70 eV，扫描范围m/z50~500，扫描速率1 000 Da/s。采集的原始质谱数据利用配套的GCMSsolution软件和代谢物质谱数据库进行处理，鉴定出132种代谢物。对于同一个色谱峰，鉴别出的代谢物浓度通过其峰面积与内标物十七烷酸峰面积的比值确定。

2 GA-Ⅰ早期检测模型

GA-Ⅰ早期检测模型的技术原理如图1所示。首先，将62例样本随机划分为训练集和测试集。然后，对划分的训练集采用自助抽样法进行有放回的抽样，生成多个样本大小与原始训练集相同的扰动数据集。基于PLS-DA分别建模，计算LW、VIP和sMC的信息向量，并根据指标的绝对值大小对变量进行排序。考虑到各变量对模型的贡献，只取排序后的变量序列中前10%的变量，作为真正表征组别间代谢差异的特征变量。最后，集成所有扰动数据集对应的排序变量序列，统计各变量跨越模型间的筛选频率，并根据筛选频率再次对变量进行排序，设置频率筛选阈值（如0．1，0．15，0．2，…，1）确定最终的特征变量序列［13］。频率筛选阈值设置得越高，最终特征变量序列的确定越缺乏偶然性，筛选越严格。对于测试集，根据确定的最终特征变量序列挑选出相应的变量组成数据集，利用PLS-DA建模，计算ŷ对样本类别进行预测。

图1 GA-Ⅰ早期检测模型流程图Fig．1 The flowchart of early detection model for GA-Ⅰ

62例样本分别按照7∶3和6∶4比例划分训练集和测试集。自助抽样法的抽样次数设为100次，每个扰动数据集对应的PLS-DA模型的隐变量个数利用10折交叉验证法确定。另外，为验证GA-Ⅰ早期检测模型的特征筛选稳健性，对62例样本随机划分50次，生成的50个训练集分别采用图1的流程筛选最终的特征变量序列，计算Kuncheva指数（KI）［17］：

其中fi表示第i个训练集筛选的最终特征变量序列；h=|fi|=|fj|，表示所有筛选的最终特征变量序列中包含元素个数的最小值；r=|fi∩fj|，表示fi和fj中共同元素的个数；N表示尿液中鉴定出的代谢物总数；修正项h2/N表示在fi和fj中选择一个共同特征其偶然性的概率。KI取值范围为-1~1，负值表明特征重叠在很大程度上是由于偶然性。KI指数越大，fi和fj的特征重叠程度越高。最后，总稳定性KItot取所有两两比较相似度值KI的平均：

为方便表述，各种方法的总稳定性用KI简化表示。对于划分50次生成的50组训练集和测试集，分别采用受试者工作特征曲线下面积（AUC）、正确率（ACC）、灵敏度（SEN）、特异性（SPE）以及马修斯相关系数（MCC）表征模型性能。

3 结果与讨论

3.1 GC-MS代谢组学图谱

图2展示了GA-Ⅰ和对照组样本的尿液代谢总离子流色谱图（TIC），在指定的保留时间范围内观察到丰富的代谢产物图谱，显示出GC-MS能够同时检测尿液中多种有机酸的能力。对鉴定的132种代谢物进行单变量统计分析，发现多种代谢物在GA-Ⅰ和对照组之间均出现显著的浓度变化（p＜0．01），揭示了GA-Ⅰ的发病机制导致复杂的组别间差异。此外，代谢物之间还存在共线性问题，如3-羟基-3-甲基戊二酸和3-羟基异戊酸的浓度相关系数为0．977 3，异枸橼酸和乌头酸的浓度相关系数为0．922 0。因此，采用多元统计分析方法PLS-DA结合自助抽样法处理变量间的高相关性，以提取GCMS尿液代谢组学图谱中更细微的代谢信息变化，选择能够很好解释疾病组与对照组间差异的稳健特征。

图2 GA-Ⅰ和对照组样本的尿液代谢总离子流色谱图Fig．2 Total ion chromatograms（TIC）of urine metabolic profiling for GA-Ⅰand controls

3.2 特征变量稳健筛选

基于LW、VIP和sMC信息向量，对每个训练集利用单个PLS-DA建立的模型分别为LW-PLSDA、VIP-PLSDA和sMC-PLSDA，结合自助抽样法生成多个扰动数据集建立的模型分别为BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA。当样本按照7∶3比例（训练集43例，测试集19例）划分时，对50个训练集分别建模后，公式（1）中的h=18，计算得到LW-PLSDA、VIP-PLSDA和sMCPLSDA的KI值均低于0．4（图3）。可见，数据集的高维小样本特点大大影响了基于单个PLS-DA建模的性能，导致LW-PLSDA、VIP-PLSDA和sMCPLSDA筛选的特征变量序列在各训练集之间差异过大，即使预测性能优异也大大降低了各领域专家对分析结果的可信度。引入自助抽样法，BS-LWPLSDA、BS-VIP-PLSDA和BS-sMC-PLSDA的KI值均显著增加，特别是BS-VIP-PLSDA在筛选特征变量个数为12时的KI值高达0．807 5，表明基于样本扰动的集成特征选择策略能够专注挑选被多个扰动数据集持续筛选的变量，有效提高特征选择的稳健性。进一步按照6∶4比例划分（训练集37例，测试集25例）增大50个训练集间的样本差异，以避免不同训练集始终筛选相同的变量。此时，公式（1）中的h=16，BS-PLSDA各模型仍然比其单独建模的PLS-DA展示出更优异的特征变量筛选稳健性（图4）。

图3 各模型筛选的特征变量序列在50个训练集之间的筛选稳健性比较Fig．3 Comparison of selection stability across 50 training sets among various techniques training set：43；test set：19

图4 各模型筛选的特征变量序列在50个训练集之间的筛选稳健性比较Fig．4 Comparison of selection stability across 50 training sets among various techniques training set：37；test set：25

3.3 稳健筛选特征变量的建模性能

为严格筛选特征变量，结合自助抽样法的3种模型BS-LW-PLSDA、BS-VIP-PLSDA和BSsMC-PLSDA对应的频率筛选阈值最高为0．3，即在100次重抽样中至少被30个扰动数据集建立的模型筛选。在频率筛选阈值设置为0．3时，3个模型由50个训练集确定的最终特征变量序列分别输入到对应的测试集中。同时，由132种代谢物组成的全变量序列也分别输入50个测试集中采用PLS-DA建模。当样本按照7∶3比例（训练集43例，测试集19例）划分时，建模结果如表1所示。全变量PLS-DA模型对50个测试集的AUC平均值为0．837 9，ACC平均值为0．821 5，MCC平均值为0．724 2。BSPLSDA 3个模型中，BS-VIP-PLSDA和BS-sMC-PLSDA对50个测试集的AUC平均值分别为0．854 8和0．847 1，ACC平均值分别为0．835 3和0．850 5，MCC平均值分别为0．783 8和0．801 3。结果表明相比全变量PLS-DA模型，BS-VIP-PLSDA和BS-sMC-PLSDA采用更少的信息变量可获得更好的建模性能，显示了特征筛选的重要性。进一步，分别统计3个模型确定的最终特征变量序列在50个训练集之间各变量的筛选频率，再次排序后的特征变量序列如表1所示。已有文献报道，在GA-Ⅰ确诊患儿的尿液中戊二酸（Glutaric acid）、3-羟基戊二酸（3-Hydroxyglutaric acid）、戊烯二酸（Glutaconic acid）的浓度增高，具有一定的诊断意义［3，18］。由表1可见，sMC信息向量的数据解释能力优于LW和VIP信息向量，上述3种具有一定诊断意义的代谢物均被筛选出，且在对应的热图分析（图5）中，GA-Ⅰ和对照组之间具有显著的浓度差异（p＜0．05）。另外，LW和VIP信息向量虽然均筛选出戊二酸和戊烯二酸，但对于LW戊烯二酸仅被50个训练集中的2个筛选，而对于VIP戊烯二酸被50个训练集中的13个筛选，体现了两者不同的特征变量搜索能力。需要注意的是，BS-LW-PLSDA的模型性能指标虽然低于全变量PLS-DA模型，但分类效果令人满意，最重要的是能够提供表征疾病组与对照组之间代谢信息差异的特征变量信息，显著提升模型的解释能力和临床应用价值。其他被筛选出的特征代谢物应进行实验验证，有助于挖掘更多与GA-Ⅰ代谢机理相关的知识。

图5 BS-sMC-PLSDA筛选的特征变量组别间差异热图分析Fig．5 Heatmap of the top-ranked significant metabolites selected by BS-sMC-PLSDA

表1 结合自助抽样法各种模型的筛选特征变量及建模结果Table 1 The top-ranked significant metabolites and classification parameters selected by various techniques combined bootstrap

当样本按照6∶4比例（训练集37例，测试集25例）划分时，BS-PLSDA 3个模型对应的频率筛选阈值最高也为0．3，建模结果如表2所示。虽然训练集样本减少，但经特征筛选后，BS-VIP-PLSDA和BS-sMC-PLSDA对50个测试集的AUC平均值分别为0．857 9和0．807 6，ACC平均值分别为0．838 4和0．813 6，MCC平均值分别为0．790 9和0．740 6，模型性能指标仍高于全变量PLS-DA模型。另外，相比于BS-VIP-PLSDA，BS-sMC-PLSDA的建模性能略低，但仍展示出最优的数据解释能力，能同时筛选出3种具有一定诊断意义的代谢物。对比两表中BS-VIP-PLSDA和BS-sMC-PLSDA分别确定的特征变量序列，发现对于不同样本划分比例，除了排序稍有不同，其组成大致相同。对于BSLW-PLSDA，相比7∶3比例样本划分，尽管在6∶4比例划分时筛选的特征变量个数降为2个，但代谢物分别为与GA-Ⅰ代谢机理密切相关的戊二酸和戊烯二酸。表明了基于自助抽样法的集成特征选择策略对提高变量筛选稳健性和建模性能的有效性。

表2 结合自助抽样法各种模型的筛选特征变量及建模结果Table 2 The top-ranked significant metabolites and classification parameters selected by various techniques combined bootstrap

为进一步展示本文GA-Ⅰ早期检测模型的可行性，在基于样本扰动的集成特征选择策略基础上，对每个扰动数据集采用支持向量机递归特征消除法（SVM-RFE）［19］进行建模。在特征变量筛选时，考虑保留其对模型贡献的信息，SVM首选线性核函数，对应模型分别为LIN-SVMRFE和BS-LINSVMRFE。另外，考虑SVM对非线性数据的优异分析能力，利用径向基核函数分别建模RBFSVMRFE和BS-RBF-SVMRFE。2种核函数对应的参数（C和gamma）利用10折交叉验证法确定。在递归循环中，SVM-RFE每次删除20%的变量个数。与BS-PLSDA 3个模型采用相同的自助抽样法配置，随着筛选变量个数的增加，BS-SVMRFE两个模型的KI值约为单独建模SVM-RFE模型KI值的2倍（图3和图4）。同样，频率筛选阈值设置为0．3，确定的最终特征变量序列分别输入到50组训练集和测试集中，BS-SVMRFE 2个模型的建模性能均优于BS-PLSDA 3个模型，特别是样本按照7∶3比例划分时BS-RBF-SVMRFE模型的性能指标平均值均超过0．900 0（表1和表2）。但在数据解释能力方面，线性核函数保留了变量对模型贡献的信息，BS-LIN-SVMRFE筛选出戊二酸和3-羟基戊二酸两个具有诊断意义的代谢物；而径向基核函数丢失了这部分信息，导致BS-RBF-SVMRFE只筛选出戊二酸代谢物，且其他特征变量与BS-LIN-SVMRFE和BS-PLSDA 3个模型也明显不同。表明GA-Ⅰ早期检测模型采用基于样本扰动的集成特征选择策略是合理的，对每个扰动数据集采用PLS-DA建模可同时兼顾建模性能和模型解释能力，符合实际临床需求；采用SVM-RFE建模虽然可获得较好的建模性能，但解释能力略低。

4 结论

为提高对高维小样本疾病数据的建模能力，基于GC-MS联用技术，结合自助抽样法和PLS-DA建立了GA-Ⅰ早期检测模型。通过对尿液的代谢组学图谱进行分析验证，相比单个模型建模，自助抽样法通过重复抽样的方式生成多个扰动数据集，结合PLS-DA的特征选择能力，建立的GA-Ⅰ早期检测模型专注于持续被筛选的变量，可有效提升特征筛选的稳健性。进一步，基于LW、VIP和sMC在特征空间的搜索能力，筛选出的稳健特征变量不仅分类效果令人满意，而且真正解释了组别间的差异与GA-Ⅰ的代谢机理密切相关，展示了丰富的生物学意义以及优异的数据解释能力。由此可见，本研究提出的模型在GA-Ⅰ的早期检测、辅助诊断以及疾病机理研究中具有一定的潜力。