鄧 魁 李貞子 侯 艷 李 康

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)

·論著·

基于二維最大重疊離散小波變換的代謝組質(zhì)譜數(shù)據(jù)的預(yù)處理方法*

鄧 魁 李貞子 侯 艷 李 康△

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)

目的通過(guò)二維最大重疊離散小波變換(maximal overlap discrete wavelet transform,MODWT)對(duì)代謝組學(xué)質(zhì)譜數(shù)據(jù)進(jìn)行預(yù)處理,去除一定的噪聲和批次效應(yīng)，提高分析方法的有效性和穩(wěn)定性。方法針對(duì)卵巢癌和卵巢囊腫的質(zhì)譜數(shù)據(jù)，選取Haar小波函數(shù)對(duì)其進(jìn)行二維MODWT的數(shù)據(jù)變換，獲得不同尺度的數(shù)據(jù)，再對(duì)其中的細(xì)節(jié)數(shù)據(jù)置0和進(jìn)行重構(gòu)。對(duì)預(yù)處理后的數(shù)據(jù)用隨機(jī)森林(RF)方法篩選差異變量和建立判別模型，評(píng)價(jià)預(yù)處理的效果。結(jié)果經(jīng)過(guò)二維MODWT處理后的質(zhì)譜數(shù)據(jù)建模判別效果明顯優(yōu)于使用原始數(shù)據(jù)得到的結(jié)果。結(jié)論針對(duì)質(zhì)譜數(shù)據(jù)，二維MODWT方法能夠更好地進(jìn)行特征提取，提高模型的判別能力，具有研究?jī)r(jià)值和應(yīng)用價(jià)值。

代謝組學(xué) 質(zhì)譜數(shù)據(jù) 數(shù)據(jù)預(yù)處理 小波變換

代謝組質(zhì)譜數(shù)據(jù)是指使用超高效液相色譜-質(zhì)譜聯(lián)用儀檢測(cè)血液、尿液或組織等生物樣品得到的數(shù)據(jù)，目的是衡量不同的內(nèi)源小分子的相對(duì)含量[1]。由于檢測(cè)結(jié)果受樣品的預(yù)處理、環(huán)境溫度、色譜柱效的改變等多種因素影響，數(shù)據(jù)中可能含有大量的噪聲和一定的批次效應(yīng)，使檢測(cè)得到的數(shù)據(jù)極不穩(wěn)定，從而影響數(shù)據(jù)的分析結(jié)果[2]。

與基因組學(xué)和蛋白質(zhì)組學(xué)數(shù)據(jù)相比，代謝組學(xué)數(shù)據(jù)由于代謝物的結(jié)構(gòu)具有更大的差異性，對(duì)其進(jìn)行標(biāo)準(zhǔn)化有更大的難度，目前尚無(wú)一種標(biāo)準(zhǔn)的公認(rèn)方法對(duì)代謝組學(xué)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化預(yù)處理[3]。針對(duì)代謝組學(xué)數(shù)據(jù)的預(yù)處理，現(xiàn)在常見的統(tǒng)計(jì)學(xué)方法有Z標(biāo)準(zhǔn)化、最大峰歸一化和總和峰歸一化等[4]。Z標(biāo)準(zhǔn)化能夠消除不同代謝物濃度數(shù)量級(jí)的差別，但無(wú)法去除噪聲的干擾。最大峰歸一化和總和峰歸一化能夠消除某些混雜因素對(duì)代謝物濃度的干擾，其局限性在于它們假定在機(jī)體紊亂時(shí)，一組代謝物濃度的上升同時(shí)伴隨著另一組代謝物濃度的下降，而這種假定在實(shí)際中往往并不成立[5]。本文給出一種新的基于二維最大重疊Haar小波變換的數(shù)據(jù)預(yù)處理方法，這種方法利用了代謝組學(xué)質(zhì)譜數(shù)據(jù)時(shí)間序列的特性，將數(shù)據(jù)分解成不同的特征，從而能夠較好地去除由于不穩(wěn)定因素引起的數(shù)據(jù)波動(dòng)，提高判別分析的效果。

原理與方法

代謝組質(zhì)譜數(shù)據(jù)可記為D={xij}(i=1,2,…,n;j=1,2,…,m)，其中n和m分別表示檢測(cè)的樣本數(shù)和物質(zhì)(變量)數(shù)目。

