Identifying biological perturbation targets through causal differential networks

Thank you for your comments! Please let us know if this helps answer your questions, and if you have any further concerns.

Theoretical guarantees

We show in Appendix C that the differential network is well-specified as a model class. Due to the complexity of the architecture, it is difficult to guarantee what the model will learn, i.e. in the identifiability sense. At the same time, it is not a given that an arbitrary Transformer can approximate any arbitrary function given limited precision, width, depth [1,2]. We are inspired by the neural network expressivity literature to show that (page 15):

• Hard interventions: given the sets of edges $E, E'$, where $E'$ represents the mutilated graph, there exists a setting of its parameters that maps $(E, E')$ to targets $I$ by detecting edges that disappeared.
• Soft interventions: under "scale" interventions and non-multiplicative noise, there exists a setting of its parameters that maps the correlation and covariance matrices (which can be derived from model inputs) to targets $I$. In the linear case, Eq. 14-15 naturally hold and are relevant to the relationships between variables (i.e. if $x$ is an intervention target, for every $y\ne x$, correlation($x,y$) does not change while cov($x,y$) does).

[1] Yun et al. Are Transformers universal approximators of sequence-to-sequence functions? ICLR 2020

[2] Sanford et al. Representational Strengths and Limitations of Transformers. NeurIPS 2023.

Alternative causal structure learners

CSIvA [1] and AVICI [2] are supervised causal discovery algorithms that could work in principle.

• CSIvA does not have a public codebase or supplemental files.

• AVICI is available and runs well on smaller datasets, but there are some practical challenges. Since the Transformer attends over raw data, its complexity scales quadratically as the number of samples. Also, AVICI was only trained on datasets of up to 50 variables and 200 samples. On graphs of 100 variables (Table 7 in [2]), the performance drops to that of classical algorithms (GES, GIES).

We chose [3] since it reported reasonable performance on larger datasets (CausalBench, 622 variables). The use of summary statistics like correlation also mirrors common practices in single cell sequencing, since the raw data are often very noisy (150-152, right).

[1] Ke et al. Learning to Induce Causal Structure. ICLR 2023.

[2] Lorch et al. Amortized inference for causal structure learning. NeurIPS 2022.

[3] Wu et al. Sample, estimate, aggregate: A recipe for causal discovery foundation models. TMLR 2025.

lines 193-198

Thank you for pointing this out, and please let us know if this makes more sense.

FCI is run $T$ times on different subsets of $k$ nodes. The "space of possible edge types" (196) refers to the edge types output by FCI, which are mapped to categorical labels (e.g. no edge to 0, -> to 1, – to 2). These elements are "embedded" similar to word embeddings in language models (201-202).

So $E'$ contains $k\times k \times T$ elements ($T$ adjacency matrices, each of size $k\times k$).

FCI noisy estimates and correlation

Our causal graph learner, finetuned from SEA [3], uses correlation and FCI estimates and as input features. Correlation does not encode CI information, but it may help the model quickly filter out edges that are not likely to exist, especially important for larger datasets.

FCI is run on subsets of $k=5$ throughout (above), rather than the full dataset [194]. These smaller subsets are sampled based on heuristics like high pairwise correlation, to prioritize CI testing for regions likely to be "interesting."

There are classic algorithms for integrating (sub)graph estimates over different variable sets [4]. SEA trains a neural network to perform this task over datasets of varying assumptions, precisely because FCI estimates may be noisy in practice [3]. In this sense, the role of the neural network is to "denoise" the graph prediction and to learn to resolve contradictions.

[4] Huang et al. Causal discovery from multiple data sets with non-identical variable sets. AISTATS 2020.

Permuting the edge embedding

This practice follows SEA [3] and loosely resembles the idea in [5], where permutations are sampled at random during training so that "in expectation" the Transformer does not rely on any particular ordering.

It's necessary to add embeddings to inform of node/edge identities across layers, since Transformers are permutation invariant. In principle, any ordering will do. However, since synthetic datasets are generated in topological order, "node 1" always corresponds to a root node and so on. As [6] shows, it's trivial for neural networks to memorize this spurious feature to recover topological order.

[5] Yang et al. XLNet: Generalized Autoregressive Pretraining for Language Understanding. NeurIPS 2019.

[6] Reisach et al. Beware of the Simulated DAG! Causal Discovery Benchmarks May Be Easy To Game. NeurIPS 2021.

Thank you for your comments! Please let us know if this response helps clarify your questions, and if you have any further concerns.

"direct application of existing CDN algorithms to the target identification task"

CDN is the name of the method we propose; are you referring to causal discovery algorithms?

Our key contribution is the idea of supervised causal discovery (predicting graphs from simulated data) as a pretraining objective, to obtain dataset representations useful for downstream tasks like inferring targets.

A direct solution to the task is to predict targets without inferring graphs. However, we empirically demonstrate that graph prediction leads to substantial improvement as an auxiliary loss (Table 3, no $L_G$) and starting from pretrained weights leads to more stable optimization dynamics (Figure 10). Furthermore, baselines that use data-mined graph priors perform substantially worse (Table 1, Figure 3).

We realize our idea with [1] as the causal structure learner, and a 2D Transformer-based architecture as the differential network, neither of which are our primary contributions. These modules were selected for their scalability and could (in principle) be replaced with alternate supervised learning algorithms [2,3] or models that operate pairs of graph [4].

[1] Wu et al. Sample, estimate, aggregate: A recipe for causal discovery foundation models. 2024.

[2] Ke et al. Learning to Induce Causal Structure. ICLR 2023.

[3] Lorch et al. Amortized inference for causal structure learning. NeurIPS 2022.

[4] Li et al. Graph Matching Networks for Learning the Similarity of Graph Structured Objects. ICML 2019.

interpretable biological analyses

This work focuses on the application of predicting perturbation targets, for efficient experimental design. The goal is not to replace experiment via simulation, but to lower the cost/scale of experiments, by reducing the search space. We do not believe that this focus diminishes from the interpretability of the results.

• Since the causal graph learner is supervised on synthetic graphs, the predicted graphs are always available for inspection. This is not a primary focus of this work, due to the lack of cell line-specific "ground truth" graphs for evaluation. However, the predicted edges represent gene-gene interactions, and the differences between the observational and interventional graphs correspond to differences in gene regulation.

• Works that do predict phenotypes often share our intended application [5], and as we write in the introduction, directly predicting the target is more efficient at inference time, especially for combinations.

• Directly predicting targets also has benefits in terms of label / evaluation quality. Perturbation target identity is one of the most reliable readouts from these single-cell data, as it is non-quantitative [6]. In contrast, due to technical and biological covariates, two datasets generated by the same lab, of the same cell line, can be largely inconsistent in terms of individual gene expression. For example, [7] writes of K562 genome-wide and essential: "The median correlation between log-fold-change point estimates [...] was only 0.16, suggesting very low replicability."

[5] Roohani et al. GEARS: Predicting transcriptional outcomes of novel multi-gene perturbations. Nat Biotech 2024.

[6] Luecken and Theis. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Mol Syst Biol (2019) 15: e8746.

[7] Nadig et al. Transcriptome-wide characterization of genetic perturbations. 2024.

Thank you for your comments! Please let us know if this response helps clarify your questions, and if you have any further concerns.

Architecture similarities and differences / novelty

Our key contribution is the idea of supervised causal discovery (predicting graphs from simulated data) as a pretraining objective, to obtain dataset representations useful for downstream tasks like inferring targets. We realize our idea with [1] as the causal structure learner, and a 2D Transformer-based architecture as the differential network, neither of which are our primary contributions. These modules were selected for their scalability and could (in principle) be replaced with alternate supervised learning algorithms [2,3] or models that operate pairs of graph [4].

For a specific comparison to [1] –

Similarities: We take the aggregator in SEA as a pretrained module, which is finetuned.

Differences:

• We finetune SEA using correlation as a summary statistic. The original paper uses the precision matrix instead (not suitable for single-cell data), but writes that it's possible to finetune weights for different statistics.
• The differential network is separate and learned from scratch.
• We generate our own datasets for training throughout.

[1] Wu et al. Sample, estimate, aggregate: A recipe for causal discovery foundation models. TMLR 2025.

[2] Ke et al. Learning to Induce Causal Structure. ICLR 2023.

[3] Lorch et al. Amortized inference for causal structure learning. NeurIPS 2022.

[4] Li et al. Graph Matching Networks for Learning the Similarity of Graph Structured Objects. ICML 2019.

Data noise and relevance of synthetic data

We agree that real data are noisy, hence the focus on summary statistics (correlation, mean, variance) as the primary inputs to the model.

Simulating mRNA expression is a task on its own, and improvements in that regard are directly applicable here. In this work, the synthetic soft interventions were loosely inspired by perturbations, e.g. in terms of fold-changes (line 796).

Summary statistics

We'll make the inputs more clear. Thank you for pointing this out.

Yes, the causal structure learner uses pairwise correlations as the summary statistic. In addition to the paired graph features, the differential network takes marginal means and variances (218-219 right, inspired by differential expression analysis).

Rank criteria

We'll be more explicit regarding the definition.

For perturbation effect predictors (GEARS, GenePT), we order genes based on the similarity of their predicted effect to the ground truth target effect (Eq 3, B.2). Methods that directly predict targets order genes based on their predicted probability of being a target. The set of candidate genes is the set of top 1000 differentially-expressed genes.

The rank (309-311) is the index of the correct answer in any ordering, divided by the number of candidates (since different perturbations have different candidate pools, B.2), where 1 is best (top of light) and 0 is worst (bottom of list).

This is equivalent to the expected recall at p (Table 3), over all p.

Rank is inspired by the application of experimental design. Since it is common to select top predictions for further testing, "rank" quantifies how many perturbations must be tested on average, before finding the ground truth.

Thank you for your comments! Please let us know if these answers are helpful, and if you have any further concerns.

Top 1000 DEGs

It's quite common to restrict analysis to subsets of genes, though the threshold and criteria may vary. Here, we believe that selecting the top 1000 genes (based on differential expression p-value) is reasonably permissive, since these genes are selected per perturbation.

For example, [1] writes: "In a genome-scale Perturb-seq screen, we find that a typical gene perturbation affects an estimated 45 genes, whereas a typical essential gene perturbation affects over 500 genes."

In addition, [2] writes regarding highly-variable gene selection (global selection across all perturbations, based on variance to mean ratio): "Depending on the task and the complexity of the dataset, typically between 1,000 and 5,000 HVGs are selected for downstream analysis [...]. Preliminary results from Klein et al (2015) suggest that downstream analysis is robust to the exact choice of the number of HVGs. While varying the number of HVGs between 200 and 2,400, the authors reported similar low‐dimensional representations in the PCA space."

Finally, a classic approach for identifying drug targets [3] argues that targets are not highly differentially expressed, but network-based approaches are useful for predicting them from gene expression. There, they also sample 1000 candidate targets as decoys.

[1] Nadig et al. Transcriptome-wide characterization of genetic perturbations. 2024.

[2] Luecken and Theis. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Mol Syst Biol (2019) 15: e8746.

[3] Isik et al. Drug target prioritization by perturbed gene expression and network information. Scientific Reports 2015.

Larger synthetic experiments

We are happy to include experiments on larger graphs ($N=100, E=200$, AUC over 10 new datasets). There is some drop in performance, especially Polynomial. This could be attributed to error propagation from SEA, which performs less well (at predicting graphs) on larger graphs, and could potentially be ameliorated by finetuning on larger synthetic data (similar to how Perturb-seq models were finetuned to be larger).

Linear + Soft remains quite good, perhaps because the theory (Appendix C) suggests that the summary statistics are sufficient for predicting targets, without the graph prediction.

Different test splits

Typically, the largest cluster (assigned to train, e.g. nearly 1/3 of K562) contains weak perturbations with minimal phenotype, so different clusterings don't introduce as much heterogeneity.

Instead, we can partition all perturbations uniformly into 5 sets for the unseen cell line setting. As seen here for K562 genome-wide, the variation across sets is minimal. We will include similar results for all cell lines.

Thank you for your comments! Please let us know if this response helps clarify your questions, and if you have any further concerns.

Intended application

We will make this more explicit in the paper. This work focuses on the application of predicting perturbation targets, for efficient experimental design. The goal is not to replace experiment via simulation, but to lower the cost/scale of experiments, by reducing the search space.

We envision for actual use cases to resemble the chemical perturbation case more, e.g.

• Inputs: before/after drug (unknown mechanism)
• Model generates hypotheses for molecular mechanisms of drugs
• Top predictions should be verified experimentally, i.e. knock down proposed genes and compare to "after drug" phenotype.

This type of method has an advantage to phenotype prediction methods when considering combinations of targets, e.g. messy drugs. Different off-targets may contribute orthogonal effects, which may make the predicted phenotype of any single target less similar to the overall profile. In addition, there are limited real training data for combinatorial perturbations, motivating methods that can leverage synthetic data.

Why Perturb-seq?

Some experiments use Perturb-seq datasets because in terms of model evaluation, the ground truth targets are relatively "clean." CRISPR is very specific compared to drugs, which may impact many proteins, known and unknown. While this is not the most "realistic" application, we chose these datasets for their size and ease of evaluation.

Single cell foundation models

We will include discussion of these works. Thank you for the suggestion!

Generalization under distribution shift

Our "unseen cell line" case resembles "zero-shot" generalization, since experiments in different cell lines are necessarily conducted separately (282). While these datasets were produced by the same lab, they are actually quite incomparable at the level of individual genes, even for the same cell line. For example, [1] writes of K562 genome-wide and essential: "The median correlation between log-fold-change point estimates [...] was only 0.16"

[1] Nadig et al. Transcriptome-wide characterization of genetic perturbations. 2024.

In controlled settings, our synthetic experiments hold out the "scale" soft intervention from training (368 right, Table 9). We include the exact formula for these interventions / causal mechanisms in D.2 (p16-17).

Interpretability of intermediate graphs

Since the causal graph learner is supervised on synthetic graphs, the predicted graphs are always available for inspection. This is not a primary focus of this work, due to the lack of cell line-specific "ground truth" graphs for evaluation. However, the predicted edges could represent gene regulatory relationships, and the differences between the observational and interventional graphs correspond to differences in gene regulation.

Theoretical analysis

We show in Appendix C that the differential network is well-specified as a model class. Due to the complexity of the architecture, it is difficult to guarantee what the model will learn, i.e. in the classical identifiability sense. At the same time, it is not a given that an arbitrary Transformer can approximate any arbitrary function given limited precision, width, depth, [2,3]. Instead, we are inspired by the neural network expressivity literature to show that (page 15):

• Hard interventions: given the sets of edges $E, E'$, where $E'$ represents the mutilated graph, there exists a setting of its parameters that maps $(E, E')$ to targets $I$ by detecting edges that disappeared.
• Soft interventions: under "scale" interventions and non-multiplicative noise, there exists a setting of its parameters that maps the correlation and covariance matrices (which can be derived from model inputs) to targets $I$. In the linear case, Eq. 14-15 naturally hold and are relevant to the relationships between variables (i.e. if $x$ is an intervention target, for every $y\ne x$, correlation($x,y$) does not change while cov($x,y$) does.

[2] Yun et al. Are Transformers universal approximators of sequence-to-sequence functions? ICLR 2020

[3] Sanford et al. Representational Strengths and Limitations of Transformers. NeurIPS 2023.

Typos

Thank you for pointing this out! We will fix it.

Thank you for your comments! Please let us know if this response helps clarify your questions, and if you have any further concerns.

Technical novelty

Our key contribution is the idea of supervised causal discovery (predicting graphs from simulated data) as a pretraining objective, to obtain dataset -> graph representations useful for downstream tasks like inferring targets. A direct solution to the task is to predict targets without inferring graphs. However, we empirically demonstrate that graph prediction leads to substantial improvement as an auxiliary loss (Table 3, no $L_G$) and starting from pretrained weights leads to more stable optimization dynamics (Figure 10). Furthermore, baselines that use data-mined graph priors perform substantially worse (Table 1, Figure 3).

We realize our idea with [1] as the causal structure learner, and a 2D Transformer-based architecture as the differential network, neither of which are our primary contributions. These modules were selected for their scalability and could (in principle) be replaced with alternate supervised learning algorithms [2,3] or models that operate pairs of graphs [4].

[1] Wu et al. Sample, estimate, aggregate: A recipe for causal discovery foundation models. 2024.

[2] Ke et al. Learning to Induce Causal Structure. ICLR 2023.

[3] Lorch et al. Amortized inference for causal structure learning. NeurIPS 2022.

[4] Li et al. Graph Matching Networks for Learning the Similarity of Graph Structured Objects. ICML 2019.

Identifiability

Our setting is slightly different from the typical causal representation learning setting, where high-dimensional observations (e.g. images) are associated with (unobserved) low-dimensional factors, on which interventions are applied.

We follow the classic causality framework (Sec 2.1), where we measure all variables (causal sufficiency), and interventions are applied directly to these variables. Many classical causal discovery algorithms apply to the observational setting, in which certain graph topologies are still identifiable, under various assumptions [5].

Supervised causal discovery algorithms do not have explicit assumptions, but they do follow implicit ones captured by the training data [6].

[5] Sprites, Glymour, and Scheines. Causation, prediction, and search. 2001.

[6] Montagna et al. Demystifying amortized causal discovery with transformers. 2024.

Value of synthetic data.

Empirically, the baselines that predict targets from real data without external information / synthetic data do not perform as well as our full method (Table 1). The synthetic soft interventions were loosely inspired by perturbations, e.g. in terms of fold-changes (796). In general, synthetic data pretraining also improves training stability (Figure 10).

Bidirected edges

The model's input and output spaces are different. The causal graph learner outputs only three edge types: "A -> B," "B -> A," "no edge." FCI graph estimates are inputs to the model, along with correlation. The input edge type embedding includes bidirected edges because FCI may estimate them.

Prior work

Thank you for the recommendation! We will include this in the discussion of related work.

L75-80 right

These lines describe how a perturbation effect "simulator" could be used to identify targets. The actual implementation for GEARS and GenePT is described in B.2. We will update this section to be more descriptive and add a reference.

• "space of perturbations $\mathcal{I}$" refers to the set of sets of perturbation targets $I$ that are possible, e.g. the powerset of all variables, if we allow interventions of arbitrary # of targets
• "an oracle $\phi : D_\text{obs}, I \to P$ refers to the simulator, which takes as input, an observational dataset (starting cells) and a set of intervention targets (genes)
• "the $I^*$ associated with $\tilde{P*}$ refers to the set of ground truth targets, associated with the target distribution (Eqs 1 and 2)
• "dissimilarity d" could be cosine distance, or mean squared error, or KL divergence between distributions
• $\arg\min$ over $\mathcal{I}$ implies an optimization over all possible sets of targets, e.g. naively we run the model $\phi$ once per $I$ and output the $I$ whose prediction is closest to the target cells.

What are k and T?

We will update the paper to clarify. FCI is run $T$ times on different subsets of $k$ nodes to produce subgraph estimates (input to the causal graph learner)

Differential network

The differential network architecture is implemented using standard multi-head attention, applied along each axis. We will expand the descriptions in the main text and include references to the appropriate details.

Conclusion is missing

We will add a conclusion to highlight the key contributions, limitations, and future directions.

Thank you for your questions and recommendations! We hope this response addresses your questions, but please let us know if you have any additional concerns.

Additional historical context

Thank you for these references and suggestions. We will incorporate them into the paper alongside additional discussion of the historical context.

Extension to dynamical systems

An empirical solution is to simulate cyclic data, and to use a classic discovery algorithm for cyclic graphs.

We can also model the dynamical system. In perturbation data, cells are sequenced after perturbations penetrate completely, e.g. 6-8 days [1], so we might assume that the data are sampled from the stationary density of a system of stochastic differential equations [2]. Here, the causal structure learner could predict the SDE parameters instead of the causal graph. Similar to the acyclic case, interventions could be "read off" the differences in the SDE parameters, e.g. $f$ and $\sigma$ in [2]. The key design choice/challenge would be to select an input featurization that allows for the targets to be identified.

For actively differentiating cells (instead of pre/post differentiation), steady state assumptions might not hold. While there exist measurements like RNA velocity (the ratio of spliced to unspliced mRNA) that can partially order cells within a population and inform of causal structure [3, 4], modeling these continuous dynamics is out of the scope of this work.

[1]. Replogle et al. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq. Cell. 2022.

[2] Lorch et al. Causal Modeling with Stationary Diffusions. AISTATS 2024.

[3] Singh et al. Causal gene regulatory analysis with RNA velocity reveals an interplay between slow and fast transcription factors. Cell Systems 2024.

[4] Atanackovic et al. DynGFN: Towards Bayesian Inference of Gene Regulatory Networks with GFlowNets. NeurIPS 2023.

Sensitivity to simulation

• CDN is not sensitive to the type of non-linearity among graphs of the same size, as performance remains about the same for novel non-linearities (Polynomial, Table 3 and Sigmoid, Table 10).
• CDN is more sensitive to the type of intervention (Table 9), so this is a key parameter to understand / simulate.
• Mismatch issues: During development, we found that finetuning on Perturb-seq was not helpful for the Sci-Plex data, presumably because CRISPR perturbations are much cleaner and stronger than drugs. We did design the synthetic datasets with less specific / soft interventions in mind [below], which may render them a better fit.

Sim-to-real gap: If a model has been trained and one would like to simulate additional finetuning data, one could be informed by the data-driven graph connectivity (e.g. based on pairwise correlation) and classic hypothesis tests (e.g. for data distribution, differences between pre/post intervention).

However, these heuristics cannot fully capture complex systems like biology. Instead, perhaps one could use a self-supervised objective to finetune the model on observational data, a la test-time training approaches [5]. Alternatively, one could incorporate biological knowledge graphs with the understanding that they are noisy priors. E.g. as an additional input feature or as a finetuning objective for the causal graph learner.

[5] Sun et al. Test-Time Training with Self-Supervision for Generalization under Distribution Shifts. ICML 2020.

Interpretability of intermediate graphs

Since the causal graph learner is supervised on synthetic graphs, the predicted graphs are always available for inspection. This is not a primary focus of this work, due to the lack of cell line-specific "ground truth" graphs for evaluation. However, the predicted edges could represent gene regulatory relationships, and the differences between the observational and interventional graphs correspond to differences in gene regulation.

Multi-target identification

Currently, CDN is trained to make an independent prediction per gene (sigmoid, not softmax, 242-244). The synthetic experiments contain up to 3 targets per graph, and the model does output multiple high scores (Table 3).

A limitation of this current approach is that the per-gene decision does not consider other genes in the proposed subset. To (efficiently) consider the entire subset, the prediction could be made autoregressively, where the order of predictions could be based on effect size, or be sampled at random during training.