哈尔滨医科大学卫生统计学教研室(150081) 李贞子 侯 艳 邓 魁 李 康

代谢组学研究方法主要是通过色谱和质谱等仪器获得组织、血液中的代谢物质，特别是尿液内源小分子代谢物的状态和含量[1]。由于质谱数据能够提供具体代谢物质化学结构的生物学信息，因此大部分的传统方法主要针对的是质谱数据。然而，大部分的传统方法在实际应用中尚存在各种不同的问题，其中最主要的问题是不易进行低浓度标志物识别[2]。质谱数据变量选取数目的多少与预处理设定的阈值有关，如果阈值设定较低，其质谱的变量个数可能达到10万以上，从而使分析不易进行。为了减少变量的个数通常将阈值设置在一个较高的水平，即对阈值以下的所有低浓度物质的含量置零。而根据生物学知识可知，肿瘤标志物更有意义往往是低浓度的化合物[3]。本文拟结合卵巢癌色谱和质谱的数据，使用小波变换的方法将原始色谱数据变换为二维图像，在此基础上使用随机森林(random forest,RF)模型筛选特征向量和相应的低浓度变量。

原理与方法

1.代谢组学数据

代谢组学数据的获取，是通过仪器检测生物体液中代谢产物的种类、含量以及状态变化，得到代谢组指纹色图谱，再通过代谢指纹色图谱获得目前常用的统计分析数据格式，如图1所示。

对于大量的小分子物质检测，色谱数据虽然不够精确，但它可以充分利用时间序列自相关信息，提取色谱峰相关性以及完整性的特征。将色谱数据分析作为整个代谢组学研究的前段工作，通过分析色谱数据，获得质谱数据分析不易识别的生物信息，色谱分析是代谢组学数据分析的前端工作，通过色谱数据差异变量筛选的结果结合保留时间对质谱数据进行定位，可以有针对性地对质谱的某一段(特征)进行重点分析。

2.小波变换

连续小波变换(continuous wavelet transform,CWT)有尺度a和b位移两个参数，其表达式可表示为

(1)

CWT(f,a,b)为信号函数f(t)在尺度a、位置b的小波变换系数；ψ(t)为满足一定条件的小波函数[4]。小波系数CWT(f,a,b)度量的是以b为中心，半径大小与a成比例的任何邻域内的信号f(t)的局部变化。在任何尺度因子a和平移因子b上，小波基函数的时-频窗面积是不变的，即时间、尺度分辨率是相互制约的。小尺度因子能够提取数据内部的局部特征，而大尺度因子显现整个数据信息的特征，细节不多，因此，需要结合研究目的选择尺度因子和平移因子的值。计算小波变化时，尺度因子和平移因子都需以小步长n增加，平移因子在尺度因子a处向右移动n个单位进行小波变换，完成时间-尺度因子相平面的采样，n决定了数据的采样点数。一维小波变换能够将需要处理的信号由某种局部变换进行不同尺度的分解与重构，建立相应尺度的一组模型，对信号的局部特征在不同尺度上进行描述和分析，不同尺度下分解只能分别分析不同尺度下的特征信号[5]。将一维代谢组数据进行多尺度小波变换，能够将微小局部特征按照尺度由第一层至最后一层逐渐增大，特征综合性愈来愈强地进行融合，获得二维小波系数图像，如图2所示。数据在任何尺度下的特征都可以通过对图像特征提取及模式识别进行分析。

图1 代谢组指纹图谱及数据格式

图2 相同保留时间段的正常对照和卵巢癌代谢组色谱转换成二维小波系数图像

3.特征提取

