基于單細胞RNA 測序數(shù)據(jù)的基因調(diào)控網(wǎng)絡(luò)推斷算法綜述

張少強，潘鏡伊

（天津師范大學計算機與信息工程學院，天津 300387）

基因調(diào)控網(wǎng)絡(luò)（gene regulatory networks，GRN）[1]由影響生物體生物過程的調(diào)控因子間的相互作用關(guān)系組成，是不同轉(zhuǎn)錄因子（transcription factors，TF）或輔因子（co-factors）與其靶基因或靶轉(zhuǎn)錄本之間相互作用的圖譜.GRN 通常表示為連接調(diào)控組件與目標組件的圖，通過網(wǎng)絡(luò)建模，能夠表示各組件如何相互作用以及它們可能產(chǎn)生的行為.雖然基因間相互作用是動態(tài)的，但當前許多生物學研究側(cè)重于基因的定性分析或基因?qū)Γ╣ene pair）間的相互作用，即研究一個基因擾動如何影響另一個基因表達.因此，通過基因表達的變化可以推斷GRN；通過分析不同細胞類型在不同生物過程中的基因表達關(guān)系可以研究基因間的動態(tài)相互作用.了解每個基因的作用及其與其他基因的關(guān)系是在分子水平上解釋生物過程的關(guān)鍵.因此，重構(gòu)和分析GRN 已成為研究基因型驅(qū)動表型潛在機制的有效方法，并有望應(yīng)用于醫(yī)學和藥物設(shè)計等領(lǐng)域[2].

GRN 具有2 個重要的特征[3]：①網(wǎng)絡(luò)拓撲結(jié)構(gòu)，它表示網(wǎng)絡(luò)連接模式.某個基因的邏輯結(jié)構(gòu)可以表示為元素是0 或1 的向量，元素的取值代表該基因與其他基因的關(guān)系.②組件之間的調(diào)控作用，即相互作用的強度和類型，該作用由一組權(quán)值來描述，權(quán)值為正表示基因的正調(diào)控（基因被激活），權(quán)值為負表示基因的負調(diào)控（基因被抑制），權(quán)值為0 表示無調(diào)控關(guān)系.

單細胞RNA 測序（single-cell RNA sequencing，scRNA-seq）為推斷細胞周期或分化等時間依賴性生物過程的GRN 提供了新的可能，并且基于scRNAseq 產(chǎn)生了大量的GRN 推斷算法[4].與批量RNA 測序（bulk RNA-seq）數(shù)據(jù)不同，scRNA-seq 沒有混合所有細胞樣本的基因表達而掩蓋特定細胞的生物信號；另外，scRNA-seq 不僅能夠提供靜態(tài)的細胞快照數(shù)據(jù)，也能夠提供細胞不同時序（time-series）的動態(tài)數(shù)據(jù).因此scRNA-seq 數(shù)據(jù)更適合于GRN 推斷[4].本文綜合分析了26 種基于scRNA-seq 數(shù)據(jù)的GRN 推斷算法，其中包括3 種早期基于批量RNA 測序數(shù)據(jù)的設(shè)計方法，但也可以用于scRNA-seq 數(shù)據(jù).根據(jù)不同的網(wǎng)絡(luò)構(gòu)造方式和模型，另外23 種GRN 推斷算法大致分為6類：①基于布爾網(wǎng)絡(luò)（Boolean network）的算法2 種[5-6]；②基于微分方程（differential equations）的算法3 種[7-9]；③基于偽時序（pseudo-time-series）基因相關(guān)性集成策略的算法5 種[10-14]；④基于共表達（co-expression）基因的算法4 種[15-18]；⑤基于細胞特異性（cell specific）的算法3 種[19-21]；⑥基于深度學習的算法6 種[22-27].詳細描述了每類算法的工作機制或方法原理、適用性、主要優(yōu)勢及局限性等，對相關(guān)的算法比較研究進行分析，并利用scRNA-seq 數(shù)據(jù)簡單驗證了這些推斷算法重建網(wǎng)絡(luò)的精確度，最后分析了當前GRN 推斷算法面臨的主要局限和突出挑戰(zhàn)，并建議新的研究思路.

1 GRN 推斷算法

1.1 早期基于批量RNA 測序數(shù)據(jù)的算法

基于批量測序數(shù)據(jù)重構(gòu)GRN 的算法主要包括3種，分別為WGCNA（weighted correlation network analysis）[28]、GENIE3（GEne network inference with ensemble of trees）[29]和GRNBoost2[30].WGCNA 使用層次聚類和動態(tài)樹分割識別高度相關(guān)（共表達）的基因簇（模塊），利用模塊特征基因或模塊內(nèi)的中心基因表征這些模塊，此外，將模塊內(nèi)基因相連和因果關(guān)系與外部信息（如基因功能語義、功能富集和蛋白質(zhì)組學等）進行關(guān)聯(lián)來尋找模塊中的關(guān)鍵驅(qū)動基因.GENIE3 是一種基于特征選擇和回歸樹集成（tree-based ensemble）的GRN推斷算法.GENIE3 算法的基本思想是將p 個基因的GRN 問題分解為p 個不同的回歸問題，目標是確定網(wǎng)絡(luò)中目標基因的調(diào)控因子，在每個回歸問題中，通過樹集成的方法（如隨機森林），基于所有其他基因的表達模式預(yù)測得到目標基因的表達模式.GRNBoost2是基于GENIE3 架構(gòu)的GRN 推斷算法，是GENIE3的快速替代方案，可擴展至數(shù)萬樣本規(guī)模的數(shù)據(jù)集.與GENIE3 類似，GRNBoost2 訓練一個回歸模型為數(shù)據(jù)集中的每個基因選擇最重要的調(diào)控因子，并利用啟發(fā)式早停正則化的隨機梯度提升機（gradient boosting machine，GBM）回歸來推斷網(wǎng)絡(luò).

