当前，对细胞动态过程的研究表明，细胞在动态过程中会发生状态变化，产生表型稳定的细胞类型以及在表型之间转变的过渡类型，这一系列转变主要由细胞内部的基因表达情况控制。在生物意义上，识别生物组织谱系不仅提供关于正常组织发育和体内平衡的信息，还提供了关于癌症等病症的相关信息。而了解组织形成的谱系，本质上是研究细胞的动态发展过程。传统的研究方式是通过细胞的可遗传标记进行追踪，然而这些标记由于数量有限，可能会掩盖标记基因在细胞亚群中的异质性，这往往会导致细胞分化研究中的错误结论。由于单细胞实验可以产生横截面群体快照，能够反映细胞在其分化或过渡状态中的变化，为解决精确描绘细胞分化问题提供了新的方式。2009年，Tang等^[1]首次发表单细胞RNA-seq技术，打破了传统研究使用细胞群体均值来表示基因表达水平的情况，能够在单细胞水平上获得细胞真实的基因表达信息。此后，相关技术迅猛发展，2012年Hashimshony等^[2]开创了cel-seq技术，采用线性扩增的测序方法，数据的错误率比较低；2013年Picelli等^[3]发表的smart-seq2技术能够获得全长的转录测序数，测得的数据具有更好的敏感性，能够在每个细胞中检测到更多基因。随着这些高通量测序技术的发展，研究者能够获得在单细胞水平下大量的基因表达数据，使得从单细胞层面了解细胞状态和分化发育提供了可能。

但是，“缺失值”问题对人们分析与解读这些数据造成了困扰。事实上，早在群体组学数据时期，已经有学者关注“缺失值”问题。2006年，Nie等^[4]提出了一个基于数据驱动的零膨胀泊松回归模型，通过转录组数据应对蛋白组数据的“缺失值”问题。而对于单细胞转录组数据，scImpute^[5]与scTSSR^[6]方法都采取了借助其他基因表达信息的方式来应对“缺失值”问题。

此外，对单细胞数据的分析与解读已经使得人们在研究包括癌症在内的内在异质性疾病上取得了许多突破。2012年，Torres-García等^[7]在单细胞水平上对细胞耗氧动态数据进行多参数分析、建模和特征提取，证实了不同类型细胞的差异。2014年，Bendall等^[8]开发了Wanderlust算法，基于k-近邻图与最短路径算法构建细胞发展轨迹，成功构建了从造血干细胞到人类B淋巴细胞的轨迹，准确预测了其发育轨迹。2016年，Grün等^[9]提出了一种从单细胞转录组数据中推导谱系树的StemID算法，通过使用k-medoid来创建细胞簇聚类中心在高维空间中进行连接，将单细胞投影到网络上的边来生成谱系轨迹网络，能够在种群内所有可检测的细胞类型中识别出干细胞。2017年，Guo等^[10]开发了Slice算法，用转录组熵作为分化的指标值，通过将细胞聚类进行简化，并以簇为单位构建最小生成树来形成分化轨迹的大框架，最终形成一个囊括所有细胞的分化轨迹。2018年，Lummertz等^[11]开发了CellRouter算法，通过降维将基于k-近邻图构建的网络图转换成流动网络图，将多状态谱系构建转化为最优化求解问题，探索复杂的细胞状态转变轨迹。2020年，Wei等^[12]开发了一种基于扩散传播的轨迹推断方法Dtflow，通过构建细胞的k-近邻图，利用高斯核函数将其转换为马尔科夫转移矩阵，进而转换为巴氏核矩阵，基于该矩阵进行降维和伪时间排序，由于减少了对降维过程的依赖性，大大增强了分化轨迹的稳定性。

