朱賢振 李雪玲

摘要：?jiǎn)渭?xì)胞 RNA 測(cè)序已成為研究生物學(xué)重要特征的強(qiáng)大高分辨率工具。然而，其測(cè)序條件苛刻，價(jià)格成本高昂。目前細(xì)胞類型反卷積能夠很好地解決這些限制問題，SMCTD（Sparse Model Cell Type Deconvolution） 使用稀疏自編碼器優(yōu)化TAPE（Tissue-AdaPtive autoEncoder） ，使其在直腸癌和PBMC模擬數(shù)據(jù)上預(yù)測(cè)細(xì)胞類型比列具有更高的靈敏度、準(zhǔn)確性和整體性能，同時(shí)在預(yù)測(cè)細(xì)胞類型特異性基因表達(dá)上表現(xiàn)更優(yōu)。

關(guān)鍵詞：?jiǎn)渭?xì)胞測(cè)序；細(xì)胞類型反卷積；深度學(xué)習(xí)；稀疏自編碼器；一致性相關(guān)系數(shù)

中圖分類號(hào)：TP311 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2024）11-0009-04

隨著二代深度測(cè)序（NGS） 、單細(xì)胞轉(zhuǎn)錄組測(cè)序（scRNA-seq） 、空間轉(zhuǎn)錄組（Spatial Transcriptomics） 技術(shù)、細(xì)胞類型反卷積算法的發(fā)展，為整合單細(xì)胞測(cè)序數(shù)據(jù)和大規(guī)模 bulk 基因表達(dá)譜，研究癌癥微環(huán)境中的細(xì)胞組成和基因表達(dá)，提供了重要手段。然而，現(xiàn)有的細(xì)胞類型反卷積方法[1]的準(zhǔn)確度和解析顆粒度有很大提升空間，發(fā)展基于深度機(jī)器學(xué)習(xí)的細(xì)胞類型反卷積算法，能夠?qū)崿F(xiàn)并加速大范圍內(nèi)高通量臨床數(shù)據(jù)的精確分析[2]。

TAPE[3]是一種連接批量RNA-seq和單細(xì)胞RNAseq[4]的深度學(xué)習(xí)方法，可在短時(shí)間內(nèi)實(shí)現(xiàn)精確的反卷積。通過構(gòu)建可解釋的解碼器并在獨(dú)特的方案下進(jìn)行訓(xùn)練，TAPE 可以自適應(yīng)地預(yù)測(cè)細(xì)胞類型分?jǐn)?shù)和細(xì)胞類型特異性基因表達(dá)。與多個(gè)數(shù)據(jù)集上的流行方法相比，TAPE 在細(xì)胞類型水平上具有更好的整體性能和相當(dāng)?shù)臏?zhǔn)確性。此外，它在不同細(xì)胞類型中更穩(wěn)健、更快、更靈敏，可以提供具有生物學(xué)意義的預(yù)測(cè)。然而TAPE是基于傳統(tǒng)自編碼器的深度學(xué)習(xí)模型，需要對(duì)所有數(shù)據(jù)進(jìn)行處理降維，容易過擬合，盡管傳統(tǒng)自編碼器可以學(xué)習(xí)數(shù)據(jù)的特征表示，但并不保證這些特征是有意義、可解釋的。模型可能會(huì)學(xué)習(xí)到捕捉數(shù)據(jù)中的一些噪聲或冗余信息，而不是真正有用的特征。

