基因組學(xué)數(shù)據(jù)整合中的批次效應(yīng)移除算法*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室（150081） 陳天成 侯 艷 李 康

基因組學(xué)數(shù)據(jù)整合中的批次效應(yīng)移除算法*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室（150081） 陳天成 侯 艷 李 康△

微陣列的基因表達(dá)譜已經(jīng)成為研究癌癥細(xì)胞轉(zhuǎn)錄的重要工具，且成功用于許多腫瘤識(shí)別和與癌癥相關(guān)的生物標(biāo)志物研究［1－2］。然而，實(shí)際中需要大量的樣本才能夠得到穩(wěn)定的分析結(jié)果，而由于實(shí)際的復(fù)雜程度及費(fèi)用問題，每個(gè)數(shù)據(jù)集的樣本量受到了限制。為此，可以將已公開的微陣列數(shù)據(jù)整合，在減少花費(fèi)同時(shí)提高結(jié)果的穩(wěn)定性。

研究者在數(shù)據(jù)整合過程中常面對(duì)由于平臺(tái)設(shè)計(jì)、合成及探針注釋不同帶來的跨平臺(tái)基因表達(dá)研究的難題［3］。另外，同一方面的微陣列實(shí)驗(yàn)研究可能在不同時(shí)間或地點(diǎn)進(jìn)行，其系統(tǒng)誤差會(huì)給數(shù)據(jù)整合帶來困難。上述系統(tǒng)誤差在文獻(xiàn)上稱之為批次效應(yīng)（batch effect）。Scherer和Andreasd將批次效應(yīng)定義為樣本在不同批次處理和測(cè)量中系統(tǒng)上的技術(shù)差異［4］。當(dāng)批次效應(yīng)引入的誤差足夠強(qiáng)以至于混雜真實(shí)的生物學(xué)差異時(shí)，未經(jīng)移除批次效應(yīng)的數(shù)據(jù)分析可能出現(xiàn)誤導(dǎo)性的結(jié)果［5］。微陣列信號(hào)強(qiáng)度的歸一化（normalization）是一種標(biāo)準(zhǔn)化分析技術(shù)，如MAS5和RMA方法，但這是一種在特定數(shù)據(jù)集內(nèi)部的標(biāo)準(zhǔn)化方法，無法移除批次效應(yīng)［6］。如何合理有效地移除批次效應(yīng)是目前基因組研究的熱點(diǎn)之一，本文將對(duì)此做一探討。

統(tǒng)計(jì)分析算法

基因組學(xué)中數(shù)據(jù)批次效應(yīng)移除的整合算法很多，本文主要介紹目前已經(jīng)明確有較好效果［7－8］及新近提出的方法。

1.經(jīng)驗(yàn)貝葉斯方法

經(jīng)驗(yàn)貝葉斯（empirical Bayes）［9］方法的思想是從基因和實(shí)驗(yàn)條件中使用“先驗(yàn)信息”，期望利用先驗(yàn)信息得到更好的估計(jì)或者更穩(wěn)定的推斷。模型的假設(shè)是基于位置和尺度的調(diào)整，形式如下:

其中，Yijg是批次 i（i＝1，2，…，b）、樣品 j（j＝1，2，…，ni）和基因 g（g＝1，2，…，mi）的表達(dá)值，αg是基因 g的平均表達(dá)值，X是樣品條件的設(shè)計(jì)矩陣，βg是對(duì)應(yīng)X的回歸系數(shù)的向量；誤差項(xiàng) εijg假設(shè)服從 N（0，σg）；γig和δig表示批次i中基因g加法和乘法的批次效應(yīng)。

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

算法分為三步:

（1）標(biāo)準(zhǔn)化數(shù)據(jù)，即通過基因表達(dá)數(shù)據(jù)，用最小二乘法和約束條件∑ini＝0得到估計(jì)值和，方差估計(jì)為

再通過下式把標(biāo)準(zhǔn)化數(shù)據(jù)Zijg計(jì)算出來，即

（2）估計(jì)批次效應(yīng)參數(shù)。假設(shè) Zijg～N（γig，δig），采用參數(shù)先驗(yàn)或者非參數(shù)先驗(yàn)方法計(jì)算出批次效應(yīng)估計(jì)值和。其中參數(shù)先驗(yàn)方法要求 γig～N（γi，）及～I(xiàn)nverseGamma（λi，θi）。