然而，这些算法需要使用除基因表达以外其他的信息，例如细胞的k-近邻图、细胞聚类和多状态谱系等，带来了额外的复杂度和不确定性。当然，也有一些学者为了避免额外的复杂度，只使用基因表达数据进行研究。2020年，Liu等^[13]基于恶性细胞较参考细胞间的差异分布提出了单细胞熵，结合高斯混合模型实现了恶性细胞与良性细胞的区分。2020年，Zhong等^[14]提出了单细胞图熵(single-cell graph entropy, SGE)，在5个胚胎分化数据集中识别出了暗基因，即不同时期在基因水平表达无差异但对SGE值敏感的基因。在此基础上，本文构建细胞特异性网络，并使其与熵相结合，提出了基因互作网络熵(gene interaction network entropy, GINE) 的方法来估算细胞的分化能力，以此来研究细胞的发展动态。通过GINE，笔者完成了从借助单细胞“不稳定”的基因表达到借助相对“稳定”的GINE值研究细胞动态的转变，通过网络的稳定性在一定程度上规避了“缺失值”问题的影响。文章将GINE方法应用于头颈部鳞状细胞癌和慢性髓细胞白血病等2个单细胞RNA-seq数据集，不仅可以将恶性细胞与良性细胞分类、显著区分不同分化时期，而且通过GINE值有效反映了疾病疗效过程。相关结果证明了GINE方法的有效性以及其研究细胞动态的潜力，可以用于细胞分化、追踪癌症发展以及疾病对药物反应过程等研究。

本文选取的两组真实数据集都来自美国国立生物技术信息中心基因表达数据库(gene expression omnibus, GEO)。第一组数据为头颈部鳞状细胞癌(head and neck squamous cell carcinoma, HNSCC) 患者的单细胞RNA-seq数据集(编号为GSE103322)，包含了来自18例患者的5 902个细胞。第二组数据是慢性髓细胞白血病(chronic myeloid leukaemia, CML) 患者的单细胞RNA-seq数据集(编号为GSE76312)，反映了经酪氨酸激酶抑制剂(tyrosine kinase inhibitor, TKI) 治疗CML过程中不同阶段的基因表达情况，该数据集包含2 055个癌症干细胞，囊括了诊断期以及经TKI治疗后不同阶段的基因表达情况。对于这两组数据，本文在使用前作了一定的质量控制(quality control, QC)，针对数据中存在多个同名基因的情况进行整合，取其均值作为其表达水平。此外，数据中存在许多基因的计数为零的情况，我们选择保留在10或更多细胞中表达的基因。

在先前的研究中，如SLICE算法，通过借助蛋白互作网络(protein-protein interaction networks, PPI) 中蛋白质i与蛋白质j间的关联强度，定义了对应编码基因之间边的权重。然而，由于不同类型细胞中的蛋白互作网络的差异或是缺失，会带来对应基因间权重的不确定性。为此，本文构建一个基因互作网络(gene-gene interaction network, GGI)，对于给定的基因i与基因j，假定两者之间的关联程度表现为它们在整体样本中的皮尔逊相关系数，并简记为r_ij。相关系数r_ij的取值范围为–1到1，非零值表示2个基因间存在相互作用关系，其绝对值则定义为相互关联强度。为了降低单细胞测序数据中的噪声干扰，首先，进行相关性检验，保留具有统计学意义(取P值小于0.05) 的相关系数作为GGI网络中的边；其次，本文认为不是关联强度非零就能表征基因相互关联，并假设小于阈值的关联强度值是噪声带来的影响。因此，在相关系数检验后进一步进行筛选，通过对已知时序数据进行分类的准确度来进行调整，选择使其分类准确度最高的值作为阈值。对于绝对值大于阈值的，两基因间存在相关性，反映为GGI网络中存在关联边。反之，则不存在关联边。对于给定的M行N列基因表达矩阵，本文先进行对数化处理，然后通过R语言Hmisc包中的rcorr函数进行计算，最终得到的M行M列的基因关联矩阵(r_ij)_M×M即可构建本文需要的GGI网络。

对于给定的M行N列的基因水平表达矩阵E_M×N，其中M代表基因维数，N代表样本维数。首先将其进行标准化处理，得到新的矩阵EN=log₂(E+1)=(e_ij)_M×N。基于先前构建的GGI网络，笔者构建了样本特异的基因互作网络熵值的矩阵H=(GINE_ij)，算法流程可参见图 1。具体而言，首先将e_ij定义为给定样本j中基因i的表达水平，将N(i)定义为GGI网络中以基因i为中心的邻近点集合(含基因i本身)。对于给定的邻近点集合k∈N(i)，本文假定基因k与基因i的关联强度近似于乘积r_ik×e_ij×e_kj。进一步地，进行归一化处理以确保，得到

