雙自編碼結(jié)合變分貝葉斯的單細(xì)胞RNA-Seq聚類

摘 要：近年來單細(xì)胞RNA測序（scRNA-seq）技術(shù)的快速發(fā)展使得在單個(gè)細(xì)胞水平上研究組織器官的異質(zhì)性成為可能。針對單細(xì)胞RNA測序數(shù)據(jù)中準(zhǔn)確鑒定細(xì)胞類型問題，提出一種新的基于雙自編碼結(jié)合變分貝葉斯高斯混合模型的聚類方法，稱之為sc-VBDAE。首先通過對抗自編碼網(wǎng)絡(luò)的編碼和解碼過程重構(gòu)數(shù)據(jù)，然后使用經(jīng)典自編碼對數(shù)據(jù)進(jìn)行降維，獲得低維且有效的數(shù)據(jù)。最后使用變分貝葉斯高斯混合模型對細(xì)胞進(jìn)行聚類，并可視化聚類結(jié)果。在10個(gè)scRNA-seq 數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果表明，該方法在6個(gè)數(shù)據(jù)集上ARI指標(biāo)均優(yōu)于其它方法，在數(shù)據(jù)集Biase和Klein上ARI指標(biāo)值達(dá)到0.90及以上。

關(guān)鍵詞：單細(xì)胞RNA測序；對抗自編碼；自編碼網(wǎng)絡(luò)；變分貝葉斯；細(xì)胞聚類

DOI：10.15938/j.jhust.2024.03.015

中圖分類號(hào)： Q811.4

文獻(xiàn)標(biāo)志碼： A

文章編號(hào)： 1007-2683（2024）03-0125-09

Single-cell RNA-Seq Clustering Based

on Dual Autoencoder with Variational Bayes

JIA Jihua， XU Yaokui， WANG Minghui

（College of Mathematics and Physics， Qingdao University of Science and Technology， Qingdao 266061， China）

Abstract：In recent years， the rapid development of single-cell RNA sequencing（scRNA-seq） technology has made it possible to research the heterogeneity of tissues and organs at the single-cell level. To accurately identify cell types in scRNA-seq data， based on dual autoencoder combined with variational Bayesian Gaussian mixture mode， a new clustering method， sc-VBDAE， is proposed. First， through the encoding and decoding process in adversarial autoencoder network， the scRNA-seq data is reconstructed. Then， the autoencoder network is used to reduce the dimensionality of the data， so as to obtain low-dimensional and effective scRNA-seq data. Finally， the variational Bayesian Gaussian mixture model is used to cluster the cells and visualize the clustering results. The experimental results on ten scRNA-seq datasets show that the ARI index of the proposed method is superior to other methods on six datasets， and the ARI index value on Biase and Klein datasets reaches 0.90 or above.

Keywords：single-cell RNA sequencing; adversarial autoencoder; autoencoder network; variational Bayes; cell clustering

0 引 言

轉(zhuǎn)錄組測序技術(shù) （RNA-seq） 是轉(zhuǎn)錄組圖譜分析的重要技術(shù)，然而傳統(tǒng)的RNA-seq技術(shù)只能從整體水平研究基因功能和基因結(jié)構(gòu)。近幾年開發(fā)的新一代高通量單細(xì)胞RNA測序 （scRNA-seq） 技術(shù)，由于可以獨(dú)立提供單個(gè)細(xì)胞的 RNA 表達(dá)譜，這允許研究人員在單個(gè)細(xì)胞水平上分析細(xì)胞異質(zhì)性和轉(zhuǎn)錄組異質(zhì)性，使其成為在單細(xì)胞規(guī)模上研究轉(zhuǎn)錄組學(xué)的有力工具。隨著STRT-Seq［1］，smart-seq2［2］，Drop-seq［3］等新的scRNA-seq技術(shù)的提出，scRNA-seq 的應(yīng)用方向越來越廣泛，例如研究癌細(xì)胞組織內(nèi)異質(zhì)性，神經(jīng)元亞型的鑒定和探索腫瘤細(xì)胞的表型狀態(tài)。

與此同時(shí)，scRNA-seq的發(fā)展給人們帶來了新的機(jī)遇，但仍然面臨著諸多挑戰(zhàn)。比如，許多現(xiàn)有的無監(jiān)督方法在細(xì)胞聚類性能的檢驗(yàn)方面存在較大的局限性。此外，還有dropout事件和維數(shù)災(zāi)難等問題。相較于傳統(tǒng)的bulk RNA-seq數(shù)據(jù)，scRNA-seq數(shù)據(jù)中包含大量的dropout event，這使得表達(dá)值為零。在測序過程中，測序技術(shù)的偏差以及轉(zhuǎn)錄組的較低讀數(shù)都會(huì)導(dǎo)致dropout event，在統(tǒng)計(jì)數(shù)據(jù)過程中，低表達(dá)值的基因缺失也會(huì)造成dropout event。這些高水平的技術(shù)噪聲以及統(tǒng)計(jì)噪聲導(dǎo)致無法捕獲準(zhǔn)確的基因表達(dá)，給scRNA-seq研究帶來巨大困難。

