腫瘤微陣列數(shù)據(jù)統(tǒng)計分析概述

張楊，張威，曹文君,2，李運明,3，陳長生

1.第四軍醫(yī)大學 軍事預防醫(yī)學院衛(wèi)生統(tǒng)計學教研室，陜西 西安 710032;2.長治醫(yī)學院 基礎部，山西 長治 046000；3.成都軍區(qū)總醫(yī)院 神經(jīng)外科，四川 成都 610083

隨著生物技術和統(tǒng)計方法的發(fā)展和改進，微陣列技術可以在一次試驗中對整個基因組進行分析，已廣泛應用于分析大規(guī)模的mRNA 表達數(shù)據(jù)。目前有4 類生物芯片平臺被廣泛應用，即Affymetrix GeneChip 芯片、寡核苷酸探針芯片、cDNA 芯片和商業(yè)化探針芯片。在腫瘤研究領域，微陣列技術對腫瘤的診斷和分型、治療和預后，以及探討腫瘤發(fā)生的分子機制和發(fā)展都有非常重要的作用?；虮磉_譜數(shù)據(jù)分析，對腫瘤患者的個性化治療和腫瘤的分子分型發(fā)揮著越來越重要的作用。

由于微陣列技術所得到的基因表達譜數(shù)據(jù)具有高維（成千上萬個基因）和樣本量小的特點，如何挖掘和解釋其中所蘊含的海量基因信息，深層次研究基因功能，選擇適當?shù)慕y(tǒng)計學方法對于芯片數(shù)據(jù)的處理至關重要。微陣列基因表達譜數(shù)據(jù)信息的提取及其統(tǒng)計分析方法的研究，已成為生物與醫(yī)學統(tǒng)計學領域中富有挑戰(zhàn)性的重要課題。

1 腫瘤研究中微陣列技術的流程

微陣列試驗的基本過程如圖1 所示?；虮磉_譜通過計算機軟件掃描圖像提取得到的原始數(shù)據(jù)首先通過標準化方法過濾掉那些低質量的探針數(shù)據(jù)，即在所有樣本中都表現(xiàn)出較低的信號強度或變化幅度，對感興趣的表型或條件不可能有貢獻[1]。芯片數(shù)據(jù)只有通過標準化處理后，才能進行下游的統(tǒng)計分析，如篩選差異表達的基因或腫瘤分型。微陣列技術通常被作為篩選工具，生物驗證和解釋應對選定的基因做進一步研究。

圖1 微陣列試驗流程圖，其中標星號的流程為統(tǒng)計分析部分

2 cRNA微陣列數(shù)據(jù)的標準化和預處理

微陣列數(shù)據(jù)的質量對下游的統(tǒng)計分析至關重要，包括RNA 的質量、探針標記、雜交條件、洗板，以及在掃描當中的信號強度和背景干擾。尤其在低豐度表達RNA 分子的研究中容易受背景系統(tǒng)的影響而導致系統(tǒng)偏差。這些偏移可導致基因表達數(shù)據(jù)研究的錯誤結論，即假陽性和假陰性的預測。但這些變異和微陣列數(shù)據(jù)的系統(tǒng)性偏差可通過科學的重復和歸一化進行控制。

在實驗過程中，當我們從實驗組和對照組樣本中獲得基因表達譜后，可以通過計算每個基因在Cy3和Cy5 通道的熒光染料強度的平均對數(shù)比值而得到該基因的表達強度[1]。通過比較在不同處理組中該基因表達強度的差異，對其進行進一步的研究。在得到該表達強度的數(shù)據(jù)后，首先要進行的是將低重復性的探針數(shù)據(jù)過濾去除。該過程可以通過控制編譯系數(shù)（小于特定的閾值）和表達密度（大于特定表達水平）進行。經(jīng)過預處理后，掃描圖像的整體亮度和試驗變化所造成的系統(tǒng)性偏差可以得到有效控制。例如，模塊和染料的影響。此步驟對微陣列數(shù)據(jù)的多重比較和下游的統(tǒng)計分析是必不可少的。這一過程統(tǒng)稱為歸一化。

