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

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

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

原理與方法

1.代謝組學(xué)數(shù)據(jù)

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

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

2.小波變換

連續(xù)小波變換(continuous wavelet transform,CWT)有尺度a和b位移兩個參數(shù)，其表達式可表示為

(1)

CWT(f,a,b)為信號函數(shù)f(t)在尺度a、位置b的小波變換系數(shù)；ψ(t)為滿足一定條件的小波函數(shù)[4]。小波系數(shù)CWT(f,a,b)度量的是以b為中心，半徑大小與a成比例的任何鄰域內(nèi)的信號f(t)的局部變化。在任何尺度因子a和平移因子b上，小波基函數(shù)的時-頻窗面積是不變的，即時間、尺度分辨率是相互制約的。小尺度因子能夠提取數(shù)據(jù)內(nèi)部的局部特征，而大尺度因子顯現(xiàn)整個數(shù)據(jù)信息的特征，細節(jié)不多，因此，需要結(jié)合研究目的選擇尺度因子和平移因子的值。計算小波變化時，尺度因子和平移因子都需以小步長n增加，平移因子在尺度因子a處向右移動n個單位進行小波變換，完成時間-尺度因子相平面的采樣，n決定了數(shù)據(jù)的采樣點數(shù)。一維小波變換能夠?qū)⑿枰幚淼男盘栍赡撤N局部變換進行不同尺度的分解與重構(gòu)，建立相應(yīng)尺度的一組模型，對信號的局部特征在不同尺度上進行描述和分析，不同尺度下分解只能分別分析不同尺度下的特征信號[5]。將一維代謝組數(shù)據(jù)進行多尺度小波變換，能夠?qū)⑽⑿【植刻卣靼凑粘叨扔傻谝粚又磷詈笠粚又饾u增大，特征綜合性愈來愈強地進行融合，獲得二維小波系數(shù)圖像，如圖2所示。數(shù)據(jù)在任何尺度下的特征都可以通過對圖像特征提取及模式識別進行分析。

圖1 代謝組指紋圖譜及數(shù)據(jù)格式

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

3.特征提取

特征提取是通過對色譜時間序列信號或圖像分析、變換來提取所需特征的一種方法[6]。對于上述二維小波系數(shù)圖像，不同尺度的小波變換會得到不同大小的矩陣，對每一個樣本的二維連續(xù)小波系數(shù)用不同子帶中小波系數(shù)的和、標準差以及最大值表示紋理特征，本算法中將子帶的統(tǒng)計量分別稱為矩形和、矩形標準差、矩形最大值3個矩形特征模板。對于圖像x×y矩陣，矩形特征篩選可以以a×b小矩陣為子矩陣對其進行分割，同時必須滿足兩個條件：x方向矩陣的邊長必須能被自然數(shù)a整除，y方向矩陣的邊長必須能被自然數(shù)b整除。由此可以獲得(x/a)×(y/b)個大小為a×b的子矩陣。對分割的小矩陣，采用矩形和、矩形標準差、矩形最大值3種矩形特征對矩陣進行一次遍歷計算，即按照從左到右和從上到下的順序，構(gòu)成新的特征向量數(shù)據(jù)集。使用二維小波系數(shù)圖像提取特征向量的原理如圖3所示。

圖3 二維小波系數(shù)圖像應(yīng)用a×b子矩陣矩形和特征提取獲得新向量的工作原理圖

針對新的特征向量數(shù)據(jù)集，使用RF(隨機森林)特征提取方法，按重要性排分篩選出對分類有作用的特征變量(矩形特征)，并利用交叉驗證的方法對模型的分類效果進行評價[7]。最后，對選出的重要特征變量，通過分析對應(yīng)的質(zhì)譜數(shù)據(jù)獲得有潛在意義的生物標志物。

實例分析

