Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Confounder-robust causal discovery and inference in Perturb-seq using proxy and instrumental variables

Kwangmoon Park, Hongzhe Li

TL;DR

This work introduces ARGEN, a confounder-robust framework for learning causal gene networks from Perturb-seq data by leveraging proxy-variable adjustments and instrumental-variable logic. By embedding a perturbation-aware SEM within a DAG and developing a two-stage descent-and-ascend estimation with online FDR, ARGEN identifiably recovers ancestors, descendants, and the full DAG even when unobserved confounders are present. Extensive simulations show superior robustness to omitted variables, and application to K562 Perturb-seq data reveals biologically coherent intra- and inter-chromosomal regulatory modules that align with 3D genome structure and epigenomic signals. The approach yields interpretable causal networks among essential genes, providing a principled path from single-cell perturbations to mechanistic regulatory insights with practical implications for understanding gene regulation.

Abstract

Emerging single-cell technologies that integrate CRISPR-based genetic perturbations with single-cell RNA sequencing, such as Perturb-seq, have substantially advanced our understanding of gene regulation and causal influence of genes. While Perturb-seq data provide valuable causal insights into gene-gene interactions, statistical concerns remain regarding unobserved confounders that may bias inference. These latent factors may arise not only from intrinsic molecular features of regulatory elements encoded in Perturb-seq experiments, but also from unobserved genes arising from cost-constrained experimental designs. Although methods for analyzing large-scale Perturb-seq data are rapidly maturing, approaches that explicitly account for such unobserved confounders in learning the causal gene networks are still lacking. Here, we propose a novel method to recover causal gene networks from Perturb-seq experiments with robustness to arbitrarily omitted confounders. Our framework leverages proxy and instrumental variable strategies to exploit the rich information embedded in perturbations, enabling unbiased estimation of the underlying directed acyclic graph (DAG) of gene expressions. Simulation studies and analyses of CRISPR interference experiments of K562 cells demonstrate that our method outperforms baseline approaches that ignore unmeasured confounding, yielding more accurate and biologically relevant recovery of the true gene causal DAGs.

Confounder-robust causal discovery and inference in Perturb-seq using proxy and instrumental variables

TL;DR

This work introduces ARGEN, a confounder-robust framework for learning causal gene networks from Perturb-seq data by leveraging proxy-variable adjustments and instrumental-variable logic. By embedding a perturbation-aware SEM within a DAG and developing a two-stage descent-and-ascend estimation with online FDR, ARGEN identifiably recovers ancestors, descendants, and the full DAG even when unobserved confounders are present. Extensive simulations show superior robustness to omitted variables, and application to K562 Perturb-seq data reveals biologically coherent intra- and inter-chromosomal regulatory modules that align with 3D genome structure and epigenomic signals. The approach yields interpretable causal networks among essential genes, providing a principled path from single-cell perturbations to mechanistic regulatory insights with practical implications for understanding gene regulation.

Abstract

Emerging single-cell technologies that integrate CRISPR-based genetic perturbations with single-cell RNA sequencing, such as Perturb-seq, have substantially advanced our understanding of gene regulation and causal influence of genes. While Perturb-seq data provide valuable causal insights into gene-gene interactions, statistical concerns remain regarding unobserved confounders that may bias inference. These latent factors may arise not only from intrinsic molecular features of regulatory elements encoded in Perturb-seq experiments, but also from unobserved genes arising from cost-constrained experimental designs. Although methods for analyzing large-scale Perturb-seq data are rapidly maturing, approaches that explicitly account for such unobserved confounders in learning the causal gene networks are still lacking. Here, we propose a novel method to recover causal gene networks from Perturb-seq experiments with robustness to arbitrarily omitted confounders. Our framework leverages proxy and instrumental variable strategies to exploit the rich information embedded in perturbations, enabling unbiased estimation of the underlying directed acyclic graph (DAG) of gene expressions. Simulation studies and analyses of CRISPR interference experiments of K562 cells demonstrate that our method outperforms baseline approaches that ignore unmeasured confounding, yielding more accurate and biologically relevant recovery of the true gene causal DAGs.
Paper Structure (15 sections, 2 theorems, 19 equations, 6 figures, 1 algorithm)

This paper contains 15 sections, 2 theorems, 19 equations, 6 figures, 1 algorithm.

Table of Contents

  1. Introduction
  2. Statistical Model for Perturb-seq Data
  3. Data observed in perturb-seq experiments and DAG
  4. Classical measurement model for scRNA-seq and its extension to Perturb-seq data
  5. Identifiability of the Causal Network
  6. DAG Estimation Based on Perturb-seq Data
  7. Estimation of the descendant sets
  8. Estimation of the parent sets
  9. A DAG search algorithm with interventions.
  10. Simulation
  11. Identification of Causal Essential Gene Networks from K562 Cells with Perturb-seq Experiments
  12. Evaluation of ARGEN intra-chromosomal causal networks from 3D genome and epigenomic datasets
  13. Modular structure of inter-chromosomal causal gene networks
  14. In- and out-degree distributions of the causal networks
  15. Discussion

Key Result

Theorem 1

Under Assumption manuscript:asm:nond and the measurement model with intervention on each node, for all $j\in V(G)$, $des(j)$ and $anc(j)$ are identifiable.