歸一化的方法很多，其中應用最廣泛的是全局歸一化（global normallization）。該方法的目的是將所有芯片探針均歸一化為具有相同的中位表達強度，這一方法可以很好地矯正模塊數(shù)據(jù)。然而，大量的統(tǒng)計學研究證明模塊數(shù)據(jù)的影響是存在的。針對這一影響，Dudoit 等經(jīng)過不斷研究，發(fā)展了對每個模塊強度數(shù)據(jù)歸一化方法，稱為“LOESS歸一化”[1]?；旧飳W假設是，在一個模塊中上調基因和下調基因的表達量是基本一致的。由于這一歸一化方法的假設條件，因此不適用于定制芯片和經(jīng)處理的特定細胞系表達數(shù)據(jù)的研究。在該方法的基礎上，Tseng等[2]使用“不變基因集”作為看家基因的代理，并僅基于該不變基因集的強度來估計模塊強度。Fan 等[3]在沒有以上方法的生物學假設條件下，利用陣列內重復，提出切片內半線性模型（semilinear in-slide）歸一化方法。該方法基于芯片內約100 次的重復探針，以避免序列特異性和噪音干擾。其基本依據(jù)是，在同一模塊中重復探針的基因表達差異基本反映了除隨機噪音之外的系統(tǒng)誤差，而這些系統(tǒng)誤差可以通過探針配對加以去除。使用切片內半線性模型，通過選取神經(jīng)瘤細胞轉移抑制因子靶向基因，以實時反轉錄PCR 反應被加以驗證，而普通的歸一化方法則容易造成一些基因的缺失。隨后，F(xiàn)an 等[4]通過芯片內重復探針和聯(lián)合其他芯片的信息，顯著擴大了這一方法的適用范圍。除了以上歸一化方法外，其他有用的歸一化方法還包括雙向半線性模型[5]和穩(wěn)?。╮obust）歸一化等[6]。

3 芯片內重復

芯片內重復不僅對歸一化的作用很大，而且對于驗證數(shù)據(jù)歸一化是否正確也非常有用。其基本的思想是，芯片內重復之間的差異是系統(tǒng)偏差除去后的純粹的隨機噪聲。當估計每個單獨基因的噪聲水平并且芯片通過歸一化，那么總的標準化方差大概服從卡方分布。這對檢測數(shù)據(jù)是否已歸一化提供了一個簡單而有用的診斷測試方法。這一檢驗統(tǒng)計量也可作為一個給定陣列選擇歸一化方法的標準——最小的檢驗統(tǒng)計量（最一致的重復）是最優(yōu)的選擇?？梢酝ㄟ^過濾方法和經(jīng)驗貝葉斯的方法估計相應的方差的差異[7]。芯片內重復也被用于改進基因集方差估計的精度，從而提高推斷方法的設計來識別差異表達基因[8]。

4 篩選差異表達基因

微陣列實驗的主要目的就是篩選實驗組和對照組或者更復雜的比較組間的差異表達基因[9]。選擇恰當?shù)慕y(tǒng)計學方法在此過程中至關重要。首先，是選擇恰當?shù)臋z驗統(tǒng)計量，通常是經(jīng)過修正的檢驗統(tǒng)計量，如在微陣列顯著性分析中，修正過的單樣本和兩樣本t檢驗、方差分析或經(jīng)典貝葉斯方法[10]。由于微陣列數(shù)據(jù)的一個顯著特點是一般只有一小部分的基因是差異表達的，因此，根據(jù)不同芯片本身的特點選擇合適的檢驗統(tǒng)計量，以提高芯片數(shù)據(jù)處理的靈敏度和特異度。

