微陣列數據分析和錯誤發現率*

王 婷，曾 平，黃水平，趙華碩

1)徐州醫學院公共衛生學院流行病學與衛生統計學教研室 徐州 221002 2)南京醫科大學公共衛生學院流行病學與衛生統計學教研室 南京 210029

現代生物和醫學技術的發展使得人們能夠收集到大量的數據，微陣列技術是其中的著名代表，為人們提供了一種從沒有過的醫學實踐方式。微陣列技術用含有成千上萬種的DNA或者蛋白質序列的微小玻璃芯片取代傳統生物醫學中的凝膠、濾器和純化柱，使得科學家們能夠在基因組規模上對基因表達水平進行快速和定量的檢測，由此產生的大規模數據也給統計領域帶來了前所未有的機遇和挑戰[1]。在微陣列實驗中研究者需要同時檢驗數以千計的基因表達水平是否與不同條件狀況之間存在關聯，由此涉及多重假設檢驗的問題。總體錯誤率(family wise error rate，FWER)是傳統多重假設中常用的錯誤控制指標，定義為至少犯一次I型錯誤的概率[2]，但是總體錯誤率控制過于嚴格并不適于微陣列數據分析。有學者[3-4]提出的錯誤發現率(false discovery rate，FDR)很好地適應了高維數據多重檢驗的錯誤控制，越來越多的應用于微陣列數據分析。作者以前列腺癌的微陣列數據為例，介紹了基于錯誤發現率控制的微陣列數據多重比較，期望為此類大規模數據分析提供恰當的方法選擇。

1 對象與方法

1.1研究對象微陣列數據通常表示為m×n的矩陣形式，m表示基因數，一般以千計，n表示生物樣本，通常只有幾個或者幾十個。表1給出了一個關于前列腺癌的微陣列數據，共收集了50(n1)個正常人和52(n2)個前列腺癌患者6 033個基因的表達數據[5]，以6 033×102矩陣排列，即m=6 033，n=102。

表1 前列腺癌微陣列數據形式

1.2t檢驗和總體錯誤率研究者的主要目的是，哪些基因在正常人和前列腺癌患者中的表達水平不同？需要識別出這些基因以進一步分析。關鍵的問題就在于如何在成千上萬的候選基因中找出特征基因，這一過程稱為特征選擇。兩個獨立樣本的t檢驗可用來檢驗兩組基因的平均表達水平是否存在差異。在零假設條件下，統計量ti(i=1，…，6 033)服從自由度為100=n-2的t分布，對應的pi為Prob(|T|≥|ti|)。表1中的最后兩列給出了每個基因的t值和P值，圖1給出了t值的直方圖和零分布曲線，0附近的t值和零分布吻合，但處于兩側的t值比零分布更加分散，暗示這些位置的基因表達可能存在差異。按照α=0.05的檢驗水準，有478個pi≤0.05，顯然這個結果值得懷疑，因為即使所有 6 033 個基因表達沒有差異，按照0.05的水準大約也會獲得302=6 033×0.05個基因，因此需要對多重假設檢驗所帶來的錯誤采取控制措施。常用的總體錯誤率控制程序包括Bonferroni、Holm、Hochberg和Hommel等方法[2]，但對前列腺癌數據都只能得到3個基因。總體錯誤率控制程序沒有發現更多基因的原因在于對大規模的微陣列數據而言總體錯誤率過于嚴格，例如，Bonferroni的檢驗水準為8.3×10-6=0.05/6 033，只有極小的pi才能被認為有差異。

1.3錯誤發現率Benjamini等[3]將錯誤發現率定義為在所有拒絕H0的次數中屬于錯誤拒絕次數的期望。表2給出了可能出現的檢驗結果[6-8]，m0和m1分別表示在m次多重檢驗中真實H0和非真實H0的個數，V表示在所有R次拒絕H0的決定中拒絕了原本真實H0的次數。除m、R和W外，表1中的其他量均為沒有觀察到的隨機變量。該研究中m即為總的基因數6 033。

