Predicting perturbation targets with causal differential networks

Thank you for your review and suggestions! We hope that this response provides clarity, especially with regards to the application, and we will incorporate your feedback into the writing. Please let us know if you have any other concerns.

Generalization

• On synthetic experiments, all testing datasets were unseen; and the Polynomial mechanism was entirely absent from pretraining/finetuning. We will emphasize this in the paper.

• We constructed our Perturb-seq train/test splits to ensure that perturbations which induce similar effects are placed in the same split (Figure 3 and 310-311). Figure 5 illustrates the correlation between log-fold change induced by each perturbation, and clear clusters emerge (which we allocate to different splits).

• Cell lines are derived from cancer tumors / other conditions (K562 myelogenous leukemia, RPE non-cancerous, HepG2 liver cancer, Jurkat acute T cell leukemia).

  Biologically they are quite distinct: different genes are expressed (the supports of each distribution); among the same genes, the marginal distributions (of expression) and gene-gene relationships may differ. So it would be reasonable to consider them different data distributions. In fact, even separate experiments on the same cell line tend to report different results at the individual gene level [1], due to differences in lab technician, equipment, experimental protocol, etc.

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

Model weights

We have several sets of model weights.

1. On synthetic data, we trained 2 models, corresponding to inverse covariance and correlation global statistics (line 239).
2. On Perturb-seq data, we finetuned the correlation model into 5 separate sets of weights. For unseen cell lines, we finetuned the model on the training sets of every cell line, minus the one in question (306-309). For seen cell lines, we trained an "overall" model on the union of each cell line's training sets.
3. For the chemical perturbations, due to limited data, we used the "overall" model without any adjustments.

Other

"reliance on accurate causal graphs and labeled data":

Since we learn the mapping from graph representations to intervention targets, the model does not require globally accurate causal graphs; only that their differences are informative. This idea has been analyzed in more theoretical detail by [2] for linear Gaussian data.

With regards to labeled data, pretraining on synthetic data actually contributes the most to performance (Table 7), so there is less reliance on labeled real data. In addition, the unseen cell line case + the Sci-Plex experiments demonstrate that it is not necessary to have labeled data from the intended distribution, though performance is higher if there is.

[2] Belyaeva et al. DCI: learning causal differences between gene regulatory networks. Bioinformatics. 2021.

DE refers to "differential expression," i.e. statistically significant change in the amount of a gene (details in B.1). We'll add a brief introduction and reference to the main paper. Thank you for pointing this out!

Thank you for your response and questions! They are very insightful for identifying which parts of our paper are confusing.

Training procedure: Thank you for pointing this out. Please let us know if this makes more sense, and if so, we will update the manuscript to reflect this explanation.

There are two modules in Figure 2, the causal featurizer and differential network. They are trained separately, in multiple stages.

1. Causal featurizer: Each set of weights corresponds to a single global statistic, i.e. one for correlation, one for inverse covariance. These steps are performed by [1], and this component is not updated afterwards.

   A. Causal featurizer is trained from scratch on synthetic data to predict causal graphs (inverse covariance).

   B. Causal featurizer is finetuned on the same synthetic data to predict causal graphs (correlation). This is necessary because the input featurization has changed. It is helpful to finetune vs. train from scratch because i) 0s still correspond to lack of edge, and ii) the "local graph estimate" track can still be used.

2. Differential network: Each set of weights corresponds to a pair of global statistic and training data.

   A. Differential network is trained from scratch on synthetic data to predict targets (inverse covariance), with frozen (inverse covariance) causal featurizer. This model is used for synthetic experiments.

   B. Differential network is trained from scratch on synthetic data to predict targets (correlation), with frozen (correlation) causal featurizer.

   C. Differential network (correlation) is finetuned on real data, 5 times for 5 data splits. These models are used for the real experiments.

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

Table 7 ablates the contribution of model components and step 2B (training differential network on synthetic data), on performance for real data.

The full procedure (on real data) is Step 1B, Step 2B, Step 2C. The ablation is Step 2C.

We observe that if we do not train the causal featurizer / differential network on synthetic data, the performance of our model drops significantly (near random guessing). Our conclusion is that the real data are too small for training the model on their own (the classification task is hard to learn), and synthetic data augmentation is crucial.

