汪 穎，張立瑩

(大連交通大學(xué)，遼寧 大連 116028)

0 引言

信息化快速發(fā)展，數(shù)據(jù)急速增加，需要更有效的聚類算法解決單細(xì)胞RNA測序數(shù)據(jù)分析問題。生物學(xué)研究中，傳統(tǒng)的批量RNA測序方法(RNA-seq)可以一次處理成千上萬個細(xì)胞，并得到變異的平均水平。但是，沒有兩個細(xì)胞是完全相同的，而單細(xì)胞測序則可以揭示出每個細(xì)胞獨(dú)特的微妙變化，甚至可以揭示全新的細(xì)胞類型。單細(xì)胞RNA測序大大提高了人們對生物系統(tǒng)的了解，而mRNA要翻譯成蛋白質(zhì)才能發(fā)揮作用。但在大多數(shù)轉(zhuǎn)錄組和蛋白組聯(lián)合組學(xué)文獻(xiàn)中，mRNA和蛋白質(zhì)的相關(guān)度只有0.4～0.6，從RNA到真正發(fā)揮功能還有很長的一段路要走，這是用RNA測序預(yù)測功能的劣勢[1]。單細(xì)胞RNA測序(單細(xì)胞測序)通過揭示具有高分辨率的單個細(xì)胞的異質(zhì)性，徹底改變了轉(zhuǎn)錄組學(xué)研究。聚類已成為識別細(xì)胞類型、描述其功能狀態(tài)和推斷潛在細(xì)胞動力學(xué)的常規(guī)分析手段。目前，關(guān)于RNA測序技術(shù)的生物信息學(xué)分析還很薄弱，僅憑KEGG通路分析、GSEA基因分析很難得到實(shí)質(zhì)有用的信息。針對RNA測序數(shù)據(jù)進(jìn)行有效分類成為一種必然趨勢，目前已存在一些聚類工具，如SC3[2]、SSCC[1]、Seurat[3]和CIDR[4]等，這些算法提高了數(shù)據(jù)聚類精度，但往往具有較高的計(jì)算復(fù)雜度。采用SCC策略，bigScale采用卷積策略通過貪婪搜索算法，將相似的單個單元合并為巨型單元。雖然這些聚類策略在某種程度上達(dá)到了一定的聚類效果，但仍有很大的改進(jìn)空間。

聚類是分析單細(xì)胞RNA測序數(shù)據(jù)的重要分析手段，但快速擴(kuò)展的數(shù)據(jù)量可能使該過程在計(jì)算上具有挑戰(zhàn)性。已有的SSCC方法是目前比較有效的方法，在保證一定聚類精確性的基礎(chǔ)上將計(jì)算復(fù)雜度從O(n2)降低到O(n)[1]。但SSCC在計(jì)算過程中使用的仍是傳統(tǒng)的相似度度量方法，其對于高維稀疏數(shù)據(jù)集的計(jì)算存在一些弊端。為了克服這些弊端，構(gòu)建了一個新的聚類公式，以識別可能屬于相同真實(shí)集群的數(shù)據(jù)點(diǎn)。與當(dāng)前的解決方案相比，此聚類方法考慮了附近數(shù)據(jù)對于聚類分類的影響，方案考慮得更全面，有更高的聚類準(zhǔn)確度。考慮到聚類算法對于大規(guī)模單細(xì)胞測序數(shù)據(jù)分析的重要性，提出了一種新的聚類框架，針對單細(xì)胞測序，構(gòu)建了二次相似性度量。已有工作表明，對于廣泛的數(shù)據(jù)分布，隨著維數(shù)的增加，傳統(tǒng)的相似性(如歐幾里德范數(shù)或余弦測度)變得不那么可靠(Beyeret al.，1999)[5]，因?yàn)橐恍?shù)據(jù)集在高維空間中變得稀疏，傳統(tǒng)的相似性度量可信度較低(Beyeret al.，1999)[5]，因此許多基于這些度量的聚類方法對于高維稀疏數(shù)據(jù)預(yù)測精度不高，效果不好。故而，利用二次相似性度量構(gòu)建相似矩陣、利用傳統(tǒng)度量構(gòu)建相似矩陣，基于共享鄰居排序，構(gòu)建二次相似性度量矩陣。