圖1 t值的直方圖和對應的零分布

實際情況不拒絕H0拒絕H0合計H0=0UVm0H0=1TSm1合計WRm

錯誤發現率定義為：FDR=E[V/(R∨1)]=E(V/R|R>0)P(R>0)。R∨1表示當R=0時FDR=0。設預先的總體錯誤率控制水平為q，q在0～1之間，Benjamini和Hochberg(BH)給出了一個基于p值的逐步向下控制程序：①將p值排序：p(1)≤p(2)…≤p(m)，H0(1)，H0(2)，…，H0(m)為對應的零假設。②定義k滿足：p(k)≤kq/m，如k存在則拒絕H0(1)，H0(2)，…，H0(k-1)，H0(k)，否則不拒絕任何一個H0。Benjamini和Hochberg證明了在BH控制下FDR=π0q≤q，π0為無差異表達基因的比例。錯誤發現率允許一定程度的錯誤拒絕，因此在假設檢驗次數很多時效能比總體錯誤率更高。

1.4統計學處理t檢驗和錯誤發現率分析在R2.13.0中完成[9]。

2 結果

2.1錯誤發現率控制選擇q=0.05，按照BH程序得到21個差異表達基因。圖2給出了BH控制的示意圖，直線過原點，斜率為0.05/6 033，參考線以左的p值被認為有差異，此時對應的p值為1.3×10-4，也即，此時的檢驗水準為1.3×10-4。q=0.05的含義在于，平均而言在所有21個差異表達基因中大約有1=21×0.05個基因屬于錯誤識別。如果選擇q=0.10或0.20，分別得到60或106個表達差異基因。同理，q=0.10或0.20的含義為，在60(106)個差異基因中約有6(21)個屬于錯誤發現。

圖2 BH控制的示意圖

2.2錯誤發現率估計BH程序將錯誤發現率控制在一個預先選擇的范圍內，一個自然的問題是，例如以|ti|≥3作為拒絕域得到105個基因，那么對應的錯誤發現率是多少呢？這一過程和上述的BH控制相反，稱為錯誤發現率估計。如果所有6 033個基因表達是無差異的，那么理論上分布在區間(-∞，-3]∪[3，∞)的t值個數約為21=6 033×2×F(-3，100)，F(x，100)表示自由度為100的t分布累計概率函數，F(-3，100)=1.7×10-3，這里取π0最保守的估計值1，當π0>0.90時這種做法不會對結果產生多大的影響，對應的錯誤發現率約為0.20=21/105。

3 討論

微陣列分析作為一種探索性的分析策略，主要目的是為后續的基因研究提供候選基因，基因數一般介于幾百到幾萬之間，因此相對于控制至少出現一次假陽性的概率，研究者更關心的是能否盡量多地識別出差異表達的基因，能夠容忍和允許在拒絕過程中發生少量的錯誤識別，只要錯誤發現的比例足夠的小即可。這就決定了在微陣列數據中需要在錯誤發現和總的拒絕次數之間尋找一種平衡，即在檢驗出盡可能多的候選基因的同時將錯誤發現控制在一個可以接受的范圍內。兩者的區別在于，總體錯誤率被定義為一個概率，錯誤發現率被定義為一個期望值。此外總體錯誤率和錯誤發現率的含義也截然不同，前者表示在H0成立的條件下錯誤拒絕H0的概率，后者表示在已經拒絕H0的次數中錯誤拒絕的比例期望。錯誤發現率為微陣列數據的多重比較提供了十分恰當的錯誤控制指標，對大規模數據分析具有十分重要的意義。