WGCNA 利用數(shù)據(jù)探索性工具或基因篩選（排序）構(gòu)造加權(quán)和未加權(quán)的無向網(wǎng)絡(luò)，GENIE3 和GRNBoost2 構(gòu)造允許存在反饋回路的有向網(wǎng)絡(luò).WGCNA 通過增加可測試的假設(shè)（如某模塊與某疾病相關(guān)），以便在獨立的數(shù)據(jù)集中進行驗證，而GENIE3 沒有對基因調(diào)控的性質(zhì)做任何假設(shè)，可以潛在地處理組合調(diào)控和非線性關(guān)系.GENIE3 和GRNBoost2 在已知大量基因表達的情況下運行良好，計算速度快，且可擴展，WGCNA 沒有考慮輸入數(shù)據(jù)是否被正確預(yù)處理或規(guī)范化，所以有可能導致輸出結(jié)果存在偏差或無效.此外，WGCNA 雖然實現(xiàn)了幾種共表達模塊的檢測方法，但沒有明確其中的最好方法.GENIE3 在基因排序中沒有考慮樹的質(zhì)量，所有樹模型的權(quán)重取值相同.這3種基于批量測序數(shù)據(jù)的GRN 推斷算法也可以用于單細胞測序數(shù)據(jù)，但在各種單細胞實驗中性能表現(xiàn)不佳，其主要原因是單細胞基因表達數(shù)據(jù)存在大量缺失（dropout）零值，零值率大大高于批量樣本數(shù)據(jù).相關(guān)實驗結(jié)果也表明，在沒有dropout 的模擬單細胞數(shù)據(jù)的實驗中，GENIE3 和GRNBoost2 的表現(xiàn)優(yōu)于很多通用方法，但在有dropout 的模擬數(shù)據(jù)的實驗中，其表現(xiàn)較差[31].

1.2 基于scRNA-seq 數(shù)據(jù)的算法

單細胞測序技術(shù)能夠在細胞分化過程中追蹤細胞譜系并識別新的細胞類型.此外，單細胞測序技術(shù)使得控制細胞分化并驅(qū)動其細胞類型轉(zhuǎn)變的GRN 成為可能.但與批量數(shù)據(jù)相比，scRNA-seq 較大的噪聲、細胞異質(zhì)性、細胞間測序深度的差異、dropout 引起的基因表達的高稀疏性和細胞周期相關(guān)效應(yīng)等都為GRN 推斷帶來新的挑戰(zhàn).

基于scRNA-seq 數(shù)據(jù)的GRN 推斷算法的總體工作流程如圖1 所示.輸入為某scRNA-seq 數(shù)據(jù)集的基因表達矩陣，行和列分別代表基因和細胞.首先進行基因過濾，將分析范圍縮小到具有高度可變性的基因或用戶感興趣的基因（預(yù)定義基因）；然后根據(jù)數(shù)據(jù)假設(shè)建模和構(gòu)造中間數(shù)據(jù)（圖1 給出了6 類建模方法的示意圖）；最后推斷網(wǎng)絡(luò)，輸出結(jié)果可以是無向的共表達網(wǎng)絡(luò)，也可以是具有基因間調(diào)控關(guān)系的有向網(wǎng)絡(luò).

圖1 基于scRNA-seq 數(shù)據(jù)的GRN 推斷算法的整體工作流程Fig.1 Overall workflow of GRN inference algorithm based on scRNA-seq data

1.2.1 基于布爾網(wǎng)絡(luò)的算法

布爾網(wǎng)絡(luò)模型是最簡單的網(wǎng)絡(luò)推斷方法之一，一個布爾網(wǎng)絡(luò)模型由n 個基因x1，…，xn和n 個基因關(guān)聯(lián)更新函數(shù)f1，…，fn：{0，1}n→{0，1}構(gòu)成.每個基因取一個布爾值x∈{0，1}表示基因表達是否存在；每個更新函數(shù)fi用布爾邏輯表示，使用布爾算子AND、OR 或NOT 指定基因x1，…，xn間的關(guān)系.布爾網(wǎng)絡(luò)的GRN推斷算法一般分為3 步：①對基因表達數(shù)據(jù)進行二值化，生成初始布爾狀態(tài)；②用布爾狀態(tài)優(yōu)化器進行模型狀態(tài)優(yōu)化；③輸出具有激活邊（edge）、抑制邊和一組布爾函數(shù)的GRN.應(yīng)用布爾網(wǎng)絡(luò)推斷基因間關(guān)系的算法主要有SCNS（single cell network synthesis）[5]和BTR（BoolTraineR）[6].

SCNS 以跨時間段（時序）單細胞基因表達數(shù)據(jù)作為輸入，通過構(gòu)建布爾網(wǎng)絡(luò)模型搜索從初始細胞狀態(tài)到后期細胞狀態(tài)的進展和轉(zhuǎn)化的邏輯規(guī)則.SCNS 產(chǎn)生的邏輯模型有助于預(yù)測基因擾動（如基因敲除或基因過表達）對特定譜系的影響.對于幾千個細胞的數(shù)據(jù)集，SCNS 在普通計算機上幾分鐘內(nèi)即可重建布爾網(wǎng)絡(luò)模型.對于較大規(guī)模的數(shù)據(jù)集，SCNS 支持云計算和并行計算，能夠輕松部署到云或高性能計算服務(wù)器上.BTR 算法是一種異步布爾模型，它使用新型布爾狀態(tài)空間評分函數(shù)訓練單細胞表達數(shù)據(jù).BTR 通過改進模型預(yù)測與表達式數(shù)據(jù)的匹配來細化現(xiàn)有的布爾模型，并重構(gòu)新的布爾模型.BTR 對細胞間表達狀態(tài)的關(guān)系不做任何假設(shè)，這有助于在條件盡可能少的情況下重構(gòu)GRN.BTR 也可用于優(yōu)化現(xiàn)有假定的調(diào)控網(wǎng)絡(luò).