(1)

图 1 基因互作网络熵计算流程示意图 Fig. 1 Schematic diagram of the GINE calculation process.

图选项

如果k∉N(i)或者k∈N(i)且e_kj=0，那么p_ik=0。接着，可以借此定义GGI网络中基因i的样本特异的基因互作网络熵，计算公式如下：

(2)

其中s_ij反映了GGI网络中以基因i为中心的局部网络的不确定性水平。同时，考虑到不同局部网络之间中心结点度数的差异，定义了归一化的局部网络熵：

(3)

其中n为集合N(i)的模。在此步骤后，本文得到基因i在样本j中对应的基因互作网络熵GINE_ij，该值不仅取决于局部网络中心基因的表达值，而且受邻近基因影响。最终，得到样本特异的基因互作网络熵矩阵H=(GINE_ij)_M×N。

为了证明基因互作网络熵的适用性，本文测试了HNSCC数据集。该数据集包含了来自不同患者共5 902个细胞，其中的癌细胞通过一种基于拷贝数差异的方法进行识别。为了减少计算量，本文随机选取HNSCC数据集中1 000个样本进行区分癌与非癌细胞的测试，为了避免取样带来的偏差，本文采用简单不重复随机取样从全部数据的矩阵中选取了维度为1 000×1 000的基因-样本的子矩阵，进行了多次测试。下文的图 2A、C、E反映了基于原始基因表达矩阵进行区分的结果，图 2B、D、F则对应以基因互作网络熵矩阵进行区分癌细胞与正常细胞的结果，这里组别1指代HNSCC，组别0指代正常细胞。事实上，图 2A、C、E达到了大致区分癌与非癌细胞的效果，但在降维形成的聚类交界处无法使两者彻底分隔，在图 2B、D、F中，癌与非癌细胞之间则存在更为明显的分界。此外，本文还比较了每一个细胞的GINE期望、方差与其T-分布邻域嵌入(T-distributed stochasticneighbor embedding, Tsne) 降维分析的关联性，发现GINE期望值与Tsne1分量呈正相关，如图 3A所示，类似地，图 3B也证实了细胞的GINE方差与Tsne2分量呈正相关。因此，本文认为细胞的GINE期望与方差可以视为GINE熵矩阵降维后的一维随机邻居嵌入，表明GINE方法是对于单细胞RNA-seq数据进行维度降低的有效方法。

图 2 HNSCC数据基于基因表达矩阵和基因互作网络熵矩阵的癌与非癌细胞区分结果图 Fig. 2 Results of differentiation between cancer and non-cancer cells based on gene expression matrix and GINE matrix for HNSCC data. (A, C, E) Visualization results based on gene expression matrix. (B, D, F) Visualization results based on GINE matrix.

图选项

图 3 HNSCC细胞的GINE统计数据与Tsne的散点图 Fig. 3 Scatter plots of GINE statistics versus Tsne for HNSCC cells. A：HNSCC细胞的GINE均值与Tsne1的散点图；B：HNSCC细胞的GINE方差与Tsne2的散点图 (A) Scatter plot of GINE mean versus Tsne1 for HNSCC cells. (B) Scatter plot of GINE variance versus Tsne2 for HNSCC cells.

图选项

Tsne是目前来说效果最好的数据降维与可视化方法，已广泛用于对单细胞数据的分析；为了阐明GINE值与基因表达值间的表现差异，本文将对真实数据集GSE76312的基因表达矩阵和基因互作网络熵矩阵进行Tsne降维可视化效果对比。图 4是降维可视化后的结果，图中紫、绿、蓝、红、黄色的分组分别代表着诊断时期以及CML患者经TKI治疗3、6、12、18个月时的细胞，图 4A是基于基因表达值矩阵进行降维可视化的结果，图 4B则是基于基因互作网络熵矩阵的降维可视化的结果。显而易见的是，基于基因互作网络熵值的降维方法成功区分出了不同时期的CML细胞，然而基于基因表达水平的降维无法做到，不同时期的细胞完全混合在一起，无法进行有效区分。通过基因互作网络熵方法对数据集的聚类分析结果表明GINE具有比原始基因表达更好的表现，证实了基于GINE值在时序聚类方面的优越性，验证了基因互作网络熵在反映生物学过程差异方面的可靠性，从而为探索细胞群体的动态信息提供了新的方法。

