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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 算法流程

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

2.1 隨機(jī)劃分

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

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

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

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

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

算法1快速劃分

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

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

2 while已選基因數(shù)

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

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ī)劃分過程L輪，根據(jù)這L輪的分類結(jié)果，對細(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ù) 測

在得到連續(xù)的分值之后，算法還希望預(yù)測出二值化的標(biāo)簽，即預(yù)測出一個細(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)過以上幾步，算法已經(jīng)得到了初步的預(yù)測結(jié)果。但在初步的結(jié)果中，往往包含一些假陽性的噪聲點。為了去除掉結(jié)果中大部分的假陽性，接下來，算法基于最近鄰算法對結(jié)果進(jìn)行調(diào)整，以此提高結(jié)果的精確率。由于稀有類型細(xì)胞并不是一些孤立的離群點，而是按類聚集的。本文認(rèn)為，如果一個細(xì)胞是稀有類型細(xì)胞，它最近的鄰居也應(yīng)該是稀有類型細(xì)胞。

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

3 實驗與結(jié)果分析

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

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

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

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

3.2 評估指標(biāo)

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

(4)

(5)

(6)

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

3.3 實驗結(jié)果

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

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

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

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

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

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

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

為了進(jìn)一步研究稀有類型細(xì)胞所占比例對算法的影響，本文調(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īng)用于單細(xì)胞RNA測序表達(dá)量數(shù)據(jù)的稀有類型細(xì)胞檢測算法。它基于多輪隨機(jī)劃分的結(jié)果，對每個細(xì)胞的稀有程度進(jìn)行打分，然后基于最近鄰對結(jié)果進(jìn)行調(diào)整，得到最終每個細(xì)胞是否為稀有類型的預(yù)測標(biāo)簽。

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