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

1)徐州醫(yī)學(xué)院公共衛(wèi)生學(xué)院流行病學(xué)與衛(wèi)生統(tǒng)計學(xué)教研室 徐州 221002 2)南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)與衛(wèi)生統(tǒng)計學(xué)教研室 南京 210029

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

1 對象與方法

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

表1 前列腺癌微陣列數(shù)據(jù)形式

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

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

圖1 t值的直方圖和對應(yīng)的零分布

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

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

1.4統(tǒng)計學(xué)處理t檢驗和錯誤發(fā)現(xiàn)率分析在R2.13.0中完成[9]。

2 結(jié)果

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

圖2 BH控制的示意圖

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

3 討論

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

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

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

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

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

[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，F(xiàn)ebbo PG，Ross K，et al.Gene expression correlates of clinical prostate cancer behavior[J].Cancer Cell，2002，1(2)：203

[6]曾平，王婷.貝葉斯錯誤發(fā)現(xiàn)率[J].山東大學(xué)學(xué)報：醫(yī)學(xué)版，2012，50(3)：120

[7]王婷，曾平，黃水平，等.錯誤發(fā)現(xiàn)率及其擴(kuò)展和應(yīng)用[J].重慶醫(yī)科大學(xué)學(xué)報，2011，36(12)：38

[8]王婷，曾平，黃水平，等.錯誤發(fā)現(xiàn)率的經(jīng)驗估計和應(yīng)用[J].鄭州大學(xué)學(xué)報：醫(yī)學(xué)版，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在假設(shè)檢驗中的應(yīng)用[J].數(shù)理統(tǒng)計與管理，2006，25(5)：616

[11]Hastie T，Tibshirani R，F(xiàn)riedman 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