二次相似性度量是基于共享最近鄰(SNN)[6]的概念，考慮了周圍鄰居數(shù)據(jù)點(diǎn)的影響。確切地說，一對樣本數(shù)據(jù)之間的相似性度量是由初步度量即傳統(tǒng)相似性度量、二次度量基于共享鄰居排序構(gòu)建相似性度量構(gòu)成。在高維情況下，SNN度量比相關(guān)的主要度量更穩(wěn)健，并產(chǎn)生更穩(wěn)定的性能(Houleet al.，2010)[7]。SNN技術(shù)已成功應(yīng)用于一些聚類問題(Ertoězet al.，2003；Guhaet al.，2000；Jarvis和Patrick，1973)[8]。受這些早期應(yīng)用的啟發(fā)，根據(jù)共享鄰居的排序，重新定義了新的相似性度量。

1 方法

使用歐幾里得距離(也可以使用其他合適的度量)來計(jì)算傳統(tǒng)相似度矩陣，基于共享鄰居排序，給出二次相似度度量公式w(xi,xj)[6]。通過公式w(xi,xj)構(gòu)建新的二次相似性度量Y(xi,xj)及公式Y(jié)(xi,xj)進(jìn)行計(jì)算，通過輪廓值來評價結(jié)果的準(zhǔn)確性。

1.1 構(gòu)建二次度量w(xi,xj)

對于第i個樣本數(shù)據(jù)xi(i=1,2,3,…,n)，使用傳統(tǒng)相似性矩陣，列出k個最近鄰居，即：

KNN(xi)={xi1,xi2,xi3,…,xik}

式中，xik表示的是xi的第k個最近鄰居。

如圖所示，圖1中，1，2，3，…，9表示xi和xj的鄰居x1，x2，x3，…，x9，取k=6，則：

KNN(xi)={x2,x6,x3,x7,x9,x4}

KNN(xj)={x1,x5,x3,x7,x8,x4},

其分別表示樣本xi和xj的6個最近鄰居。xi，xj的共享鄰居為SNN(xi,xj)={x3,x7,x4}。僅當(dāng)xi和xj具有至少1個共享的KNN時，才保留邊w(xi,xj)。

圖1 圖中數(shù)字1～9表示樣本x1～x9Fig.1 Figure 1～9 represent sample x1～x9

對于任意的xi，xj{i,j=1,2,3,…,n}，樣本xi，xj的二次度量定義如下：

k表示最近鄰居集合的大小，rank(v,xi)代表樣本v在xi的最近鄰居集合KNN(xi)中的排序。

因此，共享最近鄰居(SNN)圖捕獲了兩個樣本在鄰域中的連通性方面的相似性[6]。與其他傳統(tǒng)相似性不同，在度量中，兩個樣本之間的相似性還考慮了該樣本與鄰居樣本(公共最近鄰)之間的相似性來確認(rèn)。與傳統(tǒng)相似性相比，SNN更適用于高維數(shù)據(jù)樣本的相似性度量。

1.2 改進(jìn)的二次相似性度量Y(xi,xj)

二次相似性度量w(xi,xj)在某種程度上避免了傳統(tǒng)相似性度量對于高位稀疏數(shù)據(jù)的弊端，在此基礎(chǔ)上進(jìn)行改善，除了考慮到了最近鄰居對于樣本相似性度量的影響，還考慮了絕對相似度度量(傳統(tǒng)度量)對樣本間相似性度量的影響。利用改進(jìn)的二次相似性度量可以糾正邊緣數(shù)據(jù)樣本錯誤分類問題，改進(jìn)的二次相似性度量公式如下：

d(xi,xj)表示xi，xj的傳統(tǒng)相似性度量，d(xi,xj)與Y(xi,xj)成反比，說明在同樣的w(xi,xj)條件下，兩樣本之間傳統(tǒng)距離越大相似度相對降低，反之，傳統(tǒng)距離越小相似度相對增強(qiáng)。

1.3 聚類方法

在可用的處理大規(guī)模單細(xì)胞RNA數(shù)據(jù)方法中，一般采取的是聚類，采樣之后進(jìn)行分類，其過程為數(shù)據(jù)初步處理，數(shù)據(jù)初始分類，利用聚類方法進(jìn)行重新分類，聚類檢驗(yàn)。利用k-means[9]和k-medoids[10]對傳統(tǒng)距離構(gòu)造相似性度量矩陣，利用二次相似性度量w(xi,xj)矩陣和改進(jìn)的二次相似性度量Y(xi,xj)矩陣進(jìn)行預(yù)測和計(jì)算。為了獲得獨(dú)立于聚類算法的評估結(jié)果，使用輪廓值來檢查這些特征工程方法的全局性能，進(jìn)行評估。

