基于單細胞RNA 測序數據的基因調控網絡推斷算法綜述

張少強，潘鏡伊

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

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

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

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

1 GRN 推斷算法

1.1 早期基于批量RNA 測序數據的算法

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

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

1.2 基于scRNA-seq 數據的算法

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

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

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

1.2.1 基于布爾網絡的算法

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

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

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

1.2.2 基于微分方程的算法

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

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

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

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

由于微分方程可以模擬RNA、蛋白質等在GRN中與時間有關的各種相互作用的參數變化，并可以推斷基因間的某些因果關系，因此比其他方法更接近實際網絡和細節.然而，由于微分方程的推斷算法需要輸入大量數據來估計參數且計算量較大，因此只能應用于有限數量的基因調控關系.

1.2.3 基于偽時序基因相關性集成策略的算法

與基于微分方程的算法一樣，基于偽時序基因相關性集成策略的推斷算法也需要推斷細胞的偽時序或者需要用戶提供細胞時序，通過時序基因表達數據推斷出動態GRN[36].這類算法假設基因間的關系隨著細胞的不同發育階段而變化，并將時序數據分成較短的時間窗口，然后計算每個時間窗口的基因相關性，最后運用集成策略將多個相關矩陣聚合成單一基因鄰接矩陣.這類算法主要包括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 以偽時序數據為輸入，首先計算不同時滯的固定時長窗口讀數（read counts）的皮爾遜相關性；然后將每對基因間的相關強度設置為所有時滯中皮爾遜相關性絕對值的最大值；最后通過置換檢驗估計錯誤發現率得到最終的基因相關矩陣.由于計算的相關性不是對稱的，因此該算法可以輸出有向網絡.AR1MA1-VBEM 的輸入是單細胞時序數據或使用模塊最優葉序算法得到的偽時序，然后構造一階自回歸移動平均模型來擬合基因表達數據，并在變分貝葉斯（variational-Bayesian）框架內推斷GRN.SINCERITIES利用時序基因表達分布的時間變化，通過正則化線性回歸（嶺回歸）恢復基因間的定向調控關系，同時通過成對基因間的偏相關分析預測基因調控的激活和抑制模式.該算法的計算復雜度較低，能夠進行大規?；虻念A測.SCIMITAR 首先通過在高維曲線中使用連續參數化（時間點作為參數）的漸變高斯混合模型（morphing Gaussian mixture model）推斷靜態scRNA-seq數據的偽時序軌跡，同時考慮了數據中的異方差噪聲，并加入了檢測軌跡進展中相關基因的統計能力.SCIMITAR 從數據中標記出與時序軌跡進展相關的共表達模式基因，通過計算協方差矩陣間的距離獲得基因相關矩陣.針對偽時序的不均勻分布，SINGE 使用基于核的Granger 因果回歸來平滑不規則的偽時序和補全缺失的表達值來提高GRN 推斷的魯棒性.SINGE輸入偽時序的單細胞基因表達數據，輸出預測的調控因子和靶基因相互作用的排序列表，利用不同的超參組合的多個廣義Lasso-Granger 檢驗來控制網絡的稀疏性、核平滑和偽時序分辨率等，最后將所有廣義Lasso-Granger 檢驗得到的部分網絡聚合成整體GRN.

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

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

1.2.4 基于共表達基因的算法

基于共表達的GRN 推斷算法需要調用某種相關性度量來計算兩兩基因間的關系，其整體工作流程如圖3 所示.首先通過計算每對基因的表達相關性初始化邊的權值；然后通過假設檢驗估計每條邊的顯著性，使用預定義的閾值去除被認為不顯著的邊；最后輸出最大連通分支.這類算法主要包括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

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

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

1.2.5 基于細胞特異性的算法

傳統構建GRN 的算法都是基于所有細胞的相似性而忽略了細胞間的異質性.在單細胞的基礎上，細胞特異性網絡將數據從“不穩定”的基因表達形式轉化為“穩定”的基因關聯形式.這類算法的輸入是細胞的原始基因表達矩陣，對于每類細胞，輸出一個細胞特異性網絡，即一對基因可能在某類細胞中有關聯，但在其他類細胞中沒有關聯.這類算法主要包括CSN（cell specific network）[19]、CCSN（conditional cell specific network）[20]和CeSpGRN（cell specific GRN）[21].

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

