零膨胀校正与转录动力学整合的单细胞转录组分析方法

《PLOS Computational Biology》：Integrating zero-inflation correction and transcriptional kinetics for single-cell transcriptomic analysis

　　

单细胞转录组数据普遍存在零膨胀现象，而传统模型要么忽视该问题，要么无法捕捉由转录爆发驱动的双峰分布特征，从而阻碍了准确的基因调控研究。本研究开发了一种零膨胀电报模型（zero-inflated telegraph model），将技术性零值校正与经典电报模型（telegraph model）的随机基因状态切换动力学相整合。研究人员利用合成数据、来自狼疮和乳腺癌患者的人类单细胞RNA测序（scRNA-seq）数据，以及小鼠胚胎干细胞scRNA-seq数据进行了系统性验证。该模型展现出优越性能：能够准确拟合mRNA分布（包括双峰模式），可靠估计有效转录爆发参数并防止过拟合，从而能够校正传统模型的调控推断偏差。在差异表达基因（DEGs）检测方面，该模型也优于传统方法，在小样本条件下优势尤为显著，并能识别独特的疾病相关基因（如狼疮中的LDLR、GZMB，乳腺癌中的FAIM2、VDR）。这一具有生物学可解释性和稳健性的工具推动了单细胞转录组分析的发展。

一、研究背景与科学问题

基因转录本质上是一个随机过程，即使在遗传同质的细胞群体中，mRNA拷贝数也呈现时间性波动。量化单细胞水平的mRNA分子可生成大规模mRNA拷贝数分布数据，为估计关键系统参数和阐明随机基因调控机制提供统计框架。经典电报模型通过描述基因在激活态与非激活态之间的随机转换来刻画转录爆发行为，能够有效表征细菌、酵母和哺乳动物细胞中观察到的单峰与双峰分布，并推断两个核心参数：转录爆发频率（burst frequency，即从非激活态切换至激活态的速率）和爆发大小（burst size，即每次激活期间平均产生的mRNA拷贝数）。

然而，单细胞转录组数据普遍存在零膨胀问题，其来源于两个不同层面：生物零值（由单细胞异质性固有的随机基因状态切换导致的真实不表达）和技术零值（scRNA-seq中记录充分的伪影，如mRNA捕获效率低下或扩增失败）。这种混淆使得mRNA分布中的零概率被人为放大，严重削弱了电报模型捕捉真实分布特征的能力。虽然已有研究将零膨胀参数整合进负二项分布（negative binomial distribution）模型中，但大规模单细胞数据集已证实随机转录调控天然产生双峰分布，而负二项分布只能生成单峰模式，因此这类模型无法捕捉双峰性背后的转录爆发调控机制。基于此，研究人员开发了零膨胀电报模型，将技术性零值校正与经典电报模型的随机基因状态切换动力学相整合，以同时解决零膨胀和爆发相关表达动态的挑战。

二、研究内容与核心结论

本研究构建的零膨胀电报模型在合成数据、人类狼疮患者CD14?单核细胞、乳腺癌患者CD8?T细胞以及小鼠胚胎干细胞scRNA-seq数据上进行了系统性验证。研究结果表明：该模型在拟合mRNA分布（包括双峰模式）方面优于传统模型；能够可靠估计有效转录爆发参数且不会过拟合；在差异表达基因检测方面优于传统方法，尤其在小样本条件下优势显著；并能识别独特的疾病相关基因，如狼疮中的LDLR和GZMB，乳腺癌中的FAIM2和VDR。该模型为单细胞转录组分析提供了具有生物学可解释性和稳健性的工具。

三、关键技术方法

研究人员首先采用随机模拟算法（stochastic simulation algorithm，SSA）生成包含不同细胞样本量（N=500、1000、2000、4000）和是否含20%外源噪声的合成分布数据，并引入dropout模型添加技术性零值事件。其次，利用最大似然估计（maximum likelihood estimation，MLE）方法对三个竞争模型进行参数拟合，并通过校正的赤池信息准则（corrected Akaike information criterion，AICc）进行模型选择。再次，采用Gauss-Jacobi求积法高效计算电报模型的beta-Poisson分布表达式。最后，在人类数据验证中使用了狼疮患者外周血单核细胞（PBMCs）中的CD14?单核细胞（对照组2,931个细胞，IFN-β处理组2,765个细胞，15,701个基因）和乳腺癌患者肿瘤浸润淋巴细胞（TILs）中的CD8?组织常驻记忆T细胞（CD8?T_rm，606个细胞）与CD8?非组织常驻记忆T细胞（CD8?non-T_rm，1,097个细胞，10,061个基因），以及小鼠胚胎干细胞的Smart-seq2（9,895个基因，188个细胞）和RamDA-seq（9,182个基因，419个细胞）数据。