二維最大重疊離散小波變換(maximal overlap discrete wavelet transform，MODWT)是從二維離散小波變化改進(jìn)得來(lái)的，但又與傳統(tǒng)的二維離散小波變化有所不同[6]。二維離散小波變換的分解和重構(gòu)過(guò)程如圖1所示，即首先對(duì)原始數(shù)據(jù)的每一行樣品進(jìn)行小波變換，獲得低頻分量L和高頻分量H；然后再進(jìn)一步作小波變換的列變換，得到行列都為低頻分量(LL1)、行為低頻分量列為高頻分量(LH1)、行為高頻分量列為低頻分量(HL1)、行列皆為高頻分量(HH1)四部分?jǐn)?shù)據(jù)。第二次數(shù)據(jù)變換可以在LL1上繼續(xù)進(jìn)行得到LL2、LH2、HL2和HH2。以此類推可以變換多次。其中低頻分量在一定程度上反映了數(shù)據(jù)本質(zhì)信息，本質(zhì)特征。而高頻分量則在一定程度上反映了噪聲信息，如在實(shí)驗(yàn)過(guò)程中環(huán)境溫度、色譜柱效的改變等多種系統(tǒng)誤差的干擾以及隨機(jī)誤差的影響。最后得到的變換可以通過(guò)逆運(yùn)算，完全重構(gòu)回原始數(shù)據(jù)。實(shí)際中，可以對(duì)分解后的數(shù)據(jù)進(jìn)行適當(dāng)處理，如對(duì)某部分?jǐn)?shù)據(jù)置0或按一定閾值置0后進(jìn)行重構(gòu)，從而達(dá)到去除無(wú)用的檢測(cè)信號(hào)的目的。

圖1 二維離散小波變換分解和重構(gòu)過(guò)程示意圖

上述小波行變換的公式為：

(1)

(2)

(3)

其中i表示相應(yīng)的行，k表示變換后相應(yīng)列的數(shù)據(jù)，L(i,k)表示變換后的低頻分量，H(i,k)表示變換后的高頻分量。h(·)為尺度函數(shù)，g(·)為小波函數(shù)。

進(jìn)一步，對(duì)L和H作如下分解：

(4)

(5)

(6)

(7)

(8)

同理可以得到D3。上式中LL(k,j)表示對(duì)L作小波列變換得到的低頻分量，LH(k,j)表示對(duì)L作小波列變換得到的高頻分量，HL(k,j)表示對(duì)H作小波列變換得到的低頻分量，HH(k,j)表示對(duì)H作小波列變換得到的高頻分量。h(·)為尺度函數(shù),g(·)為小波函數(shù)。

本文使用Haar小波，其尺度函數(shù)為：

(9)

Haar小波函數(shù)定義為尺度函數(shù)的對(duì)偶函數(shù)，即g(x):=h(2x)-h(2x-1)。這是一個(gè)階躍函數(shù)，它具有非連續(xù)性，具有在時(shí)域區(qū)間內(nèi)局部化能力，其函數(shù)為

(10)

二維MODWT更適合代謝組學(xué)數(shù)據(jù)分析，因?yàn)樗试S小波系數(shù)和原始數(shù)據(jù)變量相對(duì)應(yīng)。根據(jù)小波變換理論，最有意義的信息包含在較大的系數(shù)中，噪聲一般在較小的系數(shù)中，因此，僅重構(gòu)較大的系數(shù)就可以保持原始數(shù)據(jù)的本質(zhì)特征。

本研究用上述二維最大重疊離散小波變換對(duì)質(zhì)譜數(shù)據(jù)進(jìn)行預(yù)處理，然后用隨機(jī)森林(random forest，RF)模型進(jìn)行判別分析。

實(shí)例分析

實(shí)例：2011年3月至2013年7月從哈爾濱醫(yī)科大學(xué)附屬腫瘤醫(yī)院婦科收集了140例上皮性卵巢癌(EOC)患者和158例卵巢囊腫(BOT)患者血漿樣品。得到超高效液相色譜質(zhì)譜聯(lián)用檢測(cè)結(jié)果。原始數(shù)據(jù)通過(guò)MASS-Hunter定性分析軟件轉(zhuǎn)化為mzdata格式的文件，這些文件導(dǎo)入到R中通過(guò)XCMS包進(jìn)行預(yù)處理。最后得到2106個(gè)變量。

用二維MODWT對(duì)原始數(shù)據(jù)分別進(jìn)行一層和兩層分解，得到D2和D3(參見圖1)，其中HH1和HH2部分主要包含細(xì)節(jié)噪聲信息，對(duì)其進(jìn)行置0處理，然后對(duì)數(shù)據(jù)重構(gòu)。為了驗(yàn)證預(yù)處理的效果，在140例EOC患者中隨機(jī)抽取70例患者，158例BOT患者中隨機(jī)抽取79例患者作為訓(xùn)練集，其余樣品作為測(cè)試集。在訓(xùn)練集中使用RF方法篩選出重要性評(píng)分在前50位的特征建立判別模型，對(duì)測(cè)試集樣品進(jìn)行分類和預(yù)測(cè)，通過(guò)真實(shí)分類標(biāo)簽計(jì)算ROC曲線下的面積AUC值。上述試驗(yàn)通過(guò)無(wú)重復(fù)抽樣進(jìn)行1000次，最后通過(guò)1000個(gè)AUC的平均值考核方法的有效性。小波變換和RF判別模型均通過(guò)R語(yǔ)言程序包實(shí)現(xiàn)。