自2009年9月至2011年3月，納入哈爾濱醫(yī)科大學(xué)附屬腫瘤醫(yī)院婦科采集初次發(fā)現(xiàn)的卵巢癌患者，確定76例惡性卵巢癌(EOC)和77例卵巢良性腫瘤患者(BOT)選入最終檢測數(shù)據(jù)。最后通過Waters公司UPLC/QTOF/MS系統(tǒng)處理，通過Masslynx軟件獲得代謝組指紋色譜數(shù)據(jù)，每份樣本包含1600個變量。選擇40例卵巢癌患者與40例卵巢囊腫患者進行訓(xùn)練，利用隨機森林篩選變量重要性排分在前50位的變量建立模型。

由于代謝組學(xué)色譜峰信號波形與墨西哥草帽拋面輪廓線非常相似，因此選擇Mexh小波函數(shù)對代謝組色譜數(shù)據(jù)進行變換。在本研究中，為了更好地突顯數(shù)據(jù)的不同內(nèi)在特征，選擇尺度因子1到64，平移因子初始值為1，步長為1，進行64次數(shù)據(jù)采樣。將實際卵巢癌代謝組色譜數(shù)據(jù)1600×153(1600為變量，153為樣本例數(shù))利用連續(xù)Mexh小波函數(shù)進行64個尺度變換后，獲得153個大小為1600×64的矩陣，即153個二維小波系數(shù)圖像。

進而，以10×8為子矩陣對每一個矩陣進行分割，利用3種統(tǒng)計特征提取方法(矩形和、矩形標準差、矩形最大值)進行特征提取，對矩陣進行一次遍歷計算，即從左到右從上到下的順序，構(gòu)成新的特征向量數(shù)據(jù)集，大小為1280×153(1280為新的特征)。

篩選重要性排分在前20、50、100以及200的特征建立RF模型對測試數(shù)據(jù)集進行分類判別，并分別計算不同特征數(shù)目建立RF模型判別分類后的ROC曲線下面積AUC值；將上述步驟重復(fù)1000次，得到1000個AUC值。最后取平均AUC值為模型的最終判別效果如表1所示。將色譜數(shù)據(jù)經(jīng)Mexh小波多尺度變換，以10×8為子矩陣分割矩形和特征提取的AUC值頻數(shù)分布(圖4)。從表1中可以看出色譜數(shù)據(jù)經(jīng)過小波變換后較處理前AUC值0.708都有提高，以子矩陣10×8進行分割矩形和特征提取時分類效果最佳。使用特征變量的數(shù)目在20、50、100和200時預(yù)測效果相近。

表1 不同特征提取類型對色譜數(shù)據(jù)的分類效果比較(AUC,子矩陣為10×8)

利用RF模型篩選重要性排分在前50位的特征，并計算頻數(shù)排在前20位的特征。第154位特征在1000次實驗中，出現(xiàn)了874次。由于本實驗是按照矩陣遍歷的方式特征提取，因此很容易得出，第154位特征在色譜數(shù)據(jù)中的保留時間大約在2.05～2.16分鐘，在小波系數(shù)圖像中的定位如圖5所示。

圖4 小波變換前后不同數(shù)目特征建立RF模型分類的AUC值分布(矩形和)

進而，應(yīng)用超高效液相色譜質(zhì)譜連用儀設(shè)定較低的閾值獲取質(zhì)譜數(shù)據(jù)。為了獲取更多的低濃度物質(zhì)，本文設(shè)定閾值等于2，保留時間左右擴大0.02秒，即2.03～2.18分鐘。通過這兩個參數(shù)，可以獲取包含2480個變量的質(zhì)譜數(shù)據(jù)。針對這一段保留時間內(nèi)的質(zhì)譜數(shù)據(jù)，采用RF方法篩選出分類能力在前50位的變量，重復(fù)實驗1000次，然后統(tǒng)計頻數(shù)排在前20位的變量。對質(zhì)譜數(shù)據(jù)采用RF篩選重要性排分在前50位的變量模型判別分析，重復(fù)實驗1000次，得到曲線下面積為0.761。對篩選出的代謝物變量通過保留時間、一級質(zhì)譜進行數(shù)據(jù)庫檢索，能夠初步推測出其中的一些物質(zhì)(表2)。

圖5 小波系數(shù)圖像中的保留時間定位圖

表2 區(qū)分卵巢良惡性腫瘤的血漿內(nèi)差異代謝物

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

討 論

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

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

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

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

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