1.4 輪廓值

輪廓值[11]用于測量真實(shí)標(biāo)簽與原始標(biāo)簽之間的一致性，以更好地預(yù)測數(shù)據(jù)。輪廓是基于其緊密性和分離性的比較，輪廓值顯示了具體數(shù)據(jù)在子集中的位置。針對給定具有n個樣本和聚類方案的數(shù)據(jù)集，計(jì)算每個樣本的輪廓值。對于每個樣本xi，定義一個特定的值Si，即輪廓值。

圖2Fig.2

取數(shù)據(jù)集中的任何樣本xi，圖2表示第i個樣本xi屬于集合A。當(dāng)集合A中包含除xi以外的其他對象時，可以計(jì)算出A(i)，即xi與A中的所有其他對象的平均距離。考慮到與A不同的聚類C，即xi與C的所有對象的平均距離，記為D(i，C)。計(jì)算完平均距離之后，選擇這些數(shù)字中最小的一個，表示如下：

這就像對象xi的第二個最佳選擇，如果它不能容納在集群A中，集群B將是最接近的集合。圖2中，當(dāng)A本身被丟棄時，平均而言，集合B與樣本xi最近，因此了解數(shù)據(jù)集中每個對象的鄰居是非常有用的。

式中，ai是樣本xi與其自身聚類中的樣本的平均差異，而bi是樣本xi與任何其他聚類的最低平均差異。需要注意的是，必須假設(shè)簇k的數(shù)量大于1。Si的值范圍從-1～1。接近1的值意味著樣本xi與其聚類完全匹配，而接近-1的值意味著樣本xi如果被分類到其相鄰聚類中則更合適。

對于廣泛的數(shù)據(jù)分布，傳統(tǒng)的相似性(例如歐幾里得范數(shù)或余弦度量)隨著維數(shù)的增加而變得不那么可靠。使用一個二次相似性基于共享最近鄰居(SNN)的概念，一對數(shù)據(jù)點(diǎn)之間的相似性是它們之間交點(diǎn)的函數(shù)，由主要指標(biāo)(例如歐幾里得范數(shù))確定固定大小的鄰域。在高維中，SNN度量比主要度量更有效，并具有更穩(wěn)定的性能[6]。SNN技術(shù)已成功應(yīng)用于某些聚類問題，由此基于兩個數(shù)據(jù)點(diǎn)共享鄰域的排名來定義新的相似性。

2 結(jié)果

2.1 結(jié)果比較

在圖3～8中，每個點(diǎn)代表一個單元。使用原始數(shù)據(jù)直接計(jì)算的輪廓值顯示在X軸，使用通過Y(xi,xj)計(jì)算的輪廓值顯示在Y軸。在1處的輪廓值表示標(biāo)簽和特征之間的完美匹配，在-1處的輪廓值表示該單元可能被錯誤分類。每個圖中的百分比表示對角線上方的細(xì)胞分?jǐn)?shù)，即Y(xi,xj)二次相似性度量比傳統(tǒng)相似性度量的方法更好。測試的3個數(shù)據(jù)是Kolodziejczyk數(shù)據(jù)集、Pollen數(shù)據(jù)集、Usoskin數(shù)據(jù)集。

圖3 Pollen數(shù)據(jù)集，當(dāng)k=11時使用k-medoids方法進(jìn)行計(jì)算Fig.3 Pollen data, k-mediods is used to calculate when k=11

圖4 Pollen數(shù)據(jù)集，當(dāng)k=11時使用k-means方法進(jìn)行計(jì)算Fig.4 Pollen data, k-means is used to calculate when k=11

圖5 Kolodziejczyk數(shù)據(jù)集，當(dāng)k=3時使用k-medoids方法進(jìn)行計(jì)算Fig.5 Kolodziejczyk data, utilization of k-medoids when k=3

圖6 Kolodziejczyk數(shù)據(jù)集，當(dāng)k=3時使用k-means方法進(jìn)行計(jì)算Fig.6 Kolodziejczyk data, utilization of k-means when k=3

