一種應(yīng)用于單細(xì)胞RNA測(cè)序數(shù)據(jù)的稀有類型細(xì)胞檢測(cè)算法

馬 安 琪(復(fù)旦大學(xué)計(jì)算機(jī)科學(xué)技術(shù)學(xué)院智能信息處理重點(diǎn)實(shí)驗(yàn)室 上海 200433)

0 引 言

近年來(lái)，單細(xì)胞RNA測(cè)序技術(shù)得到了迅猛的發(fā)展。在這項(xiàng)技術(shù)之前，RNA測(cè)序使用的是批量測(cè)序技術(shù)，只能夠得到基因在組織樣本中所有細(xì)胞表達(dá)量的平均值，但無(wú)法得到更細(xì)粒度的信息。與此相比，單細(xì)胞RNA測(cè)序技術(shù)是細(xì)胞層級(jí)上的測(cè)序技術(shù)，測(cè)量的是基因在每個(gè)細(xì)胞上的表達(dá)量值。因此，單細(xì)胞RNA測(cè)序技術(shù)能夠進(jìn)一步揭示細(xì)胞之間的異質(zhì)性[1]。針對(duì)單細(xì)胞RNA測(cè)序數(shù)據(jù)的研究也越來(lái)越多，并被廣泛地應(yīng)用于免疫學(xué)、神經(jīng)生物學(xué)、干細(xì)胞和腫瘤研究等多個(gè)領(lǐng)域[2]。

稀有類型細(xì)胞檢測(cè)是單細(xì)胞RNA測(cè)序數(shù)據(jù)分析中的重要課題之一。其目標(biāo)是通過(guò)分析單細(xì)胞RNA測(cè)序數(shù)據(jù)，找到樣本組織中所屬類別規(guī)模占比極少量的稀有類型細(xì)胞。盡管這些稀有類型細(xì)胞數(shù)量很少，但它們并不是可以忽略的。這些稀有類型細(xì)胞往往扮演著重要的角色，例如，抗原特異性T細(xì)胞和腫瘤干細(xì)胞都是稀有類型細(xì)胞。抗原特異性T細(xì)胞對(duì)免疫記憶的形成起關(guān)鍵作用[3]。而干細(xì)胞能夠替換受損細(xì)胞，在治療帕金森疾病、心臟病和糖尿病等方面有著重要意義[4]。因此，稀有類型細(xì)胞檢測(cè)算法在生物上有著很強(qiáng)的應(yīng)用需求。

對(duì)于生物專家而言，檢測(cè)稀有類型細(xì)胞需要通過(guò)實(shí)驗(yàn)手段或者基于某些已有的先驗(yàn)知識(shí)。例如，已知某個(gè)稀有類型細(xì)胞的標(biāo)志性基因，通過(guò)找到在這些基因上具有較高表達(dá)量的細(xì)胞來(lái)得到這一稀有類型的細(xì)胞。但并不是在所有應(yīng)用場(chǎng)景中，都具備這種領(lǐng)域性的先驗(yàn)知識(shí)。例如，樣本中可能包含沒(méi)有被發(fā)現(xiàn)過(guò)的新稀有細(xì)胞類型，先驗(yàn)知識(shí)在這種情況下就不再適用。

在不需要先驗(yàn)知識(shí)的算法中，通過(guò)K-means[5]、層次聚類[6]等算法先對(duì)所有細(xì)胞進(jìn)行聚類，再找到其中規(guī)模較小的類作為稀有細(xì)胞，是一種較為直觀的方法。例如，Giniclust[7]和RaceID[8]都采用了聚類的方法來(lái)檢測(cè)稀有細(xì)胞。但是，這些方法都依賴于底層的聚類算法，而大多數(shù)聚類算法在類分布空間大小差別懸殊且空間規(guī)則性差時(shí)都表現(xiàn)不佳。另外，聚類算法需要計(jì)算細(xì)胞兩兩之間的距離，因此，隨著數(shù)據(jù)規(guī)模的增大，基于聚類的算法具有較高的時(shí)間花費(fèi)，而現(xiàn)在單細(xì)胞測(cè)序數(shù)據(jù)的趨勢(shì)是對(duì)越來(lái)越多的細(xì)胞進(jìn)行測(cè)序。由于目標(biāo)只是檢測(cè)出稀有細(xì)胞，并不需要對(duì)所有細(xì)胞進(jìn)行類別劃分，聚類算法的代價(jià)顯得過(guò)于昂貴。因此，寄希望于尋求一種有效的非聚類算法。