BTR 將網(wǎng)絡(luò)作為圖對象輸出，而SCNS 只輸出優(yōu)化后的布爾規(guī)則.SCNS 依賴已知細胞狀態(tài)的一般軌跡，這要求單細胞表達數(shù)據(jù)至少來自2 種具有已知關(guān)系的細胞類型，而BTR 不需要細胞狀態(tài)的軌跡信息來推斷網(wǎng)絡(luò)結(jié)構(gòu)和布爾規(guī)則.布爾網(wǎng)絡(luò)對dropout 具有一定的魯棒性，但是布爾網(wǎng)絡(luò)將基因表達信息二值化會降低網(wǎng)絡(luò)預(yù)測的準確性.此外，這類算法需要通過狀態(tài)轉(zhuǎn)換優(yōu)化布爾規(guī)則，分析中包含的基因越多，構(gòu)建連接狀態(tài)轉(zhuǎn)換圖所需的細胞就越多.

1.2.2 基于微分方程的算法

基于微分方程的算法將基因表達動力學描述為與其他基因表達和環(huán)境因子相關(guān)的函數(shù).在給定一組參數(shù)的情況下，微分方程可以表達數(shù)據(jù)隨時間的變化率.使用微分方程進行網(wǎng)絡(luò)推斷需要細胞的時間順序，因此這類算法適用于推斷細胞分化期scRNA-seq 數(shù)據(jù)的GRN.這類算法主要包括Inference Snapshot[7]、SCODE[8]和SCOUP（single-cell expression data during differentiation with Ornstein-Uhlenbeck process）[9].基于微分方程算法的工作流程如圖2 所示，分為4 步：①利用外部軟件或內(nèi)嵌函數(shù)推斷細胞的偽時序；②使用微分方程描述基因與時序間的關(guān)系；③使用不同的優(yōu)化技術(shù)估計參數(shù)；④利用優(yōu)化的參數(shù)通過微分方程推斷基因間的關(guān)系，最后輸出基因親和度矩陣作為GRN.

圖2 基于微分方程的GRN 推斷算法的整體工作流程Fig.2 Overall workflow of GRN inference algorithm based on differential equations

Inference Snapshot 利用模型選擇與參數(shù)優(yōu)化相結(jié)合的方法構(gòu)建解釋偽時序的常微分方程.在構(gòu)建常微分方程模型的過程中，為減少模型數(shù)量，首先用GENIE3 生成粗糙的GRN 結(jié)構(gòu)作為先驗知識，然后重建分化通路中的基因表達動態(tài)，并推斷出關(guān)鍵的GRN 結(jié)構(gòu).SCODE 以微分方程對表達數(shù)據(jù)的最優(yōu)擬合為目標，優(yōu)化轉(zhuǎn)錄因子間調(diào)控關(guān)系的矩陣，通過整合線性微分方程變換和線性回歸實現(xiàn)矩陣優(yōu)化.SCODE 只使用少量轉(zhuǎn)錄因子即可重建表達動力學，從而顯著降低了時間復雜度.SCOUP 使用隨機微分方程直接描述細胞分化程度和細胞生命整個分化過程中的基因表達動力學.SCOUP 首先構(gòu)造從初始細胞到所有其他細胞的最小生成樹以獲取細胞的偽時序；然后將每個基因的表達隨偽時序的分化建模為一個連續(xù)的隨機擴散過程，即Ornstein-Uhlenbeck（OU）過程，其中單個基因在某一時間點的表達通過當前OU 過程的正態(tài)分布來估計；最后通過計算所有細胞的Z 值得到基因間的相關(guān)性.

這3 種基于微分方程的GRN 推斷算法存在一些技術(shù)差異.在數(shù)據(jù)降維和構(gòu)建偽時序的過程中，Inference Snapshot 通過非線性降維方法diffusion map[32]降維輸入數(shù)據(jù)后調(diào)用Wanderlust 算法[33]生成細胞偽時序；SCODE 使用軌跡推斷方法Monocle[34]降維數(shù)據(jù)并尋找最小生成樹軌跡來構(gòu)造細胞偽時序；SCOUP 使用主成分分析方法進行數(shù)據(jù)降維并采用Prim 算法[35]尋找最小生成樹的最短路徑作為細胞偽時序.另外，這3 種算法中只有SCODE 假設(shè)基因間的依賴關(guān)系是線性的.在算法輸出結(jié)果中，SCODE 和SCOUP 輸出代表基因關(guān)系的相關(guān)矩陣，并可以調(diào)用第三方網(wǎng)絡(luò)可視化應(yīng)用以實現(xiàn)可視化，而Inference Snapshot 只輸出模型的估計參數(shù).

由于微分方程可以模擬RNA、蛋白質(zhì)等在GRN中與時間有關(guān)的各種相互作用的參數(shù)變化，并可以推斷基因間的某些因果關(guān)系，因此比其他方法更接近實際網(wǎng)絡(luò)和細節(jié).然而，由于微分方程的推斷算法需要輸入大量數(shù)據(jù)來估計參數(shù)且計算量較大，因此只能應(yīng)用于有限數(shù)量的基因調(diào)控關(guān)系.

1.2.3 基于偽時序基因相關(guān)性集成策略的算法

與基于微分方程的算法一樣，基于偽時序基因相關(guān)性集成策略的推斷算法也需要推斷細胞的偽時序或者需要用戶提供細胞時序，通過時序基因表達數(shù)據(jù)推斷出動態(tài)GRN[36].這類算法假設(shè)基因間的關(guān)系隨著細胞的不同發(fā)育階段而變化，并將時序數(shù)據(jù)分成較短的時間窗口，然后計算每個時間窗口的基因相關(guān)性，最后運用集成策略將多個相關(guān)矩陣聚合成單一基因鄰接矩陣.這類算法主要包括LEAP（lag-based expression association for pseudo-time-series）[10]、AR1MA1-VBEM（first-order autoregressive moving-average）[11]、SINCERITIES[12]、SCIMITAR（single-cell inference of morphing trajectories and their associated regulation）[13]和SINGE[14].