為了解決dropout事件和維數(shù)災(zāi)難等對 scRNA-seq 研究帶來的負(fù)面影響，研究者們提出了很多基于重構(gòu)的方法。Wang等［4］利用深層自動(dòng)編碼器和貝葉斯模型，提取來自不同條件之間的基因-基因關(guān)系，以對新數(shù)據(jù)集去噪。David等［5］通過數(shù)據(jù)擴(kuò)散來共享相似細(xì)胞之間的信息，提出MAGIC算法對scRNA-seq數(shù)據(jù)矩陣的缺失值進(jìn)行處理，減輕dropout event造成的影響并提高scRNA-seq數(shù)據(jù)的分析能力。通過考慮細(xì)胞水平上的相關(guān)性，Kwak等［6］提出DrImpute技術(shù)估算dropout event。該方法在區(qū)分丟失零與真實(shí)零的方面具有更好的性能并改善聚類和可視化。Li等［7］通過借用其它相似細(xì)胞中相同基因的信息來估算細(xì)胞中基因的缺失值。這些方法處理單細(xì)胞數(shù)據(jù)時(shí)均取得不錯(cuò)的結(jié)果。Prabhakaran等［8］和Linderman等［9］都利用t-SNE作為分析scRNA-seq數(shù)據(jù)的有力工具。Becht等［10］將UMAP應(yīng)用于scRNA-seq數(shù)據(jù)研究，實(shí)現(xiàn)快速計(jì)算并具有很高的重現(xiàn)性。K-means算法［11-12］是一種基于劃分的聚類算法，把數(shù)據(jù)對象之間的距離作為相似性度量，通常對象之間距離越小越有可能在同一個(gè)簇。Yang等［13］利用迭代K-means聚類對參數(shù)進(jìn)行詳盡搜索找到最優(yōu)參數(shù)。另一種常用的聚類算法是通過計(jì)算不同數(shù)據(jù)點(diǎn)間的相似度來生成聚類簇的層次聚類。Zheng等［14］基于Spectral clustering在相似矩陣上添加低秩和非負(fù)結(jié)構(gòu)，提出了SinNLRR。Wang等［15］基于Spectral clustering提出了SIMLR算法，通過多核學(xué)習(xí)從基因表達(dá)數(shù)據(jù)中學(xué)習(xí)細(xì)胞間距離度量并構(gòu)建相似性矩陣，不僅提高聚類效果并且可以有效地適應(yīng)多個(gè)下游步驟。

為了更好地分析scRNA-seq數(shù)據(jù)中dropout events，獲得能更好表示scRNA-seq數(shù)據(jù)本質(zhì)特征的低維數(shù)據(jù)以及準(zhǔn)確將細(xì)胞聚類，本文提出了一種基于雙自編碼結(jié)合變分貝葉斯和高斯混合模型的聚類方法來分析scRNA-seq 數(shù)據(jù)的方法，稱之為 sc-VBDAE。首先使用對抗自編碼網(wǎng)絡(luò)學(xué)習(xí)數(shù)據(jù)特征，對scRNA-seq數(shù)據(jù)進(jìn)行有效重構(gòu)，去除數(shù)據(jù)中的冗余信息，提高scRNA-seq數(shù)據(jù)的基因表達(dá)能力。其次利用經(jīng)典自編碼對重構(gòu)后的數(shù)據(jù)進(jìn)行降維，獲取低維數(shù)據(jù)從而提高scRNA-seq數(shù)據(jù)分析效率。最后利用變分貝葉斯高斯混合模型揭示 scRNA-seq數(shù)據(jù)內(nèi)部結(jié)構(gòu)，更準(zhǔn)確的聚類細(xì)胞。本文在十個(gè)公開的scRNA-seq數(shù)據(jù)集上測試 sc-VBDAE的性能并與其它方法進(jìn)行比較，結(jié)果表明sc-VBDAE聚類性能略優(yōu)于其它聚類方法。

1 方 法

本文基于雙自編碼網(wǎng)絡(luò)和變分貝葉斯高斯混合聚類，構(gòu)建一個(gè)新的聚類模型 sc-VBDAE。sc-VBDAE主要包括4部分：①scRNA-seq 數(shù)據(jù)預(yù)處理；②對抗自編碼網(wǎng)絡(luò)重構(gòu)基因表達(dá)數(shù)據(jù)；③經(jīng)典自編碼網(wǎng)絡(luò)對重構(gòu)后的數(shù)據(jù)降維；④變分貝葉斯高斯混合模型聚類細(xì)胞。sc-VBDAE具體流程如圖1所示。

sc-VBDAE模型可以分為4個(gè)過程：

1）數(shù)據(jù)預(yù)處理過程。首先，通過基因篩選去除表達(dá)值中0值數(shù)量超過 95％ 的基因。然后對過濾后的數(shù)據(jù)集進(jìn)行l(wèi)og轉(zhuǎn)換。在每個(gè)數(shù)據(jù)集中，列代表細(xì)胞，行代表基因。

2）數(shù)據(jù)重構(gòu)過程。以預(yù)處理后的 scRNA-seq 表達(dá)矩陣輸入對抗自編碼器網(wǎng)絡(luò)，通過編碼器和解碼器的處理過程獲得重構(gòu)后的數(shù)據(jù)。

3）降維過程。將重構(gòu)后的數(shù)據(jù)輸入到經(jīng)典自編碼中。自編碼器由輸入層，3個(gè)隱藏層和輸出層構(gòu)成，通過無監(jiān)督訓(xùn)練有監(jiān)督調(diào)優(yōu)的兩階段方法對網(wǎng)絡(luò)參數(shù)進(jìn)行調(diào)優(yōu)，去除冗余信息后獲得降維數(shù)據(jù)。