已有的一種非聚類思路是將稀有類型細(xì)胞看作是數(shù)據(jù)中的離群點(diǎn)，并使用離群點(diǎn)檢測(cè)算法，如LOF[9]。但由于離群點(diǎn)和稀有類型細(xì)胞并不完全重合，這類方法在大多數(shù)應(yīng)用場(chǎng)景之下，并不是最佳的方案。FiRE[10]是一種專門針對(duì)測(cè)序數(shù)據(jù)的稀有類型細(xì)胞檢測(cè)的算法。它基于哈希的思想，快速辨別稀有細(xì)胞。但FiRE的精確率并不高。本文受FiRE基本思想的啟發(fā)，提出了一種改進(jìn)的稀有類型細(xì)胞檢測(cè)算法。與FiRE不同，該方法用更為簡(jiǎn)單快速的多輪隨機(jī)劃分代替了FiRE的哈希過(guò)程，同時(shí)增加了基于最近鄰調(diào)整預(yù)測(cè)結(jié)果的過(guò)程，提高了結(jié)果的精確率。

1 數(shù)據(jù)說(shuō)明

1.1 數(shù)據(jù)格式

算法應(yīng)用于單細(xì)胞RNA測(cè)序的基因表達(dá)量矩陣數(shù)據(jù)。數(shù)據(jù)的基本格式如圖1所示，為了方便說(shuō)明，這里只截取了數(shù)據(jù)的一部分。

整個(gè)數(shù)據(jù)是一個(gè)數(shù)值矩陣，在此格式的數(shù)據(jù)中，每一行對(duì)應(yīng)一個(gè)基因，每一列對(duì)應(yīng)一個(gè)細(xì)胞。行名是基因的名字，列名是每個(gè)細(xì)胞對(duì)應(yīng)的標(biāo)識(shí)碼，每個(gè)標(biāo)識(shí)碼由“A”“T”“C”和“G”四種字符組成。標(biāo)識(shí)碼具有唯一性，每個(gè)標(biāo)識(shí)碼都和某個(gè)細(xì)胞樣本相對(duì)應(yīng)。矩陣中的數(shù)值表示了對(duì)應(yīng)基因在對(duì)應(yīng)細(xì)胞中測(cè)到的轉(zhuǎn)錄子的個(gè)數(shù)。

對(duì)單細(xì)胞RNA測(cè)序而言，表達(dá)量矩陣數(shù)據(jù)通常十分稀疏，大部分的矩陣項(xiàng)為零。另外，數(shù)據(jù)往往存在著一些噪聲。

1.2 數(shù)據(jù)預(yù)處理

在執(zhí)行算法之前，需要對(duì)數(shù)據(jù)進(jìn)行一些預(yù)處理。預(yù)處理的目的主要是去除測(cè)序過(guò)程中因技術(shù)原因造成的一些異常數(shù)據(jù)，并對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理。具體操作如下：

(1) 只保留至少在3個(gè)細(xì)胞上轉(zhuǎn)錄子數(shù)量達(dá)到2的基因。

(2) 計(jì)算每個(gè)細(xì)胞對(duì)應(yīng)的文庫(kù)大小，即列和，并將每一列除以文庫(kù)大小，然后乘上文庫(kù)大小的中位數(shù)。

(3) 對(duì)所有表達(dá)量值取自然對(duì)數(shù)。

(4) 根據(jù)標(biāo)準(zhǔn)化后的Fano系數(shù)(方差除以均值)挑選出在不同細(xì)胞中表達(dá)量最多變的前1 000個(gè)基因。