LEAP 以偽時序數(shù)據(jù)為輸入，首先計算不同時滯的固定時長窗口讀數(shù)（read counts）的皮爾遜相關(guān)性；然后將每對基因間的相關(guān)強度設(shè)置為所有時滯中皮爾遜相關(guān)性絕對值的最大值；最后通過置換檢驗估計錯誤發(fā)現(xiàn)率得到最終的基因相關(guān)矩陣.由于計算的相關(guān)性不是對稱的，因此該算法可以輸出有向網(wǎng)絡(luò).AR1MA1-VBEM 的輸入是單細胞時序數(shù)據(jù)或使用模塊最優(yōu)葉序算法得到的偽時序，然后構(gòu)造一階自回歸移動平均模型來擬合基因表達數(shù)據(jù)，并在變分貝葉斯（variational-Bayesian）框架內(nèi)推斷GRN.SINCERITIES利用時序基因表達分布的時間變化，通過正則化線性回歸（嶺回歸）恢復基因間的定向調(diào)控關(guān)系，同時通過成對基因間的偏相關(guān)分析預(yù)測基因調(diào)控的激活和抑制模式.該算法的計算復雜度較低，能夠進行大規(guī)模基因的預(yù)測.SCIMITAR 首先通過在高維曲線中使用連續(xù)參數(shù)化（時間點作為參數(shù)）的漸變高斯混合模型（morphing Gaussian mixture model）推斷靜態(tài)scRNA-seq數(shù)據(jù)的偽時序軌跡，同時考慮了數(shù)據(jù)中的異方差噪聲，并加入了檢測軌跡進展中相關(guān)基因的統(tǒng)計能力.SCIMITAR 從數(shù)據(jù)中標記出與時序軌跡進展相關(guān)的共表達模式基因，通過計算協(xié)方差矩陣間的距離獲得基因相關(guān)矩陣.針對偽時序的不均勻分布，SINGE 使用基于核的Granger 因果回歸來平滑不規(guī)則的偽時序和補全缺失的表達值來提高GRN 推斷的魯棒性.SINGE輸入偽時序的單細胞基因表達數(shù)據(jù)，輸出預(yù)測的調(diào)控因子和靶基因相互作用的排序列表，利用不同的超參組合的多個廣義Lasso-Granger 檢驗來控制網(wǎng)絡(luò)的稀疏性、核平滑和偽時序分辨率等，最后將所有廣義Lasso-Granger 檢驗得到的部分網(wǎng)絡(luò)聚合成整體GRN.

在這5 種算法中，SCIMITAR 和AR1MA1-VBEM可以從基因表達數(shù)據(jù)中預(yù)測時序，其他3 種算法均需要輸入單細胞時序數(shù)據(jù).這類GRN 推斷算法將不同細胞命運決定的多分支細胞軌跡強制合并為一條線性發(fā)展軌跡，可能會影響網(wǎng)絡(luò)的準確性.此外，這些算法利用皮爾遜相關(guān)性或Granger 因果關(guān)系來推斷基因間的關(guān)系，但由于生物系統(tǒng)的復雜性，很可能產(chǎn)生非線性的基因相互作用，這也會影響網(wǎng)絡(luò)的準確性.在5 種算法中，LEAP 因計算相關(guān)性比其他算法更加簡單而具有明顯的優(yōu)勢.

與基于微分方程的算法一樣，這類算法的性能嚴重依賴于細胞時序的準確性.時序信息通常需要數(shù)據(jù)提供，或者通過偽時序軌跡推斷方法得到.大多數(shù)偽時序軌跡推斷方法需要提供某些先驗信息，如起始細胞、結(jié)束細胞和分支數(shù)，若無法獲得這些信息，則很難建立可靠的時間軌跡，所以這類推斷算法不適于處理多分支細胞軌跡.

1.2.4 基于共表達基因的算法

基于共表達的GRN 推斷算法需要調(diào)用某種相關(guān)性度量來計算兩兩基因間的關(guān)系，其整體工作流程如圖3 所示.首先通過計算每對基因的表達相關(guān)性初始化邊的權(quán)值；然后通過假設(shè)檢驗估計每條邊的顯著性，使用預(yù)定義的閾值去除被認為不顯著的邊；最后輸出最大連通分支.這類算法主要包括NLNET[15]、PIDC（PID and context）[16]、SINCERA（single-cell RNA-seq profiling analysis）[17]和SCENIC（single-cell regulatory network inference and clustering）[18].

圖3 基于基因共表達的GRN 推斷算法的整體工作流程Fig.3 Overall workflow of GRN inference algorithm based on gene co-expression

細胞之間存在著非線性的復雜關(guān)系，因此可以利用非線性依賴關(guān)系的方法推斷網(wǎng)絡(luò).NLNET 是一種基于非線性關(guān)聯(lián)的網(wǎng)絡(luò)重建算法，利用基于條件有序列表的距離[37]（distance based on conditional ordered list，DCOL）計算每對基因間關(guān)系，具有較高的靈敏度和計算效率.PIDC 使用信息論中的部分信息分解（partial information decomposition，PID）來確定基因間的關(guān)系.PIDC 將一對基因X 和Y 提供的關(guān)于目標基因Z 的信息劃分為冗余信息、協(xié)同信息和獨有信息.首先計算提供給Z 的X 獨有信息和X 與Y 的互信息（mutual information，MI）的比率；然后計算為所有其他基因提供的X 和Y 的獨有信息比率的平均值，該平均值被稱為比率獨有貢獻（proportionaluniquecontribution，PUC）.為量化PUC 值的置信度，對每個基因估計其零假設(shè)下PUC 值的分布，然后計算其置信分數(shù).由于X 和Y 的PUC 值是對稱的，因此產(chǎn)生的網(wǎng)絡(luò)是無向的.SINCERA通過給定其他基因的變量，利用一階條件相關(guān)性度量一對基因調(diào)控相關(guān)的重要性，即給定第3 個基因Z，基因X 和Y 的相關(guān)性被定義為X 和Z 的線性回歸與Y 和Z 的線性回歸所產(chǎn)生的殘差之間的相關(guān)性，然后使用最小二乘法計算回歸的目標基因和條件基因的系數(shù)，并使用單樣本t 檢驗檢測基因X 和Y 關(guān)系的顯著性.SINCERA 使用的是線性相關(guān)性方法，因此不能獲得基因間的非線性關(guān)系.SCENIC 首先使用可獲得基因間非線性關(guān)系的GENIE3 或GRNBoost2 算法推斷轉(zhuǎn)錄因子和候選靶基因之間的共表達模塊；然后通過CisTarget（https：//resources.aert-slab.org/cis-target/）識別出調(diào)控因子結(jié)合模體（binding motif）在靶基因中的顯著富集模塊，并構(gòu)建調(diào)控子；最后通過計算細胞中所有基因表達排序的恢復曲線下的面積（area under curve，AUC），對每個細胞中的每個調(diào)控子的活性進行評分，從而生成一個二值化的活性矩陣.SCENIC 最初在R/Bioconductor 中實現(xiàn)，后來被移植到python，命名為pySCENIC[38]，pySCENIC通過軟件容器和并行化大幅縮短了程序運算時間，SCENIC/pySCENIC 的主要局限是CisTarget 步驟無法預(yù)測具有未知模體的轉(zhuǎn)錄因子調(diào)控子.