四、研究结果

零膨胀电报模型。经典电报模型描述基因在激活态与非激活态之间的随机转换，其稳态mRNA分布可通过合流超几何函数精确表达。为提升计算效率，研究人员将其表示为beta-Poisson分布并通过Gauss-Jacobi求积法计算。当失活速率和mRNA合成速率均足够大时，电报模型退化为负二项分布。为处理技术性零值，研究人员引入平均技术零率参数，构建零膨胀电报模型，并证明了每个观测分布存在唯一参数集。

零膨胀电报模型稳健捕捉含dropout的mRNA分布。通过合成数据比较三个模型的拟合性能，发现当核心电报过程产生双峰或非零峰单峰分布时，零膨胀电报模型能正确输出双峰数据；ZINB模型虽可视觉上产生双峰曲线，但峰间谷值固定在mRNA计数m=1处，导致错位；电报模型则可能将双峰数据误判为单峰。在10?组参数集的模型选择测试中，零膨胀电报模型在有无外源噪声条件下均保持70%以上的正确选择率。

零膨胀电报模型可靠估计有效转录爆发频率和大小。在样本量N=4000时，电报模型低估爆发频率、高估爆发大小（PCC<0.5），ZINB模型呈现相反偏差，而零膨胀电报模型的估计紧密聚集在1:1参考线上（PCC=0.96）。在20%外源噪声条件下，零膨胀电报模型仍保持PCC>0.75的稳健性能。跨样本量比较中，零膨胀电报模型在几乎所有条件下均保持PCC>0.7，而电报模型和ZINB模型均无法同时可靠估计两个爆发参数。

零膨胀电报模型无过拟合。当拟合电报模型生成的合成数据时，零膨胀电报模型的参数估计紧密聚集在1:1线上（PCC接近1），其估计的技术零率极小，表明能准确区分dropout事件。当拟合ZINB模型生成的含dropout数据时，两个零膨胀模型均能可靠估计爆发参数，证实零膨胀电报模型在灵活捕捉dropout与避免参数化过度之间取得了良好平衡。

跨不同基因调控模型的泛化评估。利用三态模型（three-state model）和遗传双稳态开关电路（genetic toggle switch）生成合成数据，零膨胀电报模型在Hellinger距离（HD）上均取得最低值，在有无技术零值条件下均能准确区分生物零值与技术伪影，并能有效追踪整体转录动力学参数的相对变化。

零膨胀电报模型在DEGs分析中表现优越。基于狼疮患者CD14?单核细胞数据拟合的10,000个基因参数，模拟2,000个DEGs（单参数2倍扰动），在四种样本量下构建ROC曲线。零膨胀电报模型在所有样本量下均取得最高AUC，在小样本N=100时优势尤为显著。

零膨胀电报模型增强人类scRNA-seq数据DEGs检测。在PBMC和TIL数据集中，零膨胀电报模型取得最低HD值。技术零率随mRNA平均表达水平增加而降低。在DEGs数量上，零膨胀电报模型在PBMC中略多于ZINB模型，在TIL小样本中显著多于ZINB模型。GO富集分析显示其DEGs与其它模型共享14-16个前20位GO条目，同时富集独特条目如PBMC中的“U2型前催化剪接体”和TIL中的“内质网伴侣复合物”。

零膨胀电报模型量化小鼠细胞scRNA-seq数据中的技术零值。对于小鼠成纤维细胞Mbln2基因C57和CAST等位基因及大肠杆菌P_lac/ara启动子的双峰分布数据，电报模型将其误拟合为单峰，而零膨胀电报模型准确捕捉双峰模式。在两种测序策略比较中，Smart-seq2数据的技术零率中位数（0.14-0.16）显著低于RamDA-seq（0.3-0.51），与Smart-seq2对低输入RNA样本表现更优的既往报道一致。

五、结论

本研究开发的零膨胀电报模型成功整合了技术性零值校正与随机基因状态切换动力学，在合成数据、人类和小鼠单细胞转录组数据上均展现出优越性能。该模型能够准确拟合mRNA分布、可靠估计转录爆发参数、增强DEGs检测能力，并识别具有疾病相关性的独特基因。这一具有生物学可解释性的框架为单细胞转录组分析提供了有力工具，有助于从生物物理视角深入理解基因调控机制及其在疾病中的作用。