4）聚類過程。結(jié)合變分貝葉斯高斯混合模型對細(xì)胞進(jìn)行聚類，使用Bayesian Gaussian Mixture 函數(shù)，最后將聚類結(jié)果通過 t-SNE 可視化。

1.1 自編碼網(wǎng)絡(luò)

自編碼網(wǎng)絡(luò)［16］是處理scRNA-seq數(shù)據(jù)常用的深度學(xué)習(xí)方法。本文對經(jīng)典自編碼進(jìn)行訓(xùn)練，利用瓶頸層神經(jīng)元具有較少個(gè)數(shù)的特點(diǎn)，從而對高維數(shù)據(jù)進(jìn)行有效降維。經(jīng)典自編碼網(wǎng)絡(luò)包括輸入層，隱含層和輸出層。從輸入層到隱含層的過程是編碼過程，從隱含層到輸出層的過程是解碼過程。編碼是將原始scRNA-seq數(shù)據(jù)x∈Rm映射到隱含表示h（x）∈Rn的過程，可以表示為

h（x）=σh（Wx+b）（1）

其中W∈Rn×m為編碼權(quán)值矩陣；b∈Rn為編碼偏置向量；σh（x）為激活函數(shù)。

解碼是將隱含表示h（x）映射到輸出層o，對原始scRNA-seq數(shù)據(jù)x重構(gòu)的過程，可以表示為

o=σo（W′h（x）+b′）（2）

其中：W′∈Rm×n為解碼權(quán)值矩陣；b′∈Rm為解碼偏置向量；σo（x）為激活函數(shù)。

原始數(shù)據(jù)與重構(gòu)數(shù)據(jù)之間的壓縮損失函數(shù)表示為

L=12∑‖y-x‖2（3）

其中：x∈Rm為原始scRNA-seq數(shù)據(jù)；y∈Rm為重構(gòu)數(shù)據(jù)。

隱含層為具有對稱性質(zhì)的互相連接的三層神經(jīng)網(wǎng)絡(luò)。自編碼器的各層輸出函數(shù)可以表示為

h1=σh1（W1x+b1）

hk=σhk（Wkhk-1+bk），k=2，3

o=σo（W4h3+b4）（4）

其中：W1、Wk、W4為相應(yīng)的權(quán)值矩陣；b1、bk、b4為相應(yīng)的偏置。

經(jīng)過從編碼到解碼的過程后，再通過無監(jiān)督訓(xùn)練有監(jiān)督調(diào)優(yōu)的兩階段方法對網(wǎng)絡(luò)參數(shù)進(jìn)行調(diào)優(yōu)。通過兩階段方法調(diào)優(yōu)，提高了自編碼器的學(xué)習(xí)效果，且提高了學(xué)習(xí)速度和泛化性能。

1.2 對抗編碼網(wǎng)絡(luò)

對抗自編碼器［17］是一種正則化自編碼器的新方法，思想是同時(shí)訓(xùn)練兩個(gè)神經(jīng)網(wǎng)絡(luò) （生成器G和判別器D），在它們之間建立一個(gè)最小-最大對抗博弈。生成器G（z）逐步學(xué)會(huì)把樣本z從先驗(yàn)分布p（z）到數(shù)據(jù)空間，在鑒別器D（x）訓(xùn)練數(shù)據(jù)中區(qū)分?jǐn)?shù)據(jù)點(diǎn)之間的空間采樣與實(shí)際數(shù)據(jù)分布和鑒別器產(chǎn)生的數(shù)據(jù)點(diǎn)。假設(shè)訓(xùn)練G（z），利用D（x）相對于x的梯度修改其參數(shù)，使其完全混淆判別器與其生成的樣本。該算法可以形式化為如下類型的極大極小目標(biāo)，見式 （5）：

minGmaxDEX～pdata［logD（X）］+Ez～p（z）［log（1－D（G（z）））］（5）

其中：Pdata為數(shù)據(jù)分布；p（z）為模型分布。

對抗自編碼通過將聚合后驗(yàn)q（z）與任意先驗(yàn)p（z）匹配來實(shí)現(xiàn)。為了做到這一點(diǎn)，在自動(dòng)編碼器的隱藏層向量上附加了一個(gè)對抗網(wǎng)絡(luò)，與此同時(shí)，自動(dòng)編碼器試圖將重構(gòu)誤差最小化。對抗網(wǎng)絡(luò)的產(chǎn)生者也是自動(dòng)編碼器q（z|x）的編碼器。該編碼器保證聚合后驗(yàn)分布能夠欺騙判別對抗網(wǎng)絡(luò)，使其認(rèn)為聚合后驗(yàn)q（z）來自于真實(shí)的先驗(yàn)分布p（z）。且對抗網(wǎng)絡(luò)和對抗自編碼器的訓(xùn)練都是與SGD聯(lián)合進(jìn)行的，在每個(gè)小批量上分別執(zhí)行重構(gòu)階段和正則化階段［18］。在重構(gòu)階段，自動(dòng)編碼器更新編碼器和解碼器，使輸入的重構(gòu)誤差最小化。在正則化階段［18］，對抗網(wǎng)絡(luò)先更新其判別網(wǎng)絡(luò)，以區(qū)分真實(shí)樣本 （使用先驗(yàn)生成） 和生成樣本，然后對抗網(wǎng)絡(luò)更新它的生成器以混淆判別器。