這類算法需要某種度量方法計算每對基因間的關(guān)系.除SINCERA 外，其他算法均能獲得基因間的非線性關(guān)系.當樣本量不夠大時，NLNET 檢測線性關(guān)系的能力低于基于線性關(guān)系的推斷算法，因此NLNET為了獲得非線性關(guān)系可能會失去一些弱線性關(guān)系，但是可以將其與線性方法相結(jié)合，使結(jié)果更全面.此外，這類算法一般要求數(shù)據(jù)盡量同構(gòu)或盡量少的細胞類型，對于異構(gòu)數(shù)據(jù)（不同發(fā)育階段或不同細胞類型），算法可能忽略基因間的變化關(guān)系.這類算法的一個常見缺點是不能推斷出基因間的具體調(diào)控關(guān)系（如激活或抑制），但是SCENIC 和SINCERA 可以利用現(xiàn)有數(shù)據(jù)庫中的已知調(diào)控網(wǎng)絡(luò)進行富集分析來彌補這一缺點.

1.2.5 基于細胞特異性的算法

傳統(tǒng)構(gòu)建GRN 的算法都是基于所有細胞的相似性而忽略了細胞間的異質(zhì)性.在單細胞的基礎(chǔ)上，細胞特異性網(wǎng)絡(luò)將數(shù)據(jù)從“不穩(wěn)定”的基因表達形式轉(zhuǎn)化為“穩(wěn)定”的基因關(guān)聯(lián)形式.這類算法的輸入是細胞的原始基因表達矩陣，對于每類細胞，輸出一個細胞特異性網(wǎng)絡(luò)，即一對基因可能在某類細胞中有關(guān)聯(lián)，但在其他類細胞中沒有關(guān)聯(lián).這類算法主要包括CSN（cell specific network）[19]、CCSN（conditional cell specific network）[20]和CeSpGRN（cell specific GRN）[21].

CSN 設(shè)計了每個細胞中每對基因獨立性的統(tǒng)計量，以統(tǒng)計獨立性確定該對基因的關(guān)聯(lián)，該統(tǒng)計量能夠很好地區(qū)分獨立關(guān)系和非獨立關(guān)系的基因?qū)?如果一對基因相互獨立，則歸一化后的統(tǒng)計量遵循標準正態(tài)分布.CSN 作為一種無監(jiān)督算法，直接由基因表達矩陣構(gòu)建而成，不需要預(yù)先知道細胞類型或?qū)毎垲?CCSN 是為了解決CSN 對間接效應(yīng)的過高估計而提出的條件細胞特異性網(wǎng)絡(luò)算法.與CSN 不同，CCSN構(gòu)造2 個基因在第3 個基因條件下的條件獨立性統(tǒng)計量.CCSN 對每個條件基因都構(gòu)造該細胞的特異性網(wǎng)絡(luò)，然后將所有基因的特異性網(wǎng)絡(luò)整合為最終的細胞特異性網(wǎng)絡(luò).雖然CCSN 的計算復雜度高于CSN，但是它不僅可以通過消除基因間的間接關(guān)聯(lián)來測量基因間的直接關(guān)聯(lián)，還可以基于單細胞網(wǎng)絡(luò)進行細胞聚類和降維.另外，通過定義每個細胞的基因間關(guān)聯(lián)數(shù)據(jù)的網(wǎng)絡(luò)流熵（network flow entropy，NFE），CCSN可以估計單個細胞的分化潛能，從而掌握細胞分化的譜系動力學性質(zhì).CCSN 的缺點是無法構(gòu)建因果關(guān)系的基因關(guān)聯(lián).CeSpGRN 使用高斯Copula 圖模型對正調(diào)控邊和負調(diào)控邊進行推斷，而CSN 和CCSN 算法無法推斷正負調(diào)控關(guān)系.CeSpGRN 可以在任何軌跡結(jié)構(gòu)或沒有軌跡信息的數(shù)據(jù)集上工作，不需要預(yù)先知道細胞的時間信息.該算法根據(jù)基因表達相似性設(shè)計高維加權(quán)核獲得細胞軌跡平滑變化，通過計算每對細胞在K 最近鄰圖中測地距離的高斯核函數(shù)得到這對細胞的核權(quán)重，最后通過交替方向乘子法優(yōu)化模型得到每個細胞的基因鄰接矩陣.CeSpGRN 的主要缺點是沒有將GRN 建模為有向圖.

基于細胞特異性網(wǎng)絡(luò)推斷算法可以識別基因間的相互作用，并在單細胞數(shù)據(jù)上描繪網(wǎng)絡(luò)的異質(zhì)性.這類模型有助于發(fā)現(xiàn)新的細胞類型，并能夠揭示某些基因在特定細胞類型中的重要作用.

1.2.6 基于深度學習的算法