2 算法流程

本文提出的稀有類型細(xì)胞檢測(cè)算法整體流程如下。首先，該算法在預(yù)處理之后的單細(xì)胞RNA測(cè)序數(shù)據(jù)上，進(jìn)行多輪隨機(jī)劃分；然后綜合隨機(jī)劃分的結(jié)果，對(duì)每個(gè)細(xì)胞的稀有程度進(jìn)行打分，并初步預(yù)測(cè)出稀有類型細(xì)胞；最后，算法根據(jù)最近鄰對(duì)結(jié)果進(jìn)行調(diào)整，得到最終預(yù)測(cè)的稀有類型細(xì)胞結(jié)果。

2.1 隨機(jī)劃分

隨機(jī)劃分過(guò)程如下：

Step1隨機(jī)抽取一個(gè)未選擇過(guò)的基因，并在基因表達(dá)量的最小值和最大值之間以均勻概率隨機(jī)生成一個(gè)閾值。

Step2將未分類細(xì)胞中在該基因上表達(dá)量大于等于這個(gè)閾值的細(xì)胞分為一類。

Step3重復(fù)Step 1-Step 2，直到所有細(xì)胞分類完成，或者選取的基因達(dá)到了Gmax個(gè)，停止這個(gè)過(guò)程，并將剩余細(xì)胞分為一類。

隨機(jī)劃分過(guò)程的偽代碼如算法1所示。

算法1快速劃分

輸入：所有基因集合G，所有細(xì)胞集合C,表達(dá)量矩陣E。

1 已選基因數(shù)=0；

2 while已選基因數(shù)

3 gi←從G中隨機(jī)抽取一個(gè)基因；

4 已選基因數(shù)+=1；

5 G←G-{gi}；

6 threshold_i<-random(Emin,Emax)；

7 Ci←{c:E[gi][c]>=threshold_i}；

8 將Ci中的細(xì)胞標(biāo)記為一類；

9 C←C-Ci；

10 if C為空:

11 結(jié)束；

12 end if

13 end while

14 if還有未被分類的細(xì)胞:

15 將剩余細(xì)胞標(biāo)記為一類；

16 end if

輸出：細(xì)胞隨機(jī)劃分結(jié)果。

2.2 打 分

重復(fù)上一步隨機(jī)劃分過(guò)程L輪，根據(jù)這L輪的分類結(jié)果，對(duì)細(xì)胞的稀有程度進(jìn)行打分。本文認(rèn)為在隨機(jī)劃分情況下，與稀有類型細(xì)胞劃分在同一類的可能性是比較小的，用plx表示在第l輪中與細(xì)胞x分在一起的概率，表示為：

(1)

式中：n代表總細(xì)胞數(shù)，labelsl代表第l輪快速劃分后的結(jié)果類別標(biāo)簽。綜合L輪的隨機(jī)劃分結(jié)果，細(xì)胞x的稀有程度分?jǐn)?shù)為：

(2)

該分?jǐn)?shù)越高，則認(rèn)為細(xì)胞的稀有程度越高。

2.3 預(yù) 測(cè)

在得到連續(xù)的分值之后，算法還希望預(yù)測(cè)出二值化的標(biāo)簽，即預(yù)測(cè)出一個(gè)細(xì)胞是否是稀有類型細(xì)胞。算法采用和FiRE相同的判定策略，如果滿足式(3)，則把x初步標(biāo)記為稀有類型細(xì)胞，否則標(biāo)記為非稀有類型細(xì)胞。

scorex≥(75%分位數(shù))+1.5×(25%分位數(shù))

(3)

2.4 基于最近鄰算法調(diào)整

經(jīng)過(guò)以上幾步，算法已經(jīng)得到了初步的預(yù)測(cè)結(jié)果。但在初步的結(jié)果中，往往包含一些假陽(yáng)性的噪聲點(diǎn)。為了去除掉結(jié)果中大部分的假陽(yáng)性，接下來(lái)，算法基于最近鄰算法對(duì)結(jié)果進(jìn)行調(diào)整，以此提高結(jié)果的精確率。由于稀有類型細(xì)胞并不是一些孤立的離群點(diǎn)，而是按類聚集的。本文認(rèn)為，如果一個(gè)細(xì)胞是稀有類型細(xì)胞，它最近的鄰居也應(yīng)該是稀有類型細(xì)胞。

