





摘 要:單細胞轉錄組測序(scRNA-seq,single cell RNA sequencing)技術為單個細胞高通量、高分辨率的深入研究提供了機會,為在單細胞層面研究細胞功能及其背后的基因調控機制提供了重要技術手段。然而這項技術也帶來新的挑戰,單細胞數據具有規模大、噪聲高、異構性強等特點,特別是高比例的數據缺失(dropout)嚴重影響了下游分析的可靠性,甚至掩蓋了基因與基因間的重要關系。這里提出一種基于負二項分布的分治插補策略ND-Impute(Negative binomial distribution based Divide and conquer strategy for imputation)對scRNA-seq數據進行處理,該方法假設scRNA-seq數據符合負二項分布,利用包含特定損失函數的自動編碼器獲取數據的特異性參數,并使用分治策略估計潛在的基因表達值。通過聚類效果、相關性和誤差分析等比較,表明該方法可以有效地恢復缺失數據,提高了后續研究分析的準確性。
關鍵詞:單細胞轉錄組測序;數據缺失;插補策略;聚類分析
中圖分類號:O211.9 文獻標識碼:A 文章編號:2095-414X(2023)01-0014-07
0 "引言
自然界中的生物是復雜且多樣的,對生物學表征的研究,離不開對細胞的深入分析,無論是單細胞生物還是多細胞生物,單個細胞之間的差異都會產生不可估量的影響。隨著細胞分離技術以及高通量測序技術的發展,bulk RNA-seq技術逐漸被scRNA-seq技術所取代,scRNA-seq技術的出現為單個細胞的高吞吐量和高分辨率轉錄組分析提供了新途徑[1-2]。自2009年首個scRNA-seq技術發布以來,scRNA-seq被越來越廣泛地應用于基礎科學研究中,這項技術尤其在腫瘤學[3]、遺傳病學[4]、免疫學[5]等生物醫學研究中發揮了重要作用。單細胞轉錄組測序的主要步驟是分離出單個細胞、將mRNA轉為cDNA、cDNA擴增、高通量測序和數據分析等,測序所得數據可以使研究人員進一步了解異質性細胞類型、細胞發育譜系以及疾病發展狀態,對于發現新的異質細胞類型和追蹤細胞動態發育軌跡也有巨大幫助[6]。然而由于技術限制和生物因素,scRNA-seq數據比bulk RNA-seq數據更復雜,噪聲更大,目前scRNA-seq的通量已經由原來的幾個細胞發展到成千上萬個細胞,通量的增加提高了數據的復雜性。即使在現在的單細胞轉錄組測序研究中,數據噪聲問題依然普遍,尤其是缺失(dropout)事件,dropout的產生源于單細胞轉錄組的測序過程,因此難以避免[7]。
scRNA-seq數據由于其生物學特性和測序技術因素產生噪聲。首先在生物學方面,細胞的類型、大小和細胞的表達方式各不相同;其次對于測序技術方面,在使用細胞中表達量較低的RNA構建cDNA庫的過程中,導致一些轉錄組的丟失是必然的,cDNA庫的擴增會引入更多的隨機變異。同時,測序效率低會使每個細胞中的轉錄組只有很小一部分被成功測序。上述原因就導致每個細胞基因表達量中有很大比例的零值和低計數值,其中很多原本低表達的基因數據被忽略,測得結果為零,也就是數據缺失dropout,然而其中有部分零值反映的是真正的零表達量,這極大混淆了低表達基因和不表達基因,從而影響下游分析的準確性和可信度[8]。
為了解決scRNA-seq數據中存在大量數據缺失的問題,研究人員開發了許多插補方法,這些方法大致可以分為三類:統計學方法、深度學習方法和集成學習方法。早期的一些方法從scRNA-seq數據的統計學分布特性出發,利用數學的方法和統計學思維解決數據中零值過多的問題。其中最經典的統計學方法是2018年的SAVER[9],該方法利用基因與基因的關系來恢復每個細胞中每個基因的表達水平,通過具有Poisson-LASSO回歸的經驗貝葉斯方法估算先驗參數,構建負二項分布模型(NB,Negative binomial distribution),最后輸出后驗分布的均值作為插補結果。但是這種方法會改變所有基因的表達水平,包括那些不受 dropout 影響的基因表達水平,因此需要考慮基因表達數據中零值為dropout的概率。
基于對dropout概率的分析,Wei[10]等人在2018年提出scImpute,該方法基于混合模型計算每個表達值是dropout的概率,并利用相似細胞的相同基因的表達信息來估算dropout項的潛在表達值。但基于統計學的方法只能處理很小的數據量,所以目前越來越多出現和使用的插補算法都使用到了深度學習方法?;谏疃葘W習的方法不僅能夠處理大規模數據,還能利用深度學習中的神經網絡自主學習scRNA-seq數據中的一些特征,達到比基于統計學的方法更優秀的插補效果[11]。比較有代表性的方法有Eraslan G[12]等在2019年提出的基于深度自動編碼器的數據插補去噪方法DCA,該方法采用負二項噪聲模型,考慮了數據的計數分布和稀疏性,通過改進的自動編碼器學習給定數據集中的特殊參數,此方法與細胞的數量成線性比例,可以應用于數百萬細胞的數據集,但該方法著重考慮的是基因表達數據的統計分布情況,而忽略了基因表達數據的特異性。Arisdakessian[13]等提出了DeepImpute,一種基于深度神經網絡的插補算法,該方法最大的特點就是利用了分治法的思想,把scRNA-seq數據分為多個子集后,使用dropout層和loss函數來學習數據中的模式,從而實現準確的插補,與DCA相反,該方法沒有考慮到數據原本的統計分布,直接使用了深度學習方法,可能會產生更多的噪聲。2020年Xu[14]等提出了把生成對抗網絡 (GANs)用于scRNA-seq數據插補的方法scIGANs,它使用生成的細胞而不是觀察到的細胞來避免這些限制,并平衡主要和稀有細胞群之間的每個特性。該方法的主要思路首先是將基因表達矩陣轉換成圖片,具體規則為一張圖片代表一個細胞,圖片中的像素代表各個基因的基因表達水平,然后將轉換的圖片輸入GANs模型中進行訓練,最終獲得插補結果。但將生成對抗網絡運用到基因表達數據的圖片形式,數據生成的圖片與真實圖片還是有所差距。深度學習方法相對于統計學方法,同樣存在缺陷,例如模型選擇不佳、網絡訓練過擬合等問題,可能會造成更多噪聲的引入。隨著插補方法的不斷衍生,scTSSR[15]和EnTSSR[16]這樣的集成學習方法通過對比加權,實現了更全面的分析,然而要消耗大量的計算和內存資源。
在上述的工作中,各種方法都有切入面,面對基因表達數據,首先希望能找到一個合適的統計分布模型來規范數據,其次面對數據的大規模增長,需要引入深度學習方法處理巨大的計算量,而集成學習方法對資源的消耗十分巨大,因此,仍然需要一種高效的插補方法對scRNA-seq數據進行插補。自動編碼器作為一種無監督的神經網絡模型,它分為編碼器和解碼器兩個模塊,使用編碼器將原數據進行壓縮,學習數據中的特征,再使用解碼器將數據還原。自動編碼器是一個無監督的學習過程,它可以學習到輸入數據的隱含特征。
本文旨在使用深度學習方法,利用自動編碼器的特點來學習和處理scRNA-seq數據,并最終達到較好的缺失數據插補效果。這里提出一種基于負二項分布的分治插補策略(ND-Impute),以用于scRNA-seq數據的插補,該方法在兩種基于深度神經網絡的方法DCA和DeepImpute的分析思維上加以改進,既考慮了單細胞轉錄組測序數據的特殊統計學分布模型,也考慮到了基因表達數據的特異性。大量的實驗結果表明,ND-Impute在關鍵性能指標,即聚類效果、相關性和誤差分析方面優于其他相似的方法。
1 "工作介紹
1.1 "數據集介紹
本文采用五個不同大小的scRNA-seq數據集(68KPBMC、293T Cells、Jurkat Cells、CD34+ cells和CD19+ B cells)來測試ND-Impute方法對dropout事件的插補能力。其中,68KPBMC[17]數據集是外周血細胞數據集,293T cells(10X Genomic)是人體上皮細胞,Jurkat cells (10X Genomic)是某種T淋巴細胞,CD34+ cells[18]是人類骨髓造血干細胞,CD19+ B cells是免疫細胞的一種。這些公開數據集的大小如表1所示。
1.2 "負二項分布模型
首先,合理假設數據的分布情況有助于下游分析,scRNA-seq數據的統計學分布情況比較特殊,其具有高稀疏性(high sparsity),表現為數據中存在大量零值(zero inflation),且數據也是高度離散的[18-19]。因此對于基于計數讀取的scRNA-seq數據,描述連續型數據的正態分布不符合要求;其次泊松分布適合于描述單位時間(或空間)內隨機事件發生的次數(事件發生的次數只能是離散的整數),但由于均值始終與方差相等,故無法體現scRNA-seq數據的高離散度特性。這里,從方差大于均值的角度出發,負二項分布是可行的。首先,二項分布描述的是在 重伯努利試驗中,事件B恰好發生 次的概率,以二項分布 為例, 次試驗中正好得到 次成功的概率密度函數為:
(1)
負二項分布描述的也是伯努利試驗,這個分布描述的事件是當某個結果出現固定次數時,整個過程的數量分布。設試驗持續到 ( 為整數)次失敗, 表示一個事件在該伯努利試驗中每次出現的概率,那么其負二項分布的概率密度函數為:
(2)
該分布的均值和方差分別為:
(3)
(4)
進一步可以得到:
(5)
從式(5)可以看出,方差隨著均值增加呈現二次函數形式的遞增,因此符合scRNA-seq數據的特點。同時,經過學者們的研究,負二項分布被證明是一種可以合理描述scRNA-seq數據的模型,并且該分布模型已被一些經典的插補方法所使用的,如SAVER,scVI等[20]。
2 "方法
ND-Impute算法的工作流程如圖1所示,其主要包括兩個模塊:(1)數據集讀取,獲取數據集符合負二項分布模型的均值;(2)采用分治法構建訓練網絡,訓練數據集,并得到插補結果。圖中輸入基因表達矩陣,矩陣中顏色越深說明表達量越高,顏色越淺說明表達量越低。
2.1 "數據預處理
使用改進的自動編碼器實現對數據中基于負二項分布參數的學習。自動編碼器使用三個層來分析預測輸入數據:輸入層、隱藏層和輸出層。這種改進自動編碼器的特殊之處在于采用特定的噪聲模型,這種噪聲模型融合零膨脹負二項分布,構建零膨脹負二項分布損失函數,以對應scRNA-seq數據高稀疏性和高離散性的特點。零膨脹負二項分布是由負二項分量的均值和離散參數( 和 )以及表示點是dropout概率( )的混合系數來參數化的:
(6)
(7)
對于基因表達矩陣中的每一個值,與傳統的自動編碼器不同,該自動編碼器框架估算每個輸入值的三個負二項分布參數:均值( )、dropout概率 ( )和離散度( )。因此,一個輸入矩陣會得到三個輸出層,三個輸出層的大小都與輸入一致,即具有相同數量的表達量。通過該自動編碼器的處理,就得出了輸入數據符合零膨脹負二項分布的均值。其中,這一模型框架默認有三個隱藏層,包含64、32、64個神經元,32個神經元的隱藏層為瓶頸層。如圖2所示,使scRNA-seq數據的復雜性和非線性能夠被有效捕捉。
2.2 "分治法構建并訓練網絡
為了保持scRNA-seq數據中的細胞特異性,并有效考慮到細胞與基因之間的關系,這里進一步使用深度神經網絡來處理并訓練由自動編碼器獲取的均值。采用分治法的思想,將數據集中所有的基因分割成子集,并構建子網絡對數據進行訓練。
網絡構建首先要選擇“目標基因”,將數據集中的所有基因分成具有相同基因數的子集作為目標基因,這些目標基因構成了子神經網絡的輸出層。對于子集 中的每個目標基因 ,從子集 以外的基因中選擇5個相似度最高的預測基因,選擇完成后刪除在所有子神經網絡中選擇次數超過閾值的預測基因,剩余的預測基因為子神經網絡的輸入層。對于每個子集,訓練一個四層神經網絡,如圖1和圖3所示,預測基因組成的輸入層、完全連接的隱藏層、dropout層(與scRNA-Seq數據中的dropout不同,這里的dropout是按比例忽略一部分數據)和目標基因組成的輸出層。其中使用線性整流函數(ReLU)作為激活函數,使用加權MSE作為損失函數,加入dropout層是為了防止過擬合,其中dropout參數設為20%。
將95%的細胞作為訓練集進行網絡訓練,剩下的5%作為測試集,并行訓練每個子模型,在每個網絡訓練完成后測試是否產生過擬合情況。由于每個子網絡的結構較簡單,可以觀察到超參數調整導致的可變性較低。因此,將批次大小默認設置為64,學習率設置為0.0001。計算損失值時,使用加權均方誤差(MSE)損失函數,加權操作可以賦予表達值較高的基因以更高的權重,從而加強基因的表達。針對一個細胞 ,其損失值 計算如下:
(8)
其中 是細胞 的基因 的值, 為該基因訓練過程中的給定值。
3 "結果與分析
3.1 "聚類分析
在scRNA-seq數據的分析實驗中,聚類分析是必不可少的,聚類可以幫助確定數據中的細胞類型或亞型。本研究采取的數據聚類算法為Scanpy,隨機選取2000個細胞進行分析。
在聚類分析的評估中,使用本研究方法ND-Impute與其他兩種方法深度學習DCA和DeepImpute進行對比。從圖4中可以看出,ND-Impute相比其他兩種方法,聚類簇之間邊界更加清晰,能相對顯著區分出數據中的細胞類型。原始數據插補后,能在一定程度上恢復因測序技術造成的低表達量缺損,從而恢復或增強了細胞的特異性表達,使數據中的細胞類型或細胞亞群能夠被識別出來。為了評估ND-Impute在細胞聚類上的有效性,使用兩個聚類指標,歸一化互信息(NMI)和調整蘭德指數(ARI)。其中NMI用于度量兩個聚類結果的相近程度,值域為 ,數值越高說明越接近,即劃分越準確。ARI能反映兩種劃分的重疊程度,取值范圍為 ,值越大說明效果越好。
設 表示聚類(cluster)劃分, 表示實際類別劃分。
NMI定義為:
(9)
其中 表示 和 互信息:
(10)
其中, 為 的聯合分布, 分別為邊緣分布,可以理解為樣本中同時屬于兩個集合的概率,以及分別屬于 的概率。 , 分別表示 的熵:
(11)
(12)
設 為通過聚類方法正確分配到同一簇的細胞對數, 是錯誤地分配到同一簇中的單元對數, 為錯誤分配到不同簇的細胞對數, 為正確分配到不同簇的細胞對數。
ARI定義為:
(13)
三種方法在五個數據集上的NMI和ARI表現如表2所示。對比可得,ND-Impute在兩個指標的表現上要優于DCA和DeepImpute。細胞數量大的數據集相比細胞數量小的數據集可以更好地訓練網絡,所以指標結果都相對表現更佳。
3.2 "插補效果分析
對于獲取的scRNA-seq數據集,由于無法確定真實的缺失值,這里采用了一種隨機掩蔽部分scRNA-seq數據集基因表達矩陣中元素的方法,也就是用零代替表達矩陣中的部分非零值,比較插補后的數值是否能較好地還原被掩蔽的數值。選擇隨機掩碼表達式矩陣中每個細胞的非零計數值的5%。掩蔽遵循與之成比例的密度概率分布:
(14)
其中 是原始計數值。原始值可作為評價插補方法性能的參考值,使用皮爾遜相關系數和均方誤差(MSE)兩種性能指標評估插補效果,以及恢復掩蔽值的準確性。圖5和圖6中分別顯示了ND-Impute、DCA和DeepImpute的皮爾遜相關系數和均方誤差對比情況,可以看出ND-Impute獲取了更高的皮爾遜相關系數和更低的均方誤差,這說明相比DCA和DeepImpute,ND-Impute更好地恢復了所有范圍內的dropout值。
恢復掩蔽值的準確性如圖7所示,插補值與真實值越接近,圖中的點就會越接近斜線,并且分布更均勻。
4 "結論
scRNA-seq數據中的dropout事件是生物信息學分析中的重要問題,其數據的龐大以及數據的嘈雜都不利于后續的生物學研究,傳統的插補策略多使用單一思路解決這一問題。
本研究中提出了一種基于負二項分布的分治插補策略ND-Impute來處理scRNA-seq數據中的dropout。ND-Impute基于深度學習方法,在常規自
動編碼器基礎上使用符合統計學分布的特異性損失函數,并將分治思想融入神經網絡的訓練中。這種思路和結構有效地還原了缺失的數值,既考慮了數據集的分布模型,也通過學習細胞與基因之間的關系和相似性保留了細胞的特異性表達。與另外兩種深度學習方法不同,ND-Impute不會只針對數據的統計學分布或細胞特異性,該方法同時考慮這兩個方面的信息,以達到更好插補效果。通過實驗比較,ND-Impute顯示出比另外兩種深度學習方法更好的聚類效果,NMI和ARI的平均值高于DCA和DeepImpute,對于掩蔽值的恢復,不僅具有更高的皮爾遜相關系數,造成的誤差值也更低,并且插補值與真實值也更加接近。綜上所述,ND-Impute在插補效果上相比一般的深度插補模型都有所優化,是一種表現良好的插補方法。
我們下一步研究將運用其他深度學習方法來學習scRNA-seq數據的特性,以達到更好的插補效果。同時將生成模型與統計學分布相結合,進一步提高數據的特異性。
參考文獻:
[1] ShapiroE,Biezuner T, Linnarsson S. Single-cell sequencing-based technologies will revolutionize whole-organism science[J]. Nat Rev Genet, 2013, 14(09):618–630.
[2] Tang F, Barbacioru C, Wang Y, et al. mRNA-Seq whole-transcriptome analysis of a single cell[J]. Nature Methods, 2009, 6(05):377-382.
[3] Zhu D, Zhao Z, Cui G, et al. Single-cell transcriptome analysis reveals estrogen signaling coordinately augments one-carbon, polyamine, and purine synthesis in breast cancer[J]. Cell Rep, 2018, 25 (08): 2285-2298.e4.
[4] Chen, Chongyi, Xing, et al. Single-cell whole-genome analyses by Linear Amplification via Transposon Insertion (LIANTI)[J]. Science, 2017.14;356:189-194.
[5] Crinier A, Milpied P, Escalière B, et al. High-dimensional single-cell analysis identifies organ-specific signatures and conserved NK cell subsets in humans and mice[J]. Immunity, 2018, 49 (05):971-986.e5.
[6] Kolodziejczyk A, Kim J K, Svensson V, et al. The technology and biology of single-cell RNA sequencing.[J]. Molecular Cell, 2015, 58(04):610-620.
[7] Andrews T S, Hemberg M. M3Drop: dropout-based feature selection for scRNASeq[J]. Bioinformatics, 2019, 35(16):2865-2867.
[8] Zappia L, Phipson B, Oshlack A . Splatter: simulation of single-cell RNA sequencing data[J]. Genome Biology, 2017, 18(01):174.
[9] Mo H, Wang J, Torre E, et al. SAVER: gene expression recovery for single-cell RNA sequencing[J]. Nature Methods, 2018, 15(05):539-542.
[10] Wei V L, Li J J. An accurate and robust imputation method scImpute for single-cell RNA-seq data[J]. Nature Communications, 2018, 9(01):997.
[11] Ma Q, Xu D. Deep learning shapes single-cell data analysis[J]. Nature Reviews Molecular Cell Biology, 2022, 23(05): 303-304.
[12] Eraslan G, Simon L M, Mir Ce A M, et al. Single-cell RNA-seq denoising using a deep count autoencoder[J]. Nature Communications, 2019, 10:390.
[13] Arisdakessian C, Poirion O, Yunits B, et al. DeepImpute: an accurate, fast, and scalable deep neural network method to impute single-cell RNA-seq data[J]. Genome Biology, 2019, 20(01):211.
[14] Xu Y, Zhang Z, You L, et al. scIGANs: single-cell RNA-seq imputation using generative adversarial networks[J]. Nucleic Acids Research, 2020, 48(04):e85.
[15] Jin K, Le O Y, Zhao X M, et al. scTSSR: gene expression recovery for single-cell RNA sequencing using two-side sparse self-representation[J]. Bioinformatics, 2020, 36(10).
[16] Lu F, Lin Y, Yuan C, et al. EnTSSR: a weighted ensemble learning method to impute single-cell RNA sequencing data[J]. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2021, 18(06): 2781-2787.
[17] Zheng G X Y, Terry J M, Belgrader P, et al. Massively parallel digital transcriptional profiling of single cells[J]. Nature communications, 2017, 8(01): 1-12.
[18] Risso D, Perraudeau F, Gribkova S, et al. A general and flexible method for signal extraction from single-cell RNA-seq data[J]. Nature communications, 2018, 9(01): 1-17.
[19] Lopez R, Regier J, Cole M, et al. Bayesian inference for a generative model of transcriptome profiles from single-cell RNA sequencing[J]. bioRxiv, 2018: 292037.
[20] Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression[J]. Genome biology, 2019, 20(01): 1-15.
Divide and Conquer Imputation for Dropouts in Single-Cell Data based on Negative Binomial Distribution
XIONGZhen-zhena, b, ZHANG Ben-gongb, c
(a. School of Computer Science and Artificial Intelligence, b. Research Center of Nonlinear Science, c. School of Mathematical and Physical Sciences, Wuhan Textile University, Wuhan Hubei 430200, China)
Abstract:Single-cell RNA sequencing (scRNA-seq) technology provides opportunities for high-throughput, high-resolution in-depth research of single cells, and provides insights into cell functions and the underlying gene regulation mechanisms at the single-cell level important technical means. However, this technology also brings new challenges. ScRNA-seq data has the characteristics of large scale, high noise, and strong heterogeneity, especially the high proportion of data missing, which is called dropout. The problem of dropout seriously affects the reliability of the downstream analysis, and even covers up the important relationship between genes and genes. This paper proposed a divided and conquering imputation strategy based on negative binomial distribution ND-Impute to process scRNA-seq data. This method assumed that scRNA-seq data conform to the negative binomial distribution, utilized an autoencoder that incorporates a specific loss function to obtain data-specific parameters, and used a divide-and-conquer strategy to estimate potential gene expression values. The comparison of clustering effect, correlation, and error analysis showed that this method can effectively restore missing data and improve the accuracy of subsequent research and analysis.
Keywords:single-cell RNA sequencing; dropout; imputation; clustering analysis
(責任編輯:周莉)