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

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

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

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

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

統(tǒng)計分析算法

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

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

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

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

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

算法分為三步:

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

再通過下式把標準化數(shù)據(jù)Zijg計算出來，即

（2）估計批次效應參數(shù)。假設 Zijg～N（γig，δig），采用參數(shù)先驗或者非參數(shù)先驗方法計算出批次效應估計值和。其中參數(shù)先驗方法要求 γig～N（γi，）及～InverseGamma（λi，θi）。

2.改進ComBat方法

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

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

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

其中符號的意義同式（4）。

3.跨平臺歸一化

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

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

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

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

4.平臺獨立的LDA模型

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

圖1 PLIDA模型的圖解模型

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

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

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

5.基于比值的方法

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

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

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

其他的統(tǒng)計算法還包括均值中心化［14］、基因標準化［15］、分位數(shù)離散化、中位數(shù)秩分數(shù)［16］等，具體可以參閱相關文獻。

批次效應移除的評價方法

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

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

1.Mattews相關系數(shù)

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

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

2.混合分數(shù)

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

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

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

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

（1）蘭德指數(shù)

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

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

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

（2）信息變異

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

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

結(jié) 語

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

國家自然科學基金資助（81473072）

（責任編輯:郭海強）