Discovering Root Causal Genes with High Throughput Perturbations

Reviewer #1 (Public review):

Summary:

This manuscript seeks to estimate the causal effect of genes on disease. To do so, they introduce a novel algorithm, termed the Root Causal Strength using Perturbations (RCSP) algorithm. RCSP uses perturb-seq to first estimate the gene regulatory network structure among genes, and then uses bulk RNA-seq with phenotype data on the samples to estimate causal effects of genes on the phenotype conditional on the learned network structure. The authors assess the performance of RCSP in comparison to other methods via simulation. Next, they apply RCSP to two real human datasets: 513 individuals age-related macular degeneration and 137 individuals with multiple sclerosis.

Strengths:

The authors tackle an important and ambitious problem - the identification of causal contributors to disease in the context of a causal inference framework. As the authors point out, observational RNA-seq data is insufficient for this kind of causal discovery, since it is very challenging to recover the true underlying graph from observational data; interventional data are needed. However, little perturb-seq data has been generated with annotated phenotype data, and much bulk RNA-seq data has already been generated, so it is useful to propose an algorithm to integrate the two as the authors have done.

The authors also offer substantial theoretical exposition for their work, bringing to bear both the literature on causal discovery as well as literature on the genetic architecture of complex traits.

Weaknesses:

The notion of a "root" causal gene - which the authors define based on a graph theoretic notion of topologically sorting graphs - requires a graph that is directed and acyclic. It is the latter that constitutes an important weakness here - it simply is a large simplification of human biology to draw out a DAG including hundreds of genes and a phenotype Y and to claim that the true graph contains no cycles. This is briefly touched upon the discussion, but given the fundamental nature of this choice - the manuscript should devote at least some of the main results to exploring the consequence of mischaracterizing true cyclic graphs as DAGs in this framework. For example - consider the authors' analysis of T cell infiltration in multiple sclerosis (MS). CD4+ effector T cells have the interesting property that they are stimulated by IL2 as a growth factor; yet IL2 also stimulates the activation of (suppressive) regulatory T cells. What does it mean to analyze CD4+ regulation in disease with a graph that does not consider IL2 (or other cytokine) mediated feedback loops/cycles?

I also encourage the authors to consider more carefully when graph structure learned from perturb-seq can be ported over to bulk RNA-seq. Consider again the MS CD4+ example - the authors first start with a large perturb-seq experiment (Replogle et al., 2022) performed in K562 cells. To what extent are K562 cells, which are derived from a leukemia cell line, suitable for learning the regulatory structure of CD4+ cells from individuals with an MS diagnosis? Presumably this structure is not exactly correct - to what extent is the RCSP algorithm sensitive to false edges in this graph? This leap - from cell line to primary human cells - is also not modeled in the simulation. Although challenging - it would be ideal for the RCSP to model or reflect the challenges in correctly identifying the regulatory structure.

It should also be noted that in most perturb-seq experiments, the entire genome is not perturbed, and frequently important TFs (that presumably are very far "upstream" and thus candidate "root" causal genes) are not expressed highly enough to be detected with scRNA-seq. In that context - perhaps slightly modifying the language regarding RCSP's capabilities might be helpful for the manuscript - perhaps it would be better to describe it has an algorithm for causal discovery among a set of genes that were perturbed and measured, rather than a truly complete search for causal factors. Perhaps more broadly - it would also benefit the manuscript to devote slightly more text to describing the kinds of scenarios where RCSP (and similar ideas) would be most appropriately applied - perhaps a well-powered, phenotype annotated perturb-seq dataset performed in a disease relevant primary cell.

Reviewer #2 (Public review):

Summary:

This paper presents a very interesting use of a causal graph framework to identify the "root genes" of a disease phenotype. Root genes are the genes that cause a cascade of events that ultimately leads to the disease phenotype, assuming the disease progression is linear.

Strengths:

- The methodology has a solid theoretical background.
- This is a novel use of the causal graph framework to infer root causes in a graph

Weaknesses:

(1) General Comments
First, I have some general comments. I would argue that the main premise of the study might be inaccurate or incomplete. There are three major attributes of real biological systems, which are not considered in this work.