試驗(yàn)結(jié)果：經(jīng)RF判別后得到AUC值后的結(jié)果如圖2、圖3和圖4所示。由圖看出，代謝組質(zhì)譜原始數(shù)據(jù)進(jìn)行RF判別分析后得到AUC值為0.87±0.03，經(jīng)兩種不同變換再進(jìn)行RF判別得到AUC值分別為0.94±0.02和0.96±0.02。圖2和圖3給出了Haar小波變換與原始數(shù)據(jù)的AUC值頻數(shù)圖，結(jié)果顯示經(jīng)Haar小波變換后，卵巢癌質(zhì)譜數(shù)據(jù)的分類能力有了很大的提高，并且隨著分解層數(shù)從第一層增加到第二層，分類能力進(jìn)一步提高。圖4給出了原始數(shù)據(jù)，一層小波變換后數(shù)據(jù)，兩層小波變換后數(shù)據(jù)分別經(jīng)RF判別得到的判別結(jié)果的ROC曲線圖，結(jié)果顯示二維MODWT能明顯提高RF的分類能力。

圖2 一層Haar小波變換后AUC值頻數(shù)變化情況

圖3 兩層Haar小波變換后AUC值頻數(shù)變化情況

圖4 原始數(shù)據(jù)和兩種不同層小波變換的RF模型的ROC曲線

討 論

本文嘗試將二維小波變換應(yīng)用于代謝組學(xué)數(shù)據(jù)預(yù)處理，其基本思想是利用質(zhì)譜相鄰間物質(zhì)具有一定的相關(guān)性、檢測(cè)樣品間時(shí)間上的相關(guān)性，以及同類觀測(cè)樣品(如癌癥和非癌)的聚類性質(zhì)，進(jìn)行二維小波變換，適當(dāng)去除干擾信號(hào)，獲得更好的分類效果。通過(guò)小波變換后，RF模型的判別效果得到顯著的提高。

如果將本文的實(shí)例通過(guò)隨機(jī)打亂樣品標(biāo)簽后再進(jìn)行二維MODWT處理，然后進(jìn)行RF判別后得到的AUC值均在0.5附近，這進(jìn)一步證明了本研究得到結(jié)果的可信度。

在之前的研究中，我們使用常用的Z標(biāo)準(zhǔn)化、最大峰歸一化和總峰和歸一化的方法對(duì)代謝組學(xué)數(shù)據(jù)進(jìn)行預(yù)處理，發(fā)現(xiàn)RF的判別能力并不能得到提升，其AUC值分別為0.87±0.03、0.87±0.03，0.85±0.03。

不同的小波函數(shù)也有不同的特點(diǎn)。因此在二維離散小波變換中，需要選擇合適的小波函數(shù)。研究表明，Haar小波函數(shù)對(duì)代謝組質(zhì)譜峰數(shù)據(jù)進(jìn)行質(zhì)譜峰提取能夠獲得更好的結(jié)果[8]，因此本文選用Haar小波。然而，實(shí)際中究竟選用何種小波函數(shù)好，與數(shù)據(jù)結(jié)構(gòu)和研究目的有關(guān)，對(duì)此需要進(jìn)一步研究。

本研究中分解層數(shù)為一層和二層時(shí)，達(dá)到了較好的預(yù)測(cè)能力。實(shí)際我們?cè)谶M(jìn)行更多層次處理后發(fā)現(xiàn)，分類能力雖然越來(lái)越強(qiáng)，甚至接近1，但隨機(jī)打亂分類標(biāo)簽再進(jìn)行處理后，RF的判別效果依然很高，甚至接近0.9，這說(shuō)明這種分類能力提高是虛假的提高，可能出現(xiàn)過(guò)擬合的現(xiàn)象。在實(shí)際研究中，要根據(jù)數(shù)據(jù)結(jié)構(gòu)和特點(diǎn)選擇最佳的分解層數(shù)。