在選擇好檢驗統(tǒng)計量之后，下一個步驟是計算檢驗統(tǒng)計量并由此得出顯著性P值。芯片數(shù)據(jù)常常要同時處理成千上萬個基因，相應的P值的計算也較為復雜。例如，當樣本量較大，且服從正態(tài)分布時，可采用student't 分布以計算P值；而在小樣本非正態(tài)分布時，permutation 或bootsrapping 方法是較為合適的選擇。

判斷出P值之后，接下來將是篩選差異表達的基因。由于微陣列試驗同時檢驗成千上萬個基因，很可能出現(xiàn)較高的陽性結果錯誤發(fā)現(xiàn)率。因此，控制錯誤發(fā)現(xiàn)率對結果的生物學解釋至關重要[11-12]。目前已發(fā)表了很多關于控制錯誤發(fā)現(xiàn)率方法的論文。Storey等[11]和Dudoit等[12]對關于基因集中控制錯誤發(fā)現(xiàn)率的方法進行了綜述。這些方法要求非常精確地P值，通常在10-6數(shù)量級。此外，有些方法可以通過判斷檢驗統(tǒng)計量是否超過某一界值或相應的P值（小于閥值）來篩選差異表達基因，從而控制錯誤發(fā)現(xiàn)率[2,10]。例如，假設在15 000 個基因中有100 個基因相應的P值小于0.001，則期望的被錯誤發(fā)現(xiàn)的基因數(shù)不超過0.001×15 000=15，那么我們可以估計錯誤發(fā)現(xiàn)率為15/100=15%[2]。

5 微陣列技術在腫瘤分類和聚類中的應用

微陣列技術腫瘤臨床研究的一個重要應用是發(fā)現(xiàn)腫瘤生物學標志物和對腫瘤進行病理分類。Inamura 等[13]在其研究中對這一領域進行了較為詳細的闡述。他們在肺癌和正常肺組織標本的研究中，通過分層聚類和非負矩陣因子化的方法將肺鱗狀細胞癌分為2個不同的亞型，2個亞型具有完全不同的分子特征和臨床結局。

分類也被稱為有監(jiān)督學習。許多統(tǒng)計分析方法被應用于聚類和分類，這些方法包括決策樹分類方法、線性鑒別分析、支持矢量機法，以及神經(jīng)網(wǎng)絡特征分析[14]。使用歸一化芯片數(shù)據(jù)作為輸入向量，可以建立分類規(guī)則。Svrakic 等[15]對基因集中所應用的聚類方法做了全面綜述。

通常我們在對腫瘤的分類研究中，希望在篩選基因或腫瘤標志物的過程中得到具有較高判別效能和低誤判率的差異表達基因。這不僅提高了基因在腫瘤中功能的理解，而且降低了錯誤分類的幾率?；虮磉_譜收縮因子方法在腫瘤分類研究中被作為一種重要的分類方法[16]。統(tǒng)計變量選擇法，如誤判的出發(fā)散度，也可被應用于篩選重要的差異表達基因和腫瘤標志物[17]。

聚類也稱為無監(jiān)督學習方法，常常用于對基因表達譜中具有相似表達特征的基因進行歸類[15]。這有利于我們發(fā)現(xiàn)共表達或表達特性相似的基因群。聚類算法的一個重要步驟是如何在輸入空間中定義合適矩陣指標，如分層和K最近鄰分類法[15-16]。這些輸入向量既可以是不相關的基因集中具有相似表達特性的基因，也可以是不同樣本的同一表達基因。在聚類過程中，常用的計算指標包括歐氏距離和Pearson 相關系數(shù)。相似表達基因通常通過系統(tǒng)樹圖或彩色編碼表示。

6 時間序列和調控網(wǎng)絡