在求最近鄰前，算法先通過(guò)PCA[12]降維算法，把數(shù)據(jù)降至20維。接下來(lái)，對(duì)上一步中每一個(gè)稀有細(xì)胞，基于歐氏距離計(jì)算降維后最近的k個(gè)鄰居。算法將檢查它的k個(gè)最近鄰居是否都在上一步標(biāo)記成稀有類型細(xì)胞，如果有至少一個(gè)不是，則重新標(biāo)記為非稀有類型細(xì)胞。

3 實(shí)驗(yàn)與結(jié)果分析

3.1 數(shù)據(jù)集

實(shí)驗(yàn)使用了兩個(gè)來(lái)自10X Genomics測(cè)序平臺(tái)的單細(xì)胞RNA測(cè)序數(shù)據(jù)集。這兩個(gè)數(shù)據(jù)集都已經(jīng)根據(jù)已有的生物知識(shí)打好了細(xì)胞類型的標(biāo)簽，作為實(shí)驗(yàn)中計(jì)算各個(gè)指標(biāo)的標(biāo)準(zhǔn)答案。

第一個(gè)數(shù)據(jù)集PBMC[11]是中測(cè)試檢測(cè)稀有類型細(xì)胞的一個(gè)生成數(shù)據(jù)，采樣自外周血單核細(xì)胞測(cè)序數(shù)據(jù)，包含CD14+Monocyte、CD56+NK和CD19+B三種類型的細(xì)胞，其中CD14+Monocyte是該采樣數(shù)據(jù)中占比少數(shù)的稀有類型細(xì)胞。第二個(gè)數(shù)據(jù)集是FiRE工具包中的樣例數(shù)據(jù)，包含293T、Jurkat兩種類型的細(xì)胞，Jurkat是該數(shù)據(jù)集中的稀有類型細(xì)胞。

各數(shù)據(jù)集中稀有類型細(xì)胞具體情況如表1所示。

3.2 評(píng)估指標(biāo)

實(shí)驗(yàn)中使用了三種評(píng)估指標(biāo)，分別是精確率P(Precision)、召回率R(Recall)和綜合前兩項(xiàng)指標(biāo)的F1值。具體計(jì)算公式如下：

(4)

(5)

(6)

式中：TP、FP和FN分別表示預(yù)測(cè)結(jié)果中真陽(yáng)性、假陽(yáng)性和假陰性的個(gè)數(shù)。

3.3 實(shí)驗(yàn)結(jié)果

本文實(shí)驗(yàn)統(tǒng)一取參數(shù)Gmax=50，L=100，k=2。對(duì)于對(duì)比算法，均使用官方推薦的默認(rèn)參數(shù)。

圖2是本文方法在Jurkat-293T數(shù)據(jù)集上各細(xì)胞稀有程度分?jǐn)?shù)的灰度圖，該圖中細(xì)胞的分布根據(jù)原數(shù)據(jù)經(jīng)過(guò)t-SNE[13]降維后得到，稀有程度分?jǐn)?shù)越高，圖中對(duì)應(yīng)亮度越低。圖2中下方真實(shí)稀有類型細(xì)胞所在簇的顏色都較暗，而非稀有細(xì)胞在灰度圖上較亮。另外，盡管有一些非稀有類型的細(xì)胞也顯示出了較暗的顏色，它們中的大部分會(huì)在基于最近鄰調(diào)整的階段被從預(yù)測(cè)結(jié)果中剔除。

圖2 稀有程度分?jǐn)?shù)分布的灰度圖

