裴晶晶，余彩裙，佘玉梅

( 云南民族大學 數(shù)學與計算機科學學院，云南 昆明 650500 )

腫瘤不僅僅是惡性細胞群，而且是由不同類型細胞組成的復雜生態(tài)系統(tǒng).腫瘤細胞的子代具有無限增殖遺傳特性，這種性質(zhì)不顧正常免疫系統(tǒng)約束，給個體健康造成了致命的危險.腫瘤按是否轉(zhuǎn)移可分為是良性腫瘤和惡性腫瘤.前者可通過手術(shù)切除達到治療的目的，而后者則會在生物體的其它部位形成繼發(fā)性(轉(zhuǎn)移)腫瘤，散布在身體周圍以及其他組織中，以至難以通過簡單切除來進行有效的治療.在針對惡性腫瘤治療過程中，腫瘤浸潤性免疫細胞在腫瘤控制和對治療的反應中起重要作用[1-3]，不同類型的腫瘤細胞的定量可以揭示抗癌應答的潛在機制，并有助于評估抗癌治療過程中的整體效果，對后續(xù)深入治療起到重要的指導作用.目前，利用流式細胞計數(shù)技術(shù)可以準確的定量腫瘤組織中免疫細胞的類型及比例，但是該方法需要大量的人力和實驗成本.然而，利用計算的方法可以直接推斷出復雜組織中的細胞類型及其比例，該策略具有快速、準確的特點，對腫瘤診斷、治療以及機制的研究具有重要的意義.

目前，針對復雜組織細胞的鑒定和含量的計算方法均是圍繞著反卷積策略來進行的.在過去幾年中，一些已發(fā)表的反卷積算法試圖解決逆向解決免疫細胞基因表達譜的任務. 如2010年Shen Orr[4]等提出了一種稱為“CSSAM”的算法(微陣列的細胞類型特異性顯著性分析)，該算法的開發(fā)是基于傳統(tǒng)的微陣列分析方法而忽略了樣本細胞類型的組成.以至于無法準確區(qū)分基因表達與不同的細胞類型之間的關(guān)系.Buettner等在2015中提出了一種稱為“SCLVM”的計算方法(單細胞潛變量模型)，它使用潛變量模型來解釋尋找新的細胞亞群.該技術(shù)允許在未分化的T細胞分化為T輔助細胞的過程中，識別與不同階段對應的細胞的不可檢測的亞群.Renaud Gaujoux已經(jīng)出版了一個R包，CellMix，其包含了一些已發(fā)表的計算反卷積方法[4].紐曼等在2015年發(fā)表了一篇論文中展示了一種名為CIBERSORT的新開發(fā)工具，該工具通過使用適用于免疫細胞譜系中廣泛細胞類型的支持向量機(SVM)算法的變體，其性能優(yōu)于所有其他方法[5].與早期的算法(通常是線性方法)相比，這是一種全新的反卷積方法.盡管一些已發(fā)表的反卷積方法顯示出不同表型細胞的一些有望的結(jié)果.但是，這些文章中僅限于對特定動物或人體特定部位的測試，如來自肝臟的腫瘤，或者甚至試圖對遠處相關(guān)的免疫細胞進行反卷積[6].這使得對彼此不同的方法進行基準測試變得非常困難.我們需要一種工具能夠在性能和結(jié)果方面比較相同數(shù)據(jù)集上的各種技術(shù).這樣的數(shù)據(jù)集應該來自真實的組織，其中存在不同細胞和基因的確切數(shù)量，并且還可以摻入腫瘤和噪聲以模擬來自腫瘤的真實樣品.

1 樣本數(shù)據(jù)準備及數(shù)據(jù)預處理

我們獲取的數(shù)據(jù)包含2個方面：①真實的組織樣本微陣列實驗；②真實的組織樣本RNA-Seq測序數(shù)據(jù)(見表1).我們利用上述數(shù)據(jù)集對本文反卷積算法與其他3種主流算法進行比較.其中對于真實組織來源的數(shù)據(jù)集，組織樣本中細胞混合比例是已知的.另外，我們對采集到的數(shù)據(jù)作了如下預處理：①對于真實來源的微陣列實驗數(shù)據(jù)，我們以2為底數(shù)對探針的信號進行對數(shù)轉(zhuǎn)換，并采用R語言中的bioMart包將探針映射到對應的基因上.②對于真實來源的RNA-seq數(shù)據(jù)，我們將每個樣品采用TPM進行標準化，然后將觀察到的每個基因的read數(shù)目加1，并以2為底數(shù)進行對數(shù)轉(zhuǎn)化.所有數(shù)據(jù)(真實組織來源和模擬數(shù)據(jù))均采用分位數(shù)標準化策略進行歸一化，以此消除批次或文庫大小帶來的差異.