图 4 CML数据的Tsne降维效果图 Fig. 4 Plot of the effect of Tsne dimensionality reduction for CML data. A：基于基因表达矩阵的降维可视化图；B：基于基因互作网络熵值矩阵的降维可视化图 (A) Dimensionality reduction visualization plot based on the gene expression matrix. (B) Dimensionality reduction visualization plot based on the GINE entropy matrix.

图选项

本文又研究了基于GINE值不同阶段的CML细胞异质性情况，如图 5所示。值得注意的是，GINE均值可以反映细胞的多能性程度，即细胞分化成所有主要细胞系的能力。细胞多能性越高，越不会表现出对特定细胞系信号路径的偏好，因此所有信号路径会具有相似的活性，对应的GINE均值越高。对于一个分化的细胞，或者对于一个致力于特定谱系的细胞，其不确定性是降低的，因为这需要激活反映该谱系选择的特定的信号途径，对应的GINE均值越低；由于终端分化细胞的多能性是低于肿瘤干细胞的，因此，GINE均值下降，表明细胞群体越靠近正常细胞状态。按照不同时期计算了转化后细胞的GINE均值，如图 5A所示，CML患者在TKI治疗后的3–12个月，细胞的GINE均值整体呈下降趋势，反应了TKI治疗出现好转的情况。而TKI治疗18个月后，GINE均值的提高表明整体的状态又在偏离正常细胞状态，因此代表细胞产生了耐药性。图 5A显示，细胞在诊断时期的GINE均值最低，这表明细胞不进行基因筛选就计算GINE均值是不合适的。事实上，由于基因的选择性表达，基因会处于开启或者关闭的状态，而细胞的状态主要是由处于开启状态的高表达基因所决定，这需要选取反映细胞真正状态的高表达基因。为确保正确反映细胞状态，更好地呈现不同时期细胞对TKI的反应状态，本文又对不同时期细胞分别取其前5%和1%的GINE均值并考察其变化情况，如图 5B和5C所示，变化曲线与整体情况大致保持一致，且曲线变化更加明显。由于前1%的GINE均值更能凸显TKI治疗后不同时期的状态变化，在取前1%的GINE均值作为细胞状态的基础上，计算同一时期所有细胞的GINE均值(mean of stage, MOS)，以此作为该状态的指标值。图 5D是以MOS曲线来衡量CML患者经TKI治疗后不同时期的状态变化，横轴上的数字1–5分别代表诊断时期，以及CML患者经TKI治疗3、6、12、18个月等5个时期。如图 5D所示，整体曲线呈下降趋势，尤其是在第2个和第3个时期，MOS值大幅度下降，而在第5个时期MOS值开始有回升迹象。因此，有理由相信CML患者经TKI治疗后3–6个月是药物反应最有效的时期，而到18个月将可能产生耐药性。这个发现为研究疾病对药物的反应，进而研究细胞的动态提供了方法上的支持。

图 5 基于基因互作网络值不同阶段的CML细胞异质性情况 Fig. 5 Heterogeneity of CML cells at different stages based on GINE values. (A) Change in mean GINE value of CML cells. (B) GINE mean change of the top 5% genes in CML cells. (C) GINE mean change of the top 1% genes in CML cells. (D) Change in mean GINE values across time. Numbers 1–5 on the horizontal axis represent the 5 periods of CML, representing the period of diagnosis, and the periods of 3, 6, 12 and 18 months after treatment with TKI, respectively.

图选项

