图1
随着生物数据的数量和类型持续以指数速度增长，生物网络的数量和类型也在增长，包括蛋白质-蛋白质相互作用网络、代谢网络、遗传相互作用网络、基因/转录调节网络(GRNs)、细胞信号网络。虽然不同的网络模型根据其基本假设具有固有的优势和局限性，但它们都具有图形化模型的共同特征，即描述生物系统中的信息流，以帮助理解和解释基本的生物过程。
在过去的几十年里，网络建模被广泛应用于帮助理解关键的生物过程和健康和疾病的调节。特别是，人类生理和病理生理学的巨大复杂性要求在系统水平上理解生物分子如何在单个细胞和组织内相互作用，细胞和组织之间如何相互作用以维持体内平衡，以及这些相互作用的干扰如何导致疾病。omnigenic disease模型指出，网络中所有基因的相互作用都可能导致复杂疾病，该模型正日益被人们所认可和接受。这些概念框架完全符合网络生物学，因此，在生物学的所有领域中，网络建模方法的使用越来越多也就不足为奇了。
例如，许多遗传变异可以影响疾病，每一种变异都通过很小的影响使生物学解释变得困难。这些复杂的遗传效应可以通过它们在转录、信号网络和生物途径中的关系来更好地理解。我们的团队和其他人已经利用网络模型来解释复杂疾病的遗传原因。类似地，网络可以用来了解与各种环境引起的疾病有关的分子级联。例如，Chella Krishnan等人通过整合组织特异性GRNs的遗传关联，发现与非酒精性脂肪肝相关的大量遗传变异影响多种生物途径，包括脂质代谢、免疫系统、细胞周期、转录调节、胰岛素信号、Notch信号和氧化磷酸化，这些途径在肝脏和脂肪组织的GRNs中相互作用。
基于网络拓扑结构，他们确定了疾病通路和亚网络中心参与线粒体功能的关键调节因子。在另一项研究中，利用组织特异性GRNs对心血管疾病和2型糖尿病的遗传风险进行网络建模，揭示了共享的和疾病特异性的网络和调节因子。格林等建造144组织/特异性网络和使用这些网络来预测和理解lineage-specific IL1B刺激的反应。
虽然基于网络的方法促进了我们对复杂疾病的理解，但需要注意的是，大多数网络方法和应用主要依赖于从bulk组织中获得的组学数据。在组织水平上，已经开发了许多用于网络建模的方法和算法，主要关注于预测组织内和组织间的GRNs，并具有合理的准确性。然而，肝脏的非实质细胞等组织是由包括库普弗细胞、窦状内皮细胞和肝细胞卫星细胞在内的异质细胞群组成的，它们都具有与独特的基因调控谱相关的独特功能。考虑到组织的异质性，组织网络主要代表了所有细胞群的平均活动，这些细胞群可以由最丰富的细胞类型所控制。因此，组织网络无法捕捉单个细胞群的独特行为，以及细胞如何相互作用来执行更高层次的组织功能。
最近的高通量单细胞组学技术的爆炸带来了令人兴奋的可能性，包括但不限于模型的动态，内部和细胞间的基因网络，以阐明背后的过程，细胞发育，功能状态和细胞-细胞通讯，而这些不可能在bulk水平获得。这些单细胞组学技术给了我们前所未有的能力以检查转录，蛋白质和表观基因组的概况在单细胞解决，在调节和功能关系的生物分子在单个细胞或细胞类型以及细胞群之间。在理论上，类似的框架和方法已经被用于组织网络建模，可以扩展到单细胞数据，以揭示调控细胞内部和细胞之间的功能的调节机制。然而，正如Chen和Mar在他们最近的研究中所指出的，bulk组织模型可能不太适合克服单细胞数据带来的独特挑战。
在此，我们将讨论为bulk组织组学数据开发的现有网络建模方法，使用单细胞组学数据进行网络建模所面临的独特挑战，利用单细胞数据建立网络模型的方法的最新发展及其关键的底层算法优缺点。最后，我们讨论了有待克服的问题以及我们认为该领域将在哪些方面实现基于单细胞组学数据的更高效、更准确的基因调控网络建模。
bulk组织数据GRN建模方法
针对bulk组织数据开发和优化的常见GRN方法通常基于相关、回归、常微分方程(ode)、互信息、高斯图形模型和贝叶斯方法。例如：
基于相关关系的加权基因共表达网络分析(WGCNA)是最常用的方法。WGCNA用于发现高度相关的基因簇(或模块)，这些基因通常代表参与类似生物途径或功能的受到严格调控的基因。虽然基于共表达的方法计算效率高，且较少依赖于假设，但这些方法主要对功能相似或调控相似的基因进行分组，但不能推断出方向性或直接的调控关系，需要整合其他信息以促进可解释性。
基于回归的方法，如GENIE3，通过基于回归模型确定每个网络基因的最具预测性的基因子集来解析网络。这些方法适用于线性级联，但不适用于前馈回路。
对于基于相互信息的方法，如ARACNE和CLR，网络结构是由基因对之间的依赖程度决定的。这些相互信息的网络方法可以推断方向性和潜在的因果关系，可以更准确地预测前馈回路，但线性级联的性能有限。
贝叶斯网络(BN)建模方法提供了灵活的框架来合并和整合多组数据作为先验信息，以推断因果性和方向性基因-基因相互作用。BN编码基因之间的条件依赖性，其中每个基因由其父节点的值决定。为了提高精确度，BNs通过可能图的多元空间进行搜索，这样做的代价是较高的计算成本，并且不能保证可以检测到最优拓扑。
常用的GRN推理算法各有优缺点，多方法的集成可以弥补每种方法固有的缺点，更好地解释数据。需要注意的是，这些方法是针对bulk组织级数据进行优化的，这些数据通常符合标准数据分布，并且几乎没有缺失值。
单细胞技术和数据结构
利用最近开发的单细胞技术，我们现在能够检测转录(DropSeq ， inDrop， 10X，SmartSeq v4，Marsseq ， Seq-Well ， SPLiT-seq ， sci-RNA-seq)，蛋白质(CITE-seq)，以及表观基因组如开放染色质(scacc -seq)和甲基化景观。这些单细胞技术带来了令人兴奋的可能性，以前所未有的分辨率和规模探索生物学。单细胞RNA测序(scRNAseq)是目前最流行和广泛应用的高通量检测单细胞的技术。通常，这些高通量单细胞转录组技术是基于从3 '端计数转录片段，然后与参考基因组对齐。由此产生的数据结构集合了每个单细胞的基因计数，称为数字基因表达谱(digital gene expression ，DGE)。
对于其他数据类型，相似的细胞标记(如蛋白质，染色质位置，和甲基化位点)矩阵形成主要的数据结构。虽然已经将单细胞表观基因组投射到单细胞转录组上，但就我们所知，整合多组数据用于GRN建模还没有尝试过，这是方法学未来发展的方向。多组数据可以以多种方式合并，包括构建一个跨组学层外推的具有边缘置信度的单一网络，以及从单个组学层构建多个网络，通过关联关系或已知功能相关性在各层之间进行交互。例如，位于特定基因启动子或增强子区域的开放染色质将允许在scac -seq和scRNAseq层之间绘制有向边;蛋白质组数据可以帮助推断蛋白质之间的相互作用，并提供有关调控蛋白质的信息，如转录因子(TFs)和调控转录组和表观基因组的表观基因组调控。在这篇综述中，我们将重点关注scRNAseq数据，因为它们是GRN建模中研究的最丰富的单细胞数据类型。
现有GRN方法在单细胞网络建模中的性能
最近，Chen和Mar在利用经验和模拟单细胞数据进行网络构建时，评估了五种常用的用于bulk组织数据的广义网络重建方法的能力。在他们的分析中使用的方法包括部分相关、BN、GENIE3、ARACNE和CLR。利用精确召回和接受者工作特性曲线（ precision-recall and receiver operating characteristic curves ）来评估每种方法是否能准确再现参考网络，发现在模拟和实验单细胞数据集中，所有方法都没有明显优于随机生成方法。此外，在网络预测中，不同方法之间也只存在有限的重叠。这表明，现有的基于单细胞数据的网络建设方法缺乏通用性和适用性。然而，在解释这种比较结果时需要谨慎，因为使用的金标准参考网络的有效性和质量评价指标会显著影响比较结果。
scRNAseq数据网络建模中特有挑战与机遇
现有方法的潜在性能不足可能是由于与数据稀疏性、分布以及数据维数和容量增加相关的独特挑战造成的。
首先，对于使用最近的高通量平台的scRNAseq来说，由于单细胞中存在的mRNA数量非常少，以及目前技术限制导致DGE矩阵中的大多数条目都是零，这导致矩阵非常稀疏，使得为bulk组织数据设计的方法的直接扩展非常困难。重要的是，尽管这些零可能是个体细胞中随机基因表达的结果(生物学上的零)，但它们并不一定意味着mRNA分子的缺失，而是对中至低表达基因的低技术敏感性的结果，称为缺失。【著名的双零问题】值得注意的是，基于读取计数的scRNAseq是零膨胀的，而包含独特分子标识符(UMI)计数的scRNAseq被发现具有“非零膨胀”特征，导致与基于读取计数的技术相比，其分布不同。reads计数和基于uml的scRNAseq中底层数据分布的差异要求在未来实现新方法时考虑这些不同技术的数据特征。
在试图为缺失值赋值时，许多单细胞归算方法，例如MAGIC、scImpute、DrImpute、SAVER、BISCUIT、ScUnif、PBLR、deepImpute等得到了开发和应用。但是，这些方法的性能差异很大。在基准测试中，scImpute和DrImpute在模拟数据上成功，但在面对非共线经验数据时失败，而SAVER和BISCUIT只能持续地将dropouts归为接近零的值。此外，用于测量性能的主要指标(例如rand指数或相互信息)对这些方法定义细胞集群的能力进行基准测试；目前还不清楚这些估算值如何影响网络结构。由于对这些数值插入方法的结果没有一致的意见和实验验证，在使用数值插入数据进行网络构建时需要谨慎。简单和直观的方法由汉等使用相同的细胞亚群的细胞类型，平均每个基因在细胞的非零值从每个子集获得一个超细胞（supercell）的基因表达矩阵，与零值和更少的膨胀可能更多的生物有关。值得注意的是，这种做法会减少细胞数量并牺牲统计能力。
第二个挑战是与细胞数据中的dropout 问题相关的非标准数据分布模式。大量的缺失值显著地使数据分布从单峰分布(如高斯分布)向多模态分布倾斜，这违反了大多数经典GRN建模方法的统计假设。对于单细胞网络的构建，需要仔细评估数据分布模式和适当的统计方法。有几种统计方法，如零膨胀因子分析(ZIFA)和ZINB-WaVE(基于零膨胀负二项式的期望变异提取)已经被开发来专门建模零膨胀的单细胞数据分布。ZIFA是一种降维方法，其假设是低表达的基因比高表达的基因更有可能导致缺失。ZIFA扩展了因子分析，在非零均值表达式的基础上，将dropout率模型作为指数衰减。然而，ZIFA也有限制，因为它严格地对零测量进行建模，不能解释接近零的值。此外，子发还有一个底层的线性转换框架;然而，非线性降维技术，如t-SNE和UMAP已经被证明在解释单细胞数据中是有用的，所以零膨胀模型的扩展到这些非线性方法可能是有用的。ZINB-WaVE是另一种降维技术，它独特地模拟了scRNAseq数据的计数性质，并使用样本级截获（sample-level intercept）和灵活的基因级和样本级协变合并标准化，以解决批次效应和序列组成效应(如基因长度或GC内容)。为了解决数据的零膨胀和过度分散问题，ZINB-WaVE修改了一个不适合数据的标准负二项分布，用一个术语给出观察到0的概率，而不是实际的计数。虽然ZINB-WaVE主要被证明是单细胞数据的降维技术，但作者认为低维表示可以用于下游分析，如聚类或伪时间。
最近，Townes等发现多项式方法在特征选择和降维方面优于现有的其他方法。考虑这些替代的统计方法在GRN推断可能被证明是有用的。应该指出的是，这些统计方法是为读计数数据开发的，可能不适合基于UMI的单细胞数据集，因为它们有不同的底层数据分布，而这些数据分布不是零膨胀的。
第三，该领域必须掌握校正混杂因素的能力，并将从多个实验获得的数据推断为一个共同的图谱。挑战出现在各批次数据由变量构成和研究，甚至批次还包含相同的细胞类型，个别细胞类型的细胞数量和转录状态由于程序上可以有很大的不同噪声(组织分离、排序和试剂批次)，scRNAseq平台(例如10×与Dropseq)，和试剂版本(10×的2和3版本)。就像在bulk组织设置中使用批次校正来调整混杂因素一样，不同实验甚至实验室产生的数据集的集成是重要的，因为它增强了统计强度和重现性。最初用于bulk组织校正的方法，如limma和ComBat，已应用于单细胞数据的批量校正；而，已有研究表明，将这些为bulk数据开发的方法同时应用于模拟和真实单细胞数据存在局限性。最近，该领域取得了重大进展，产生了专门用于单细胞批次校正的方法，如典型相关分析(CCA)和mnnCorrect，以及基于带标记的参考数据集的细胞类型识别方法，如scmap和singleR。然而，在对单细胞数据应用批次修正方法后，谨慎地进行下游分析(如GRN构建)是很重要的，而且有必要了解底层算法和假设。
像CCA和mnnCorrect这样的方法只利用跨数据集共享的高度可变的基因进行集成，并返回一个校正后的基因表达矩阵，其中只包含用于集成的可变基因。这些基因主要定义细胞类型特异性标记，而CCA的过程固有地引入了基因之间的依赖关系，并违反了用于下游分析(如差异表达)的统计测试的假设，因此CCA的作者警告说，不要使用CCA进行跨数据集的保守细胞类型鉴定。一般来说，为批次数据开发的批次校正方法在批次校正中执行得更差，而为单细胞数据开发的方法在从不同批次聚集细胞类型方面更准确，但可能不能扩展到下游分析。因此，有必要开发能两者兼顾的方法。
最后，与通常由试验组id、样本id和特征测量组成的bulk组织数据相比，单细胞数据通过从每个样本中添加数十种细胞类型和数千个细胞，也呈现出维数和数据量的增加。这种维数和数据量的增加，不仅使网络建模更加复杂，计算成本更高，而且从生物学角度带来了现有方法无法承受的新可能性。除了基因在网络中是如何组织和相互作用这一典型问题之外，人们还可以提出许多新的令人激动的问题。例如：
什么定义了细胞类型?
基因在每种细胞类型中是如何组织的?
细胞类型之间的网络架构有何不同?
细胞之间的关系是什么?
它们来自相同的还是不同的血统?
这些血统是如何进化的?
同一细胞类型是否有不同的状态?
什么基因调控通路决定细胞状态?
细胞如何从一种状态过渡到另一种状态?
哪些细胞相互沟通以确定更高层次的功能，以及它们通过哪些基因和途径进行沟通?
许多这些新问题在bulk时代并没有被考虑或容易解决。除了提供回答这些重要问题的机会外，在每个样本中测量的 细胞中细胞间的差异性或异质性也提供了足够的信息来构建样本内或特定于轮廓的网络。这样的网络描述单个生物样本的GRN，这在bulk时代是不可能的。换句话说，利用大细胞数维的能力允许为每个样本构建基于其组成细胞剖面的grn，这可以用于样本间的共识网络，以提高准确性。
最近的scRNAseq GRN建模方法
认识到对单细胞数据的新的GRN建模方法的需要，最近开发了许多主要基于scRNAseq数据的方法。我们根据基本的生物学问题对这些方法进行分类(动态建模、细胞内网络和细胞-细胞相互作用网络。然后是具体的生物学假设(例如TF目标相互作用，配体-受体相互作用)和算法(例如共表达，回归，ode，贝叶斯和布尔型)，如表1所示：
| Category | Example methods | Underlying biological assumption | Algorithmic basis | Advantages | Limitations | |
|---|---|---|---|---|---|---|
| Dynamic network **(extensively reviewed in refs | SCNS | Single-gene changes between cell transition states can inform on gene regulatory relations | Boolean | Does not rely on prior knowledge. Has a web UI. Resulting models are executable and can be used to make predictions | Need data discretization; limit to small numbers of genes; regulatory relations need to follow Boolean rules | |
| Dynamic network (extensively reviewed in refs) | SCODE [82] | TF expression dynamics (pseudo-time) and TF regulatory relations (GENEI3) | ODE; Bayesian model selection | Estimate relational expression efficiently using linear regression; reduction of time complexity; fast algorithm | Need dimension reduction first for computing speed and memory feasibility; assumes that all cells are on the same trajectory; optimization is computationally intractable | |
| Dynamic network (extensively reviewed in refs) | GRISLI [83] | Variability in scRNAseq data caused by cell cycle, states, etc. allows the inference of pseudo-time associated with each individual cell | ODE | Makes no restrictive assumption on the gene network structure; can consider multiple trajectories; fast algorithm | Has to estimate the velocity of each individual cell using information from neighbors | |
| Dynamic network (extensively reviewed in refs) | SINCERITIES [84] | Changes in the expression of a TF will alter the expression of target genes | Ridge regression and partial correlation analysis | Low computational complexity and able to handle large-scale data | Requires scRNAseq data at multiple time points. Restricted to TFs and their targets to infer edges | |
| Dynamic network (extensively reviewed in refs) | Scribe [85] | Cell ordering can be improved with time-series or cell velocity estimations | RDI | Outperforms other pseudo-time methods given time-series data. Can be applied to any data type if the data structure is appropriate | Requires time-ordered gene expression profiles or velocity estimation from introns and exons | |
| Dynamic network (extensively reviewed in refs) | AR1MA1-VBEM [40] | The cell differentiation process or response to external stimulus reveals the hierarchical structure of the transcriptome | First-order autoregressive moving-average and variational Bayesian expectation-maximization | Weighted interactions between genes along psuedotime. Model used accounts for noisy data | Data are expressed as fold changes between timepoints/conditions or scaled by housekeeping genes | |
| Dynamic network (extensively reviewed in refs) | SCINGE [86] | Learned target regulator genes can be used to assign each cell to their progress along a trajectory | Granger causality | Smooths irregular pseudo-times and missing expression values | Near random performance for predicting targets of individual regulators | |
| Dynamic network (extensively reviewed in refs) | SoptSC [87] | Similarities between whole transcriptomes of single cells can be used to order them | Cells ordered by minimum paths on weighted cluster-to-cluster graph derived from cell similarity matrix | Includes comprehensive single-cell workflow; leverages information from other parts of the workflow to improve performance | Cannot be run with other tools, have run the full workflow to get pseudo-time inference | |
| Within-cell or cell population network | SCENIC [88] | TF target-based regulation | Combining TF regulatory relations (GENIE3) with TF-binding motif analysis | Robust against dropouts, get a TF score for individual cells (no averaging of cells). | Limited to TF-based relations | |
| Within-cell or cell population network | Pina et al. [89] | TFs drive lineage commitment | Odds ratio for on/off gene associations and spearmen correlation for expression levels associations | Robust to dropouts | Based on single-cell multiplex qRT-PCR, may be difficult to extend the method to sparse single-cell data (selected 44 genes to test) | |
| Within-cell or cell population network | Iacono et al. [90] | Coexpression is regulated by TFs, cofactors, and signaling molecules which can be captured with gene–gene correlations | Pearson correlation using z-score-transformed counts | Can compute correlations at the single-cell level and it is robust to dropouts and noise inherent to single-cell data | Networks are very dense (some have millions of significant edges) | |
| Within-cell or cell population network | PIDC [39,91] | Gene regulatory information reflected in dependencies in the expression patterns of genes | Partial information decomposition using gene trios | Compared with correlation, captures more complicated gene dependencies | Networks are influenced by data discretization, choice of mutual information estimator, method developed for sc-qPCR data, may not be extendable to higher throughput and sparser scRNAseq data | |
| Within-cell or cell population network | Jackson et al. [92] | Deletion of TFs combined with experimental conditions allows for the inference of gene relationships | MTL to leverage cross-dataset commonalities and incorporate prior knowledge | Does not require sophisticated normalization of single-cell data or imputation. Able to combine multiple conditions/datasets for more accurate inference. TF deletions give strong causal link to affected genes | Requires single-cell data with TF deletions and/or environmental perturbations | |
| Within-cell or cell population network | Wang et al. [93] | Gene perturbations allow for inference of causal relationships | Scoring of conditional independence test to identify optimal DAG | Gives causal relationships between genes | Requires interventional data. No loops allowed in DAG | |
| Within-cell or cell population network | ACTION [94] | Functional identity of cells is determined by a weak, but specifically expressed set of genes which are mediated by TFs | Kernel-based cell similarity and geometric approach to identify primary functions | Robust to dropout and does not require averaging. Identifies functions unique to cell types | Requires TFs and their targets. Only provides TF-driven networks | |
| Within-cell or cell population network | SINCERA [95] | TF target-based regulation | First-order conditional dependence on gene expression to construct a DAG | Key TFs identified using multiple importance metrics | Only considers TFs and their targets. Requires genes/TFs to be DEGs or expressed in >80% of cells | |
| Cell–cell communication network | iTALK [96] | Ligand–receptor interactions | Threshold ranked list of genes from two cell types for ligand–receptor pairs | Allows for the inference of directionality of interaction | Requires curation of ligand–receptor interactions (not all interactions are known). Average expression at the cell-type level (no longer single cell). Cannot reveal novel interactions beyond known ligand–receptor knowledge | |
| Cell–cell communication network | Zhou et al. [97] | Ligand–receptor interactions | Expression of ligand and corresponding receptor more than three standard deviations greater than the mean | Allows for the inference of directionality of interaction | Requires curation of ligand–receptor interactions (not all interactions are known). Average expression at the cell-type level (no longer single cell) | |
| Cell–cell communication network | Kumar et al. [98] | Ligand–receptor interactions | Product of the average expression of ligand and corresponding receptor | Allows for the inference of directionality of interaction. Interaction score gives the strength of interaction (rather than just significance) | Requires curation of ligand–receptor interactions (not all interactions are known). Average expression at the cell-type level (no longer single cell) | |
| Cell–cell communication network | Arneson et al. [99] | Ligand to downstream signaling | Coexpression of ligand genes in source cells with other genes in target cells | Use secreted ligands as a guidance for directional inference between cell populations | Gene expression is summarized to the cell population level and coexpression is at the sample level, requiring large sample sizes | |
| Cell–cell communication network | SoptSC [87] | Ligand–receptor interactions | Likelihood estimate of the interaction | Cell–cell communication networkbetween two cells based on expression of the ligand, receptor, and downstream pathway target genes (including expression direction). Consensus signaling network derived from all cells in each cluster | Incorporates target genes of pathways and their directionality. Computes interaction likelihood at the single-cell level and summarizes across all cells in the cluster for higher confidence | Requires curation of ligand–receptor interactions and their downstream pathways |
| Cell–cell communication network | scTensor [100] | Ligand–receptor interactions | Tensor decomposition with cell–cell interactions as hypergraphs | Allows L–R pairs to function across multiple cell-type pairs (not restricted to a single-cell-type pair), which is more reflective of underlying biology | Requires curation of ligand–receptor interactions. Averages single cells to the cell-type level |
最直接的算法是共表达，即一个基因与另一个基因相互作用的可能性取决于它们成对相关系数的强度。虽然在计算上易于处理，但这些方法中的大多数不提供方向性，而且可能推断函数相关性而不是直接调节。更复杂的方法包括ode、布尔网络和BNs，如前所述，每种方法都有其优点和局限性。
布尔网络要求离散化基因表达值，并应用布尔函数来描述调控相互作用，这可能导致过度简化。
基于ode的方法使用线性、非线性或分段微分方程以连续而不是离散的方式对mRNA含量的动态特性建模。
BN是一个有向无环图(DAG)，它整合了先验信息来指导其基因-基因相互作用的预测，本质上是概率性的。
最后，信息论度量描述了生物实体之间的统计相关性，包括熵(熵是基于随机变量的不确定性对信息进行量化的概念)和互信息(互信息是对一个随机变量的观察可以告知或减少另一个随机变量的不确定性)。这种方法产生了更一般的关联，允许捕获非线性依赖，并在网络推理中被普遍使用。
值得注意的是，由于新的方法正在迅速发展，不可能详尽地记录所有现有的方法。在这里，我们强调了单细胞GRN建模的广泛类别，并讨论了示例方法来说明这些概念，并注意到它们的优点和潜在的局限性。我们还排除了基于旧的低通量单细胞平台(如单细胞qPCR)数据开发的方法，这些方法与稀疏高通量的scRNAseq没有相同的挑战。
动力学网络
到目前为止，大多数基于scRNAseq的GRN建模方法被设计用于处理动态细胞状态转换(图1B)，因为scRNAseq数据包含来自显示时间动态的异步细胞种群的信息，允许在拟（伪）时间（pseudo-time）尺度上映射细胞转换。表达动力学或伪时间估计的常用模型假设细胞变化(即发育、激活和失活)沿着连续曲线或理想化树进行，每个中间阶段都很短，并通过对大量细胞的测序得到。在这些假设下，计算建模可以推断细胞动力学的轨迹，可以根据已知的调控关系，如TF靶标信息、基因表达的相似性、以未成熟和成熟mRNA含量为代表的RNA速度来推导。但是，需要注意的是，在给定快照上同时出现的各种细胞状态并不代表序列或沿袭信息推断的实时过程。因此，加入伪时间不一定能改善GRN的构造。
到目前为止，已经开发了超过50种方法用于轨迹推断来推导伪时间信息，这些方法之前已经被回顾和比较过。伪时间排序为动态GRN建模提供了方向性和交互性信息。这种伪时间信息与上述常用的网络构造算法如correlation ， ODE ， Boolean ， BN，信息论和其他方法相结合。许多动态GRN方法已经被其他人广泛地回顾过，我们在这里只讨论不同类别中的几个例子。
一种布尔网络方法，即SCNS，是基于有序细胞之间的单基因变化，细胞被离散到一个开/关的状态。
另一种方法SCODE使用线性ODE(一种假设所有细胞都在同一轨迹上的伪时间估计)和基于TF的框架来建模TF动力学，捕获基因间的调控关系。
在此基础上，GRISLI最近被开发出来，它使用了与SCODE类似的方法，但考虑了多个细胞轨迹，不采用网络结构，计算时间更快。GRISLI首先估计每个细胞的速度，然后解决一个稀疏回归问题，将细胞的基因表达与其速度分布联系起来，以估计GRN。
以信息论为基础的方法，SINCERITIES，利用Granger因果关系获取方向性信息，量化每个基因表达在两个后续(伪)时间点之间的时间变化。通过岭回归分析，利用TF表达的变化预测相应基因在下一个时间窗口的变化，通过对每个基因对的表达进行偏相关分析，推断边缘方向和符号。SCINGE还对有序单细胞数据使用基于核函数的因果回归来预测调控因子与靶基因的相互作用，然后对回归结果进行聚合，对预测的相互作用进行排序。
另一种方法是PIPER，它使用局部泊松图形建模来更有效地捕获细胞分化过程中的网络变化，并突出显示驱动这些变化的关键TFs。
NB推理方法，AR1MA1-VBEM(变分贝叶斯采用)，应用一阶自回归移动平均(AR1MA1)模型适合代表观测时间序列的线性模型的组合数据前面的计算和噪声项，并使用一个问题的框架，利用变分法来优化网络模型的边际似然和后验分布。
Scribe是另一种最近发展起来的方法，它使用限制性定向信息(RDI)，通过引用相关的时间序列数据或从内含子(指示未成熟RNA)和外显子读子中推断细胞速度来推断因果grn。作者证明当有真实的时间序列数据时，Scribe优于其他伪时间方法；然而，当测量的时间信息丢失时，所有方法的性能都会受到极大影响。有趣的是，Deshpande等人最近对各种方法进行了比较，发现加入伪时间并不一定会带来更好的性能，但在某些情况下会损害网络重建。如前所述，这可能是由于伪时间方法的假设存在问题造成的。
细胞内网络
第二类方法侧重于在不考虑细胞轨迹或动力学的情况下对细胞内群体的grn进行建模。这些方法包括共表达和基于tf的方法、共表达和不依赖tf的方法以及信息论方法(表1和图1B)。这符合组织基因-基因相互作用GRN建模的基本概念，除了这里为特定细胞群建模的单细胞数据。
与动态网络建模相似，对细胞内群体建模的最简单方法是基于共表达。在这里，共表示方法分为两组：利用TFs形式的先验信息的方法和不利用TFs形式的方法。对于与tf无关的方法，一个基因与另一个基因相互作用的可能性取决于它们成对相关系数的强度，并考虑了所有可能的基因对。在基于tf的方法中，根据与不同TFs的成对相关系数最强的基因分组到模块中，或者根据之前的文献或motif证据分离到潜在的相互作用。定义细胞内grn的一种更复杂的方法是部分信息分解，它可以捕获非线性的基因依赖性。在这里，由一对基因提供的信息被用来量化所有三组基因中关于第三个基因的独特的、共享的和协同的信息，从而推断出一个网络结构。
几种基于相关的方法已经被开发出来，用来比较已知或预测的TFs与靶基因或所有基因之间的基因表达模式。例如，
通过对共表达基因模块进行带有TF结合基序的SCENIC couples基因共表达分析，以识别GRN模块，预测TF调节因子，并识别假定的TF靶标(称为调节因子)的单细胞水平活性。这些调节因子的活性可用于群集细胞类型，比较网络保护，并确定参与疾病的重要细胞状态和grn。
另一种方法是使用完整的分析管道来处理scRNAseq数据。它首先识别每个细胞类型的候选TFs和它们的靶标。然后利用基因表达的一阶条件依赖性来确定两个或一个转录因子与目标基因之间的相互作用，并通过整合六个不同节点重要性指标来确定每个GRN中的关键转录因子。
其他的coexpression-based GRN方法，使用一个细胞类型特异的GRN正交化方法构建基于细胞的功能标识的关键假设是由一组弱，但具体表达基因介导的一组TFs。
ACTION将每个细胞描述为高维空间中的一组“细胞函数”，这些函数的数量使用非参数方法确定。使用正交化法确定每种细胞功能所特有的基因，并评估TFs在控制这些细胞功能基因中的作用。细胞内的TF和相关靶基因构成了这个网络。
Pina等和最近的Iacono等也利用共表达构建了不限于TF目标关系的全局GRNs。前者计算一个细胞类型内所有细胞间的Spearman等级（rank）相关性，以推断造血过程中的细胞型GRNs，并利用线性转化表达数据的比值比确定显著的成对关联。Iacono等人使用了一种基于皮尔逊相关的方法，该方法首先使用bigSCale转换表达值，使用概率模型推导出每个基因的z分数，以考虑单细胞数据固有的噪音和变异性。利用z得分的两两相关关系来构建grn。z分数的使用增加了显著的基因对基因的相关性。
为了揭示简单的相关策略所不能提供的复杂的基因依赖关系，GRN推理方法采用了信息论的技术。具体来说，PIDC使用部分信息分解，在所有其他可能的基因中找到任何一对两个基因所提供的唯一信息。这种多元信息的方法利用第三个基因之间的依赖关系识别非线性双基因关系。
细胞通信网络
一个给定的异质组织的基本功能不仅由组织内不同细胞类型的活动决定，而且由细胞群体之间密切的沟通和协调决定。例如，神经元和星形胶质细胞相互作用以保证大脑的基本功能，免疫细胞与脂肪组织中的脂肪细胞相互作用以调节能量代谢和产热。因此，细胞-细胞间的通讯是一个关键的生物学问题，但由于之前缺乏高通量、高分辨率的单细胞数据而尚未得到全面的解决。单细胞方法同时捕获多种细胞类型的独特能力，使得建立细胞-细胞通信网络模型成为可能。建立这种网络模型的基本假设是，细胞之间的通信可以通过测量单个细胞群体的分子模式来捕获。例如，一对相互联系的细胞可能以协调的方式表达参与特定功能的基因和蛋白质(例如，一个细胞表达配体，另一个细胞表达相应的受体，以触发信号通路)。
早期对细胞-细胞通信网络模型的尝试主要是基于基因共同表达的概念，无论是否考虑配体-受体的相互作用信息。潜在的假设是细胞之间的基因相关模式反映了真实的生物相互作用。在组织-组织相互作用的水平上，已有证据支持这一假设的有效性。例如，脑区域间的基因共表达可以概括出小鼠脑连接体功能衍生的相互作用，5种不同小鼠组织间的基因共表达揭示了介导沟通的新内分泌因子，这些新因子随后被实验验证。
当Han等人基于不同细胞类型的基因表达谱的相似性建立细胞-细胞连接时，共表达方法迅速适用于单细胞数据。然而，这些网络更有可能反映细胞类型之间的相似性，而不是相互作用或通信。为了修改经典的共表达框架，基于配体和受体的方法已经被提出，它依赖于这样的假设:细胞间通信的重要部分是通过释放化学分子从一个细胞结合到另一个细胞的受体。利用这个假设允许基于配体受体的方法来构建可靠的基于生物的定向网络。然而，这样做的代价是在固有的稀疏数据模式中严重限制了潜在基因的集合。值得注意的是，基于共表达的分析通常使用皮尔逊相关系数，由于零膨胀的性质和独特的分布模式，皮尔逊相关系数可能不适用于基于读取的单细胞数据集。在对单细胞数据使用基于共表达的分析时，重要的是要考虑数据转换和适当的统计。
有几种方法说明了通过配体-受体相互作用的细胞-细胞通信。
Zhou等人编制了一份>25000对已知配体受体的清单，以检测它们在关于4000个黑色素瘤细胞的转录组中的变化。为了确定一对细胞是否在交流，配体和相应的受体必须在这两种细胞中表达超过一定的可调阈值。
类似地，Kumar等人关注的是基于约1800文献的配体-受体对，但采用了不同的评分方案，考虑了各自被测细胞类型中平均受体表达和平均配体表达的产物。
iTALK是另一种新的基于配体-受体交互的网络构建方法，它被移植为带有数据可视化工具的R包。对于每一个细胞类型和iTALK数据库中的配体对，iTALK标识中的配体-受体对(> 2600对)两个细胞类型之间通过询问排名列表基因来源于平均差异表达基因(单一的计算/条件)或(多个时间点/条件)。此外，iTALK还能够使用元数据(例如，时间点、组和群组)通过识别不同表达的配体-受体对来发现细胞-细胞相互作用的变化。
类似地，Smillie等人使用了FANTOM5数据库中数千篇支持文献的受体-配体相互作用来识别细胞-细胞相互作用，要求基因是细胞标记基因或差异表达基因来表示细胞间的重要相互作用。
在大多数配体-受体方法中，配体-受体对仅限于细胞类型；然而，在scTensor中，Tsuyuzaki等人采用了更灵活的方法，不存在此类限制。在scTensor中，细胞-细胞相互作用被表示为超图，超图描述了用张量分解确定的配体-受体对的有向边。
Vento-Tormo等人最近提出的一种方法也考虑了分泌分子和细胞表面分子，并使用基于置换的方法来寻找细胞类型之间丰富的配体-受体对。为了实现这一点，作者开发了CellPhoneDB，一个配体-受体相互作用的公共知识库，由蛋白质-蛋白质相互作用的公共资源管理，其中包括配体和受体的亚基组成，以充分表达它们的相互作用。对于由多个亚基组成的蛋白质，需要表达所有亚基才能推断出准确的相互作用。
上述方法都只专注于配体-受体对，这依赖于假定的基因，使其局限于一组可通知细胞-细胞通讯的基因对。此前，一种限制较少的建模方法已经被开发出来，该方法基于编码源组织分泌肽和靶组织中所有基因的共同表达，来解剖组织-组织通信网络。Arneson等人采用这一概念，在假小鼠和脑外伤小鼠的海马中构建细胞-细胞通信网络图，揭示了脑损伤中广泛的网络重新布线。这种方法推断出细胞之间的联系，其基础是假设一个细胞通过分泌信号分子与另一个细胞通讯，这些信号分子与靶细胞上的受体结合，从而触发靶细胞的下游分子事件。因此，可能在源细胞类型中编码分泌信号分子(即配体)的基因与受体以及靶细胞类型的下游通路基因之间存在共表达。通过考虑细胞类型之间所有表达基因的模式，其他方法可以将细胞-细胞相互作用的范围扩展到基于配体-受体的关系之外，尽管对这种方法的生物学解释不是那么直接。
综合的方法
wang et al提出了SoptSC，一个统一的框架来进行单细胞分析从基因表达矩阵基本分析工作流(如标准化、集群、降维,并确定细胞标记基因)，随后推断信息交流网络和pseudotemporal。SoptSC的关键前提是结构化的细胞间相似度矩阵有助于改进网络推理步骤。相似矩阵也被用于伪时间排序，在加权簇到簇图中寻找细胞之间的最短路径。为了推断细胞-细胞信号网络，根据配体-受体对的表达和下游通路靶基因的方向来计算两个细胞间相互作用的可能性估计值。通过总结任意两种细胞类型的所有细胞之间的信号转导概率，就可以形成一个集群/细胞类型之间的一致网络。
基因扰动网络
上述方法都是利用TF级联、配体-受体关系等信息流假设，没有直接的因果信息。含有基因扰动信息的单细胞数据对于提供GRN构建的因果信息极为有用，因为一个基因的靶向扰动是其他基因下游反应的来源或触发器。Jackson等人提出了利用基因缺失突变体的方法。具体地说，他们汇集了横跨12种不同基因型(TF缺失)和11种不同条件的72种不同酵母株，生成了38000个细胞的scRNAseq数据。除了表达数据，该方法使用来自TF目标和生物物理参数(如TF活性和mRNA衰减率)的先验信息，使用多任务学习(MTL)框架构建GRN。这允许在不同的条件和实验中整合信息，以解释TF扰动和观察到的基因表达变化之间的关系。通过直接删除TFs，作者创建了一个有价值的数据集，可以作为其他单细胞网络推理方法的有用基准。Wang等人提出了一种推断因果DAGs的算法。将CRISPR/ cas9介导的基因扰动与单细胞测序相结合，生成高通量的介入基因表达数据。该算法基于贪心SP来限制基于置换的DAG搜索空间，利用贪心干涉等价搜索来评估潜在的网络分数。为了进一步扩展因果网络推理的研究，Wang等人引入了一种方法，可以识别从不同数据集推导出的DAGs之间的差异。同一组也表明,软干预用于Perturb-seq，比如那些导致局部破坏的基因相关性(例如RNAi或CRISPR-mediated基因激活)，提供相同数量的因果信息困难的干预(例如CRISPR / Cas9-mediated基因删除)，导致完全中断，尽管只是轻微扰动。
单细胞GRN建模方法的性能评估
Chen和Mar最近将一些单细胞网络建模方法，包括SCENIC、SCODE和PIDC，应用于模拟和经验的单细胞数据集，以评估其捕获已知网络交互的能力。他们发现这些方法之间的一致性很低。然而，由于每一种方法都有独特的假设，并且可能不被设计来捕捉类似的交互作用，因此方法之间的一致并不一定适合于评估性能。另一项比较研究考察了包含伪时间信息的多种网络推理方法(如SCINGE、SCODE和SINCERITIES)的性能，也表明许多调节目标预测对于每一种被测试的方法都可以接近随机。
这些发现需要对单细胞网络建模方法进行改进，并对现有的单细胞GRN方法的性能进行全面评估。另一方面，由生物学假设和数据驱动的基因共同表达驱动的配体-受体框架似乎在细胞-细胞通信网络建模中很有前景。例如，用这种方法对scRNAseq数据进行建模，再现了海马体内已知的细胞-细胞相互作用。
理想与现实之间
单细胞多组学分析技术正在迅速发展，带来革命性的力量以提高我们对生命的基本单位----细胞-----以及在生理和病理条件下细胞之间的综合了解。在更准确地分类细胞类型、纠正混杂因素、描述细胞谱系和细胞状态转变等方面取得了重大进展。然而，这些进展还不足以使我们完全理解个体细胞群功能的调节机制，以及决定更高水平组织功能的细胞-细胞相互作用。现有方法模型基因网络优化的主要组织数据为单细胞数据表现不佳或不能适应新的生物单细胞数据,带来的问题和方法，有效地和精确地模型流出单细胞数据到全面的RGNS的图谱还在初级阶段。特别是，目前仍急需新的网络方法来解决单细胞数据的独特挑战，如数据稀疏性、多模态分布和更高维数。数据稀疏性问题可以通过改进单细胞技术来增强信号捕获，或者通过更精确的注入方法来解决，这些方法得到了强有力的实验验证数据的支持。这些努力将有助于缓解与非标准数据分布相关的问题，这些问题限制了现有网络方法的使用。另外，建立在更合适的统计数据和算法上的方法可以更好地适应dropout值和独特的数据分布，这是有必要的。
在单细胞数据的网络建模中，另一个重要但不太突出的缺陷是缺少空间信息来约束建模空间。目前许多高通量的单细胞测序方法缺乏保持单个细胞的空间身份的能力，这降低了准确解析细胞网络的能力，特别是在开发阶段。各种高通量荧光原位杂交(FISH)方法已被开发为解决空间信息的工具。假设细胞间的距离越近，就越有可能进行通信，可以利用成对单细胞间的空间距离作为建立更复杂、更准确的网络模型的先验。最近发现，产生配体的细胞与表达相应受体的靶细胞直接相邻，这一假设得到了支持。基于单分子鱼的方法的另一个关键优势是，它们是非常定量的，并且不会出现中断，而这种中断会困扰基于高通量单细胞测序方法。在空间单细胞方法中，也可以将表型(即行为)与细胞激活(即cFos)结合，在假设特定表型或刺激中活跃的细胞更有可能进行交流的前提下，整合到模型中。此前，Moffitt等人曾使用这种方法来识别在养育过程中激活的神经元。因此，将单细胞测序方法与高通量单分子成像相结合，在提高单细胞分辨率下的网络建模方面具有巨大潜力。尽管有潜力，但使用空间数据构建GRNs仍存在局限性和复杂性。首先，基于单分子fish方法的细胞分割是非平凡的，没有它GRN的构建是不可能的。此外，单个图像承载有限的动态细胞景观表示。事实上，许多这些技术只能实现单个细胞的成像深度，因此它本质上是一个给定时间的二维快照，可能无法捕捉到在成像平面和时间框架之外的细胞动态。
目前，大多数方法都是为scRNAseq设计的，需要结合其他单细胞组学指标(遗传、表观遗传和蛋白质)的方法。这与大组织GRN推断所面临的挑战相同，而多组学集成和建模的最新进展可能为单细胞多组学建模提供指导。
最后，从经验数据中预测的网络的准确性很难评估，因为通过在体内单个细胞中干扰预测的调节因子进行高通量验证比进行全身敲除或敲除更具挑战性。值得肯定的是，新的高通量基因干扰技术，如Perturb-seq 与scRNAseq的结合，有可能深入了解基因和细胞之间的真正关系。来自这些平台的数据可以作为更合适的基准数据集，通过测试每种方法检索扰动-响应实验中已知的真实调节关系或交互关系的效果，来评估现有网络方法的预测。同样，利用文献中已知的、经过实验验证的基因-基因、细胞-细胞通路可以作为这些方法的基准。即使在没有经过验证的网络连接的情况下，也可以采用基于社区的方法，通过结合多种方法推断出的多个网络来获得一致的网络，从而提高网络性能。这种方法已经被证明对提高预测网络的质量非常有价值。
总而言之，我们正在进入一个黄金时代，在这个时代，生物发现可以以前所未有的分辨率和通量进行。单细胞多组学数据的网络建模是解开病理生理学背后复杂的分子机制和指导精准医疗的关键工具之一。尽管面临诸多挑战，但该领域正在迅速发展，方法创新的大量机会等待着更准确地描绘健康和疾病细胞的分子图谱。