Separate models for each statistic: (Inverse) covariance maintains data scale, and [1] reports that inverse covariance performs better on synthetic data, so theoretically, these two would indeed be preferable.

However, correlation is numerically more stable in high dimensions (e.g. number of genes), which goes up to $N=1000$ here. Computing covariance can lead to large values in the input (suboptimal for subsequent neural network optimization), which is why we opted for correlation on larger, biological datasets.

OOD-ness: We believe there are two ways to interpret unseen cell lines (yes and no). With an analogy to language: suppose genes are "words," each cell is a "document," and each cell line is a category of news, e.g. "sports" vs. "popular science."

1. There are domain-specific words (genes), e.g. "football" vs. "DNA," which only appear in their respective categories (cell lines). The style of reporting might also be different, e.g. "sports" focuses on statistics of the latest games, while "popular science" focuses on telling a story.

   Cell lines may express different genes and run different "gene programs." If we compute the Jaccard index (# intersection / # union) of their genes (vocabulary), the highest is 0.766 (K562 essential and Jurkat), while the lowest is 0.702 (Jurkat and HepG2). In this sense, the news categories / cell lines are out-of-distribution from each other.

2. There are also aspects that are similar, on a distributional level. For example, if the goal is to determine whether the news article is positive / negative towards its subject, we can ignore common words like "the" "and" "it", regardless of the genre. Domain-specific words for "goodness" may exist, e.g. "star player" vs. "groundbreaking discovery," but these phrases should have "sharper" distributions, e.g. high tf-idf.

   In cells, we can ignore all genes whose relationships do not change. Across cell lines, the "effect size" of perturbations varies, e.g. how many genes are affected on average. However, the true target may have a "sharper" change, even if "by how much" differs. In this latter sense, there is information that can be shared, albeit softly.

Based on these considerations, our design choice when finetuning the differential network on real data (Step 2C), was not to maintain the mapping of genes across datasets, i.e. the model does not know that gene X in K562 is the same as gene X in Jurkat. This was to enable generalization across cell lines.

what a “target” may not make much sense in the context of biology

We agree that in biology, the concept of a "target" can be hard to define. In this work, we loosely consider targets as molecules which are physically impacted by perturbation, rather than any requirement on conditional independencies.

• Since CRISPRi targets promoters of expressed genes, we treat the intended gene as the "target" (line 291-292).
• In the case of drugs, the "target" is less clear in the context of transcriptomics, since drugs bind physically to protein targets instead of genes. However, prior work [4] has shown that a drug's mechanism of action can be inferred from transcriptomic readouts, presumably due to feedback mechanisms, which is why we consider the protein target's gene counterpart as the "target" (293).

[4] Iorio et al. Discovery of drug mode of action and drug repositioning from transcriptional responses. PNAS 2010.

"I also don’t agree with the statement that 'it is common to assume' that the data was generated by a structural causal model"

We agree that in the mechanistic modeling literature, dynamical systems are more common + correct representations of the relationships between biological variables.

However, we intended to convey that the SCM formalism is a simplifying assumption, often made to reduce the space of plausible explanations. For example, [5] is an early work that demonstrated that causal methodologies can automatically recapitulate known signaling pathways, and [6] is a more recent work that used causal models to narrow down the space of drug repurposing candidates. Thus, following the classic statistics aphorism, we acknowledge that this formulation may be wrong, but as our experiments show, we hope that the resultant model may still be useful.

[5] Sachs et al. Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data. Science. 2005. [6] Belyaeva et al. Causal network models of SARS-CoV-2 expression and aging to identify candidates for drug repurposing. Nat Commun. 2021.

perturbation targets

• Label quality: Due to standard quality control procedures, the identity of each cell's perturbation target is one of the most reliable readouts from Perturb-seq experiments, as it is non-quantitative. When we refer to "noise" in transcriptomic data, we typically refer to the variability in the amount of each gene, rather than sequencing error (modern sequencing machines have error rates of 0.1-0.5% [7], and library design takes this into account).

  More specifically [8]: In Perturb-seq (Figure 3), single cells are annotated with a unique barcode after they are isolated into individual droplets. Thus, each "read," i.e. individual instance of a gene / intended perturbation, is associated with one particular cell. Cells that cannot be confidently mapped to a single perturbation (doublets, unviable cells, empty droplets) are discarded prior to any analysis.

• Strong perturbations: Since the goal of our framework is to achieve a particular phenotypic profile, we focus on strong perturbations to ensure that the desired phenotype is indeed different from the control. In this regard, it appears that "label quality" is a poor word choice, and we will update the paper to reflect this – thank you for pointing this out.

[7] Stoler and Nekrutenko. Sequencing error profiles of Illumina sequencing instruments. NAR Genomics and Bioinformatics, Volume 3, Issue 1, March 2021, lqab019.

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

"forward" vs. "backward" modeling

To the best of our knowledge, there does not currently exist a genome-scale experiment to combinatorially screen genetic perturbations, nor to parallelize across different cell lines (they may have different growth requirements, characteristics, etc.).

Given the experiments that we do have, we introduced the concepts of "forward" and "backward" modeling as two different strategies to extract information from perturbation screens, to inform iterative experimental design or biological understanding. In both cases, the available training data are identical. In the "forward" case, we train a model to map X_obs, Y to X_int. In the "backward" case, we map X_obs, X_int to Y.

Given these trained models, there are two applications that could be of interest.

Thank you for your review and questions. We hope this response clarifies your confusion, and we look forward to discussing with you!

"the paper does not compare to any prior specialist methods that are specifically used to predict intervention targets in biology applications" As we mention in Section 4.1 (338-361), the majority of baselines on our biological experiments are specialist methods and the current state-of-the-art in perturbation modeling:

• PDGrapher is designed for this exact task, i.e. predicting perturbation targets.
• GEARS is an established baseline for unseen perturbation effect prediction, whose goal is to "guide the design of perturbational experiments"
• GenePT is a recent gene/cell "foundation model" that reports state-of-the-art results on a variety of related tasks.

algorithmic novelty

We will improve the writing to highlight the novelty of our architecture.

Ho et al. (2020) is commonly cited for the idea of axial attention, e.g. in MSA Transformer for protein sequence modeling [1], but the architecture we use is quite different.

• The original "Axial Transformer" was designed for images, with considerations for image patch size, spatial locality, color channels, etc.
• In our case, the "axial attention based model" operates over paired graph adjacency matrices, so it is invariant to node labeling order, but does consider the ordering of observational vs. interventional graphs. Axial attention was a natural choice for scaling to large adjacency matrices with $O(N^2)$ entries.
• Pooling over incoming edges was also inspired by the intuition that exogenous interventions sever or alter the relationships between parents and target nodes.

[1] Rao et al. MSA Transformer. ICML 2021.

nothing about the conception or design of the method is specifically motivated by the (biological) problem at hand

From a causality perspective, the task of unknown intervention prediction has been "solved" in many theoretical settings. However, existing algorithms 1) impose strict assumptions regarding the data generation process, 2) are analyzed in the infinite data limit, and 3) scale poorly to >100 variables (lines 141-152).