非深度學習的推斷算法大多是無監(jiān)督的，不能在標記好的數(shù)據(jù)上進行訓練，屬于特定假設(shè)下的特定方法.這些算法雖然在某些情況下表現(xiàn)較好，但是容易導致更多的假陽性或假陰性結(jié)果.而深度學習模型一般不需要對輸入數(shù)據(jù)進行模型或分布假設(shè)，因此能夠模擬細胞異質(zhì)性下的基因相互作用并克服噪聲和減少錯誤.這類算法主要包括CNNC（convolutional neural network for co-expression）[22]、DeepDRIM（deep learning-based direct regulatory interaction model）[23]、Deep-SEM[24]、scGeneRAI（single-cell gene regulatory network prediction by explainable AI）[25]、STGRNS[26]和GENELink（gene regulatory network via link prediction）[27].

CNNC 運用卷積神經(jīng)網(wǎng)絡(luò)（convolutional neural network，CNN）進行有監(jiān)督的基因關(guān)系推斷.該算法首先將所有細胞的每對基因匯總成共現(xiàn)（co-occurrence）直方圖，然后將直方圖轉(zhuǎn)換為類圖像形式（經(jīng)典概率函數(shù)矩陣）以便于進行神經(jīng)網(wǎng)絡(luò)訓練.CNNC 以少量帶標記的基因?qū)献鳛橛柧毤梢詫W習區(qū)分各類基因關(guān)系.CNNC 的輸出層可以根據(jù)應(yīng)用需要輸出多種類型的關(guān)系，如對于因果推斷關(guān)系，CNNC 可以輸出兩基因無交互的概率和一個基因調(diào)控另一個基因的概率等.DeepDRIM 是構(gòu)造細胞類型特異性（cell type specific）GRN 的有監(jiān)督深度神經(jīng)網(wǎng)絡(luò)模型，與CNNC不同，DeepDRIM 主要針對轉(zhuǎn)錄因子基因?qū)Γ═F-gene pair）進行預(yù)測.DeepDRIM 模型的結(jié)構(gòu)如圖4 所示.

圖4 DeepDRIM 模型的結(jié)構(gòu)Fig.4 Structure of DeepDRIM

DeepDRIM 首先將轉(zhuǎn)錄因子基因?qū)Φ穆?lián)合表達二維直方圖轉(zhuǎn)化為主圖像，并將與該基因?qū)哂袕娬齾f(xié)方差的基因?qū)D像作為主圖像的最近鄰，以消除傳遞交互引起的錯誤.然后將2 個堆疊卷積層嵌入結(jié)構(gòu)網(wǎng)絡(luò)A 和B，分別用于處理主圖像和最近鄰圖像.網(wǎng)絡(luò)A 類似于VGGnet[39]，包含堆疊卷積層和最大池化層.網(wǎng)絡(luò)B 的結(jié)構(gòu)與A 相似，是A 的類孿生網(wǎng)絡(luò)，用于處理多個最近鄰圖像，所有子網(wǎng)絡(luò)共享權(quán)值.網(wǎng)絡(luò)B 由已知轉(zhuǎn)錄因子的基因?qū)ο嗷プ饔糜柧殻@些相互作用來自公開可用的細胞類型特異性染色質(zhì)免疫共沉淀測序（ChIP-seq）數(shù)據(jù).最后DeepDRIM 將每個圖像對應(yīng)的向量串聯(lián)成高維向量，經(jīng)過2 個堆疊全連接層壓縮，并使用Sigmoid 函數(shù)進行二分類，得到有或無交互2 種結(jié)果.雖然DeepDRIM 使用的神經(jīng)網(wǎng)絡(luò)模型比CNNC 復雜，但是CNNC 能夠得到一個基因調(diào)控另一基因的概率，而DeepDRIM 只能得到簡單的二分類.

DeepSEM 的算法架構(gòu)是一個beta 變分自動編碼器（beta variational autoencoder，beta-VAE），包括編碼器、GRN 層、用于矩陣逆運算的逆GRN 層和解碼器.DeepSEM 的編碼器每次輸入一個基因的表達，不同基因共享beta-VAE 的權(quán)重.GRN 層和逆GRN 層均為基因交互矩陣，不同基因之間條件依賴關(guān)系的鄰接矩陣使用結(jié)構(gòu)方程模型表示，GRN 層用來擬合基因表達矩陣和鄰接矩陣之間的非線性關(guān)系.由于矩陣逆運算的復雜性，該算法的運行時間會隨基因數(shù)的增加而增加.scGeneRAI 使用分層相關(guān)傳播（layer-wise relevance propagation，LRP）估計每個基因的相關(guān)性，從靜態(tài)scRNA-seq 數(shù)據(jù)中推斷GRN，以LRP 得分衡量目標基因與所有預(yù)測基因之間相互作用的強度.STGRNS的模型結(jié)構(gòu)包括4 個模塊，分別為GEM（gene expression motif）模塊、位置編碼層、Transformer 編碼器和分類層.GEM 模塊將基因?qū)D(zhuǎn)化為可輸入Transformer編碼器的格式；位置編碼層用于捕獲位置或時間信息；Transformer 編碼器用于計算不同子向量（特別是關(guān)鍵子向量）的相關(guān)性；分類層輸出最終的分類結(jié)果.對于靜態(tài)數(shù)據(jù)、偽時序數(shù)據(jù)或時序數(shù)據(jù)，STGRNS 都可以根據(jù)已知的基因間關(guān)系推斷GRN.GENELink 是基于圖注意力網(wǎng)絡(luò)（graph attention network，GAT）的有監(jiān)督模型，該算法首先聚合條件特異性基因表達譜和基于知識的基因相互作用矩陣，用以學習基因的低維向量化表示，然后通過優(yōu)化嵌入空間，并學習特定的基因特征，用于下游相似性測量或基因調(diào)控的因果推理.