為了監(jiān)測基因的隊列表達模式，基因隨著疾病的進展時間，或在不同治療中的表達情況，我們可以在不同的時間點取樣獲得基因的表達數(shù)據(jù)。統(tǒng)計學上一個重要的問題是，在某個特定時間點經(jīng)過處理后的基因表達是否有差異。Hotelling T2檢驗可用來驗證隨著時間的推移，基因的表達譜是否發(fā)生變化。隨著時間的進程，某些基因的表達上調或下調，某些基因的表達保持不變，可以發(fā)現(xiàn)隨著時間及條件的變化表現(xiàn)出不同表達水平的標識基因。基因表達模式的時間序列也有利于理解與疾病相關的表達通路及其功能。單樣本t檢驗可用于評估在每個特定的時間點基因的表達是否上調和下調或保持不變。那么，基因表達模式的時間進程可通過一類分類技術來加以分析，這對理解基因間的調控過程和生物通路提供了非常有用的工具[15]。Schulte 等[18]研究了不同的基因表達模式，包括即早基因、“延遲”基因和效應基因在神經(jīng)母細胞瘤中TrkA和TrkB 受體中的表達情況，進而發(fā)現(xiàn)在誘導即早基因和下游的靶點調控中的分子機理。

目前，在時間序列微陣列研究中已發(fā)展了許多方法，這些方法都是針對提取時間過程中的差異表達基因。Storey等[19]利用基于spline的方法發(fā)現(xiàn)了在時間過程中基因表達的改變。Yuan 等[20-21]分別在2006和2008 年用基于隱馬爾科夫算法的模型分析了多種生物條件下的時間序列微陣列數(shù)據(jù)。Tai和Speed[22]提出了一種多變量經(jīng)驗貝葉斯的統(tǒng)計方法鑒別差異表達基因。Ma 等[23]在2009 年通過功能性ANOVA 混合效應模型將時間序列基因表達值分類，并且鑒別差異表達基因。Zhou 等[24]在2010 年根據(jù)時間序列基因表達數(shù)據(jù)發(fā)展了一種對生存結果的預測模型。Tibshirani 等[25]在2013 年提出了一種根據(jù)時間序列基因表達譜對病人分類的方法，這種方法的提出為腫瘤的個性化治療提供了新的手段。

基于通路和GSEA 的方法在功能基因組學研究中已經(jīng)發(fā)展了十余年。由于不完整的信息和通路數(shù)據(jù)注釋不佳，研究人員開始結合基因集富集分析方法和基于網(wǎng)絡模塊的方法，以鑒別較大幅度的分子機制。第三代基因表達譜分析方法（包括基因集/通路/網(wǎng)絡分析）可被定義為一個以知識為導向的數(shù)據(jù)驅動的方法，這不僅是基于先驗基因集的知識，而且利用了基因集內部或基因集之間的通路/網(wǎng)絡的拓撲結構。2007 年，Vidal[26]的研究小組在哈佛大學利用各種生物信息學數(shù)據(jù)集對乳腺癌的易感性構建蛋白質相互作用網(wǎng)絡，并確定HMMR 為新的疾病易感位點。隨后，TreyIdeker[27]的研究小組在加州大學圣地亞哥分校綜合蛋白質網(wǎng)絡和基因表達數(shù)據(jù)，以提高乳腺癌患者轉移形成的預測。這2 項研究是一個新的里程碑，標志著網(wǎng)絡和通路的激動人心的開始，雖然容易出錯，而且不完整，但可作為一種導向，引導今后的微陣列數(shù)據(jù)分析。

基因芯片技術已被廣泛應用于遺傳變異，基因網(wǎng)絡、調控過程中的相互作用，以及生物通路等方面的研究，已成為理解基因的相互作用、協(xié)同、網(wǎng)絡調控等的有力工具[28]。

7 結語

微陣列技術正深入到人類腫瘤疾病研究的各個方面。與其他研究方法相比，該技術更關注腫瘤在不同條件下基因表達的變化。可以通過微陣列技術對基因組進行分析來確定新的潛在的治療途徑或研發(fā)新的診斷試劑，即所謂的生物標志物發(fā)現(xiàn)研究[29]。