(1)

其中，ni表示基因i的read數(shù)目；lj表示基因i的長度.

表1 樣本信息表

2 混合樣本反卷積方法

2.1 反卷積總體策略

從Venet[9]等開始，許多研究者提供了如何從全部基因表達量估計細胞類型、組織特異性特征以及細胞類型比例的方法.概括來說，根據(jù)所需輸入數(shù)據(jù)的不同，這些方法可歸納為兩種不同的類型，具體如下：

1) 部分基因表達反卷積方法 該方法需要提供細胞類型特異性特征[5,10-13]或不同細胞類型在特定組織樣品中混合比例統(tǒng)計特征[14-15].

2) 完整的去卷積方法 該方法直接從異質(zhì)樣本的全部基因表達數(shù)據(jù)中估計細胞、組織特征和比例[9,16-18].

本文提出的反卷積算法是基于部分基因表達反卷積策略，是一種半監(jiān)督的卷積方法.需要提供參考細胞的表達信息，通過所有基因在不同細胞類型中的表達情況推斷出具有細胞特異性表達的基因.并以特異性表達基因為基礎，對混合樣本進行反卷積.推斷混合樣本的細胞類型和比例信息.具體流程策略如圖1所示：

2.2 細胞特異性表達基因的鑒定

構(gòu)建具有細胞特異性表達的基因表達譜矩陣是后續(xù)去卷積算法的基礎，即在去卷積之前過濾掉非特異性或者特異性較弱的基因[19-22].具體而言，使用細胞特異性基因表達譜矩陣具有如下優(yōu)點：①減少內(nèi)存和計算層面上的消耗，加快了運行時間；②高信噪比-篩選出具有高信噪比的細胞特異性表達基因，提升了算法的有效性和靈敏度.目前半監(jiān)督卷積策略鑒定細胞特異性表達基因的主流策略是等方差或者異方差的t-test. 本文提出一種新的鑒定細胞特異性表達基因的策略，具體過程如下：

1) 計算出每個基因在不同樣本中的平均表達量

(2)

2) 去除冗余的背景信號

(3)

3) 擬合高斯分布，估計出高斯分布的參數(shù)μ,σ

①高斯分布的概率密度：

(4)

②對數(shù)似然函數(shù)：

(5)

③被估計的參數(shù)為

(6)

(7)

④鑒定細胞特異性表達的基因

(8)

其中，K表示參考細胞樣本的數(shù)目，N表示基因的數(shù)目，gij表示第i個基因在第j個樣本中的表達量，Si表示第i個基因在所有樣本中最大的信號，μ和σ是高斯分布的2個參數(shù).

2.3 細胞特異性表達基因可靠性分析

為了驗證細胞特異性表達基因鑒定算法的可靠性，以GSE19830數(shù)據(jù)集作為實例，該數(shù)據(jù)中包括來自組織Lung、Liver、Brain 3種不同的細胞類型.利用數(shù)據(jù)集提供的參考樣本推斷出標簽矩陣，并對矩陣中包含的細胞特異性基因進行了分析.結(jié)果顯示所鑒定的基因在不同的組織中具有明顯不同的表達模式，表明我們給出的鑒定標簽矩陣的策略是有效的(見圖2(a)).另外，進一步分析了不同算法鑒定出標簽基因的重疊程度.結(jié)果顯示文中的算法鑒定出的標簽基因與CIBERSORT[7]和dtangle[8]鑒定具有高重疊度(見圖2(b)).

2.4 反卷積算法

針對混合樣本去卷積問題常被建模成帶有約束條件的二次規(guī)劃或者LARSOR回歸問題，盡管CIBERSORT[7]提出了一種基于SVM的全新去卷積算法，但是其本質(zhì)仍然是回歸的問題.在此，提出了一種基于逐步回歸的去卷積策略，這一策略的優(yōu)勢在于逐步剔除不顯著的變量，該方法在現(xiàn)有文獻中未見報道.另外，在文章中所涉及到的加粗字母均表示向量.

1) 逐步回歸的主要思路：逐步回歸是以常規(guī)的線性回歸策略為基礎，考慮的全部自變量對響應變量的作用大小，將作用不顯著的變量剔除.以此保證預測方程的顯著性.