設(shè)x為帶有深度編碼器和解碼器的自動(dòng)編碼器的輸入，z為潛在的代碼向量 （隱藏單位）。設(shè)p（z）為希望施加在碼上的先驗(yàn)分布，q（z|x）為編碼分布，p（x|z）為解碼分布。設(shè)pd（x）為數(shù)據(jù)分布，p（x）為模型分布。對抗自編碼器q（z|x）的編碼函數(shù)定義了q（z）在對抗自編碼器隱藏層向量上的后驗(yàn)聚集分布，如式（6）所示：

q（z）=∫Xq（z|X）pd（X）dX（6）

生成器G和判別器D都可以被搭建成完全連接的神經(jīng)網(wǎng)絡(luò)，然后用一個(gè)合適的優(yōu)化器進(jìn)行反向傳播訓(xùn)練。本文使用了自適應(yīng)矩估計(jì)算法（Adam），這是對隨機(jī)梯度下降的擴(kuò)展。一旦訓(xùn)練過程完成，自動(dòng)編碼器的解碼器將定義生成模型，將施加的先驗(yàn)p（z）映射到數(shù)據(jù)分布。

1.3 高斯混合模型的變分貝葉斯

變分貝葉斯可以看做是期望最大化算法 （EM）的擴(kuò)展，因?yàn)樗彩遣捎脴O大后驗(yàn)估計(jì) （MAP）。另外，變分貝葉斯也通過一組相互依賴 （mutually dependent） 的等式進(jìn)行不斷的迭代來獲得最優(yōu)解。這類實(shí)現(xiàn)了兩種類型的權(quán)重分布的先驗(yàn)：有限混合模型的Dirichlet分布和無限混合模型的Dirichlet過程。在實(shí)踐中，Dirichlet過程推理算法是近似的，并使用具有固定最大組件數(shù)量的截?cái)喾植?。?shí)際使用的聚類數(shù)量幾乎總是取決于數(shù)據(jù)。對于變分貝葉斯高斯混合模型，本文使用 sklearn.mixture 模塊的 Bayesian Gaussian Mixture 函數(shù)，并將參數(shù) n_components 設(shè)置為數(shù)據(jù)集中已知細(xì)胞類型的個(gè)數(shù)，其它參數(shù)默認(rèn)。

1.4 實(shí)驗(yàn)評(píng)價(jià)指標(biāo)

為評(píng)估聚類方法的性能，本研究選擇4個(gè)常用的聚類評(píng)價(jià)指標(biāo)：標(biāo)準(zhǔn)化互信息 （NMI）［19］，調(diào)整后的蘭德指數(shù) （ARI）［20］，Homogeneity［21］和 Completeness［21］。4個(gè)指標(biāo)均是根據(jù)聚類方法得到的預(yù)測標(biāo)簽與數(shù)據(jù)集中提供的真實(shí)標(biāo)簽進(jìn)行計(jì)算得到。

NMI：互信息 （MI） 是通過聚類標(biāo)簽和預(yù)測標(biāo)簽的熵來度量兩個(gè)集合之間相關(guān)程度的指標(biāo)，標(biāo)準(zhǔn)化互信息 （NMI） 通過計(jì)算聚類結(jié)果與真實(shí)劃分之間的差異比率，用于檢測聚類結(jié)果的準(zhǔn)確性。NMI是衡量聚類結(jié)果好壞的常用指標(biāo)之一。通過聚類方法預(yù)測的標(biāo)簽和真實(shí)標(biāo)簽的信息熵以及互信息來計(jì)算 NMI。NMI的范圍是 （0，1），NMI越接1近說明聚類結(jié)果越準(zhǔn)確。假設(shè)共有N個(gè)樣本，U，V分別是預(yù)測結(jié)果標(biāo)簽和真實(shí)標(biāo)簽，NMI 可以表示為

NMI（U，V）=MI（U，V）H（U）H（V）（7）

H（U）=－∑|U|i=1P（i）log（P（i））（8）

H（V）=－∑|V|j=1P（j）log（P（j））（9）

其中H（U）和H（V）分別是U和V的熵。U和V的互信息如式（10）所示：

MI（U，V）=∑|U|i=1∑|V|j=1P（i，j）log（P（i，j）P（i）P（j））（10）

其中：P（i）為樣本屬于Ui的概率；P（j）為樣本屬于Vj的概率；P（i，j）為樣本屬于Ui和Vj的概率。

ARI：蘭德指數(shù) （RI） 通過預(yù)測結(jié)果和真實(shí)聚類中分配在相同或不同簇中的標(biāo)簽對來計(jì)算兩個(gè)聚類之間的相似性。調(diào)整后的蘭德指數(shù) （ARI） 是RI調(diào)整后的指標(biāo)，ARI比RI具有更高的區(qū)分度。ARI也是度量聚類結(jié)果的一個(gè)重要指標(biāo)，ARI與NMI不同之處在于，ARI是比較兩種聚類結(jié)果之間的吻合程度。兩者計(jì)算方式也不同，ARI只需真實(shí)標(biāo)簽和預(yù)測標(biāo)簽。ARI的范圍是 （-1，1），ARI值越大意味著聚類結(jié)果與真實(shí)結(jié)果越吻合。

假設(shè)N是樣本數(shù)量，U，V分別為預(yù)測結(jié)果標(biāo)簽和真實(shí)標(biāo)簽，Nij表示在U中i類和V中j類重疊的個(gè)數(shù)，ai表示出現(xiàn)在U中i類的數(shù)量，bj表示出現(xiàn)在V中j類的數(shù)量。ARI定義如式（11）所示：

ARI=∑ijNij2－∑iai2∑jbj2N2

12∑iai2+∑jbj2－∑iai2∑jbj2N2（11）

