基于小波變換特征提取的代謝組低濃度標志物的篩選*

哈爾濱醫(yī)科大學衛(wèi)生統計學教研室(150081) 李貞子 侯 艷 鄧 魁 李 康

代謝組學研究方法主要是通過色譜和質譜等儀器獲得組織、血液中的代謝物質，特別是尿液內源小分子代謝物的狀態(tài)和含量[1]。由于質譜數據能夠提供具體代謝物質化學結構的生物學信息，因此大部分的傳統方法主要針對的是質譜數據。然而，大部分的傳統方法在實際應用中尚存在各種不同的問題，其中最主要的問題是不易進行低濃度標志物識別[2]。質譜數據變量選取數目的多少與預處理設定的閾值有關，如果閾值設定較低，其質譜的變量個數可能達到10萬以上，從而使分析不易進行。為了減少變量的個數通常將閾值設置在一個較高的水平，即對閾值以下的所有低濃度物質的含量置零。而根據生物學知識可知，腫瘤標志物更有意義往往是低濃度的化合物[3]。本文擬結合卵巢癌色譜和質譜的數據，使用小波變換的方法將原始色譜數據變換為二維圖像，在此基礎上使用隨機森林(random forest,RF)模型篩選特征向量和相應的低濃度變量。

原理與方法

1.代謝組學數據

代謝組學數據的獲取，是通過儀器檢測生物體液中代謝產物的種類、含量以及狀態(tài)變化，得到代謝組指紋色圖譜，再通過代謝指紋色圖譜獲得目前常用的統計分析數據格式，如圖1所示。

對于大量的小分子物質檢測，色譜數據雖然不夠精確，但它可以充分利用時間序列自相關信息，提取色譜峰相關性以及完整性的特征。將色譜數據分析作為整個代謝組學研究的前段工作，通過分析色譜數據，獲得質譜數據分析不易識別的生物信息，色譜分析是代謝組學數據分析的前端工作，通過色譜數據差異變量篩選的結果結合保留時間對質譜數據進行定位，可以有針對性地對質譜的某一段(特征)進行重點分析。

2.小波變換

連續(xù)小波變換(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上，小波基函數的時-頻窗面積是不變的，即時間、尺度分辨率是相互制約的。小尺度因子能夠提取數據內部的局部特征，而大尺度因子顯現整個數據信息的特征，細節(jié)不多，因此，需要結合研究目的選擇尺度因子和平移因子的值。計算小波變化時，尺度因子和平移因子都需以小步長n增加，平移因子在尺度因子a處向右移動n個單位進行小波變換，完成時間-尺度因子相平面的采樣，n決定了數據的采樣點數。一維小波變換能夠將需要處理的信號由某種局部變換進行不同尺度的分解與重構，建立相應尺度的一組模型，對信號的局部特征在不同尺度上進行描述和分析，不同尺度下分解只能分別分析不同尺度下的特征信號[5]。將一維代謝組數據進行多尺度小波變換，能夠將微小局部特征按照尺度由第一層至最后一層逐漸增大，特征綜合性愈來愈強地進行融合，獲得二維小波系數圖像，如圖2所示。數據在任何尺度下的特征都可以通過對圖像特征提取及模式識別進行分析。

圖1 代謝組指紋圖譜及數據格式

圖2 相同保留時間段的正常對照和卵巢癌代謝組色譜轉換成二維小波系數圖像

3.特征提取

特征提取是通過對色譜時間序列信號或圖像分析、變換來提取所需特征的一種方法[6]。對于上述二維小波系數圖像，不同尺度的小波變換會得到不同大小的矩陣，對每一個樣本的二維連續(xù)小波系數用不同子帶中小波系數的和、標準差以及最大值表示紋理特征，本算法中將子帶的統計量分別稱為矩形和、矩形標準差、矩形最大值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月，納入哈爾濱醫(yī)科大學附屬腫瘤醫(yī)院婦科采集初次發(fā)現的卵巢癌患者，確定76例惡性卵巢癌(EOC)和77例卵巢良性腫瘤患者(BOT)選入最終檢測數據。最后通過Waters公司UPLC/QTOF/MS系統處理，通過Masslynx軟件獲得代謝組指紋色譜數據，每份樣本包含1600個變量。選擇40例卵巢癌患者與40例卵巢囊腫患者進行訓練，利用隨機森林篩選變量重要性排分在前50位的變量建立模型。

由于代謝組學色譜峰信號波形與墨西哥草帽拋面輪廓線非常相似，因此選擇Mexh小波函數對代謝組色譜數據進行變換。在本研究中，為了更好地突顯數據的不同內在特征，選擇尺度因子1到64，平移因子初始值為1，步長為1，進行64次數據采樣。將實際卵巢癌代謝組色譜數據1600×153(1600為變量，153為樣本例數)利用連續(xù)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 區(qū)分卵巢良惡性腫瘤的血漿內差異代謝物

經數據庫與二級質譜的標準品比對，確定V140為2-哌啶酮，這是我們新發(fā)現的卵巢癌生物標志物。代謝組數據經小波變換后的小波系數圖像，原本特征差別不大的低濃度代謝組經過多尺度小波變換，通過上述二維小波系數圖像能夠鑒別出有意義的低濃度“差異特征”，如圖6所示。其他鑒定的代謝物有V622為人體必需氨基酸，V679為不飽和脂肪酸氧化中間產物，V549為神經遞質和激素調節(jié)劑，V746為色氨酸代謝物，V298為氨基酸代謝產物，V705為半胱氨酸和蛋氨酸代謝，V456為尼古丁代謝物。通過與代謝組數據庫與相關文獻查閱，有些物質已被認定與卵巢癌、淋巴癌、結腸癌等疾病有關。

討 論

1.本文利用連續(xù)小波變換具有時間序列性的特點，將其應用于代謝組卵巢癌色譜峰指紋圖譜進行多尺度小波變換。隨著小波函數尺度和位置的不斷變化，色譜峰局部微小特征逐漸增大，RF模型的分類判別效果較原始數據AUC值得到顯著的提高。

2.連續(xù)小波函數相同尺度變換后的數據，應選擇適合的子矩陣與特征類型對其進行分析。實例分析表明，特征和對色譜數據進行特征提取能夠獲得更好的結果。不同小波函數、不同尺度變換后的數據及結果會有所不同，但考慮連續(xù)Mexh小波函數變換圖形與色譜峰有很大的相似之處，因此本文選用Mexh函數。

*箭頭為新發(fā)現低濃度物質標志物2-哌啶酮所在位置

3.針對Mexh連續(xù)小波變換后的特征進行分析，通過差異特征在色譜數據中的位置對質譜數據定位，然后重點分析這一段保留時間內的質譜數據。通過分析推測出若干潛在的低濃度生物標志物。實例表明，將連續(xù)小波變換應用于色譜指紋圖譜分析中提取有差異的特征，特別是篩選低濃度的生物標志物具有研究與應用價值。

[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.