稀疏自編碼器相比自編碼器特征的學(xué)習(xí)與選擇能力更好，稀疏自編碼器通過強(qiáng)制編碼層的神經(jīng)元保持較低的激活率，鼓勵(lì)網(wǎng)絡(luò)只激活最重要的特征神經(jīng)元[5]。這可以使更有意義、更具區(qū)分性的特征被學(xué)習(xí)和保留，從而提高模型對(duì)數(shù)據(jù)的表征能力。也能降低過擬合的風(fēng)險(xiǎn)： 通過在編碼層施加稀疏性約束，稀疏自編碼器可以降低過擬合的風(fēng)險(xiǎn)。這是因?yàn)槟Ｐ捅黄葍H僅選擇最重要的特征，而不會(huì)過度適應(yīng)訓(xùn)練數(shù)據(jù)中的噪聲或不重要的變化。更好的泛化能力： 由于稀疏自編碼器傾向于學(xué)習(xí)更有意義的特征，它們通常能夠更好地泛化到新的、未見過的數(shù)據(jù)，從而提高模型的泛化能力。還擁有更高效的特征表示：稀疏自編碼器可以學(xué)習(xí)數(shù)據(jù)的更緊湊、更高級(jí)的表示。這種表示可以更好地捕捉數(shù)據(jù)中的關(guān)鍵模式和結(jié)構(gòu)，從而在后續(xù)的任務(wù)（如分類、聚類等）中表現(xiàn)更好。

1 細(xì)胞類型反卷積模型構(gòu)建

細(xì)胞類型反卷積是能夠?qū)⒋笠?guī)模bulk基因表達(dá)譜反卷積得到bulk中細(xì)胞類型比例和細(xì)胞類型特異性基因的一種方法，傳統(tǒng)細(xì)胞類型反卷積需要制作簽名矩陣，簽名矩陣是基于一組特定基因的表達(dá)模式構(gòu)建的。這些特定基因通常與某個(gè)生物學(xué)特性、狀態(tài)或功能相關(guān)聯(lián)。通過分析大量單細(xì)胞數(shù)據(jù)，可以從中提取出這些特定基因的表達(dá)模式，并將這些模式整合到一個(gè)矩陣中，即簽名矩陣[1]。非常煩瑣，最近幾年興起的深度學(xué)習(xí)細(xì)胞類型反卷積能夠自主學(xué)習(xí)bulk基因表達(dá)譜中特征，無須構(gòu)建簽名矩陣即可進(jìn)行反卷積。這些方法中，TAPE是反卷積性能最出色的之一，比同樣使用深度學(xué)習(xí)的Scaden [2]更穩(wěn)健，與傳統(tǒng)細(xì)胞類型反卷積方法CIBERSORTx [1]性能旗鼓相當(dāng)。而TAPE 使用自編碼器的缺點(diǎn)也讓我們思考用更優(yōu)秀的網(wǎng)絡(luò)構(gòu)建模型來提升性能。

1.1 模型設(shè)置

如圖1所示，我們引入稀疏自編碼器，與傳統(tǒng)自編碼器相比，引入稀疏性目標(biāo)值和稀疏性權(quán)重，稀疏性目標(biāo)值是一個(gè)預(yù)設(shè)的小數(shù)，表示隱藏層神經(jīng)元的期望平均激活度，通常是0.05或者更小。稀疏性權(quán)重是一個(gè)常數(shù)，表示稀疏性懲罰項(xiàng)在總代價(jià)函數(shù)中的權(quán)重，用來控制稀疏性的強(qiáng)度。圖1中B 表示輸入的bulk基因表達(dá)數(shù)據(jù)，C 表示經(jīng)過稀疏自編碼器編碼再解碼重現(xiàn)的bulk基因表達(dá)數(shù)據(jù)，X 表示通過編碼器得到的細(xì)胞類型比例。圖1左半部分表示為編碼器，是一個(gè)回歸模型，負(fù)責(zé)將高維bulk基因表達(dá)數(shù)據(jù)映射到低維的細(xì)胞類型比例數(shù)據(jù)。相反，圖1右半部分可以根據(jù)細(xì)胞類型比例數(shù)據(jù)重建bulk基因表達(dá)數(shù)據(jù)。

1.2 模型訓(xùn)練

