張偉偉， 吳志強

(東華理工大學 理學院，江西 南昌 330013)

癌癥是一類由于細胞分裂和凋亡機制失常而導致的疾病，通常表現為惡性腫瘤(Hanahan et al., 2011)。差異表達的基因可能直接或間接關系到癌癥的發(fā)生和易感性。差異表達基因的檢測能夠從基因層次揭示癌癥的發(fā)生機理并且能夠更好地探索癌癥治療的新方法(Wu et al., 2013)。

腫瘤樣本的差異表達分析通常是將腫瘤樣本和正常樣本直接進行比對。然而，臨床上獲得的腫瘤組織往往包含不同比例的正常細胞，如果不考慮正常細胞的混入，直接進行差異甲基化/基因表達或樣本聚類等下游分析，勢必會導致結果出現偏差甚至錯誤。在表觀遺傳層面上，針對腫瘤純度導致差異甲基化分析出現偏差的問題，Zheng等(2017)提出了腫瘤純度校正的差異甲基化分析模型 InfiniumDMC。Zhang等(2017)提出了一種考慮腫瘤細胞純度的腫瘤樣本亞型聚類方法 InfiniumClust。同樣，腫瘤樣本的“不純”也會使腫瘤樣本在基因上的均值估計產生偏差，進而導致差異表達分析結果出現偏差。為此Aran等(2015)將腫瘤純度作為線性“協(xié)變量”加入回歸模型，盡管這樣做得到了更多差異表達的基因。但是，通過分析發(fā)現，腫瘤純度與基因差異之間的關系是乘性的而非原來認為的線性關系。忽視腫瘤純度，或者將腫瘤純度作為“協(xié)變量”加入回歸模型都會使差異分析結果出現偏差。

為解決上述問題，本文在正態(tài)分布假設下構建了腫瘤樣本基因表達量的混合概率模型，通過廣義最小二乘法和Wald檢驗來確定每個基因在癌癥和正常細胞之間的差異性，并通過TCGA數據分析來說明所提模型的優(yōu)越性。

1 考慮腫瘤細胞純度的差異表達基因模型

基因的原始數據不服從正態(tài)分布，經過log2變換后數據的正態(tài)性會變強。圖1以肺腺癌為例，圖1a顯示肺腺癌正常樣本變換前后基因正態(tài)分布擬合檢驗的統(tǒng)計量，每個點代表一個基因，顯然變換后的統(tǒng)計量變得更加小，特別是對于那么違反正態(tài)假設的基因，因此log2變換有助于數據正態(tài)化。其次，對變換后的正常樣本進行正態(tài)分布擬合檢驗，圖1b顯示大部分p值均分分布，這表明正態(tài)分布假設對于絕大多數基因是成立的。雖然有一小部分基因具有較小的p值，但這些基因可能是正常樣本中不同細胞的差異表達基因。因此，在以下模型假設中可以將每個基因在正常樣本中的表達量假設成一個單一的正態(tài)分布。

圖1 正常樣本上所有基因正態(tài)分布檢驗的p值和統(tǒng)計量Fig.1 P-values and statistics for the normal distribution of all genes in normal samples

Yis′ =(1-λs)Xis+λsYis

=(1-λs)Xis+λs(Xis+δis)

=Xis+λsδis

令Zi=[Xi1,Xi2,…,Xin0,Yi1′,Yi2′…,Yin1′]。其中前n0項為正常樣本的基因表達量，后n1項為腫瘤樣本的基因表達量。上述問題變?yōu)榫€性模型Zis=mi+asμi+s,s=1,2,…,n0+n1。模型的參數可以通過廣義最小二乘法來求解，令

把腫瘤純度作為線性模型的一個實驗設計因子引入，會降低腫瘤樣本的內部方差，因此，從理論上來說對于差異表達基因的估計會有所提升。選取 InfiniumPurify估計的腫瘤細胞純度數據，并將所提方法命名為LinReg，通過與limma (Ritchie et al., 2015)在TCGA的14種癌癥類型(正常樣本數量超過10個的癌癥類型)上的差異分析比較來說明所提方法的優(yōu)越性。