2) 去卷積算法過程：

①建立混合樣本基因表達量與細胞特異性表達基因表達量之間的線性關(guān)系.

E(Y|X)=β0+X·βT,

(11)

其中，Y表示細胞特異性表達基因在混合樣本中的表達量，E(Y|X)為目標期望值；β0,β為待估參數(shù).X為參考細胞樣本對應的細胞特異性基因的表達量.

②向后剔除：即將所有變量均放入模型，然后嘗試將其中一個自變量從模型中剔除，看整個模型解釋因變量的變異是否有顯著變化，之后將使解釋量減少最少的變量剔除；此過程不斷迭代，直到?jīng)]有自變量符合剔除的條件.

③計算混合細胞可能的細胞類型的相對豐度.另外，需要說明的是，如果βi<0,則令βi=0

(12)

2.5 結(jié)果比較與分析

文中收集了4個具有真實比例的數(shù)據(jù)集(表1)，其中包含93個芯片樣本數(shù)據(jù)和54個RNA-Seq測序數(shù)據(jù)的樣本.這些混合樣本對應的不同細胞真實比例的信息是已知的.利用這些數(shù)據(jù)來測試本文算法，并與CIBERSORT進行比較，結(jié)果顯示本文算法的測試結(jié)果與真實細胞的比例具有較好的一致性.其中圖中的點越靠近對角線，表明與真實結(jié)果越接近.

通過圖3可以看出，Our算法在GSE5350與PRJEB8231數(shù)據(jù)集中的預測結(jié)果與真實比例之間的相關(guān)系數(shù)高于CIBERSORT(圖3(b)和(d))，盡管在GSE19830和GSE64098中相關(guān)系數(shù)低于CIBERSORT, 但是仍表現(xiàn)出不錯的預測性能(圖3(a)和(c)).為進一步分析Our算法與其他算法的比較情況，引入了DSA[23]和dtangle[8]進行了更進一步的比較分析.通過表2可以看出，Our算法在GSE5350和PRJEB8231 2個測試集中表現(xiàn)性能排第1，在GSE19830和GSE64098中盡管僅排第3，但仍優(yōu)于DSA算法的性能.

表2 多個算法解卷積結(jié)果比較

數(shù)據(jù)集CIBERSORTdtangleDSAOurRankGSE198300.9920.9910.9620.9763GSE53500.9700.9650.9120.9731GSE640980.9940.9890.9750.9883PRJEB82310.8890.8980.7520.9021

3 腫瘤免疫浸潤細胞組分應用

為了進一步探討本文去卷積算法在腫瘤浸潤性免疫細胞混合樣本的有效性，我們從 https://github.com/gjhunt/dtangle 下載了20個腫瘤浸潤性細胞混合樣本，該數(shù)據(jù)集包含了多個不同的免疫細胞類型，分別是：B、Dendritic、Eosinophils、Macrophages、Mast、Monocytes、Neutrophils、NK、Plasma、T. 且這些樣本具有真實的不同細胞類型的混合比例信息.我們利用本文提出的算法對這些樣本進行解卷積，并將去卷積后預測比例與真實比例進行比較.

通過圖4可以看出，我們提出的算法在分析腫瘤免疫浸潤性細胞比例中具有不錯的預測效果，混合樣本去卷積預測結(jié)果與真實細胞的混合比例具有高一致性，顯示了該算法的有效性和可靠性.

4 結(jié)語

提出了一種基于逐步回歸模型對混合細胞樣本去卷積新算法.該算法主要包含如下2個方面：①鑒定具有細胞特異性的基因，該過程首先計算每個基因在不同參考細胞中的前景信號，然后將每個基因在不同樣本中的最強信號擬合高斯分布并通過最大似然方法估計相應的參數(shù)，最后將落在Si>μ+3σ的基因作為具有細胞特異性表達的基因.②利用具有細胞特異性表達的基因?qū)旌蠘颖救ゾ矸e，該過程采用逐步回歸的策略，過濾掉與模型擬合不顯著的協(xié)變量.我們利用4個數(shù)據(jù)集共147個樣本對該算法進行測試，并與CIBEROSRT、DSA、dtangle進行比較，結(jié)果顯示我們的算法具有不錯的解卷積能力.最后我們將算法應用于20個腫瘤浸潤性混合樣本數(shù)據(jù)，結(jié)果顯示我們的算法在對免疫浸潤定量的過程中，具有不錯的準確性能.