2.改進(jìn)ComBat方法

改進(jìn)ComBat［10］是前述經(jīng)驗(yàn)貝葉斯方法的改進(jìn)方法，該法通過改變參數(shù)的估計(jì)值及為批次水平的估計(jì)值及，實(shí)現(xiàn)把樣本轉(zhuǎn)換到“金標(biāo)準(zhǔn)”參考批次的均數(shù)和方差，而不是總均值和合并方差。標(biāo)準(zhǔn)化數(shù)據(jù)表示如下:

最后批次效應(yīng)調(diào)整數(shù)據(jù)中使用的均值和方差估計(jì)值通過參考批次的參數(shù)估計(jì)值進(jìn)行估計(jì)。

M-ComBat調(diào)整的數(shù)據(jù)可由下式計(jì)算得到，即

其中符號(hào)的意義同式（4）。

3.跨平臺(tái)歸一化

跨平臺(tái)歸一化（cross-platform normalization，XPN）方法［11］，其基本思想是識(shí)別在兩個(gè)研究中具有相似表達(dá)特征的基因和樣本的同質(zhì)性區(qū)組。數(shù)據(jù)模型構(gòu)建基于簡(jiǎn)單區(qū)組線性模型，即對(duì)于每個(gè)基因，其表達(dá)式如下:

其中Yijg為批次 i、樣品 j和基因 g的觀測(cè)值，α*:｛1，…，m｝→｛1，…，K｝和:｛1，…，ni｝→｛1，…，L｝，i＝1，2分別為聚類后的基因和樣本的組別，Aα*（g），β*i（j），i是區(qū)組的均數(shù)，big和cig分別表示不同批次和特定基因的靈敏度及偏移參數(shù)，噪聲變量εijg是相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布。

算法首先對(duì)研究數(shù)據(jù)進(jìn)行樣本標(biāo)準(zhǔn)化和中心化，然后采用k-均值聚類法分別對(duì)樣品和基因聚類，給出分配函數(shù)α和β。進(jìn)而，對(duì)函數(shù)αg和βi（j）使用極大似然方法得到模型參數(shù)的估計(jì)值，各次聚類后基因表達(dá)值為

上述過程需要進(jìn)行30次，算法的輸出結(jié)果是重復(fù)運(yùn)行獲得的歸一化平均值。

4.平臺(tái)獨(dú)立的LDA模型

平臺(tái)獨(dú)立的LDA（platform independent LDA，PLIDA）模型［12］，假設(shè)不存在批次引起的差異，不同批次的基因表達(dá)芯片獲得的數(shù)據(jù)能夠合并在一起進(jìn)行分析?；诖思僭O(shè)，該方法通過同時(shí)學(xué)習(xí)主題模型分解實(shí)現(xiàn)跨批次歸一化，該分解過程描述了來自所有批次數(shù)據(jù)的歸一化及每個(gè)批次使用調(diào)整因子進(jìn)行校正。其建模思想與潛在狄利克雷模型（latent Dirichle tallocation，LDA）［13］相似。該模型可用圖1表示，其中圖中的陰影圈表示可觀測(cè)變量，非陰影圈表示潛在變量。這里，tm表示基因表達(dá)水平（拷貝數(shù)），θ表示基因表達(dá)水平不同的概率分布是批次和特定基因的調(diào)整因子，為基因g在離散化表達(dá)水平z時(shí)的概率，即＝P（g＝z）。

圖1 PLIDA模型的圖解模型

上述模型參數(shù)可以用模型（9）進(jìn)行估計(jì)，即對(duì)于第i個(gè)批次的第j個(gè)樣品，在離散化基因表達(dá)水平z（z∈｛1，2，…，K｝）情況下有概率模型:

在該等式中，左側(cè)表示所有基因增加1個(gè)拷貝數(shù)屬于基因g的概率，右側(cè)中的fig是批次和特定基因的調(diào)整因子，該因子調(diào)整在離散化基因表達(dá)水平πz下基因g拷貝數(shù)的概率。假定模型中的先驗(yàn)概率關(guān)系如下:

