一種應用于單細胞RNA測序數據的稀有類型細胞檢測算法

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

0 引 言

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

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

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

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

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

1 數據說明

1.1 數據格式

算法應用于單細胞RNA測序的基因表達量矩陣數據。數據的基本格式如圖1所示，為了方便說明，這里只截取了數據的一部分。

整個數據是一個數值矩陣，在此格式的數據中，每一行對應一個基因，每一列對應一個細胞。行名是基因的名字，列名是每個細胞對應的標識碼，每個標識碼由“A”“T”“C”和“G”四種字符組成。標識碼具有唯一性，每個標識碼都和某個細胞樣本相對應。矩陣中的數值表示了對應基因在對應細胞中測到的轉錄子的個數。

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

1.2 數據預處理

在執行算法之前，需要對數據進行一些預處理。預處理的目的主要是去除測序過程中因技術原因造成的一些異常數據，并對數據進行標準化處理。具體操作如下：

(1) 只保留至少在3個細胞上轉錄子數量達到2的基因。

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

(3) 對所有表達量值取自然對數。

(4) 根據標準化后的Fano系數(方差除以均值)挑選出在不同細胞中表達量最多變的前1 000個基因。

2 算法流程

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

2.1 隨機劃分

隨機劃分過程如下：

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

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

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

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

算法1快速劃分

輸入：所有基因集合G，所有細胞集合C,表達量矩陣E。

1 已選基因數=0；

2 while已選基因數

3 gi←從G中隨機抽取一個基因；

4 已選基因數+=1；

5 G←G-{gi}；

6 threshold_i<-random(Emin,Emax)；

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

8 將Ci中的細胞標記為一類；

9 C←C-Ci；

10 if C為空:

11 結束；

12 end if

13 end while

14 if還有未被分類的細胞:

15 將剩余細胞標記為一類；

16 end if

輸出：細胞隨機劃分結果。

2.2 打 分

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

(1)

式中：n代表總細胞數，labelsl代表第l輪快速劃分后的結果類別標簽。綜合L輪的隨機劃分結果，細胞x的稀有程度分數為：

(2)

該分數越高，則認為細胞的稀有程度越高。

2.3 預 測

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

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

(3)

2.4 基于最近鄰算法調整

經過以上幾步，算法已經得到了初步的預測結果。但在初步的結果中，往往包含一些假陽性的噪聲點。為了去除掉結果中大部分的假陽性，接下來，算法基于最近鄰算法對結果進行調整，以此提高結果的精確率。由于稀有類型細胞并不是一些孤立的離群點，而是按類聚集的。本文認為，如果一個細胞是稀有類型細胞，它最近的鄰居也應該是稀有類型細胞。

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

3 實驗與結果分析

3.1 數據集

實驗使用了兩個來自10X Genomics測序平臺的單細胞RNA測序數據集。這兩個數據集都已經根據已有的生物知識打好了細胞類型的標簽，作為實驗中計算各個指標的標準答案。

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

各數據集中稀有類型細胞具體情況如表1所示。

3.2 評估指標

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

(4)

(5)

(6)

式中：TP、FP和FN分別表示預測結果中真陽性、假陽性和假陰性的個數。

3.3 實驗結果

本文實驗統一取參數Gmax=50，L=100，k=2。對于對比算法，均使用官方推薦的默認參數。

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

圖2 稀有程度分數分布的灰度圖

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

表2是該方法與FiRE、LOF在兩個單細胞RNA測序數據集上的對比結果。從整體表現來看，LOF最差，FiRE稍好一些，而本文方法是最好的。

表2 在各數據集上的對比結果

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

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

圖4 三種算法在不同稀有類型細胞比例數據集上的F1值

4 結 語

本文提出了一種應用于單細胞RNA測序表達量數據的稀有類型細胞檢測算法。它基于多輪隨機劃分的結果，對每個細胞的稀有程度進行打分，然后基于最近鄰對結果進行調整，得到最終每個細胞是否為稀有類型的預測標簽。

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