錯誤發現率和假陽性率(false positive rate，FPR)之間有著本質的區別：例如，如果前列腺癌數據中所有6 033個基因的表達水平原本是沒有差異的，5%的FPR表示將大約有302=5%×6 033的基因會被錯誤地認為差異表達；5%的FDR則表示，如果在所有6 033個多重假設檢驗中，如果有100次拒絕H0(即R=100)，那么大約有5=5%×100基因是被錯誤識別的。錯誤發現率突破了傳統多重檢驗的范疇，在高維數據分析中扮演著十分重要的角色，錯誤發現率不但具有十分現實的應用性，還具有潛在的理論意義。雖然錯誤發現率是在頻率統計下發展起來的，但同時也具有可解釋的貝葉斯和經驗貝葉斯含義，曾平等[6]闡述了錯誤發現率的貝葉斯解釋和經驗估計。王婷等[5]研究了基于錯誤發現率發展的q值和局部錯誤發現率。

作者對2組微陣列實驗選擇了t檢驗執行統計分析，并假設統計量服從t分布從而獲得P值。當不能確定t統計量的零分布形式時，可采用基于重抽樣技術得到P值[10]。

作者介紹了基于BH程序的錯誤發現率控制和估計，兩者得到了相近的結果，Efron[4]證明了錯誤發現率估計和控制的等價性。雖然有其他控制程序被提出，但是BH程序的一個優點在于簡單易行，并且對相關具有一定的穩健性。錯誤發現率估計和控制相反，即估計一個具體拒絕域的錯誤發現率，這對研究者根據得到有差異表達基因的個數和對應的錯誤發現率大小做出決策有重要意義，作者選擇了對稱的拒絕域，生物學家和其他專家的信息可作為選擇合適拒絕域的重要參考，如選擇非對稱的拒絕域[11]。

此外，無差異表達基因的比例π0往往也是研究者關心的參數，在大規模多重比較中具有十分重要的作用。該文中估計錯誤發現率時采用了π0最保守的估計值，即π0=1，Benjamini和Hochberg的BH控制程序中同樣采用π0=1，這一選擇在大多時候并不會帶來嚴重的影響，因為很多微陣列數據的π0一般都在0.9以上。有學者[12]報道了更加精細的方法用來估計π0，以提高和改進錯誤發現率的控制和估計。

[1]Johnstone IM，Titterington DM.Statistical challenges of high-dimensional data[J].Philos Transact A Math Phys Eng Sci，2009，367(1906)：4237

[2]Bretz F，Hothorn T，Westfall P.Multiple comparisons using R[M].London：Chapman & Hall，2010：11

[3]Benjamini Y，Hochberg Y.Controlling the false discovery rate：A practical and powerful approach to multiple testing[J].J Royal Statist Soc：Series B，1995，57(1)：28

[4]Efron B.Large-scale inference：empirical Bayes methods for estimation，testing，and prediction[M].New York：Cambridge University Press，2010：46

[5]Singh D，Febbo PG，Ross K，et al.Gene expression correlates of clinical prostate cancer behavior[J].Cancer Cell，2002，1(2)：203

[6]曾平，王婷.貝葉斯錯誤發現率[J].山東大學學報：醫學版，2012，50(3)：120

[7]王婷，曾平，黃水平，等.錯誤發現率及其擴展和應用[J].重慶醫科大學學報，2011，36(12)：38

[8]王婷，曾平，黃水平，等.錯誤發現率的經驗估計和應用[J].鄭州大學學報：醫學版，2012，47(5)：636

[9]R Development Core Team.R：A language and environment for statistical computing[EB/OL].R Foundation for Statistical Computing，Vienna，Austria，2007.URL http：//www.R-project.org.

[10]荀鵬程，趙楊，易洪剛，等.Permutation Test在假設檢驗中的應用[J].數理統計與管理，2006，25(5)：616

[11]Hastie T，Tibshirani R，Friedman J.The elements of statistical learning：data mining，inference，and prediction，second edition[M].New York：Springer-Verlag，2009.

[12]Benjamini Y.Discovering the false discovery rate[J].J Royal Statist Soc：Series B，2010，72(4)：405