其中c＝1/a。通過塊共軛梯度下降的方法使得后驗(yàn)對(duì)數(shù)似然比最大化，從而估計(jì)出模型參數(shù)｛｝，…，｛｝，｛θ1｝，…，｛θK｝，｛α｝和｝，…，｛，最后得到調(diào)整后的基因表達(dá)值

5.基于比值的方法

基于比值的方法（ratio based method）［5］以各自批次的單個(gè)（或一組）參考樣品為基礎(chǔ)調(diào)節(jié)樣品的基因表達(dá)值。該方法假定參考樣品中的基因與剩下樣品的具有相同的批次效應(yīng)，因此通過減去對(duì)應(yīng)批次的參考樣品的每個(gè)基因平均值來移除批次效應(yīng)。這里把第i個(gè)批次中第j個(gè)參考樣品中的第g個(gè)基因表達(dá)值用表示，第i個(gè)批次的參考樣品數(shù)量為ni。如果參考樣品數(shù)ni≥1，則可以使用參考樣品中的算數(shù)均數(shù)或者幾何均數(shù)的基因表達(dá)值。兩種方法如下:

（1）基于算術(shù)均數(shù)比值方法（Ratio-A）

（2）基于幾何均數(shù)比值方法（Ratio-G）

其他的統(tǒng)計(jì)算法還包括均值中心化［14］、基因標(biāo)準(zhǔn)化［15］、分位數(shù)離散化、中位數(shù)秩分?jǐn)?shù)［16］等，具體可以參閱相關(guān)文獻(xiàn)。

批次效應(yīng)移除的評(píng)價(jià)方法

目前文獻(xiàn)中將移除批次效應(yīng)評(píng)價(jià)方法歸納為兩種:可視化方法及定量度量方法［17］?？梢暬椒僭O(shè)不同批次中感興趣的生物學(xué)變量的樣本分布相同，該假設(shè)對(duì)于數(shù)據(jù)整合過程至關(guān)重要這也意味著，合并的數(shù)據(jù)集在病例及對(duì)照組樣本上應(yīng)當(dāng)包含相似的比例，由于估計(jì)的參數(shù)具有不可比性，整個(gè)分析可能傾向于誤導(dǎo)性的結(jié)果。基于這個(gè)假設(shè)，若無批次效應(yīng)的存在，不同研究間相同基因的表達(dá)水平應(yīng)當(dāng)有相似的分布，如果相同的批次聚在一起則表明批次效應(yīng)的存在。根據(jù)上述兩種理想情況，可視化方法主要分為基因水平及綜合水平。基因水平的可視化方法從單個(gè)基因水平上呈現(xiàn)出批次效應(yīng)的局部可視化，常見方式為基因表達(dá)數(shù)據(jù)的箱式圖和概率密度曲線圖。而綜合可視化工具提供樣本水平批次效應(yīng)存在的圖形，常見方式為樹狀圖，主成分圖及相對(duì)對(duì)數(shù)表達(dá)圖［18］。

盡管可視化手段只對(duì)批次效應(yīng)移除算法的有效性進(jìn)行粗略估計(jì)，但是能夠?qū)λ惴ㄒ瞥涡?yīng)結(jié)果進(jìn)行快速檢查，因此這種方法是實(shí)際中最常見和最直接的方法。如果想要準(zhǔn)確評(píng)價(jià)批次效應(yīng)移除的效果，應(yīng)當(dāng)使用定量度量的方法，這里簡(jiǎn)要介紹3種:

1.Mattews相關(guān)系數(shù)

Mattews相關(guān)系數(shù)（Matthews correlation coefficient）是從診斷目的出發(fā)，通過跨批次的預(yù)測(cè)準(zhǔn)確性評(píng)價(jià)批次效應(yīng)移除效果的評(píng)價(jià)指標(biāo)。其基本思想是，如果已經(jīng)移除批次效應(yīng)，則分類器的錯(cuò)誤率理應(yīng)下降。需要注意的是，預(yù)測(cè)準(zhǔn)確性是高度依賴于分類比率的指標(biāo)，當(dāng)終點(diǎn)指標(biāo)在研究中高度不平衡，結(jié)果可能使人誤解。為此可以使用 Mattews相關(guān)系數(shù)（MCC）［5］，即