While not limited to biological applications, certain challenges posed by transcriptomics include: 1) it is unclear what assumptions should hold, if they can even be enumerated; 2) each perturbation has limited data; and 3) measurements are over hundreds or thousands of variables (lines 84-86).

Amortized causal discovery algorithms are designed to address (1) by training over large numbers of datasets, which represent diverse data generating processes (Section 2.3). To address data sparsity and scalability, [2] only requires ~200 data points for high causal discovery performance, and runs on up to 1000 variables in seconds.

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

"no algorithmic component underlying the approach is rooted in a mechanistic model, intervention modeling, or a causal graph"

The causal featurizer is an amortized causal discovery algorithm that infers graph representations from each cell population. While biology largely lacks "ground truth" causal graphs for evaluation, the predicted graphs from [2] have been validated on simulated mRNA data (directed graphs) and on physical protein-protein interaction graphs (undirected graphs), which inspire its application here.

With regards to interventions, many classical causal discovery algorithms are dedicated towards the observational setting, in which certain graph topologies are still identifiable, under many conditions, see [3]. Here, the outputs of the FCI algorithm are input to our model, which was trained to reconcile errors and inconsistencies, as real data are noisy and often misspecified. FCI estimates are asymmetric and represent rich information, e.g. ancestry or confoundedness.

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