本研究通過(guò)二維最大重疊離散小波變換，將數(shù)據(jù)分解為行列都為低頻分量(LL)、行為低頻分量列為高頻分量(LH)、行為高頻分量列為低頻分量(HL)、行列皆為高頻分量(HH)四部分?jǐn)?shù)據(jù)，考慮到細(xì)節(jié)系數(shù)(HH)主要包含噪聲信息[10]，因此將其置0處理，再進(jìn)行重構(gòu)，從而達(dá)到去除一定的隨機(jī)波動(dòng)，保留原始數(shù)據(jù)主要特征信息的目的，這只是一種簡(jiǎn)單的處理方法。更好的方法應(yīng)該根據(jù)實(shí)際數(shù)據(jù)，對(duì)行為低頻分量列為高頻分量(LH)、行為高頻分量列為低頻分量(HL)這兩部分?jǐn)?shù)據(jù)有選擇地使用硬閾值[11]或軟閾值[12]的方法進(jìn)行處理，既能去除干擾信號(hào)，又能夠保留有效的檢測(cè)數(shù)據(jù)信息。對(duì)此有待于更進(jìn)一步的研究。

[1] Wang Y,Hu H,Su Y,et al.Potential of monitoring isotopologues by quantitative gas chromatography with time-of-flight mass spectrometry for metabolomic assay.Journal of Separation Science,2016,39(6):1137-1143.

[2] De Livera AM,Sysi-Aho M,Jacob L,et al.Statistical methods for handling unwanted variation in metabolomics data.Analytical Chemistry,2015,87(7):3606-3621.

[3] Wu Y,Li L.Sample normalization methods in quantitative metabolomics.Journal of Chromatography A,2016,1430:80-95.

[4] 柯朝甫,張濤,武曉巖,等.代謝組學(xué)數(shù)據(jù)分析的統(tǒng)計(jì)學(xué)方法.中國(guó)衛(wèi)生統(tǒng)計(jì),2014,31(2):357-359.

[5] Sysi-Aho M,Katajamaa M,Yetukuri L,et al.Normalization method for metabolomics data using optimal selection of multiple internal standards.BMC Bioinformatics,2007,8(1):1-17.

[6] P Liò.Wavelets in bioinformatics and computational biology:state of art and perspectives.Bioinformatics,2003,19(1):2-9.

[7] 彭娟,李川.基于最大重疊離散小波變換的油中顆粒污染物特征信號(hào)提取.重慶工商大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,30(6):24-28.

[8] Davis RA,Chariton AJ,Godward J,et al.Adaptive binning:An improved binning method for metabolomics data using the undecimated wavelet transform.Chemometrics & Intelligent Laboratory Systems,2007.85(1):144-154.

[9] Yang Y,He Y,Cheng J,et al.A gear fault diagnosis using Hilbert spectrum based on MODWPT and a comparison with EMD approach.Measurement,2009,42(4):542-551.

[10]Lewis AS,Knowles G.Image compression using the 2-D wavelet transform.IEEE Transactions on Image Processing,1992,1(2):244-50.

[11]Chen C,Zhou N.A new wavelet hard threshold to process image with strong Gaussian Noise.IEEE Fifth International Conference on Advanced Computational Intelligence,2012:558-561.

[12]Fang Y.A Method of Wavelet Image Enhancement Based on Soft Threshold.Computer Engineering & Applications,2002,38(23):16-19.

ThePreprocessingMethodofMetabolomicMassSpectrumDataBasedontheTwo-dimensionalMaximalOverlapDiscreteWaveletTransform

Deng Kui,Li Zhenzi,Hou Yan,et al

(DepartmentofMedicalStatistics,HarbinMedicalUniversity(150081),Harbin)

ObjectiveTo preprocess metabolomic mass spectrum data through using the two-dimensional maximal overlap discrete wavelet transform(MODWT)with the purpose of removing noise and batch effects to some extent and improving the effectiveness and stability of the analytical methods.MethodsWe conducted the two-dimensional MODWT with the wavelet function of Haar to metabolomic mass spectrum data of ovarian cancer and ovarian cyst and obtained the data of different scales.And then we set the detail data to zero and reconstructed the data.After that,the random forest method was applied to the preprocessed data to screen variables,establish the discrimination model and evaluate the effects of preprocessing.ResultsThe classification performance of mass spectrum data using the two-dimensional MODWT is obviously better than the original data.ConclusionFor the mass spectrum data,the two-dimensional MODWT can well conduct feature extraction and improve the discriminant ability of the model and it has the research value and application value.

Metabolomics;Mass spectrometric data;Data preprocessing;Wavelet transform

國(guó)家自然科學(xué)基金(81302511,81473072)；哈爾濱醫(yī)科大學(xué)創(chuàng)新科學(xué)研究基金(2016JCZX13)；哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)

△通信作者：李康，E-mail:likang@ems.hrbmu.edu.cn

郭海強(qiáng))