在上一小节中，本文识别了CML患者经TKI治疗3–6个月为药物反应最有效的时期，故将该时期前5%差异表达基因作为GINE关键基因集，当然这些基因是本文选取P值为0.05得到的。为了进一步研究验证这些基因的生物功能，使用在线生物信息学数据库KOBAS，对这前5%的差异表达基因进行京都基因与基因组百科全书(kyoto encyclopedia of genes and genomes, KEGG) 通路分析，表 1为通路分析的结果中P值小于0.01的部分，图 6则是GINE基因集富集分析。根据KEGG通路分析结果，所选基因大多与代谢通路(hsa01100)、PI3K-Akt信号通路(hsa04151)、病毒蛋白与细胞因子和细胞因子受体的相互作用通路(hsa04061)、MAPK信号通路(hsa04010)、癌症途径通路(hsa05200) 等紧密相关。值得注意的是，PI3K-Akt信号通路的P值小于0.01，MAPK通路的P值也不超过0.02。事实上，PI3K-Akt信号通路会被多种类型的细胞刺激或毒性刺激所激活，并调节细胞的基本功能，如转录、翻译、增殖、生长和存活^[15]。进一步的研究表明，PI3K-Akt信号通路还参与了CML的发病机制^[16]。而MAPK通路是细胞增殖、应激、炎症、分化、功能同步化、转化、凋亡等信号转导通路的共同交汇通路之一，把胞外信号经受体、G蛋白/小G、蛋白激酶、转录因子等组成的信号网络，传递到胞内，参与细胞增殖、分化、癌变、转移、凋亡等。

图 6 GINE基因集富集分析 Fig. 6 GINE gene set enrichment analysis.

图选项

此外，经过文献挖掘，本文筛选出的诸多基因已经显示和CML或者其他肿瘤相关。例如，趋化因子受体CCR5在急性髓细胞白血病(acute myeloid leukaemia, AML) 的恶性特征确定和髓外浸润的定位方面发挥作用^[17]。血红素加氧酶1 (heme oxygenase-1, HO-1) 在多种实体瘤和慢性白血病具有细胞保护作用，它的过表达可抑制伊马替尼诱导的细胞凋亡^[18]。在加入HO-1诱导剂(cobalt protoporphyrin, CoPP) 后，NTRK1基因表达水平显著增加。基因PDGFB更是从22号染色体转位到9号染色体的证据，是判定慢性粒细胞白血病的重要依据^[19]。

在研究细胞分化与疾病发展方面，研究生命过程中细胞动态是生物学和临床重要的任务，了解这种细胞命运的变化有助于构建个体特异的疾病模型，有助于针对与细胞分化相关的复杂疾病设计具有高度特异性的疗法。在基于单细胞数据研究分化、致癌转化等的方法中，大多都需要大量的信息，例如，在蛋白作用强度与相应编码基因相关性的假设下引入PPI网络，带来了不确定性和复杂度。当前用于分析单细胞数据的大多数现有方法都是基于基因的表达，然而由于基因表达的不稳定性，对其表征生物过程带来了困难。本文通过构建基因互作网络，依据基因在细胞中的整体表达确定相互关联的基因，借助相关基因信息来应对基因不表达或者低表达的情况，即“缺失值”问题。事实上，在不同的时间、不同的条件下，网络具有更好的稳定性，与单个基因相比，通过对基因进行集聚，网络水平的基因表达更加稳定，一定程度上规避了“缺失值”事件的影响。本文开发的GINE方法，通过将不稳定的基因表达矩阵转换为相对稳定的GINE熵矩阵，可以可靠地表征动态生物学过程的状态，探索来自单细胞数据的基因-基因关联的动态信息，从而研究细胞动态。GINE为单细胞分析提供了新的计算见解，并有助于发现细胞命运的潜在信号。本文的分析结果表明，GINE在表征时间序列数据的细胞聚类方面的效果比基因表达数据要好。其次，GINE值的变化还可以识别疾病对药物反应过程中的差异，展示了时间序列数据中样本的动态变化，依据细胞GINE值的分布特征，可以识别疾病演变过程中的关键时期，进一步可以识别其关键分化途径，这有助于追踪细胞异质性并阐明细胞分化的分子机制。