基因调控网络 (GRNs) 对于理解基因如何协调细胞过程至关重要。如今大规模的单细胞扰动研究为GRNs推断提供了机遇，但目前许多方法未能充分利用干预信息。因此，作者提出了PSGRN，该方法通过自训练框架结合合成金标准，整合了干预和观测的单细胞数据。
PSGRN推断基因调控网络的核心原理：利用“自训练”框架，将基因表达的“相关性”转化为可迭代优化的“伪标签”，并通过机器学习模型来捕捉超越简单相关性的高阶调控模式。其原理可以拆解为以下四个关键步骤：
-
生成初始伪标签
对于所有可能的基因对 (Gene A → Gene B)，PSGRN首先计算它们在表达数据中的相关性。具体计算时，会将干预数据（A基因被扰动时的B表达）和观测数据（未扰动时的B表达）合并。
- 若某基因对的绝对相关系数大于阈值，则将其初始标注为“正例”（存在调控关系），否则标注为“负例”。
- 这一步提供的只是粗略的“伪黄金标准”，并非真实标签。
-
提取基因对特征
PSGRN并不直接把表达值输入模型，而是为每一对基因<g_i, g_j>提取四个统计特征：
-
g_i和g_j在观测数据中的平均表达水平。 -
g_i和g_j在干预数据 (即人为扰动某个基因后) 中的平均表达水平。 - 这些特征包含了基因的基础表达量和对外界扰动的响应信息。
-
分类模型训练与自训练
利用步骤1生成的“伪标签”作为监督信号，步骤2提取的特征作为输入，训练一个LightGBM分类模型。
- 关键机制：模型在训练过程中，不仅看到了相关性，还学习到了基因表达水平和扰动响应模式。训练完成后，模型会对全部基因对 (包括之前作为训练集的同一批数据) 重新进行预测打分。
- 这被称为“自训练”或“自监督”——模型试图利用额外的表达水平信息，自动“纠正”第一步中相关性产生的噪声标签。
-
输出最终调控网络
模型预测出所有基因对的概率分数后，按分数从高到低排序。用户可以选择输出排名前topN的基因对，作为最终推断的基因调控网络。
PSGRN之所以有效，是因为它跳出了“直接预测因果关系”的复杂数学约束，转而采用“相关性粗略排序 + 机器学习精细矫正”的策略。当有更多的干预数据加入时，模型提取的“扰动前后表达变化”特征会更加丰富，从而使自训练过程能更有效地提升真实调控关系的排名 (文献中的秩位移分析证明了这一点)。这种方法使其在处理含干预数据的单细胞数据时，显著优于纯因果发现算法和纯相关性方法。
下面的命令行代码用来跑文章的测试数据没有问题，并不能直接用来跑自有数据：
export PYTHONPATH="."
python causalscbench/apps/main_app.py \
--dataset_name weissmann_k562 \
--output_directory output \
--exp_id psgrn_rpe1_1_1 \
--data_directory data \
--training_regime "partial_interventional" \
--partial_intervention_seed 0 \
--fraction_partial_intervention 1.0 \
--model_name "custom" \
--inference_function_file_path "./src/main.py" \
--subset_data 1.0 \
--model_seed 0 \
--omission_estimation_size 2000 \
--do_filter
因为软件里面很多代码都是针对weissmann_k562、weissmann_rpe1这两个数据集而写成了硬编码。想要自己分析使用只能动一动流程来适应自有数据了：
import warnings
warnings.filterwarnings('ignore')
import sys
sys.path.insert(0, '.')
from causalscbench.apps import main_app
ds = main_app.CreateDataset('data', filter=True)
dataset = ds.preprocess_and_save(path_k562, summary_stats_k562, "dataset_k562")
app = main_app.MainApp(output_directory='output', data_directory='data', model_name='custom', inference_function_file_path='src/main.py', dataset_name='dataset_k562', training_regime=main_app.training_regimes.TrainingRegime('partial_interventional'), exp_id='k562_result', omission_estimation_size=2000, filter=Ture)
app.dataset_splitter = main_app.DatasetSplitter(dataset, app.subset_data)
app.load_evaluators()
app.load_model()
metrics = app.train_and_evaluate()
INFO [1785481659.1561604]: Loading and preprocessing data.
INFO [1785481661.3930368]: Loading and processing evaluators.
INFO [1785482246.3828545]: Loading model.
INFO [1785482247.5626926]: Starting model training.
Cache observations: 100%|████████████████████████████| 1158/1158 [00:27<00:00, 42.76it/s]
Cache self interventions: 100%|██████████████████████| 1158/1158 [00:13<00:00, 84.16it/s]
26%|████████ | 306/1158 [1:22:19<3:43:37, 15.75s/it]
如果filter为真时需要提前准备好一个统计表格，然后流程会用其来过滤数据，否则为None。MainApp类的所有参数解释如下：
Args:
output_directory (str): Directory for output results
data_directory (str): Directory to store the datasets
model_name (str, optional): Which method to run. Defaults to METHODS[0].
inference_function_file_path (str, optional): Path to file for custom inference function. Default to empty string.
dataset_name (List[str], optional): Which dataset to use. Defaults to DATASET_NAMES[0].
model_seed (int, optional): Seed for model reproducibility. Defaults to 0.
training_regime (training_regimes.TrainingRegime, optional): Choice of training regime. Defaults to training_regimes.Interventional.
partial_intervention_seed (int, optional): If training_regime is partial intervention, seed for random selection of perturbed genes. Defaults to 0.
fraction_partial_intervention (float, optional): If training_regime is partial intervention, fraction of genes which should have interventional data. Defaults to 1.0.
subset_data (float, optional): Option to subset the whole dataset for easier training. Defaults to 1.0.
exp_id (str, optional): Unique experiment id (6 digit number). Default to randomly generated.
max_path_length (int, optional): Maximum length of path to consider for statistical evaluation. Default to -1 (all paths).
omission_estimation_size (int, optional): Number of negative samples to draw to estimate the false omission rate. If 0, the FOR is not checked.
filter (bool, optional): controls whether to select only the strong perturbations
其中model_name参数的可选值如下：
METHODS = [
"random100",
"random1000",
"random10000",
"fully-connected",
"lasso",
"random_forest",
"grnboost",
"genie",
"ges",
"gies",
"pc",
"mvpc",
"gsp",
"igsp",
"notears-lin",
"notears-lin-sparse",
"notears-mlp",
"notears-mlp-sparse",
"DCDI-G",
"DCDI-DSF",
"DCDFG-LIN",
"DCDFG-MLP",
"corum",
"lr",
"string_network",
"string_physical",
"custom",
"chipseq",
"pooled_biological_networks",
"sortnregress",
]
文章实现的方法在main.py文件中，里面有一个N变量用来控制获取的调控基因对数量，默认为1000。文章的模型底层依托于LightGBM包来实现，程序的运行速度还可以，在拥有24核的环境里面5个小时跑完测试集k562 (310385 × 8563)。