"The CDN is a simple classifier that maps Xs to Ys (interaction embeddings to whether or not variables were targeted)"

On a high level, mapping datasets X to graphs G and/or intervention targets Y is the goal of causal discovery. The key lies in how this mapping is learned. Our insight is that modeling putative causal graphs provides useful inductive biases towards this predictive task.

In Table 7, we ablated removing the pretrained amortized causal discovery weights and found that the resultant model (which is a "simple classifier that maps Xs to Ys") performs near random.

Therefore, we concluded that causal structure prediction is a useful pretraining objective for this task.

Thank you for your replies, I remain with my current rating, since the answers do not address my concerns. Your description of the algorithmic novelty of the approach does not provide any additional evidence to me that there is more contribution. I also still remain with my view that nothing about this approach is causal, or mechanistic, so I find it misleading to use the term for a supervised learning method.

Below, I reply to some other comments:

GEARS is an established baseline for unseen perturbation effect prediction,

To my knowledge, this is not true. GEARS is primarily designed to predict perturbation effects, not perturbation targets.

Due to standard quality control procedures, the identity of each cell's perturbation target is one of the most reliable readouts from Perturb-seq experiments, as it is non-quantitative.

My question concerned the drug perturbation data, not the perturb-seq (gene perturbation) data. I am unsure about how confident we can be in the targets there.

We agree that many architectures can approximate arbitrary functions, which is why we wanted to illustrate that our axial attention also has sufficient capacity to perform this task. The AdamW optimizer provably converges [9]

This is confusing in two ways. First, illustrating that axial attention has sufficient capacity is not an informative claim since, as stated above, “many architectures can approximate arbitrary functions”. Second, it is irrelevant whether AdamW converges, because the question is whether it converges to the right solution (this is a non-convex optimization problem)

Thank you for your extensive comments and recommendations! We hope that this and the general response provide clarity on your concerns, and we look forward to discussing with you!

Markov boundary + CI

We implemented a proof of concept of this approach to better illustrate its advantages/drawbacks. We first perform graph lasso with scikit-learn using default parameters. To ensure comparability to existing metrics, instead of discrete CI decisions, we compute the mutual information of MB variables with the domain variable.

Advantages: Overall, the performance is reasonable, and the method is generally fast on small data.

On the synthetic data, we report mAP / AUC with ours and the next best baseline (from Table 4). MB+CI performs similarly to CDN on Polynomial, while CDN is much better on Linear. This baseline is indeed fast on the synthetic data (averaging ~20s / dataset, vs. CDN's 1s for MLP, 39s for Axial).

On the real Sci-Plex data, this simple baseline performs well on 2/6 cases, and underperforms in the rest. A larger evaluation set might reveal more about the failure modes of both.

Disadvantages:

• For theoretical investigation, it makes complete sense to pre-specify the scope of each algorithm or model. In terms of practical application, however, these assumptions are not known in advance.
• Practically speaking, mismatched assumptions lead to difficult / unstable numerical optimization. Graphical lasso converges within seconds on linear datasets, but when given a budget of 1000 iterations, graphical lasso fails to converge on ~15% of our synthetic test sets, and after 10,000 iterations, only ~2% more converge. Numerical stability (NaN values) also became an issue in ~20 synthetic cases.
• This method was very slow on Perturb-seq data, sometimes requiring the same amount of time per perturbation (minutes), as CDN takes on the whole dataset (cell line). For this reason, we chose not to benchmark MB+CI on the full Perturb-seq data during the rebuttal period.
• The two stage procedure can cause error propagation, and the outputs of the first stage are also subject to hyperparameters like the sparsity regularization weight.

Identifiability and training data

Recent work [1] has argued that amortized causal discovery algorithms trade explicit assumptions for implicit ones, encoded by the training data (while still abiding by classical identifiability theory, e.g. linear Gaussian and invertible data are non-identifiable).

• Our real finetuning data (training sets of Perturb-seq data) reflect the assumptions that may be found in similar (biological) data. We constructed our Perturb-seq train/test splits to ensure that perturbations which induce similar effects are placed in the same split (Figure 3 and 310-311). Figure 5 illustrates the correlation between log-fold change induced by each perturbation, and clear clusters emerge (which we allocate to different splits).
• Our synthetic pretraining data contain a variety of linear and non-linear data, as well as both Erdos-Renyi and scale-free graphs. All testing datasets were unseen; and the Polynomial mechanism was entirely absent from pretraining/finetuning. We will emphasize this in the paper.

Therefore, while we do not provide formal guarantees in the classical sense, we expect that the diverse synthetic data + real finetuning data would convey implicit assumptions relevant to the end task.

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

Thank the authors for the response. My concerns are still not addressed: just to predict the intervention target, 1) to learn the Markov blanket is theoretically sufficient, and 2) to use supervised method to learn two graphs is too overkill.