基于細胞特異性網絡推斷算法可以識別基因間的相互作用，并在單細胞數據上描繪網絡的異質性.這類模型有助于發現新的細胞類型，并能夠揭示某些基因在特定細胞類型中的重要作用.

1.2.6 基于深度學習的算法

非深度學習的推斷算法大多是無監督的，不能在標記好的數據上進行訓練，屬于特定假設下的特定方法.這些算法雖然在某些情況下表現較好，但是容易導致更多的假陽性或假陰性結果.而深度學習模型一般不需要對輸入數據進行模型或分布假設，因此能夠模擬細胞異質性下的基因相互作用并克服噪聲和減少錯誤.這類算法主要包括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 運用卷積神經網絡（convolutional neural network，CNN）進行有監督的基因關系推斷.該算法首先將所有細胞的每對基因匯總成共現（co-occurrence）直方圖，然后將直方圖轉換為類圖像形式（經典概率函數矩陣）以便于進行神經網絡訓練.CNNC 以少量帶標記的基因對集合作為訓練集，可以學習區分各類基因關系.CNNC 的輸出層可以根據應用需要輸出多種類型的關系，如對于因果推斷關系，CNNC 可以輸出兩基因無交互的概率和一個基因調控另一個基因的概率等.DeepDRIM 是構造細胞類型特異性（cell type specific）GRN 的有監督深度神經網絡模型，與CNNC不同，DeepDRIM 主要針對轉錄因子基因對（TF-gene pair）進行預測.DeepDRIM 模型的結構如圖4 所示.

圖4 DeepDRIM 模型的結構Fig.4 Structure of DeepDRIM

DeepDRIM 首先將轉錄因子基因對的聯合表達二維直方圖轉化為主圖像，并將與該基因對具有強正協方差的基因對圖像作為主圖像的最近鄰，以消除傳遞交互引起的錯誤.然后將2 個堆疊卷積層嵌入結構網絡A 和B，分別用于處理主圖像和最近鄰圖像.網絡A 類似于VGGnet[39]，包含堆疊卷積層和最大池化層.網絡B 的結構與A 相似，是A 的類孿生網絡，用于處理多個最近鄰圖像，所有子網絡共享權值.網絡B 由已知轉錄因子的基因對相互作用訓練，這些相互作用來自公開可用的細胞類型特異性染色質免疫共沉淀測序（ChIP-seq）數據.最后DeepDRIM 將每個圖像對應的向量串聯成高維向量，經過2 個堆疊全連接層壓縮，并使用Sigmoid 函數進行二分類，得到有或無交互2 種結果.雖然DeepDRIM 使用的神經網絡模型比CNNC 復雜，但是CNNC 能夠得到一個基因調控另一基因的概率，而DeepDRIM 只能得到簡單的二分類.

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

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

2 GRN 推斷算法綜合比較

2.1 各算法的相關信息和應用場景

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

表1 GRN 推斷算法的相關信息Tab.1 Relevant information of GRN inference algorithms

表2 GRN 推斷算法的應用場景、優勢和局限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 能夠反映真假陽性在不同閾值下的分數變化，而AUROC 比ROC 能更直觀地反映算法的性能，因此大部分GRN 推斷算法使用AUROC 代替ROC，AUROC的數值越大，說明算法性能越好.AUPRC 也常用來評價GRN 推斷算法性能，該指標定義為精確率-召回率（precision-recall，PR）曲線下面積，同AUROC 一樣，AUPRC 的數值越大，說明算法性能越好.EP 為預測網絡前k 條邊中真陽性的比例，其中k 為基準真實（ground true）網絡的邊數.EPR 為模型隨機預測的早期精度與模型早期精度的比值，其中隨機預測的精度為基準真實網絡的邊密度.

2.3 推斷算法比較研究

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

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

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

2023 年，Xu 等[26]在48 個涉及不同物種、譜系、網絡類型和網絡規模的基準數據集上進行實驗，發現STGRNS 優于有監督算法CNNC 以及4 種無監督算法PIDC、LEAP、SINGE 和DeepSEM，并且STGRNS 在“TF 基因對”預測上具有一定的可移植性，且有更好的超參數魯棒性.

本文在人類胚胎干細胞（hESC）數據集上簡單評估以上26 種算法，并計算其AUROC 值，結果如圖5所示，可以發現DeepDRIM、GENELink 和STGRNS 的表現優于其他算法.

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

3 GRN 推斷算法面臨的機遇與挑戰

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

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

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