特征提取是通过对色谱时间序列信号或图像分析、变换来提取所需特征的一种方法[6]。对于上述二维小波系数图像，不同尺度的小波变换会得到不同大小的矩阵，对每一个样本的二维连续小波系数用不同子带中小波系数的和、标准差以及最大值表示纹理特征，本算法中将子带的统计量分别称为矩形和、矩形标准差、矩形最大值3个矩形特征模板。对于图像x×y矩阵，矩形特征筛选可以以a×b小矩阵为子矩阵对其进行分割，同时必须满足两个条件：x方向矩阵的边长必须能被自然数a整除，y方向矩阵的边长必须能被自然数b整除。由此可以获得(x/a)×(y/b)个大小为a×b的子矩阵。对分割的小矩阵，采用矩形和、矩形标准差、矩形最大值3种矩形特征对矩阵进行一次遍历计算，即按照从左到右和从上到下的顺序，构成新的特征向量数据集。使用二维小波系数图像提取特征向量的原理如图3所示。

图3 二维小波系数图像应用a×b子矩阵矩形和特征提取获得新向量的工作原理图

针对新的特征向量数据集，使用RF(随机森林)特征提取方法，按重要性排分筛选出对分类有作用的特征变量(矩形特征)，并利用交叉验证的方法对模型的分类效果进行评价[7]。最后，对选出的重要特征变量，通过分析对应的质谱数据获得有潜在意义的生物标志物。

实例分析

自2009年9月至2011年3月，纳入哈尔滨医科大学附属肿瘤医院妇科采集初次发现的卵巢癌患者，确定76例恶性卵巢癌(EOC)和77例卵巢良性肿瘤患者(BOT)选入最终检测数据。最后通过Waters公司UPLC/QTOF/MS系统处理，通过Masslynx软件获得代谢组指纹色谱数据，每份样本包含1600个变量。选择40例卵巢癌患者与40例卵巢囊肿患者进行训练，利用随机森林筛选变量重要性排分在前50位的变量建立模型。

由于代谢组学色谱峰信号波形与墨西哥草帽抛面轮廓线非常相似，因此选择Mexh小波函数对代谢组色谱数据进行变换。在本研究中，为了更好地突显数据的不同内在特征，选择尺度因子1到64，平移因子初始值为1，步长为1，进行64次数据采样。将实际卵巢癌代谢组色谱数据1600×153(1600为变量，153为样本例数)利用连续Mexh小波函数进行64个尺度变换后，获得153个大小为1600×64的矩阵，即153个二维小波系数图像。

进而，以10×8为子矩阵对每一个矩阵进行分割，利用3种统计特征提取方法(矩形和、矩形标准差、矩形最大值)进行特征提取，对矩阵进行一次遍历计算，即从左到右从上到下的顺序，构成新的特征向量数据集，大小为1280×153(1280为新的特征)。

筛选重要性排分在前20、50、100以及200的特征建立RF模型对测试数据集进行分类判别，并分别计算不同特征数目建立RF模型判别分类后的ROC曲线下面积AUC值；将上述步骤重复1000次，得到1000个AUC值。最后取平均AUC值为模型的最终判别效果如表1所示。将色谱数据经Mexh小波多尺度变换，以10×8为子矩阵分割矩形和特征提取的AUC值频数分布(图4)。从表1中可以看出色谱数据经过小波变换后较处理前AUC值0.708都有提高，以子矩阵10×8进行分割矩形和特征提取时分类效果最佳。使用特征变量的数目在20、50、100和200时预测效果相近。

表1 不同特征提取类型对色谱数据的分类效果比较(AUC,子矩阵为10×8)

利用RF模型筛选重要性排分在前50位的特征，并计算频数排在前20位的特征。第154位特征在1000次实验中，出现了874次。由于本实验是按照矩阵遍历的方式特征提取，因此很容易得出，第154位特征在色谱数据中的保留时间大约在2.05～2.16分钟，在小波系数图像中的定位如图5所示。

图4 小波变换前后不同数目特征建立RF模型分类的AUC值分布(矩形和)