由于腫瘤受多重因素的影響，因此我們所獲得的差異表達基因，哪怕是一個簡單的比較試驗，也可能受到其他信號的干擾。癌癥樣本的來源各異，包括直接手術活檢的樣本、很有限的針刺活檢樣本、尸檢樣本、特定癌癥所建立的細胞系，甚至是石蠟固定樣本的切片。當比較腫瘤樣本及其正常對照時，必須確保樣本的匹配度，如乳腺癌樣本必須與相對應的正常乳腺細胞相比較。但這可能很難做到，如活檢樣本不可能是均質的，含有多種細胞類型，是正常和惡性細胞或不同階段腫瘤細胞的混合物。此外，腫瘤細胞的正常對照也可能并不確定，因此需要分析多個不同的正常樣本。同樣，不同患者樣本的遺傳差異也可能會影響結果，需要增加足夠的對照以減少這些因素的干擾。顯然，正常樣本和腫瘤樣本的匹配程度也取決于實驗的目標。如探索性分析實驗的目的是了解某一系統(tǒng)的基本生物學特征，那么樣本匹配度的要求可能不像生物標志物篩選那么嚴格，因為生物標志物篩選的目標是確定可靠的診斷工具。因此，實驗設計、原始數(shù)據(jù)分析及統(tǒng)計方法的選擇，是腫瘤研究中至為重要的步驟。

將一個基因表達譜或特殊基因表達信號轉換成生物學上可以理解的概念，仍然是一個需要很大努力和充滿挑戰(zhàn)的任務。在這方面，對人類基因組及其他模式生物基因組功能越來越多的了解，將為基因表達研究提供大量的補充信息。另外，近年出現(xiàn)的系統(tǒng)生物學，在轉錄組水平的目的是能夠描述支持個體基因表達狀態(tài)的基因調節(jié)網(wǎng)絡，也將及時提供細胞轉錄水平的可預見的詳細圖譜。通過表達譜分析，可以預測疾病狀態(tài)的細胞和組織中受影響的特殊生化途徑和生物學過程。

腫瘤基因表達譜數(shù)據(jù)挖掘不僅對認識腫瘤發(fā)生發(fā)展的機理具有重要意義，而且也會為腫瘤的分子診斷和防治開辟全新的途徑，并有助于腫瘤個性化治療的實現(xiàn)。利用基因表達譜對腫瘤樣本進行準確診斷，構建腫瘤基因調控網(wǎng)絡，是一項具有重要意義的大課題。

[1]Dudoit S，Yang Y,Callow M J,et al.Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments[J].Stat Sin,2002,12:111-139.

[2]Tseng G C,Oh M K,Rohlin L,et al.Issues in cDNA microarray analysis:quality filtering,channel normalization,models of variations and assessmentof gene effects[J].Nucleic Acids Res,2001,29(12):2549-2557.

[3]Fan J,Tam P,Vande Woude G,et al.Normalization and analysis of cDNA micro-arrays using within-array replications applied to neuroblastoma cell response to a cytokine[J].Proc Natl Acad Sci USA,2004,101(5):1135-1140.

[4]Fan J,Peng H,Huang T.Semilinear high-dimensional model for normalization of microarray data:a theoretical analysis and partial consistency[J].J Am Stat Assoc,2005,100(471):781-813.

[5]Huang J,Wang D,Zhang C H.A Two-way semi-linear model for normalization and analysis of cDNA microarray data[J].J Am Stat Assoc,2005,100:814-829.

[6]Ma S,Kosorok M R,Huang J,et al.Robust semiparametric cDNA microarray normalization and significance analysis[J].Biometrics,2006,62(2):555-561.

[7]Cui X,Hwang J T G,Qiu J,et al.Improved statistical tests for differential gene expression by shrinking variance components estimates[J].Biostatistics,2005,6(1):59-75.

[8]Smyth G K,Michaud J,Scott H S.Use of within-array replicate spots for assessing differential expression in microarray experiments[J].Bioinformatics,2005,21(9):2067-2075.