The authors experiments above further support my concern. MB learning can always be tuned more accurate (there are ones for nonlinear relations instead of just graph LASSO or graphical LASSO). In terms of efficiency, one don't need to learn the whole MB; there are variants to learn only the MB of a local variable. In this sense, I have no excuse to believe that separately learning two graphs is more efficient than a local MB learning.

Thank you for your review and comments! We hope that this response clarifies some of your confusion.

Algorithmic contribution: We will improve the writing to highlight the novelty of our architecture. Ho et al. (2020) is commonly cited for the idea of axial attention, e.g. in MSA Transformer for protein sequence modeling [1], but the architecture we use is quite different.

• The original "Axial Transformer" was designed for images, with considerations for image patch size, spatial locality, color channels, etc.
• In our case, the "axial attention based model" operates over paired graph adjacency matrices, so it is invariant to node labeling order, but does consider the ordering of observational vs. interventional graphs. Axial attention was a natural choice for scaling to large adjacency matrices with $O(N^2)$ entries.
• Pooling over incoming edges was also inspired by the intuition that exogenous interventions sever or alter the relationships between parents and target nodes.

Off-target effects: In terms of modeling, we output a separate binary score per gene, so in principle, the model is able to predict multiple targets per condition. In the synthetic data, the number of targets varies from 1-3, and we don't impose any additional inference-time constraints.

In terms of the biological data:

• Many of the Sci-Plex drugs are annotated with off-targets. However, these drugs did not result in differential expression of any annotated on/off target, so we omitted them from the evaluation. These drugs with off-target effects also tended to have smaller effect sizes, median 59 differentially expressed genes vs. 373.
• For the Perturb-seq datasets, CRISPRi tends to be quite specific. Replogle et al 2022. reports 92.5% on target rate, and "neighbor gene knockdown was not enriched in perturbations with a negligible growth defect that produced a transcriptional phenotype," i.e. in our subset.

Table 6 and trivial baseline: On the Sci-Plex data, we can compare to classic differential expression analysis, where we rank targets by their adjusted p-value (Wilcoxon ranked sum test, BH correction). Interestingly, the targets that DE / our method predict better seem to be mutually exclusive, so a larger evaluation set might reveal more about the failure modes of both. This also suggests that these two approaches could be complementary, as one focuses primarily on pointwise differences, while the other focuses more on pairwise changes.

Re: large numbers of "trivial" perturbations – Our subset of Perturb-seq data are somewhat "artificial" in that the intended target is often that with the largest log-fold change, due to selection bias and CRISPRi's specificity (above). Since the goal of our framework is to achieve a particular phenotypic profile, we subsetted to strong perturbations, to ensure that the desired phenotype is indeed different from the control. For example, the raw K562 genome-wide data actually had 9,851 perturbations, but the vast majority had little effect. After filtering, we were only left with 1,767 perturbations (678 test), which by design, tend to report the intended target as the gene with the largest change.

Thank you for your time and insights!

Regarding the Sci-Plex results, we suspect that DE might have performed particularly well on doxorubicin, since its target TOP2A is directly involved in transcription (the only case out of the six). As a result, we might expect to see high transcriptomic response as a result of its inhibition.

On a limited screening budget, DE does miss the target entirely in half the cases, so perhaps it would be preferable to test half of the top predictions from each approach, to provide better coverage.