在基于深度學習的算法中，CNNC 能夠區(qū)分相互作用、因果對、負對，或任何可以定義的基因關(guān)系；DeepDRIM 只可以識別TF 和基因的相互作用及因果關(guān)系；DeepSEM 不僅可以利用GRN 作為“橋梁”，整合不同的單細胞形態(tài)，構(gòu)建共同的潛在空間，還可以集成其他分子相互作用網(wǎng)絡(luò)，如蛋白和蛋白相互作用網(wǎng)絡(luò)，開放染色質(zhì)數(shù)據(jù)、DNA 的結(jié)合模體和遺傳相互作用網(wǎng)絡(luò)；STGRNS 可以處理靜態(tài)和時序數(shù)據(jù)，根據(jù)已知基因關(guān)系訓練Transformer 進行GRN 推斷，基于Transformer 的方法在適用性和性能方面優(yōu)于其他深度學習模型的算法[26].基于深度學習的算法相較于傳統(tǒng)算法有3 個顯著優(yōu)勢：①scRNA-seq 數(shù)據(jù)隱含著生物實體間和過程間復雜的相互依賴關(guān)系，這些實體和過程通常具有內(nèi)在的噪聲，并且發(fā)生在多個尺度上，深度學習算法能夠檢測生物信息數(shù)據(jù)中的復雜模式；②神經(jīng)網(wǎng)絡(luò)可以自然地擴展到包括表觀遺傳和序列信息在內(nèi)的互補數(shù)據(jù)；③深度學習算法適用于高通量的測序數(shù)據(jù)，隨著測序深度和規(guī)模的提升以及訓練樣本的增加，深度學習算法的性能會得到提升.

2 GRN 推斷算法綜合比較

2.1 各算法的相關(guān)信息和應(yīng)用場景

本文綜述的GRN 推斷算法都有自身的優(yōu)勢和局限性.表1 給出了這些算法的相關(guān)信息，包括算法工具、使用語言和網(wǎng)絡(luò)類型等，表2 給出了這些算法的主要優(yōu)勢和局限性.在實際應(yīng)用中，可根據(jù)不同的應(yīng)用場景和具體研究問題選擇GRN 推斷算法.基于單細胞數(shù)據(jù)研究細胞發(fā)育、細胞分化、疾病發(fā)展或疾病療程，則需要考慮基于時序數(shù)據(jù)或者基于偽時序軌跡推斷的算法，如SCOUP、SCODE、SCENIC、LEAP、PIDC和SINGE 等.對于需要有向網(wǎng)絡(luò)的研究，可以選擇SCOUP、SCODE、LEAP 和SINGE 等算法.非線性關(guān)系的基因網(wǎng)絡(luò)比線性關(guān)系更貼合實際的基因調(diào)控關(guān)系，因此SCOUP 和SINGE 更適合于時序數(shù)據(jù)的非線性有向網(wǎng)絡(luò)的預(yù)測.對于僅需了解與某些基因相關(guān)的關(guān)鍵調(diào)控因子，則應(yīng)選擇能夠嚴格界定基因關(guān)系的二值化網(wǎng)絡(luò)，如SCNS、BTR 和DeepDRIM.對于已得到足夠多的已知基因調(diào)控關(guān)系的數(shù)據(jù)（包括ChIP-seq 或ACTC-seq 等），則可以選擇基于深度學習的算法進行訓練和預(yù)測.

表1 GRN 推斷算法的相關(guān)信息Tab.1 Relevant information of GRN inference algorithms

表2 GRN 推斷算法的應(yīng)用場景、優(yōu)勢和局限Tab.2 Application scenarios，advantages and limitations of GRN inference algorithms

2.2 算法主要評價指標

GRN 推斷算法的主要評價指標包括AUROC（area under the ROC）[41]、AUPRC（area under the precisionrecall curve）[42]、早期精確度（early precision，EP）[31]和早期精確比（early precision ratio，EPR）[31].

ROC（receiver operating characteristic）曲線的橫軸和縱軸分別為假陽性率和真陽性率，AUROC 定義為ROC 曲線下面積.相對于準確率和召回率等指標，ROC 能夠反映真假陽性在不同閾值下的分數(shù)變化，而AUROC 比ROC 能更直觀地反映算法的性能，因此大部分GRN 推斷算法使用AUROC 代替ROC，AUROC的數(shù)值越大，說明算法性能越好.AUPRC 也常用來評價GRN 推斷算法性能，該指標定義為精確率-召回率（precision-recall，PR）曲線下面積，同AUROC 一樣，AUPRC 的數(shù)值越大，說明算法性能越好.EP 為預(yù)測網(wǎng)絡(luò)前k 條邊中真陽性的比例，其中k 為基準真實（ground true）網(wǎng)絡(luò)的邊數(shù).EPR 為模型隨機預(yù)測的早期精度與模型早期精度的比值，其中隨機預(yù)測的精度為基準真實網(wǎng)絡(luò)的邊密度.

2.3 推斷算法比較研究

2018 年，Chen 等[43]使用GeneNetWeaver[44]生成的兩組大腸桿菌單細胞仿真數(shù)據(jù)和2 組分別來自老鼠胚胎干細胞（ESC）和造血干細胞（HSC）樣本的真實單細胞數(shù)據(jù)評估了8 種早期的GRN 推斷算法.該研究分別繪制了算法在仿真數(shù)據(jù)集和真實數(shù)據(jù)集上的PR 曲線和ROC 曲線，發(fā)現(xiàn)8 種算法的性能都不高（SCODE在仿真數(shù)據(jù)集上的AUROC 值最高，但沒超過0.65），分析表明，基于scRNA-seq 數(shù)據(jù)的GRN 推斷算法并不一定比基于批量RNA 測序的推斷算法性能更好.

2020 年，Pratapa 等[31]分析了GeneNetWeaver 的缺點，設(shè)計了一個基于布爾模型的仿真數(shù)據(jù)生成工具BoolODE，使用BoolODE 生成了400 多個包含6 個合成網(wǎng)絡(luò)和4 個精選布爾模型的仿真數(shù)據(jù)集，然后在生成的仿真數(shù)據(jù)和5 個實驗性人類和小鼠scRNA-seq數(shù)據(jù)集上評估了12 種非監(jiān)督、無附加信息和非細胞特異性的GRN 推斷算法，構(gòu)建了基準數(shù)據(jù)集綜合評估框架BEELINE.該研究結(jié)果表明，PIDC、GENIE3 和GRNBoost2 在精確度方面領(lǐng)先其他算法；GENIE3 和PIDC 在靜態(tài)細胞數(shù)據(jù)集上多次運行可表現(xiàn)出更好的穩(wěn)定性，而GRNBoost2 對dropout 的靈敏度較低；在給定較好偽時序數(shù)據(jù)的情況下，SINCERITIES 具有較好的性能和穩(wěn)定性.同年，Dai 等[45]簡單評估了19 種GRN推斷算法，發(fā)現(xiàn)非線性關(guān)系的推斷算法通常比線性關(guān)系的算法能更準確地預(yù)測基因調(diào)控關(guān)系.