Homogeneity：如果聚類結(jié)果中所有的簇都只包含屬于單個(gè)簇的細(xì)胞，則聚類結(jié)果滿足同質(zhì)性。假設(shè)H（V|U）是簇U分配到簇V的細(xì)胞類型條件熵，H（V）是簇V的熵，homogeneity的定義如式（12）所示：

homogeneity=1－H（V|U）H（V）（12）

H（V|U）=－∑|V|v=1∑|U|u=1Nv，uNlog（Nv，uNu）（13）

H（V）=－∑|V|v=1NvNlog（NvN）（14）

其中：N為樣本總數(shù);Nu為屬于簇U的樣本數(shù);Nv為屬于簇V的樣本數(shù);Nv，u為從簇V分配到簇U的樣本數(shù)。

Completeness：如果聚類結(jié)果的簇中所有細(xì)胞都是屬于同一簇，則聚類結(jié)果滿足完整性。completeness 的定義如式（15）所示：

completeness=1－H（U|V）H（U）（15）

同質(zhì)性和完整性都是基于條件熵的互信息分?jǐn)?shù)來衡量簇向量間的相似度，兩者的范圍都是 （0，1），并且數(shù)值越大說明聚類效果越好。需要注意簇標(biāo)簽值的排列不會(huì)更改分?jǐn)?shù)值。

2 實(shí) 驗(yàn)

2.1 數(shù)據(jù)與預(yù)處理

為了評(píng)估sc-VBDAE方法的性能，本文使用了10個(gè)公開的scRNA-seq數(shù)據(jù)集，數(shù)據(jù)集分別來自人類和小鼠的細(xì)胞。所有數(shù)據(jù)集均提供每個(gè)樣本細(xì)胞所屬細(xì)胞類型的高度可信的標(biāo)簽，它們被用來與聚類的預(yù)測標(biāo)簽作對比。數(shù)據(jù)集先通過基因篩選，去除表達(dá)值中0值超過95%的基因。然后對篩選后的數(shù)據(jù)進(jìn)行l(wèi)og轉(zhuǎn)換處理。每個(gè)數(shù)據(jù)集中列代表細(xì)胞，行代表基因。前9個(gè)數(shù)據(jù)集均來自 https：//hemberg-lab.github.io/scRNA.seq.datasets/ 網(wǎng)站。第10個(gè)數(shù)據(jù)集是來自3名COVID-19患者和3名相關(guān)對照者的pbmc的數(shù)據(jù)，該數(shù)據(jù)集可以從BIG data Centre的GSA下載，登錄號(hào)為CRA002390；10個(gè) scRNA-seq 數(shù)據(jù)集的具體信息如表 1 所示。

2.2 對抗自編碼網(wǎng)絡(luò)性能分析

為評(píng)估sc-VBDAE中對抗自編碼網(wǎng)絡(luò)的性能，本文探究了去除sc-VBDAE中對抗自編碼網(wǎng)絡(luò)重構(gòu)數(shù)據(jù) （No auto） 對模型性能的影響，即直接對scRNA-seq數(shù)據(jù)進(jìn)行自編碼降維和變分貝葉斯高斯混合模型對細(xì)胞進(jìn)行聚類。

sc-VBDAE 和No auto 模型輸出得到的 ARI 聚類指標(biāo)如圖2所示。從圖中可以明顯看出，使用對抗自編碼重構(gòu)后的數(shù)據(jù)在8個(gè)數(shù)據(jù)集上可以明顯提升聚類性能，在 Goolam 數(shù)據(jù)集和 Darman 數(shù)據(jù)集上略微提升了聚類性能。

為了進(jìn)一步測試 sc-VBDAE中對抗自編碼網(wǎng)絡(luò)的性能，將有無對抗自編碼網(wǎng)絡(luò)的兩種模型在10個(gè)scRNA-seq 數(shù)據(jù)集上的聚類結(jié)果繪制成基因表達(dá)熱圖。從熱圖可以清晰地看到每種細(xì)胞類型的聚類情況，而且可以得到每種細(xì)胞類型相應(yīng)的標(biāo)記基因。圖3為對比兩種聚類模型的基因表達(dá)熱圖，從圖上可得，雖然 progenitor1 和 neuron2 的標(biāo)記基因基本類似，但是sc-VBDAE 模型的progenitor12 和 mesenchyme 的標(biāo)記基因和 No auto 差距很大。正是不同的標(biāo)記基因?qū)е聝煞N模型的聚類結(jié)果不同，根據(jù)兩種模型的聚類結(jié)果對比顯然 sc-VBDAE 結(jié)果更優(yōu)。這是由于對抗自編碼網(wǎng)絡(luò)重構(gòu) scRNA-seq 數(shù)據(jù)后，減輕了數(shù)據(jù)中的 dropout 事件并提高了標(biāo)記基因的表達(dá)值 （由熱圖可知），更容易根據(jù)標(biāo)記基因聚類細(xì)胞以及進(jìn)行其它下游分析。

通過以上結(jié)果分析可知，對抗自編碼網(wǎng)絡(luò)通過對scRNA-seq數(shù)據(jù)表達(dá)矩陣進(jìn)行重構(gòu)，提高了基因的表達(dá)值，不僅更容易識(shí)別標(biāo)記基因，而且提高聚類的準(zhǔn)確性。sc-VBDAE通過對抗自編碼網(wǎng)絡(luò)對 scRNA-seq 數(shù)據(jù)進(jìn)行重構(gòu)，使數(shù)據(jù)具有更強(qiáng)的特征學(xué)習(xí)能力，進(jìn)而提高了 sc-VBDAE 的數(shù)據(jù)分析能力?？梢?，對抗自編碼網(wǎng)絡(luò)可以準(zhǔn)確分析scRNA-seq 數(shù)據(jù)，對scRNA-seq數(shù)據(jù)集的研究具有重要意義。

