基于基因表達(dá)譜預(yù)測(cè)腫瘤浸潤(rùn)免疫細(xì)胞類(lèi)型及比例的解卷積算法

裴晶晶，余彩裙，佘玉梅

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

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

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

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

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

(1)

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

表1 樣本信息表

2 混合樣本反卷積方法

2.1 反卷積總體策略

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

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

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

本文提出的反卷積算法是基于部分基因表達(dá)反卷積策略，是一種半監(jiān)督的卷積方法.需要提供參考細(xì)胞的表達(dá)信息，通過(guò)所有基因在不同細(xì)胞類(lèi)型中的表達(dá)情況推斷出具有細(xì)胞特異性表達(dá)的基因.并以特異性表達(dá)基因?yàn)榛A(chǔ)，對(duì)混合樣本進(jìn)行反卷積.推斷混合樣本的細(xì)胞類(lèi)型和比例信息.具體流程策略如圖1所示：

2.2 細(xì)胞特異性表達(dá)基因的鑒定

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

1) 計(jì)算出每個(gè)基因在不同樣本中的平均表達(dá)量

(2)

2) 去除冗余的背景信號(hào)

(3)

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

①高斯分布的概率密度：

(4)

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

(5)

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

(6)

(7)

④鑒定細(xì)胞特異性表達(dá)的基因

(8)

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

2.3 細(xì)胞特異性表達(dá)基因可靠性分析

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

2.4 反卷積算法

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

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

2) 去卷積算法過(guò)程：

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

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

(11)

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

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

③計(jì)算混合細(xì)胞可能的細(xì)胞類(lèi)型的相對(duì)豐度.另外，需要說(shuō)明的是，如果βi<0,則令βi=0

(12)

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

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

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

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

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

3 腫瘤免疫浸潤(rùn)細(xì)胞組分應(yīng)用

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

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

4 結(jié)語(yǔ)

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