我們先預(yù)設(shè)稀疏性目標(biāo)值和稀疏性權(quán)重，稀疏性目標(biāo)值代表期望的神經(jīng)元平均激活度，稀疏性權(quán)重用來控制稀疏性強(qiáng)度。然后在代碼中定義了KL散度函數(shù)（Kullback-Leibler Divergence） ，使用KL散度函數(shù)和稀疏性目標(biāo)值和權(quán)重計(jì)算稀疏性懲罰損失[5]。

我們使用大約5 000個(gè)bulk樣本進(jìn)行訓(xùn)練。使用預(yù)測(cè)細(xì)胞類型比例和真實(shí)細(xì)胞類型比例之間的MAE （平均絕對(duì)誤差）與稀疏性懲罰損失的和來優(yōu)化編碼器的參數(shù)，并使用重構(gòu)bulk數(shù)據(jù)和原始bulk數(shù)據(jù)之間的MAE 與稀疏性懲罰損失的和來優(yōu)化解碼器和編碼器。

1.3 模型預(yù)測(cè)

使用反卷積模型進(jìn)行預(yù)測(cè)需要預(yù)先準(zhǔn)備單細(xì)胞參考數(shù)據(jù)，行為細(xì)胞類型，列為基因名稱，文件為TXT 格式。為了使反卷積結(jié)果更精準(zhǔn)更具有生物學(xué)意義，單細(xì)胞參考數(shù)據(jù)要與需要預(yù)測(cè)的bulk 數(shù)據(jù)為同一組織的，并且擁有相同的細(xì)胞類型。

需要預(yù)測(cè)的bulk 數(shù)據(jù)需要指定分隔符，行為樣本名稱，列為基因名稱，數(shù)據(jù)類型最好為“counts”，若使用“TPM”或者“FPKM”格式須自備基因長(zhǎng)度文件使數(shù)據(jù)最終以counts格式運(yùn)行在程序上。

我們把準(zhǔn)備好的單細(xì)胞參考數(shù)據(jù)和bulk數(shù)據(jù)輸入程序，將模式選擇為“overall”，然后選擇合適的數(shù)據(jù)類型及基因長(zhǎng)度文件，自適應(yīng)參數(shù)選擇為“True”或者“False”。如果是“ True”，那么它將會(huì)預(yù)測(cè)輸出簽名矩陣，反之，則返回空值。等待程序運(yùn)行完成，會(huì)得到預(yù)測(cè)的細(xì)胞類型比例數(shù)據(jù)（行為樣本，列為細(xì)胞類型）及選擇可得到的簽名矩陣。

2 細(xì)胞類型反卷積模型性能比較

由于公共數(shù)據(jù)庫(kù)中同一樣本中既測(cè)bulk數(shù)據(jù)又測(cè)單細(xì)胞數(shù)據(jù)的少之又少，所以為了精準(zhǔn)測(cè)出反卷積模型的性能，因此有必要進(jìn)行偽bulk數(shù)據(jù)測(cè)試進(jìn)行估計(jì)。偽bulk數(shù)據(jù)是通過具有基本事實(shí)（預(yù)定義的細(xì)胞類型比例）的單細(xì) 胞基因表達(dá)數(shù)據(jù)在計(jì)算機(jī)中生成的。也就是說，偽bulk數(shù)據(jù)是許多單細(xì)胞基因表達(dá)數(shù)據(jù)的總和。我們將使用TAPE中的偽bulk模擬程序模擬bulk。

2.1 GSE176078單細(xì)胞數(shù)據(jù)模擬預(yù)測(cè)比較