圖7 Usoskin數(shù)據(jù)集，當(dāng)k=4時使用k-medoids方法進(jìn)行計(jì)算Fig.7 Usoskin data, utilization of k-medoids when k=4

圖8 Usoskin數(shù)據(jù)集，當(dāng)k=4時使用k-means方法進(jìn)行計(jì)算Fig.8 Usoskin data, utilization of k-means when k=4

從6個圖形中能發(fā)現(xiàn)，除了利用k-medoids方法針對Usoskin數(shù)據(jù)集的計(jì)算結(jié)果(圖7)不太理想之外，其他的輪廓值圖形結(jié)果均顯示利用本研究提出的二次度量Y(xi,xj)比傳統(tǒng)相似性度量總體來說輪廓值體更高一些，效果更好。3個數(shù)據(jù)庫的計(jì)算結(jié)果顯示，Y(xi,xj)二次相似性度量的計(jì)算結(jié)果要比傳統(tǒng)距離相似性度量的結(jié)果更為準(zhǔn)確，準(zhǔn)確率是研究進(jìn)行的先決條件。

2.2 使用的單細(xì)胞測序數(shù)據(jù)集

使用4個單細(xì)胞測序數(shù)據(jù)集來評估特征空間中的聚類性能，包括Kolodziejczyk數(shù)據(jù)集[12]、Pollen數(shù)據(jù)集[13]、Usoskin數(shù)據(jù)集[14]。下面列出了這些數(shù)據(jù)集的詳細(xì)描述。

Kolodziejczyk數(shù)據(jù)集[12]包含704個具有3個簇的細(xì)胞，這些簇是在不同培養(yǎng)條件下從小鼠胚胎干細(xì)胞獲得的。通過使用Fluidigm C1系統(tǒng)并應(yīng)用SMARTer試劑盒獲得cDNA和大約10 000個基因以高測序深度(平均每個細(xì)胞9 000 000個讀數(shù)，>80%讀數(shù)映射到Mus musculus基因組GRCm38，外顯子>60%)進(jìn)行分析。 用于Illumina文庫制備的Nextera XT試劑盒。

Pollen數(shù)據(jù)集[13]包含249個細(xì)胞，其中11個簇是從皮膚細(xì)胞、多能干細(xì)胞、血細(xì)胞、神經(jīng)細(xì)胞等中獲得的?；贑1單細(xì)胞自動準(zhǔn)備集成流體回路的低或高測序深度，SMARTer超低RNA試劑盒和Nextera XT DNA樣品制備試劑盒用于描述單個細(xì)胞的基因表達(dá)譜(每個細(xì)胞約50000個讀數(shù))。

Usoskin數(shù)據(jù)集[14]包含622個具有4個簇的小鼠神經(jīng)元細(xì)胞，即具有肽能的傷害性感受器官，含有非肽能的傷害性感受器，含有神經(jīng)絲的細(xì)胞和含有酪氨酸羥化酶的細(xì)胞。用機(jī)器人細(xì)胞采集裝置挑選神經(jīng)元細(xì)胞，并在RNA-seq(114 000個讀數(shù)和每個細(xì)胞3574個基因)之前將其定位于96孔板的孔中。

3 討論

單細(xì)胞轉(zhuǎn)錄組分析中，通常需要根據(jù)單個細(xì)胞的基因表達(dá)水平對其進(jìn)行分組，以便每組對應(yīng)具有特定功能的細(xì)胞類型。這種分析有助于描述組織中的細(xì)胞組成及區(qū)分發(fā)育階段，從而更好地理解組織的生理學(xué)和病理學(xué)及發(fā)育過程。由于存在生物和技術(shù)差異，一種理想的全基因組單細(xì)胞數(shù)據(jù)聚類方法能夠有效區(qū)分細(xì)胞類型和高基因表達(dá)水平。針對一些高維較稀疏數(shù)據(jù)集，提出了二次相似性度量Y(xi,xj)，該相似性度量在一定程度上克服了傳統(tǒng)相似性度量對高維稀疏數(shù)據(jù)可信度較低的弊端。二次相似性度量Y(xi,xj)將為廣泛存在的高維稀疏數(shù)據(jù)集提供一種較為有效的相似性度量，為進(jìn)一步研究和分析提供有利條件。