2.3 自動(dòng)編碼網(wǎng)絡(luò)的性能分析

為評(píng)估Autoencoder network對模型性能的影響，本文將重構(gòu)后的數(shù)據(jù)集作為輸入，將 sc-VBDAE 中的 Autoencoder network （AE） 分別替換為PCA［22］，t-SNE［23］，UMAP［10］和ZIFA［24］，再對降維后的數(shù)據(jù)進(jìn)行聚類。其中PCA 和 t-SNE 使用 sklearn［25］包中的函數(shù)，PCA 的n_components 參數(shù)與 Autoencoder network 參數(shù)相同，t-SNE 的 perplexity 參數(shù)設(shè)置為樣本中細(xì)胞數(shù)量的 0.2 倍。對于 ZIFA，維數(shù)參數(shù)k與Autoencoder network的參數(shù) n_components 一致。UMAP 使用模型默認(rèn)參數(shù)。5種方法在10個(gè)數(shù)據(jù)集上得到的 ARI 如表2所示。

由表2可以看出，在10個(gè)scRNA-seq數(shù)據(jù)集上，Autoencoder network的 ARI值幾乎均高于其他4種降維方法的指標(biāo)值，這說明 Autoencoder network比另外四種降維方法更有效地捕獲這10個(gè)scRNA-seq數(shù)據(jù)中重要的獨(dú)立特征。sc-VBDAE利用Autoencoder network得到scRNA-seq數(shù)據(jù)中的關(guān)鍵獨(dú)立信息，降低數(shù)據(jù)維數(shù)并減少數(shù)據(jù)冗余。不僅為scRNA-seq數(shù)據(jù)分析提高了效率，而且使聚類結(jié)果更加準(zhǔn)確?？梢妔c-VBDAE結(jié)合Autoencoder network降維scRNA-seq數(shù)據(jù)的性能優(yōu)于另外4種方法。

2.4 聚類性能分析

為測試 sc-VBDAE 的聚類性能，本文將scScope［26］，SIMLR［27］，SNN-cliq［28］，Seurat［29］，scGMAI［30］和sc-VBDAE分別在10個(gè)scRNA-seq數(shù)據(jù)集上運(yùn)行，并分析比較6種模型的4個(gè)度量指標(biāo)。對于sc-VBDAE的變分貝葉斯高斯混合模型，本研究使用sklearn包中的函數(shù)，其中的n_components參數(shù)選擇數(shù)據(jù)集提供的細(xì)胞類型個(gè)數(shù)，其它參數(shù)默認(rèn)。SIMLR，scScope，scGMAI，SNN-cliq和Seurat均使用模型默認(rèn)參數(shù)。6種聚類模型在10個(gè)數(shù)據(jù)集上的ARI結(jié)果如表3所示。

由表3可知，sc-VBDAE在十個(gè)scRNA-seq數(shù)據(jù)集上的ARI指標(biāo)均超過了scScope，SNN-cliq，SiMLR和Seurat 4種聚類方法。除了Biase數(shù)據(jù)集，Camp2數(shù)據(jù)集，Goolam數(shù)據(jù)集和Klein數(shù)據(jù)集外，其余6個(gè)數(shù)據(jù)集的ARI指標(biāo)也均超過scGMAI方法，特別的，在Darmanis數(shù)據(jù)集ARI指標(biāo)值比其他模型高0.49～0.06，在Deng數(shù)據(jù)集ARI指標(biāo)值比其他模型高0.61～0.07，在Baron1數(shù)據(jù)集ARI指標(biāo)值比其他模型高0.38～0.09，這說明sc-VBDAE的聚類性能優(yōu)于現(xiàn)存的聚類模型。而且sc-VBDAE的NMI，homogeneity和completeness也優(yōu)于其它聚類方法。通過四種性能度量指標(biāo)說明sc-VBDAE的聚類更加準(zhǔn)確，與真實(shí)情況更加吻合。

相比于其它聚類方法，變分貝葉斯高斯混合聚類是一種基于統(tǒng)計(jì)的聚類模型，而且采用后驗(yàn)概率可以比先驗(yàn)概率更準(zhǔn)確地模擬數(shù)據(jù)的分布。通過統(tǒng)計(jì)方法計(jì)算細(xì)胞簇之間的分布結(jié)構(gòu)，揭示scRNA-seq數(shù)據(jù)內(nèi)部性質(zhì)及規(guī)律并更準(zhǔn)確的聚類細(xì)胞?？梢?，sc-VBDAE是從scRNA-seq數(shù)據(jù)中準(zhǔn)確聚類細(xì)胞并識(shí)別細(xì)胞類型的有力工具。

2.5 細(xì)胞軌跡的推斷

除了通過聚類描述細(xì)胞外，scRNA-seq還可以通過時(shí)間進(jìn)程或發(fā)育階段 （即細(xì)胞軌跡） 來描述細(xì)胞。一般來說，軌跡分析首先降低scRNA-seq數(shù)據(jù)集的維數(shù)，然后推測細(xì)胞的分化軌跡，最后將每個(gè)細(xì)胞投射到該軌跡的適當(dāng)位置。盡管單細(xì)胞實(shí)驗(yàn)可以闡明各種生物環(huán)境中的軌跡，但沒有一種單細(xì)胞軌跡推斷方法可以解釋dropout事件。人們推測，在重構(gòu)后的scRNA-seq數(shù)據(jù)上推斷細(xì)胞軌跡可以提高偽時(shí)間分析的準(zhǔn)確性。