首先，我們預(yù)設(shè)真實(shí)細(xì)胞類型比例，再將從GEO 中下載的乳腺癌數(shù)據(jù)集GSE176078 [6]的樣本作為參考的單細(xì)胞數(shù)據(jù)生成模擬bulk對(duì)五種細(xì)胞類型進(jìn)行反卷積性能預(yù)測(cè)。設(shè)定兩個(gè)參數(shù)指標(biāo)：MAE（平均絕對(duì)誤差）和CCC（一致性相關(guān)系數(shù)）[7]，MAE是對(duì)每個(gè)細(xì)胞類型的預(yù)測(cè)值與其對(duì)應(yīng)的實(shí)際值之間的絕對(duì)差值進(jìn)行求和，然后對(duì)所有數(shù)據(jù)點(diǎn)的絕對(duì)差值求平均值，數(shù)值越小性能越好。CCC是評(píng)價(jià)細(xì)胞類型比例預(yù)測(cè)值與真實(shí)值之間的一致性的指標(biāo)，CCC值越接近1代表性能越好。

最終結(jié)果如圖2所示，我們的模型SMCTD（稀疏自編碼器）在誤差方面要比TAPE和Scaden都低，在CCC方面也要比TAPE更出色，略遜于Scaden，綜合兩方面來看，SMCTD是三者中反卷積性能最出色的。

接下來測(cè)試模型在細(xì)胞類型增加的情況下的性能，我們繼續(xù)用GSE176078模擬bulk，這次bulk包含14 種細(xì)胞亞型，例如：Monocyte、Fibroblasts、NK cells 等。評(píng)價(jià)指標(biāo)同上。

如圖3所示，結(jié)果表明，在這種情況下，所有方法都出現(xiàn)性能下降的情況，但這些方法的MAE與預(yù)測(cè)五種細(xì)胞類型情形中的 MAE相當(dāng)，這表明這些方法可以預(yù)測(cè)接近真實(shí)值的值。同時(shí)SMCTD也是這種情況下誤差最低的算法，其CCC值也為三種方法第二高，說明其性能是三者間最出色的。

2.2 PBMC 單細(xì)胞數(shù)據(jù)模擬預(yù)測(cè)比較

接下來我們使用10X Genomics官網(wǎng)的PBMC（外周血單個(gè)核細(xì)胞）單細(xì)胞數(shù)據(jù)[8]模擬bulk（其中包含七種細(xì)胞類型）進(jìn)行反卷積，評(píng)價(jià)指標(biāo)同上。

最終結(jié)果如圖4 所示，可以發(fā)現(xiàn)在PBMC 數(shù)據(jù)上Scaden 表現(xiàn)是最出色的，我們的模型SMCTD 雖然比Scaden 略差，但要比TAPE在MAE、CCC值兩個(gè)方面都要更優(yōu)秀。

綜上所述，我們發(fā)現(xiàn)模型SMCTD在預(yù)測(cè)細(xì)胞類型方面有著不遜于TAPE、Scaden的性能，雖然這三種方法在細(xì)胞類型增多的情況下都會(huì)出現(xiàn)性能下降，這種情況是可以預(yù)見的，所以解決這一問題也是未來研究的方向之一。

2.3 組織適應(yīng)性細(xì)胞類型特異性基因表達(dá)預(yù)測(cè)

SMCTD不僅可以預(yù)測(cè)細(xì)胞類型分?jǐn)?shù)，同樣可以自適應(yīng)地預(yù)測(cè)細(xì)胞類型特異性基因表達(dá)。也就是說，SMCTD 只需要模擬數(shù)據(jù)來訓(xùn)練，如果給出相應(yīng)的bulk RNAseq數(shù)據(jù)，它可以預(yù)測(cè)細(xì)胞類型特異性的基因表達(dá)。此功能使SMCTD 能夠剖析不同細(xì)胞類型中的bulk基因表達(dá)，并發(fā)現(xiàn)不同細(xì)胞類型中一些潛在差異表達(dá)的基因。