One is that the process from health-to-disease is not linear most of the time with many checks along the way that aim to prevent the disease phenotype. This leads to a non-deterministic nature of the path from health-to-disease. In other words, with the same root gene perturbations, and depending on other factors outside of gene expression, someone may develop a phenotype in a year, another in 10 years and someone else never. Claiming that this information is included in the error terms might not be sufficient to address this issue. The authors should discuss this limitation.

Two, the paper assumes that the network connectivity will remain the same after perturbation. This is not always true due to backup mechanisms in the cells. For example, suppose that a cell wants to create product P and it can do it through two alternative paths:
Path #1: A -> B -> P Path #2: A -> C -> P
Now suppose that path #1 is more efficient, so when B can be produced, path #2 is inactive. Once the perturbation blocks element B from being produced, the graph connectivity changes by activation of path #2. I did not see the authors taking this into consideration, which seems to be a major limitation in using perturb-seq results to infer connectivities.

Three, there is substantial system heterogeneity that may cause the same phenotype. This goes beyond the authors claim that although the initial gene causes of a disease may differ from person to person, at some point they will all converge to changes in the same set of "root genes". This is not true for many diseases, which are defined based on symptoms and lab tests at the patient level. You may have two completely different molecular pathologies that lead to the development of the same symptoms and test results. Breast cancer with its subtypes is a prime example of that. In theory, this issue could be addressed if there is infinite sample size. However, this assumption is largely violated in all existing biological datasets.

All the above limit the usefulness of this method for most chronic diseases, although it might still lead to interesting discoveries in cancer (in which the association between genes' dysregulation and development of cancer is more direct and occurs in less amount of time).

With these in mind, the theoretical and algorithmic advances this paper offers are interesting. And the theoretical proofs are solid.

(2) Specific comments.
I am curious how the simulated data were generated and processed. Specifically, were the values of the synthetic variables Z-scored? If not, then I would expect that the variance of every variable will increase from the roots of the graph to the leaves. That will give an advantage in any algorithm aiming to identify causal relations based on error terms. For fairness and completeness, the authors should Z-score the values in the synthetic data and compare the results.

The algorithm seems to require both RNA-seq and Perturb-seq data (Algorithm 1, page 14). Can it function with RNA-seq data only? What will be different in this case?

(3) Additional comments:
Although the manuscript is generally written clearly, some parts are not clear and others have missing details that make the narrative difficult to follow up. Some specific examples:
- Synthetic data generation: how many different graphs (SEMs) did they start from? (30?) How many samples per graph? Did they test different sample sizes?
- The presentation of comparative results (Suppl fig 4 and 7) is not clear. No details are given on how these results were generated. (what does it mean "The first column denotes the standard deviation of the outputs for each algorithm"?) Why all other methods have higher SD differences than RCSP? Is it a matter of scaling? Shouldn't they have at least some values near zero since the authors "added the minimum value so that all histograms begin at zero"? also, why RCSP results are more like a negative binomial distribution and every other is kind of normal?
- What is the significance of genes changing expression "from left to right" in a UMAP plot? (eg Fig. 3h and 3g)

The authors somewhat overstate the novelty of their algorithm. Representation of GRNs as causal graphs dates back in 2000 with the work of Nir Friedman in yeast. Other methods were developed more recently that look on regulatory network changes at the single sample level which the authors do not seem to be aware (e.g., Ellington et al, NeurIPS 2023 workshop GenBio and Bushur et al, 2019, Bioinformatics are two such examples). The methods they mention are for single cell data and they are not designed to connect single sample-level changes to a person's phenotype. The RCS method needs to be put in the right background context in order to bring up what is really novel about it.

Reviewer #3 (Public review):

Summary:

The authors provide an interesting and novel approach, RCSP, to determining what they call the "root causal genes" for a disease, i.e. the most upstream, initial causes of disease. RCSP leverages perturbation (e.g. Perturb-seq) and observational RNA-seq data, the latter from patients. They show using both theory and simulations that if their assumptions hold then the method performs remarkably well, compared to both simple and available state-of-the-art baselines. Whether the required assumptions hold for real diseases is questionable. They show superficially reasonable results on AMD and MS.

Strengths:

The idea of integrating perturbation and observational RNA-seq dataset to better understand the causal basis of disease is powerful and timely. We are just beginning to see genome-wide perturbation assay, albeit in limited cell-types currently. For many diseases, research cohorts have at least bulk observational RNA-seq from a/the disease-relevant tissue(s). Given this, RCSP's strategy of learning the required causal structure from perturbations and applying this knowledge in the observational context is pragmatic and will likely become widely applicable as Perturb-seq data in more cell-types/contexts becomes available.

The causal inference reasoning is another strength. A more obvious approach would be to attempt to learn the causal network structure from the perturbation data and leverage this in the observational data. However, structure learning in high-dimensions is notoriously difficult, despite recent innovations such as differentiable approaches. The authors notice that to estimate the root causal effect for a gene X, one only needs access to a (superset of) the causal ancestors of X: much easier relationships to detect than the full network.

The applications are also reasonably well chosen, being some of the few cases where genome-scale perturb-seq is available in a roughly appropriate (see below) cell-type, and observational RNA-seq is available at a reasonable sample size.

Weaknesses:

Several assumptions of the method are problematic. The most concerning is that the observational expression changes are all causally upstream of disease. There is work using Mendelian randomization (MR) showing that the _opposite_ is more likely to be true: most differential expression in disease cohorts is a consequence rather than a cause of disease (https://www.nature.com/articles/s41467-021-25805-y). Indeed, the oxidative stress of AMD has known cellular responses including the upregulation of p53. The authors need to think carefully about how this impacts their framework. Can the theory say anything in this light? Simulations could also be designed to address robustness.

A closely related issue is the DAG assumption of no cycles. This assumption is brought to bear because it required for much classical causal machinery, but is unrealistic in biology where feedback is pervasive. How robust is RCSP to (mild) violations of this assumption? Simulations would be a straightforward way to address this.

The authors spend considerable effort arguing that technical sampling noise in X can effectively be ignored (at least in bulk). While the mathematical arguments here are reasonable, they miss the bigger picture point that the measured gene expression X can only ever be a noisy/biased proxy for the expression changes that caused disease: 1) Those events happened before the disease manifested, possibly early in development for some conditions like neurodevelopmental disorders. 2) bulk RNA-seq gives only an average across cell-types, whereas specific cell-types are likely "causal". 3) only a small sample, at a single time point, is typically available. Expression in other parts of the tissue and at different times will be variable.