其中TP和FP表示真陽性和假陽性，TN和FN表示真陰性和假陰性。－1≤MCC≤1，1表示完全預(yù)測(cè)正確，0表示隨機(jī)預(yù)測(cè)，－1表示預(yù)測(cè)結(jié)果完全相反。

2.混合分?jǐn)?shù)

混合分?jǐn)?shù)（mixture score）［19］的原理是，計(jì)算數(shù)據(jù)A的樣品在數(shù)據(jù)B的k-近鄰范圍內(nèi)的數(shù)目與數(shù)據(jù)B中樣品數(shù)目k的比值。指標(biāo)計(jì)算如下:

其中x是數(shù)據(jù)A的樣品在數(shù)據(jù)B的k近鄰范圍內(nèi)的數(shù)目，0≤Mixturescore≤1。實(shí)際中，分?jǐn)?shù)愈接近0.5說明兩個(gè)批次的數(shù)據(jù)整合愈好，而分?jǐn)?shù)接近0或1則意味兩個(gè)數(shù)據(jù)集整合不好。

3.校正蘭德指數(shù)和信息變異

Saman指出目前定量評(píng)價(jià)方法通常不能同時(shí)滿足以下兩個(gè)標(biāo)準(zhǔn):①技術(shù)因素引起的混雜信號(hào)應(yīng)當(dāng)從數(shù)據(jù)中移除（批次效應(yīng)移除）；②感興趣的生物學(xué)信號(hào)在批次校正算法后應(yīng)該保持不變（生物學(xué)信號(hào)保留）。為此可以使用校正蘭德指數(shù)和信息變異［20］。

（1）蘭德指數(shù)

蘭德指數(shù)（Rand index，RI）基本思想是假設(shè)k個(gè)批次數(shù)據(jù)中共有n個(gè)樣品，用某種聚類算法將樣品分成k類，計(jì)算批次和聚類中樣品配對(duì)組合完全一致的數(shù)目n11與完全不一致的數(shù)目n00之和與樣品中隨機(jī)配對(duì)組合數(shù)的比值［21］，即

當(dāng)樣品隨機(jī)分配到聚類的組別中，RI取值接近0。但是，當(dāng)批次與樣品是相同數(shù)量級(jí)排序時(shí)，RI無法為0。為此，可以使用校正蘭德指數(shù)（CRI），即

E（RI）和max（RI）有專門的計(jì)算公式。當(dāng)對(duì)象隨機(jī)分配到聚類的組別中，CRI取值接近0，即批次和聚類結(jié)果無關(guān)。

（2）信息變異

信息變異（VI）與互信息概念有關(guān)。假設(shè)有n個(gè)樣品的 k個(gè)批次 A＝｛A1，A2，…，Ak｝，聚類為 k個(gè)組別B＝｛B1，B2，…，Bk｝，則 VI指標(biāo)用下式計(jì)算，即

其中H（·）為熵，I（·，·）是兩種分類（批次和聚類結(jié)果）的互信息，0≤VI（A，B）≤min（log（n），2log（k）），VI＝0表示兩個(gè)批次完全分開，實(shí)際中越大越好。

結(jié) 語

基因組學(xué)數(shù)據(jù)的快速增長(zhǎng)為腫瘤生物學(xué)、尋找生物標(biāo)志物及藥物治療靶點(diǎn)提供了更好的研究基礎(chǔ)。盡管研究者需要面對(duì)批次效應(yīng)的問題，但近年數(shù)據(jù)整合算法已經(jīng)有了實(shí)質(zhì)性的發(fā)展，從而能夠?yàn)槲覀兲峁└煽康慕y(tǒng)計(jì)處理結(jié)果。Jason Rudy和FaramarzValafar的研究表明，EB、XPN和DWD在對(duì)數(shù)轉(zhuǎn)化的基因組學(xué)微陣列數(shù)據(jù)中表現(xiàn)最好，PLIDA則要求表達(dá)值在線性范圍內(nèi)變化。在樣本量方面，EB算法可用于單批次小樣本的基因組學(xué)數(shù)據(jù)（檢測(cè)的n＜10），而其他算法通常要求單批次樣本量n≥25。可以預(yù)測(cè)，今后在不同平臺(tái)或批次檢測(cè)的數(shù)據(jù)整合技術(shù)上會(huì)有進(jìn)一步的發(fā)展。