我們測(cè)試了預(yù)測(cè)的細(xì)胞類型特異性 GEP 的正確性。為了測(cè)試這一點(diǎn)，我們測(cè)量了每種細(xì)胞類型的預(yù)測(cè)基因表達(dá)值與從單細(xì)胞RNA-seq獲得的原始基因表達(dá)值之間的一致性（圖5） 。這里，bulk 數(shù)據(jù)用GSE176078單細(xì)胞數(shù)據(jù)模擬生成，而單細(xì)胞數(shù)據(jù)是乳腺癌癌單數(shù)據(jù)集GSE176078。由于在訓(xùn)練階段使用Log2 和 MinMaxScaler（） 將輸入的 RNA-seq 數(shù)據(jù)轉(zhuǎn)換為 0-1 值，因此按細(xì)胞類型分組的基因表達(dá)值的總和也以這種方式轉(zhuǎn)換以進(jìn)行比較與預(yù)測(cè)的相對(duì)基因表達(dá)值。

由圖5、圖6和表1可得，SMCTD在四種細(xì)胞類型上的基因表達(dá)預(yù)測(cè)一致性相關(guān)系數(shù)都要高于TAPE，僅在免疫細(xì)胞中表現(xiàn)不佳，考慮到其他四種細(xì)胞類型的良好一致性，這種失真可能是由個(gè)體差異引起的。圖表中顯示的一致性證明SMCTD正確預(yù)測(cè)了細(xì)胞類型特異性基因表達(dá)，為進(jìn)一步的基因表達(dá)分析奠定了基礎(chǔ)。

3 結(jié)束語(yǔ)

本文優(yōu)化了一個(gè)細(xì)胞類型反卷積模型TAPE，使用了稀疏自編碼器作為模型基礎(chǔ)，增強(qiáng)了模型的稀疏性，提高了模型性能，使細(xì)胞類型反卷積模型的誤差降低，且提高了相關(guān)性，在我們測(cè)試的兩類數(shù)據(jù)和多種細(xì)胞類型上均能體現(xiàn)出。此外，在預(yù)測(cè)細(xì)胞類型特異性基因表達(dá)上SMCTD也比TAPE在大多數(shù)細(xì)胞類型上的結(jié)果更準(zhǔn)確。

但模型在細(xì)胞類型過多的情況下表現(xiàn)下降，這是未來要攻克的方向之一。在預(yù)測(cè)細(xì)胞類型特異性基因表達(dá)方面，SMCTD同樣有一定的優(yōu)化上升空間。

參考文獻(xiàn)：

[1] NEWMAN A M， STEEN C B， LIU C L， et al. Determining cell type abundance and expression from bulk tissues with digital cytometry [J]. Nat Biotechnol， 2019， 37（7）： 773.

[2] MENDEN K， MAROUF M， OLLER S， et al. Deep learningbased cell composition analysis from tissue expression profiles[J]. Science Advances，2020，6（30）：eaba2619.

[3] CHEN Y S，WANG Y X，CHEN Y L，et al. Deep autoencoder for interpretable tissue-adaptive deconvolution and cell-typespecific gene analysis[J]. Nature Communications，2022，13（1）：6735.

[4] WANG Z，GERSTEIN M，SNYDER M. RNA-Seq：a revolutionary tool for transcriptomics[J]. Nature Reviews Genetics，2009，10（1）：57-63.

[6] WU S Z，AL-ERYANI G，RODEN D L，et al. A single-cell and spatially resolved atlas of human breast cancers[J]. Nature Genetics，2021，53（9）：1334-1347.

[5] NG A. Sparse autoencoder[R]. CS294A Lecture Notes， 2011：72.

[6] WU S Z，AL-ERYANI G，RODEN D L，et al. A single-cell and spatially resolved atlas of human breast cancers[J]. Nature Genetics，2021，53（9）：1334-1347.

[7] LIN L I. A concordance correlation coefficient to evaluate reproducibility[J]. Biometrics，1989，45（1）：255-268.

[8] CHEMISTRY） K P F A H D V. single cell gene expression dataset by cell ranger 2. 1. 0 [DS]. 10X Genomics， 2017，

【通聯(lián)編輯：李雅琪】