My remaining concerns are more minor.

While there are connections to the omnigenic model, the latter is somewhat misrepresented. 1) The authors refer to the "core genes" of the omnigenic model as being at the end (longitudinally) of pathogenesis. The omnigenic model makes no statements about temporally ordering: in causal inference terminology the core genes are simply the direct cause of disease. 2) "Complex diseases often have an overwhelming number of causes, but the root causal genes may only represent a small subset implicating a more omnigenic than polygenic model" A key observation underlying the omnigenic model is that genetic heritability is spread throughout the genome (and somewhat concentrated near genes expressed in disease relevant cell types). This implies that (almost) all expressed genes, or their associated (e)SNPs, are "root causes".

The claim that root causal genes would be good therapeutic targets feels unfounded. If these are highly variable across individuals then the choice of treatment becomes challenging. By contrast the causal effects may converge on core genes before impacting disease, so that intervening on the core genes might be preferable. The jury is still out on these questions, so the claim should at least be made hypothetical.

The closest thing to a gold standard I believe we have for "root causal genes" is integration of molecular QTLs and GWAS, specifically coloc/MR. Here the "E" of RCSP are explicitly represented as SNPs. I don't know if there is good data for AMD but there certainly is for MS. The authors should assess the overlap with their results. Another orthogonal avenue would be to check whether the root causal genes change early in disease progression.

The available perturb-seq datasets have limitations beyond on the control of the authors. 1) The set of genes that are perturbed. The authors address this by simply sub-setting their analysis to the intersection of genes represented in the perturbation and observational data. However, this may mean that a true ancestor of X is not modeled/perturbed, limiting the formal claims that can be made. Additionally, some proportion of genes that are nominally perturbed show little to no actual perturbation effect (for example, due to poor guide RNA choice) which will also lead to missing ancestors.

The authors provide no mechanism for statistical inference/significance for their results at either the individual or aggregated level. While I am a proponent of using effect sizes more than p-values, there is still value in understanding how much signal is present relative to a reasonable null.

I agree with the authors that age coming out of a "root cause" is potentially encouraging. However, it is also quite different in nature to expression, including being "measured" exactly. Will RCSP be biased towards variables that have lower measurement error?

Finally, it's a stretch to call K562 cells "lymphoblasts". They are more myeloid than lymphoid.