[9]Kerr M K,Churchill G A.Experimental design for gene expression microarrays[J].Biostatistics,2001,2(2):183-201.

[10]Tusher V G,Tibshirani R,Chu G.Significance analysis of microarrays applied to the ionizing radiation response[J].Proc Natl Acad Sci USA,2001,98(9):5116-5121.

[11]Storey J D,Tibshirani R.Statistical significance for genomewide studies[J].Proc Natl Aca Sci USA,2003,100(16):9440-9445.

[12]Dudoit S,Shaffer J P,Boldrick J C.Multiple hypothesis testing in microarray experiments[J].Stat Sci,2003,18(1):71-103.

[13]Inamura K,Fujiwara T,Hoshida Y,et al.Two subclasses of lung squamous cell carcinoma with different gene expression profiles and prognosis identified by hierarchical clustering and non-negative matrixfactorization[J].Oncogene,2005,24:7105-7113.

[14]Hastie T J,Tibshirani R,Friedman J.The elements of statistical learning:data mining,inference and prediction[M].2nd ed.New York:Springer,2005:83-85.

[15]Svrakic N M,Nesic O,Dasu M R K,et al.Statistical approach to DNA chip analysis[J].Recent Prog Horm Res,2003,58:75-93.

[16]Tibshirani R,Hastie T,Narasimhan B,et al.Diagnosis of multiple cancer types by shrunken centroids of gene expression[J].Proc Natl Acad Sci USA,2002,99(10):6567-6572.

[17]Fan J,Li R.Statistical challenges with high dimensionality:feature selection in knowledge discovery[J].Proc Madrid Intl Congress Math,2006,111:595-622.

[18]Schulte J,Schramm A,Klein-Hitpass L,et al.Microarray analysis reveals differential gene expression patterns and regulation of single target genes contributing to the opposing phenotype of TrkA-and TrkB-expressing neuroblastomas[J].Oncogene,2005,24(1):165-177.

[19]Storey J D,Xiao W J,Tompkins R,et al.Significance analy-sis of time course microarray experiments[J].Proc Natl Acad Sci USA,2005,102(36)：12837-12842.

[20]Yuan M,Kendziorski C.Hidden markov models for microarray time course data under multiple biological conditions[J].J Am Stat Assoc,2006,101(476):1323-1332.

[21]Yuan Y,Li C T,Wilson R.Partial mixture model for tight clustering of gene expression time-course[J].BMC Bioinf,2008,9:287.

[22]Tai Y C,Speed T P.A multivariate empirica bayes statistic for replicated microarray time course data[J].Ann Stat,2006,34(5):2387-2412.

[23]Ma P,Zhong W,Liu J S.Identifying differentially expressed genes in time course microarray data[J].Stat Biosci,2009,1:144-159.

[24]Zhou B,Xu W,Herndon D,et al.Analysis of factorial timecourse microarrays with application to a clinical study of burn injury[J].Proc Natl Acad Sci USA,2010,107(22):9923-9928.

[25]Zhang Y,Tibshirant R,Davis R.Classification of patients from time-course gene expression[J].Biostatistics,2013,14(1):87-98.

[26]Pujana M A,Han J D J,Starita L M,et al.Network modeling links breast cancer susceptibility and centrosome dysfunction[J].Nat Genet,2007,39(11):1338-1349.

[27]Chuang H Y,Lee E,Liu Y T,et al.Network-based classification of breast cancer metastasis[J].Mol Syst Biol,2007,3:140-149.

[28]Akiyoshi T,Kobunai T,Watanabe T.Predicting the response to preoperative radiation or chemoradiation by a microarray analysis of the gene expression profiles in rectal cancer[J].Surgery Today,2012,42(8):713-719.

[29]Qian Z,Qingshan C,Chun J,et al.High expression of TNFSF13 in tumor cells and fibroblasts is associated with poor prognosis in non-small cell lung cancer[J].Am J Clin Pathol,2014,141(2):226-233.