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

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

原理與方法

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

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

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

2.小波變換

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

(1)

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

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

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

3.特征提取

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

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

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

實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

討 論

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

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

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

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

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