2 TCGA腫瘤樣本差異表達基因分析結果

圖2 兩種方法在不同癌癥類型中各指標比對圖Fig.2 Comparison diagram of the two methods in different cancer types

首先，把每一種癌癥類型的腫瘤樣本隨機平均分成兩份，將這兩份腫瘤樣本分別和正常樣本利用limma和LinReg進行差異表達基因分析(圖2)。在同樣的顯著性指標下(FDR < 0.001)，對每一種方法將這兩組探測到的差異表達基因取交集，并且來比較一下交集中所包含的基因的數量。圖2a顯示對于絕大多數TCGA的癌癥類型，線性模型比limma得到了更多的交集基因，表明線性方法相比較limma具有更好的探測出診斷生物標志物的性能。接下來比較這兩種方法在同樣的顯著性指標下(FDR < 0.001)得到的差異表達的基因數量。如圖2b所示，對于絕大多數TCGA的癌癥類型，線性模型比limma得到了更多的差異表達的基因，表明線性方法敏感性更強。雖然不同癌癥類型具有不同的病理學基礎和差異表達的基因，但從整體上看，不同癌癥類型應具有一些相同點，這些共有的性質會帶來癌癥類型之間較高的統(tǒng)計量相關性。圖2c中列舉了不同癌癥類型與其他癌癥類型的平均統(tǒng)計量相關性，結果同樣表明線性方法比limma的相關性更強，說明其更具一致性。最后，來分析一下差異表達基因的生物學意義。選取FDR最小的前2 000個基因，將這2 000個基因與KEGG的“PATHWAY_IN_CANCER”上的328個基因取交集，來看一下這14種癌癥上兩種方法的交集基因個數。圖2d顯示線性方法在BLCA，BRCA，ESCA，KICH，KIRC，LIHC，LUSC，PRAD，STAD和UCEC這10種癌癥上具有更多的交集中的基因，其中LinReg在UCEC上的交集基因的個數幾乎是limma的5倍。這些結果說明線性方法得到的基因更具有生物學意義，這些基因對于癌癥診斷和治療具有重要意義。

圖3 只被線性方法探測到的差異表達的基因與腫瘤純度的關系Fig.3 The relationship between differentially expressed genes only detected by linear method and tumor purities

為了更好地理解腫瘤純度在差異表達基因分析中所起的作用，本文任意選取了一個只能被LinReg探測到差異的基因來研究該問題。圖3左面面板顯示該基因log2變換后在腫瘤樣本和正常樣本的表達量分布，該基因的均值在兩組之間是沒有差異的，limma檢驗的p值為0.872；中間的面板顯示該基因在腫瘤樣本上log2變換后的表達量和腫瘤樣本純度的相關系數，相關系數達-0.67；右面面板顯示校正腫瘤純度后，兩組之間均值有明顯差異，利用LinReg進行檢驗發(fā)現p=9.41×10-33。這些結果說明腫瘤樣本間較大的組內方差是由于腫瘤樣本的純度引起的，因此利用癌旁的正常樣本進行線性提純后，發(fā)現該基因在兩組之間具有明顯差異，這個結果進一步說明了差異分析時校正腫瘤純度的必要性。

3 結論

本文針對腫瘤樣本“不純”導致差異分析結果出現偏差的問題，構建了腫瘤樣本基因表達量的混合概率模型。該模型是把腫瘤純度作為影響差異表達基因分析的因素引入差異分析，相比直接使用腫瘤樣本作為癌癥細胞樣本與正常樣本直接比對的方法更加合理，有效避免了因為腫瘤樣本“不純”而導致差異表達分析出現的偏差。同樣腫瘤純度也影響著基因共表達網絡和共甲基化網絡，計劃將

在進一步工作中構建腫瘤純度校正的基因共表達網絡和共甲基化網絡，并分析其與傳統(tǒng)的共表達/共甲基化網絡在功能模塊分析中的差異。