进而，应用超高效液相色谱质谱连用仪设定较低的阈值获取质谱数据。为了获取更多的低浓度物质，本文设定阈值等于2，保留时间左右扩大0.02秒，即2.03～2.18分钟。通过这两个参数，可以获取包含2480个变量的质谱数据。针对这一段保留时间内的质谱数据，采用RF方法筛选出分类能力在前50位的变量，重复实验1000次，然后统计频数排在前20位的变量。对质谱数据采用RF筛选重要性排分在前50位的变量模型判别分析，重复实验1000次，得到曲线下面积为0.761。对筛选出的代谢物变量通过保留时间、一级质谱进行数据库检索，能够初步推测出其中的一些物质(表2)。

图5 小波系数图像中的保留时间定位图

表2 区分卵巢良恶性肿瘤的血浆内差异代谢物

经数据库与二级质谱的标准品比对，确定V140为2-哌啶酮，这是我们新发现的卵巢癌生物标志物。代谢组数据经小波变换后的小波系数图像，原本特征差别不大的低浓度代谢组经过多尺度小波变换，通过上述二维小波系数图像能够鉴别出有意义的低浓度“差异特征”，如图6所示。其他鉴定的代谢物有V622为人体必需氨基酸，V679为不饱和脂肪酸氧化中间产物，V549为神经递质和激素调节剂，V746为色氨酸代谢物，V298为氨基酸代谢产物，V705为半胱氨酸和蛋氨酸代谢，V456为尼古丁代谢物。通过与代谢组数据库与相关文献查阅，有些物质已被认定与卵巢癌、淋巴癌、结肠癌等疾病有关。

讨 论

1.本文利用连续小波变换具有时间序列性的特点，将其应用于代谢组卵巢癌色谱峰指纹图谱进行多尺度小波变换。随着小波函数尺度和位置的不断变化，色谱峰局部微小特征逐渐增大，RF模型的分类判别效果较原始数据AUC值得到显著的提高。

2.连续小波函数相同尺度变换后的数据，应选择适合的子矩阵与特征类型对其进行分析。实例分析表明，特征和对色谱数据进行特征提取能够获得更好的结果。不同小波函数、不同尺度变换后的数据及结果会有所不同，但考虑连续Mexh小波函数变换图形与色谱峰有很大的相似之处，因此本文选用Mexh函数。

*箭头为新发现低浓度物质标志物2-哌啶酮所在位置

3.针对Mexh连续小波变换后的特征进行分析，通过差异特征在色谱数据中的位置对质谱数据定位，然后重点分析这一段保留时间内的质谱数据。通过分析推测出若干潜在的低浓度生物标志物。实例表明，将连续小波变换应用于色谱指纹图谱分析中提取有差异的特征，特别是筛选低浓度的生物标志物具有研究与应用价值。

[1] Xia J,Sinelnikov I,Han Beomsoo,et al.MetaboAnalyst 3.0-making metabolomics more meaningful.Nucleic Acids Research,2015,43(1):251-257.

[2] Gruhl F,Länge K.Surface modification of an acoustic biosensor allowing the detection of low concentrations of cancer markers.Analytical Biochemistry,2012,420(2):188-190.

[3] Gao J,Lv F,Wang J,et al.Postoperative dynamic changes in the concentration of CK19-2G2 in lung cancer patients and the clinical value of this marker.Tumor Biology,2015,36(11):8295-8299.

[4] Zheng Y,Fan R,Qiu C,et al.An improved algorithm for peak detection in mass spectra based on continuous wavelet transform.International Journal of Mass Spectrometry,2016,409:53-58.

[5] Pandey J N,Jha N K,Singh O P.The continuous wavelet transform in n-dimensions.International Journal of Wavelets Multiresolution and Information Processing,2016,14(5):1650037.

[6] Zhang T,Chen W,Li M.AR based quadratic feature extraction in the VMD domain for the automated seizure detection of EEG using random forest classifier.Biomedical Signal Processing and Control,2017,31:550-559.

[7] Wang C,Zhang Y.Improving scoring-docking-screening powers of protein-ligand scoring functions using random forest.Journal of Computational Chemistry,2017,38(3):169-177.