圖4利用了細(xì)胞從NPC分化到GW21+3時(shí)間過程的scRNA-seq數(shù)據(jù)，然后使用Mococle 包重建細(xì)胞分化軌跡。從圖中可以明顯看出和聚類結(jié)果對應(yīng)的分化軌跡。結(jié)果表明sc-VBDAE可以很好地捕獲scRNA-seq數(shù)據(jù)的主要特征并準(zhǔn)確聚類細(xì)胞，有助于模擬分析單細(xì)胞數(shù)據(jù)分化軌跡和恢復(fù)基因表達(dá)的時(shí)間動(dòng)態(tài)。

3 結(jié) 論

隨著越來越多單細(xì)胞RNA測序技術(shù) （scRNA-seq） 的研發(fā)，允許從稀有細(xì)胞或者難以獲得的細(xì)胞中獲取基因表達(dá)信息，在單個(gè)細(xì)胞水平上揭示更多未知疾病來源以及其它生物學(xué)問題。但是很多挑戰(zhàn)同時(shí)存在，比如維數(shù)災(zāi)難，dropout event以及更準(zhǔn)確的細(xì)胞聚類。

針對 scRNA-seq 研究中遇到的挑戰(zhàn)，本文基于對抗自編碼網(wǎng)絡(luò)和經(jīng)典自編碼并結(jié)合變分貝葉斯高斯混合聚類，提出了一種新的聚類 scRNA-seq數(shù)據(jù)的模型，稱為sc-VBDAE。對抗自編碼網(wǎng)絡(luò)的編碼和解碼過程可以去除冗余信息，學(xué)習(xí)scRNA-seq數(shù)據(jù)特征并減輕 dropout events，提高scRNA-seq數(shù)據(jù)的分析效率。sc-VBDAE 利用對抗自編碼網(wǎng)絡(luò)重構(gòu)數(shù)據(jù)，不僅提高聚類結(jié)果和可視化性能，而且更準(zhǔn)確的識(shí)別細(xì)胞簇的標(biāo)記基因，為scRNA-seq下游分析提供強(qiáng)有力幫助。scRNA-seq數(shù)據(jù)具有既相互獨(dú)立又相互聯(lián)系的特點(diǎn)，自編碼網(wǎng)絡(luò)可以根據(jù)數(shù)據(jù)之間的相互聯(lián)系選擇scRNA-seq數(shù)據(jù)中具有關(guān)鍵獨(dú)立性的特征，形成代表數(shù)據(jù)本質(zhì)特征的潛在低維空間。并且自編碼網(wǎng)絡(luò)具有快速計(jì)算和處理海量數(shù)據(jù)的能力，從而提高了scRNA-seq分析效率。變分貝葉斯高斯混合模型使用概率模型描述聚類原型，可以很好地模擬scRNA-seq數(shù)據(jù)分布，更準(zhǔn)確的揭示scRNA-seq數(shù)據(jù)內(nèi)在性質(zhì)及規(guī)律。通過比較6種聚類方法在10個(gè)scRNA-seq數(shù)據(jù)集上的性能度量指標(biāo)。結(jié)果顯示sc-VBDAE的聚類性能優(yōu)于其它5種scRNA-seq聚類方法。

特別的，本文首次使用對抗自編碼網(wǎng)絡(luò)對scRNA-seq數(shù)據(jù)進(jìn)行重構(gòu)，不僅提高了模型聚類scRNA-seq數(shù)據(jù)的精度，而且為 scRNA-seq 數(shù)據(jù)得到以及其它生物領(lǐng)域研究提供新的方法。盡管sc-VBDAE可以有效聚類并鑒定 scRNA-seq 數(shù)據(jù)中細(xì)胞類型，但仍然存在一定的提升空間。下一步我們將會(huì)使用更高效的深度學(xué)習(xí)方法得到 scRNA-seq 數(shù)據(jù)信息，進(jìn)一步提高 scRNA-seq 數(shù)據(jù)的聚類精度。

參 考 文 獻(xiàn)：

［1］ ISLAM S， KJLLQUIST U， MOLINER A， et al. Highly Multiplexed and Strand-specific Single-cell RNA 5′ End Sequencing［J］. Nature Protocols， 2012， 7（5）：813.

［2］ PICELLI S， BJORKLUND K， FARIDANI O R， et al. Smart-seq2 for Sensitive Full-length Transcriptome Profiling in Single Cells［J］. Nature Methods， 2013， 10（11）：1096.

［3］ MACOSKO E Z， BASU A， SATIJA R， et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets［J］. Cell， 2015， 161（5）：1202.

［4］ WANG Jingshu， AGARWAL D， Huang Mo， et al. Data Denoising with Transfer Learning in Single-cell Transcriptomics［J］. Nature Methods， 2019， 16（9）：875.

［5］ DAVID V D， ROSHAN S， JUOZAS N， et al. Recovering Gene Interactions from Single-Cell Data Using Data Diffusion［J］. Social Science Electronic Publishing， 2018：S0092867418307244.