［1］Golub TR，Slonim DK，Tamayo P，et al.Molecular classification of cancer:class discovery and class prediction by gene expression monitoring.Science，1999，286（5439）:531-537.

［2］Bittner M，Meltzer P，Chen Y，et al.Molecular classification of cutaneous malignant melanoma by gene expression profiling.Nature，2000，406（6795）:536-540.

［3］Tinker AV，Boussioutas A，Bow tell DD.The challenges of gene expression m icroarrays for the study of human cancer.Cancer Cell，2006，9（5）:333-339.

［4］Scherer A.Variation，Variability，Batches and Bias in Microarray Experiments:An Introduction.John Wiley＆Sons，Ltd，2009:1-4.

［5］Luo J，Schumacher M，Scherer A，et al.A comparison of batch effect removal methods for enhancement of prediction performance using MAQC-II microarray gene expression data.Pharmacogenomics J，2010，10（4）:278-291.

［6］Leek JT，Scharpf RB，Bravo HC，et al.Tackling the w idespread and critical impact of batch effects in high-throughput data.Nat Rev Genet，2010，11（10）:733-739.

［7］Rudy J，Valafar F.Empirical comparison of cross-platform normalization methods for gene expression data.BMC Bioinformatics，2011，12:467.

［8］Chen C，Grennan K，Badner J，et al.Removing batch effects in analysis of expression microarray data:an evaluation of six batch adjustment methods.PLoSOne，2011，6（2）:e17238.

［9］Johnson WE，Li C，Rabinovic A.Adjusting batch effects in microarray expression data using empirical Bayes methods.Biostatistics，2007，8（1）:118-127.

［10］Stein CK，Qu P，Epstein J，et al.Removing batch effects from purified plasma cell gene expression microarrays with modified Com Bat.BMC Bioinformatics，2015，16（1）:63.

［11］Shabalin AA，Tjelmel and H，F(xiàn)an C，et al.Merging two gene-expression studies via cross-platform normalization.Bioinformatics，2008，24（9）:1154-1160.

［12］Deshwar AG，Morris Q.PLIDA:cross-platform gene expression normalization using perturbed topic models.Bioinformatics，2014，30（7）:956-961.

［13］Heinrich G.Parameter estimation for text analysis.Technical Report，2004.

［14］Sims AH，Smethurst GJ，Hey Y，et al.The removal of multiplicative，systematic bias allows integration of breast cancer gene expression datasets-improving meta-analysis and prediction of prognosis.BMC Med Genomics，2008，1:42.

［15］Li C，Wong WH.Model-based analysis of oligonucleotide arrays:expression index computation and outlier detection.Proc Natl Acad Sci USA，2001，98（1）:31-36.

［16］Jiang H，Deng Y，Chen HS，et al.Joint analysis of two m icroarray gene-expression data sets to select lung adenocarcinoma marker genes.BMC Bioinformatics，2004，5:81.

［17］Lazar C，Meganck S，Taminau J，et al.Batch effect removal methods form icroarray gene expression data integration:a survey.Brief Bioinform，2013，14（4）:469-490.

［18］Gagnon-Bartsch JA，Speed TP.Using control genes to correct for unwanted variation in microarray data.Biostatistics，2012，13（3）:539-552.

［19］Kim KY，Kim SH，Ki DH，et al.An attempt for combining microarray data sets by adjusting gene expressions.Cancer Res Treat，2007，39（2）:74-81.

［20］Vaisipour S.Detecting，correcting，and preventing the batch effects in multi-site data，with a focus on gene expression Microarrays.University of Alberta，2014.

［21］Vinh NX，Epps J，Bailey J.Information Theoretic Measures for Clusterings Comparison:Is a Correction for Chance Necessary？Proceedings of the 26th International Conference on Machine Learning（ICML-09），2009.

國(guó)家自然科學(xué)基金資助（81473072）

（責(zé)任編輯:郭海強(qiáng)）