圖3是在Jurkat-293T數(shù)據(jù)集上基于最近鄰調(diào)整前后預(yù)測(cè)結(jié)果的對(duì)比，其中：深黑色的點(diǎn)代表本文算法預(yù)測(cè)出的稀有類型細(xì)胞，淺灰色的點(diǎn)代表非稀有類型細(xì)胞。算法初步預(yù)測(cè)出的結(jié)果中，雖然大部分都是真實(shí)的稀有類型細(xì)胞，但也包含了一些非稀有類型細(xì)胞。而在基于最近鄰調(diào)整后，結(jié)果中的大部分假陽(yáng)性都被剔除了。盡管右上角仍然存在一些被誤判為稀有類型細(xì)胞的點(diǎn)，但總體上調(diào)整后的結(jié)果和數(shù)據(jù)集真實(shí)答案基本一致。

表2是該方法與FiRE、LOF在兩個(gè)單細(xì)胞RNA測(cè)序數(shù)據(jù)集上的對(duì)比結(jié)果。從整體表現(xiàn)來(lái)看，LOF最差，F(xiàn)iRE稍好一些，而本文方法是最好的。

表2 在各數(shù)據(jù)集上的對(duì)比結(jié)果

其中，F(xiàn)iRE和本文方法在兩個(gè)數(shù)據(jù)集中的Recall都是1，也就是說(shuō)，這兩種方法都能找出兩個(gè)數(shù)據(jù)集中所有的稀有類型細(xì)胞。而LOF會(huì)遺漏掉部分稀有類型細(xì)胞。另外，LOF和FiRE的精確率都不夠高，因此導(dǎo)致最后的F1值也比本文方法低。在這兩個(gè)單細(xì)胞RNA測(cè)序數(shù)據(jù)集的結(jié)果表明，本文方法在檢測(cè)稀有細(xì)胞類型問(wèn)題上比其他方法更加精確有效。

為了進(jìn)一步研究稀有類型細(xì)胞所占比例對(duì)算法的影響，本文調(diào)整了Jurkat-293T數(shù)據(jù)中稀有細(xì)胞的比例，記錄三種算法在不同稀有類型細(xì)胞比例數(shù)據(jù)集F1值的結(jié)果。圖4表明，三種算法的效果都隨著稀有細(xì)胞所占比例的下降而變差，但本文方法始終是表現(xiàn)最好的，且下降幅度比其他兩種算法要低。即使在最低比例(0.5%)的數(shù)據(jù)集上，本文方法依然有著較高的F1值，這表明本文方法在不同稀有類型比例的數(shù)據(jù)集上，具有一定的穩(wěn)定性。

圖4 三種算法在不同稀有類型細(xì)胞比例數(shù)據(jù)集上的F1值

4 結(jié) 語(yǔ)

本文提出了一種應(yīng)用于單細(xì)胞RNA測(cè)序表達(dá)量數(shù)據(jù)的稀有類型細(xì)胞檢測(cè)算法。它基于多輪隨機(jī)劃分的結(jié)果，對(duì)每個(gè)細(xì)胞的稀有程度進(jìn)行打分，然后基于最近鄰對(duì)結(jié)果進(jìn)行調(diào)整，得到最終每個(gè)細(xì)胞是否為稀有類型的預(yù)測(cè)標(biāo)簽。

這一方法不依賴于生物的領(lǐng)域知識(shí)，同時(shí)也避免了聚類算法中較高的時(shí)間復(fù)雜度。從單細(xì)胞RNA測(cè)序數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果來(lái)看，它具有較高的精確率和召回率，表現(xiàn)優(yōu)于其他非聚類方法。盡管當(dāng)稀有類型細(xì)胞比例下降時(shí)，本文方法表現(xiàn)也有所下降，但仍然要比其他方法在同比例的情況下要好。未來(lái)工作將考慮把本文方法和聚類相結(jié)合，先用本文提出的稀有類型檢測(cè)方法預(yù)測(cè)出稀有類型細(xì)胞，然后對(duì)稀有類型細(xì)胞和非稀有類型細(xì)胞采用不同的聚類策略，以在單細(xì)胞RNA測(cè)序數(shù)據(jù)上得到更加準(zhǔn)確的聚類結(jié)果。