［6］ KWAK I Y， GONG Wuming， KOYANO-NAKAGAWA N， et al. DrImpute： Imputing Dropout Events in Single Cell RNA Sequencing Data［J］. Cold Spring Harbor Laboratory， 2017， 19（1）：220.

［7］ LI W V， LI J J. An Accurate and Robust Imputation Method ScImpute for Single-cell RNA-seq Data［J］. Nature Communications， 2018， 9（1）：1.

［8］ PRABHAKARAN S， AZIZI E， CARR A， et al. Dirichlet Process Mixture Model for Correcting Technical Variation in Single-cell Gene Expression Data［C］. International Conference on International Conference on Machine Learning. JMLR.org， 2016： 1070.

［9］ LINDERMAN G C， RACHHM， HOSKINS J G， et al. Fast Interpolation-based t-SNE for Improved Visualization of Single-cell RNA-seq Data［J］. Nature Methods， 2019， 16（3）：243.

［10］BECHT E， MCINNES L， HEALY J， et al. Dimensionality Reduction for Visualizing Single-cell Data Using UMAP［J］. Nature Biotechnology， 2018， 37（1）：38.

［11］HARTIGAN J A， WONG M A. A k-means Clustering Algorithm. Applied Statistics［J］. Algorithms， 1978， 1326（28）：100.

［12］王寧， 陳晨， 陳德運(yùn)， 等. 哼唱檢索中旋律特征的聚類與優(yōu)化方法［J］. 哈爾濱理工大學(xué)學(xué)報(bào)，2022，27（1）：61.

WANG Ning， CHEN Chen， CHEN Deyuan， et al. Melody Feature Clustering and Optimization for Query-by-humming［J］. Journal of Harbin University of Science and Technology，2022，27（1）：61.

［13］YANG Lu， LIU Jiancheng， LU Qiang， et al. SAIC： An Iterative Clustering Approach for Analysis of Single Cell RNA-seq Data［J］. BMC Genomics， 2017， 18（S6）：689.

［14］ZHENG Ruiqing， LI Min， LIANG Zhenlan， et al. SinNLRR： A Robust Subspace Clustering Method for Cell Type Detection by Non-negative and Low-rank Representation［J］. Bioinformatics， 2019， 35（19）：3642.

［15］WANG Bo， ZHU Junjie， PIERSON E， et al. Visualization and Analysis of Single-cell RNA-seq Data by Kernel-based Similarity Learning［J］. Nature Methods， 2017， 14（4）：414.

［16］ERASLAN G， SIMON L M， MIRCEA M， et al. Single-cell RNA-seq Denoising Using a Deep Count Autoencoder［J］. Nature Communications， 2019， 10（1）：1.

［17］DING Jiarui， CONDON A， SHAH S P. Interpretable Dimensionality Reduction of Single Cell Transcriptome Data with Deep Generative Models［J］. Cold Spring Harbor Laboratory， 2017， 9（1）：1.

［18］TIAN Yingjie， ZHANG Quqi. A Comprehensive Survey on Regularization Strategies in Machine Learning［J］. Information Fusion， 2022， 80：146.

［19］ARJOVSKY M， CHINTALA S， BOTTON L. Wasserstein GAN［J］. arXiv， 2017， doi：1701.07875.

［20］GULRAJANI I， AHMED F， ARJOVSKY M， et al. Improved Training of Wasserstein Gans［J］. Machine Learning， 2017： 5767.

［21］RAO Jiahua， ZHOU Xiang， LU Yutong， et al. Imputing Single-cell RNA-seq Data by Combining Graph Convolution and Autoencoder Neural Networks［J］. iScience， 2021： 102393.

［22］WOLD S， ESBENSEN K， GELADI P， et al. Principal Component Analysis［J］. Chemometrics & Intelligent Laboratory Systems， 1987， 2（1/3）：37.

［23］MAATEN L， HINTON G. Visualizing Data Using t-SNE［J］. Journal of Machine Learning Research， 2008， 9（2605）：2579.

［24］PIERSON E， YAU C. ZIFA： Dimensionality Reduction for Zero-inflated Single-cell Gene Expression Analysis［J］. Genome Biology， 2015， 16（1）：241.

［25］LIN Peijie， TROUP M， HO J W. CIDR： Ultrafast and Accurate Clustering Through Imputation for Single-cell RNA-seq Data［J］. Genome Biology， 2017， 18（1）： 59.

［26］DENG Yue， BAO Feng， DAI Qionghai， et al. Scalable Analysis of Cell-type Composition from Single-cell Transcriptomics Using Deep Recurrent Learning［J］. Nature Methods， 2019， 16（4）： 311.

［27］WANG Bo， ZHU Junjie， PIERSION E， et al. Visualization and Analysis of Single-cell RNA-seq Data by Kernel-based Similarity Learning［J］. Nature Methods， 2017， 14（4）： 414.

［28］XU Chen， SU Zhengchang. Identification of Cell Types From Single-cell Transcriptomes Using a Novel Clustering Method［J］. Bioinformatics， 2015， 31（12）：1974.

［29］BULER A， HOFFMAN P， Smibert P， et al. Integrating Single-cell Transcriptomic Data Across Different Conditions， Technologies， and Species［J］. Nature Biotechnology， 2018， 36（5）： 411.

［30］YU Bin， CHEN Chen， QI Ren， et al. scGMAI： a Gaussian Mixture Model for Clustering Single-cell RNA-Seq Data Based on Deep Autoencoder［J］. Briefings in Bioinformatics， 2021， 22（4）：1.

（編輯：溫澤宇）