Figures (6)

  • Figure 1: Schematic overview of ARGEN.a. ARGEN takes as input single-cell--level UMI counts $Y=(Y_1,\ldots,Y_p)\in\mathbb{R}^p$, gRNA binding indicators $D=(D_1,\ldots,D_p)\in\{0,1\}^p$, and technical confounders $X=(X_1,\ldots,X_J)\in\mathbb{R}^J$, all of which are observable. Each cell receives a perturbation on a single gene. ARGEN relies on a classical measurement model sarkar2021separating that generates UMI counts $(Y_1,\ldots,Y_p)$ from latent true expression levels $(\mu_1,\ldots,\mu_p)$. Each $\mu_j$ directly receives the perturbation treatment and is influenced by both observed ($X$) and unobserved confounders $(U=(U_1,\ldots,U_m))$. The gRNA binding indicator is also associated with the observed technical factors $X$. b. A simplified DAG with $p=4$, representing the structural equation model of ARGEN. The goal of ARGEN is to learn the (red) edges among the unobserved true expression variables. For simplicity, $X$ is omitted from the visualization. c. Overview of ARGEN's DAG search procedure. Left: ARGEN first estimates the descendant sets of each $\mu_j$ by testing \ref{['manuscript:condind_null']} across $k\neq j$. The estimated $\{\mathrm{des}(j)\mid j\in[p]\}$ are then passed to Algorithm \ref{['manuscript:alg:DAG_search']}. Right: ARGEN subsequently identifies the parent set $pa(j)$ for each $j\in[p]$, proceeding from the bottom layer to the top within Algorithm \ref{['manuscript:alg:DAG_search']}. Parents are estimated via the following steps: (i) construct proxy measurements $\eta_k$ for the true expression $\mu_k$ for all $k\in anc(j)$; (ii) solve a QMLE to estimate $\theta_{jk}$ for each gene $j$; and (iii) adjust the p-values of $\theta_{jk}$, denoted by $p_{\theta_{jk}}$, using Online FDR zrnic2020power across $j\in[p]$ to obtain the final estimate of $pa(j)$.
  • Figure 2: Simulation experiment results for omitted gene analysis.a. Visualization of the true data-generating DAG (True DAG) and the DAGs recovered by each of the three methods when genes $Y_7$ and $Y_8$ (yellow nodes) were omitted from the data. The true DAG is defined on the latent true expressions $\mu_j$, but for simplicity the nodes are labeled using the corresponding observed variables $Y_j$. Edges are colored by the sign of the coefficient values (blue: positive; red: negative). Gray edges denote the edges that are omitted in actual applications. b. Boxplots of the estimated coefficients for six edges $\theta_{jk}$ across 100 replicates. Each panel corresponds to one of the three methods, and the blue and red dashed lines indicate the true values 0.5 and -0.5, respectively. Within each panel, the x-axis denotes the six edges, and the boxplot colors distinguish the Full and Omitted analyses. Red stars above each edge indicate the significance of the Wilcoxon rank-sum test comparing the Full and Omitted analyses for that edge.
  • Figure 3: Directional relationships between genes learned from the intra-chromosomal application of ARGEN.a. Visualization of DAG edges from all parent genes to MCM3. Edge colors correspond to the sign of each coefficient (red: negative; blue: positive). b. Comparison of normalized and adjusted MCM3 expression levels between control cells and cells with perturbations of each of the eight parent genes shown in a. Red stars denote the significance of the Wilcoxon rank-sum test comparing the two groups. c. As in b, but showing the expression levels of each parent gene in control cells and in cells with MCM3 perturbation. d. Comparison between ARGEN edge coefficients (x-axis) and the treatment effect multiplied by $-1$ (y-axis) when each parent gene is intervened. All edges across 23 chromosomes are shown. e. As in d, but with interventions applied to the child genes of each target gene rather than to the parent genes.
  • Figure 4: Biological evaluation of intra-chromosomal DAG edges based on 3D genomics and epigenomics.a. Relationship between genomic distance (x-axis) and edge coefficient $p$-values (y-axis). Each dot represents an edge from the 23 DAGs, and the dahsed red line shows a LOESS fit to the data. The y-axis is displayed with a break between 2.5 and 10, compressing larger values while preserving all observations. b. Box plots comparing the p-values of the edges called as significant by ARGEN. The edges whose nodes are contained in the A compartment (blue) are compared against those contained in the B compartment (orange). c. Proportion of parent genes located within the same Topologically Associating Domain (TAD) as the target gene among the parent genes separated from the target gene by less than 2Mb (y-axis). Each data point represents a target gene with at least one parent gene within 2Mb distance from the TSS. The proportion of ARGEN-significant edges (red) is compared with that of non-significant edges (blue), edges identified by INSPRE (purple), Naive GLM (yellow) and randomly selected edges (gray). d. Enrichment of K562-specific transcription factors and the H3K27ac histone mark in promoter regions (1Kb upstream of the transcription start site) of genes involved in each edge. For each protein, we first computed the proportion of random edges with significant peaks at both nodes (genes). We then report $\log_2(\mathrm{odds\ ratio})$ enrichment values relative to random edges for ARGEN-significant edges (red), non-significant edges (blue), Naive GLM edges (yellow), and INSPRE edges (purple). Each panel corresponds to a protein, and panels with a yellow background highlight proteins for which ARGEN edges show a higher proportion than random edges.ARGEN edges exhibit a higher proportion than random edges.
  • Figure 5: Global structure of the inter-chromosomal causal gene network of K562 essential genes.a. Visualization of the learned DAG. Nodes are colored according to Louvain module clustering, and for each module, the top three Gene Ontology terms are listed using the corresponding colors. b. Histogram of the topological layers of genes in the DAG. c. As in b, but showing the distribution of Katz centrality values.
  • ...and 1 more figures

Theorems & Definitions (2)

  • Theorem 1: Identifiability of descendant and ancestor nodes
  • Theorem 2: Identifiability of the parents