2021 年，Nguyen 等[46]使用GeneNetWeaver 工具生成不同基因規(guī)模、不同細胞數(shù)量和不同dropout 率的仿真scRNA-seq 數(shù)據(jù)集，評估了15 種GRN 推斷算法，考察其重構(gòu)參考網(wǎng)絡(luò)的AUROC 值、不同dropout率和稀疏水平的敏感度以及算法的時間復雜性.該研究發(fā)現(xiàn)，大部分算法只能推斷小規(guī)模基因的網(wǎng)絡(luò)，如BTR 要求基因數(shù)小于30，SCNS 要求基因數(shù)小于64，SCINCERA 要求基因數(shù)不超過500，而SCODE、NLNET、SCENIC、LEAP 等算法能夠分析所有基因規(guī)模的數(shù)據(jù)；SCENIC 在大部分模擬數(shù)據(jù)上AUROC 的中位值最高；LEAP 和NLNET 是運行最快的算法，即使是3 000 個基因、1 000 個細胞的數(shù)據(jù)集也可以在數(shù)分鐘內(nèi)完成一次分析；在不同dropout 率和稀疏水平下，SCODE、SCOUP 和SCENIC 的AUROC 值始終高于0.5，其中SCOUP 的表現(xiàn)最穩(wěn)定.

2023 年，Xu 等[26]在48 個涉及不同物種、譜系、網(wǎng)絡(luò)類型和網(wǎng)絡(luò)規(guī)模的基準數(shù)據(jù)集上進行實驗，發(fā)現(xiàn)STGRNS 優(yōu)于有監(jiān)督算法CNNC 以及4 種無監(jiān)督算法PIDC、LEAP、SINGE 和DeepSEM，并且STGRNS 在“TF 基因?qū)Α鳖A(yù)測上具有一定的可移植性，且有更好的超參數(shù)魯棒性.

本文在人類胚胎干細胞（hESC）數(shù)據(jù)集上簡單評估以上26 種算法，并計算其AUROC 值，結(jié)果如圖5所示，可以發(fā)現(xiàn)DeepDRIM、GENELink 和STGRNS 的表現(xiàn)優(yōu)于其他算法.

圖5 GRN 推斷算法的性能評估Fig.5 Performance evaluation of GRN inference algorithms

3 GRN 推斷算法面臨的機遇與挑戰(zhàn)

由于scRNA-seq 數(shù)據(jù)可能無法為可靠的網(wǎng)絡(luò)推斷提供足夠的分辨率，因此，當前的推斷算法難以預(yù)測真實的GRN 或推測性能較低.一般來說，GRN 推斷算法包括3 個關(guān)鍵組成部分：數(shù)據(jù)預(yù)處理、特征提取和潛在基因調(diào)控模式的推斷.數(shù)據(jù)預(yù)處理的目的是轉(zhuǎn)換數(shù)據(jù)，使其更易于分析，如對細胞進行偽時序排序和降維等.即使單細胞數(shù)據(jù)不構(gòu)成實際的時間序列，也可以根據(jù)偽時序索引對其進行排序.特征提取旨在選取特定基因和提取基因調(diào)控信息，如基因表達均值和基因互信息等，但是在很多情況下，統(tǒng)計量沒有完全刻畫基因間的關(guān)系，表達模式間的統(tǒng)計關(guān)系對應(yīng)于調(diào)控相互作用的假設(shè)也可能存在固有缺陷.對潛在基因調(diào)控模式的推斷是在所有可能的GRN 空間中找到與數(shù)據(jù)統(tǒng)計最匹配的GRN，但是基因調(diào)控的隨機性引起的數(shù)據(jù)內(nèi)在噪聲會降低預(yù)測的準確性.

單細胞數(shù)據(jù)本身的生物異質(zhì)性和技術(shù)異質(zhì)性也為GRN 推斷帶來挑戰(zhàn).生物異質(zhì)性主要由基因隨機突變或內(nèi)部細胞相互作用和環(huán)境變化引起的表型變化導致.即使是批次相同的細胞類型，由于噪聲和化學批次不同等因素，單個細胞類型的細胞數(shù)量和轉(zhuǎn)錄狀也可能會有顯著差異.GRN 主要代表了所有細胞群的平均活動，這些細胞群可以由最豐富的細胞類型所控制.因此，組織網(wǎng)絡(luò)無法捕捉單個細胞群的獨特行為，以及細胞如何相互作用來執(zhí)行更高層次的組織功能.盡管一些GRN 推斷算法聲稱能夠模擬基因表達的隨機部分，但是區(qū)分不同異質(zhì)性仍是一個重大挑戰(zhàn).此外，GRN 推斷算法的驗證和性能評估是最具挑戰(zhàn)性的問題，目前仍缺乏高標準的指標評估這些數(shù)學模型的性能.針對不同的應(yīng)用場景，整合其他類型的數(shù)據(jù)用于訓練和驗證非常重要，如運用已知的轉(zhuǎn)錄因子結(jié)合位點信息的ChIP-seq 數(shù)據(jù)，或者依靠染色質(zhì)可及性和基因甲基化等多組學測序數(shù)據(jù).

隨著單細胞測序技術(shù)的發(fā)展，單細胞數(shù)據(jù)GRN推斷算法的進步將推動細胞功能的研究，在細胞水平上發(fā)現(xiàn)缺失的GRN 成分，可應(yīng)用于識別疾病的生物標志物和信號通路，從而為疾病精準醫(yī)療和設(shè)計藥物提供支持.