pyscenic做转录因子分析时，GRN是限速步骤，数据大点就很耽误时间。此时，RegDiffusion就是一个加速选择。
以往从表达数据中推断GRN的方法常常难以应对数据的高维性和固有的噪声问题。单细胞RNA测序中普遍存在的dropout和背景噪声进一步加剧了这些挑战。作者提出RegDiffusion，一种受去噪扩散概率模型启发、专注于特征变量间调控效应的新型神经网络结构。与现有用于单细胞RNA测序数据的GRN方法不同，RegDiffusion通过扩散调度表在输入数据中引入高斯噪声，随后使用参数化的邻接矩阵神经网络训练以预测添加的噪声。
RegDiffusion的核心工作原理可以概括为：利用扩散模型的去噪过程，通过一个参数化的可训练邻接矩阵，学习基因表达数据在加噪与去噪过程中的线性关系，从而逆向推导出基因间的调控网络（GRN）。可以从下面三个层面来理解原理：
核心数学假设
RegDiffusion的成功依赖于一个关键生物学假设——线性加性假设。即假设一个基因的表达量是其所有调控基因表达量的线性组合：X = X * A（其中A就是待求的基因调控邻接矩阵）。
在此基础上，作者将扩散过程也做了线性化扩展：当原始数据X被加入高斯噪声Z变成Xt时，这种线性关系依然近似存在，即Xt ≈ Xt * A + λ_t * Z。经过简单变换可得：Z ∝ Xt * (I - A)。这意味着，模型预测噪声Z的过程，本质上就是在反向推算矩阵A。训练与推理流程
整个工作流程分为三步：
正向加噪：遵循DDPM（去噪扩散概率模型）的固定调度表，向原始的基因表达数据
X0逐步添加高斯噪声，生成一系列越来越乱的样本X1...Xt。训练时，模型只需随机抽取一个时间步t，利用重参数化技巧一次性生成加噪后的Xt，无需逐步迭代。反向去噪与参数学习：设计一个神经网络来预测步骤1中加入的噪声
Z。网络结构的关键是包含一个可训练的、参数化的邻接矩阵A（对角线强制为0以防自调控，并施加L1正则化和Dropout以增加稀疏性和鲁棒性）。网络输入Xt和t，经过嵌入和MLP（多层感知机）提取特征后，直接通过矩阵乘法h_x * (I - A)计算预测噪声，再经线性层输出。损失函数就是预测噪声与实际噪声的MSE。网络提取：训练收敛后，无需进行去噪采样，直接取出模型中学习到的邻接矩阵
A。矩阵中绝对值最大的权重值，即为预测的基因调控边（权重越大，调控关系越强）。
-
设计优势与速度提升
相比同类型的VAE（变分自编码器）方法（如DeepSEM、DAZZLE），RegDiffusion有两个显著特点：
- 省去矩阵求逆：VAE方法在重构时需要昂贵的矩阵求逆运算，而RegDiffusion通过预测噪声直接计算，避开了这一步，因此速度快了数十甚至数百倍（如论文中mDC数据集从45分钟降至3.2秒）。
- 稳定性更高：扩散过程强制数据经历一条“退化为正态分布”的轨迹，比VAE仅约束单一隐空间的KL散度更容易优化，因此在10次重复运行中标准差极小，性能稳定。
总而言之，RegDiffusion不像传统方法那样直接拟合表达值，而是利用“预测噪声”作为代理任务，迫使模型学习基因间的线性依赖，最终将该依赖关系固化在可提取的邻接矩阵中，从而实现高效且准确的GRN推断。
import regdiffusion as rd
import numpy as np
adata.layers["logcnt"] = np.log(adata.X.toarray() + 1)
rd_trainer = rd.RegDiffusionTrainer(adata.layers["logcnt"])
rd_trainer.train()
- 删除方差为0的基因
-
regdiffusion的输入数据是经过log转换的原始count - 表达值归一化：(表达值 - 该细胞最小值) / 该细胞极差
grn = rd_trainer.get_grn(adata.var_names, tf_names, top_gene_percentile=80)
# 或者
grn.remove_weak_edges(threshold=0.448)
edgelist = grn.extract_edgelist(k=50, workers=4)
edgelist.columns = ['TF', 'target', 'importance']
- gene_names和tf_names参数提供基因顺序需要与表达矩阵保持一致
- 根据权重过滤，保留top-k%的边。或者，手动
remove_weak_edges方法删除低于阈值的边 - 每个基因保留top-k个regulator/targets
RegDiffusion速度优势来自于GPU加速，如果在CPU环境中运行速度也将受限。比如下面的数据：
adata
AnnData object with n_obs × n_vars = 4257 × 21114
var: 'n_cells'
layers: 'logcnt'
rd_trainer = rd.RegDiffusionTrainer(adata.layers["logcnt"])
You specified cuda as your computing device but apprently it's not available. Setting device to cpu for now.
rd_trainer.train()
Training loss: 0.472, Change on Adj: -0.013: 10%|██ | 102/1000 [20:01<3:01:57, 12.16s/it]
从上面的运行过程可知，4257 × 21114大小的表达谱在拥有24个CPU的环境里完整运行下来预估时间为3小时，可以看到时间优势荡然无存。
