Gene Regulatory Network Inference in the Presence of Selection Bias and Latent Confounders
Authors: Gongxu Luo
Haoyue Dai
Loka Li
Chengqian Gao
Boyang Sun
Kun Zhang
Venue: NeurIPS 2025 poster ()
Keywords: Selection bias; Causal discovery; latent confounders; Gene regulatory newwork inference
Primary area: machine_learning_for_sciences

Abstract: Gene regulatory network inference (GRNI) aims to discover how genes causally regulate each other from gene expression data. It is well-known that statistical dependencies in observed data do not necessarily imply causation, as spurious dependencies may arise from *latent confounders*, such as non-coding RNAs. Numerous GRNI methods have thus been proposed to address this confounding issue. However, dependencies may also result from *selection*&mdash;only cells satisfying certain survival or inclusion criteria are observed&mdash;while these selection-induced spurious dependencies are frequently overlooked in gene expression data analyses. In this work, we show that such selection is ubiquitous and, when ignored or conflated with true regulations, can lead to flawed causal interpretation and misguided intervention recommendations. To address this challenge, a fundamental question arises: can we distinguish dependencies due to regulation, confounding, and crucially, selection? We show that gene perturbations offer a simple yet effective answer: selection-induced dependencies are *symmetric* under perturbation, while those from regulation or confounding are not. Building on this motivation, we propose GISL (Gene regulatory network Inference in the presence of Selection bias and Latent confounders), a principled algorithm that leverages perturbation data to uncover both true gene regulatory relations and non-regulatory mechanisms of selection and confounding up to the equivalence class. Experiments on synthetic and real-world gene expression data demonstrate the effectiveness of our method.

Gene Regulatory Network Inference in the Presence of Selection Bias and Latent Confounders
Gongxu Luo', Haoyue Dai?, Loka Li', Chengqian Gao', Boyang Sun', Kun Zhang 1.2
Mohamed bin Zayed University of Artificial Intelligence, Carnegie Mellon University (gongxu.luo, kun.zhang/ @mbbuuii.c..e
Abstract
Gene regulatory network inference (GRNI) aims to discover how genes causally regulate each other from gene expression data. is well-known that statistical dependencies in observed data do not necessarily imply causation, spurious dependencies may arise from latent confounders, such as non-coding RNAs. Nu- merous GRNI methods have thus been proposed to address this confounding issue. However, dependencies may also result from selection-one only cells satisfying certain survival inclusion criteria are observed-while these celection-induced spurious dependencies are frequently overlooked in gene expression data analyses. In this work, show that such selection is ubiquitous and, when ignored conflated with true regulations, can lead to flawed causal interpretation and misguided inter- vention recommendations. To address this challenge, fundamental question arises: can we distinguish dependencies due to regulation, confounding, and crucially, selection? We show that gene perturbations offer simple yet effective answer: selection-induee dependencies are symmetric under perturbation, while those from regulation or confounding are not. Building on this motivation, we propose GISL (Gene regulatory network Inference in the presence of Selection bias and Latent confounders). a principled algorithm that leverages perturbation data to uncover both true gene regulatory relations and non-regulatory mechanisms of selection and confounding up to the equivalence class. Experiments on synthetic and real-world gene expression data demonstrate the effectiveness of our method.
Introduction
Gene regulatory network inference (GRNI 212) is fundamentally problem of causal discovery, that is, identifying causal regulatory relationships from observational and experimental gene expression
data [3 (4). Existing GRNI studies include dependence-beee
methods using correlation [5] 6),
regression 7 8 91. and mutual information [10 I[: dynamic modeling using pseudotemporal
15 16); and perturbation modeling using differential analysis [77 18 19 20 Central to these efforts is fundamental question:
trajectories 12 13
statistical dependence observed in gene expression data be interpreted as a regulatory relationship? is knoll that the answer is not always yes yes-yeesseeaee dependencies do not imply causation, spurious dependencies may arise from latent confounding such as non-coding RNAS or environmental stimuli 21
22
23 Numerous GRNI
methods have thus been proposed to address this confounding issue 24
25
26
27
Yet another important source of spurious dependencies remains underexplored: selection bias-the preferential inclusion of data points based on specific criteria (28). In gene expression data, this arises when only cells meeting certain survival inclusion criteria are observed. As result, two genes may appear statistically associated not because one regulates the other they share common regulator, but because only cells in which both genes satisfy the survival criteria persist and are sequenced. Dynamic proliferation studies support this mechanism: perturbations selectively eliminate cells that fail to meet the survival criteria, yielding the differential proliferation rates before and after We show that such selection is pervasive and, if ignored conflated with
perturbation 29 30
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

genuine regulatory interactions, can severely bias GRNI results. Conversely, explicitly modeling selection bias can reveal non-regulatory dependencies and yield novel biological insights.
fundamental question then arises towards GRNI under selection bias: can we distinguish depen- dencies due regulation, confounding, and selection and if so, how? This is challenging. because despite the very different biological nature between regulatory and selection processes, both of them occur upstream of the data collection process (i.e., gene screening), and may thus leave indistin- guishable statistical patterns in observational data. Fortunately, with gene perturbation experiments becoming increasingly feasible in practice BI], this challenge can be effectively addressed. We show that gene perturbations offer simple yet powerful answer to the question: selection-induced dependencies are synmettic under perturbation, while those from regulation or confounding are not.
To illustrate how perturbations help identify selection for GRNI, we examine a case study in leukemia cells 32. In this dataset, the genes AURKA and YOR/AIPI exhibit strong statistical dependence that cannot be explained by any other genes, yet no known regulatory relationship between them is documented in existing databases [33]. Could this be hidden confounding? Perturbation suggests otherwise: perturbing AURKA produces a shift in the marginal distribution of TORIAIPI expression- an outcome inconsistent with the confounding hypothesis, as perturbing gene should not affect its upstream confounders. Furthermore, perturbing TORIAIPI also leads to notable change in AURKA, contradicting the asymmetrie nature expected from a pure causal relationship. Together, these symmetric dependencies under perturbation point to an alternative explanation: a selection process between the two, which aligns with P53 pathway coupling analyses in cancer cells [34 35 36).
Building on the motivation above, we develop a flexible causal framework that models both observa- tional and perturbation gene expression data, and allows for the presence of both latent confounders and selection bias. We characterize the information provided by perturbations through conditional independence (CI) relations in data, and show that regulatory relations, latent confounders, and selection processes typically exhibit distinct CI patterns. Based on these findings, we propose GISL (Gene regulatory network Inference in the presence of Selection bias and Latent confounders), general nonparametric algorithm that not only identifies regulatory relations from potentially biased data, but also detects the underlying confounding and selection processes themselves, shedding lights on non-regulatory relations that, while often overlooked, also play important roles in cellular systems.
The contribution of this work is threefold. 1. This is the first, to the best of our knowledge, to identify and to address the issue of selection bias in gene expression data and its impact on GRNI. 2. We propose novel algorithm for identifying regulatory relationships, as well latent confounders and selection processes up to the equivalence class, Our algorithm is general, without requiring any parametric or graphical assumptions except for the standard ones. 3. We validate our claims and demonstrate the effectiveness of our proposed GISL on both synthetic and real-world experimental single-cell gene expression data, showing its superiority over canonical causal discovery methods and computational GRNI baselines.
Preliminaries
2.1 Causal formulation of gene regulatory networks and gene perturbations
Gene regulatory networks (GRNs) represent the causal relationships governing gene activities in cells (37] Since regulatory interactions fundamentally correspond to causal relationships, we refer to them such throughout. We represent the whole data generating process by a Directed Acyelic Graph (DAG) G whose vertex set can be partitioned by V = [X, SJ. X' = { correspond to observed variables where each X; represents the expression of an individual gene. L = (L,) f accounts for the latent factors that regulate gene expression like non-coding RNAs and environmental constraints. S = [S;]] J,M are selection variables that capture the underlying selection processes. Each S, is a binary indicator variable for an independent selection criterion, with value 1 indicating that criterion being satisfied in a cell, and 0 otherwise [38] 39]. Only cells with all the 5; = are harvested in the dataset. Observed genes influenced by latent confounders and selection variables designated as confounder pairs and selection pairs, respectively, and are formally defined as follows:
(
Definition 2.1 (Confounded pair). A pair (X,,X)) is referred to as confounded pair, denoted (X,, X, J),, if and only if there exist an L3 E L such that the structure X; + LL + X, exists in G.
2

---

x
Z
X
x
Z
mut (G, 0): - Y
mur (V]):X Y Z
T (V.
(0, (X])
(b)
(c)
Figure 1: Alternative graphical representations of interventions. (a) Mutilated DAGs depicting hard
intervention [40]. (b) Generalized intervention representation using the augmented DAG [41]. (c)
Augmented DAG for confounded pairs, where L denotes a latent confounder 442
Definition 2.2 (Selection pair). pair (X;, Xj) is referred as selection pair, denoted Xy)s. if and only if there exist an Sy E S such that the structure X; - S₁ + X, exists in g.
Given gene expression and perturbation data, each intervention (perturbation) target X E X denotes variable X3 is intervened on. Let T {T.) T, S X represent the collection of intervention tar- gets. To model the "action of do perturbation" perturbation indicators I = (T.) are incorporated into the DAG as exogenous variables with directed edges pointing to corresponding intervention tar- gets represented by an augmented DAG 43). . Ls 0 indicates the observational data Do. 4a indi- cates the interventional data D with X being intervened on. Other basic concepts are in Appendix
Established CRISPR-bace gene perturbation methodologies encompass gene knockout (CRISPR- Cas9) and transcriptional modulation (CRISPRa/i), which can be mathematically formalized within causal inference frameworks as hard and soft interventions, respectively. For hard interventions, [40] consider each T3 as factoring a mutilated DAG over [N]. denoted by mul Xx), where the edges incoming to the target X₂ are removed and others remain shown in Figure [[a).
For soft interventions, the mutilated DAG representation does apply, as soft interventions do not remove incoming edges, and all settings may factor in the same G. Thus, the augmented DAG is utilized as generalized framework for representing interventions 433 41, illustrated in Figure where T denotes the intervention target set. Intervening on a cause changes the marginal p(cause) and p(effect), but the conditional p(effect|cauu cause) remains invariant. Conversely, intervening on an effect leaves p(cause) unchanged, p(cause jeffect) changes 44 45
This invariance has been leveraged by numerous interventional causal discovery methods 46 47 predominantly implemented through parametric regression analysis. More related works are discussed in Appendix B
2.2 Understanding selection bias: principles and key characteristics
Selection bias arises when only samples satisfying underlying criteria are systematically included- excluding all others-and thereby induces spurious dependencies (28). The criteria can be categorized in two distinct ways: (1) inherent constraints ('survival'), which arise prior to treatment and are always operative, and (2) sampling bias, which results from non-randomly sampling, both which can be conceptualized within the framework of exogenous selection in causal graphs (28) 44). In
gene expression data, inherent selection constraints predominate, notably competition for shared cellular resources 50, 51 and differential regulatory activity across promoters and enhancers (2].
Example To illustrate how the lection process makes two indepen- dent variables statistically dependent, consider two measurements from a tumor growth study: X (Concentra- tion of inflammatory cytokine IL-6) and (Tumor growth rate). Both variables are independent with no causal relations or confounders (both observed and latent). We simulate this independence in Figure by sampling X and Y independently from a uniform distribution 21(0,2 (sample size = 3,000). However,
L.5
15.
1.0
0.8
0.5
Survired
survised
Figure 2: (a) Scatterplot of X and Y showing selected patients ('e") and excluded individuals ('x"). (b) and (c) Distributions after two distinct interventions variables X and Y, respec- tively, in the selected population (* in (a)).

---

the study protocol restricts analysis to patients with "favorable clinical outcomes" (i.e., those who survived and subsequently presented for hospital care). We model this selection process by retaining only cases where X +
> 2, resulting in a subset of n = 1,485 patients (marked as in (a) Within this selected subset, X and Y appear negatively correlated-crettig the illusion of statistical peeendence-despiin being truly independent in the full population.
To elucidate how selection processes interact with perturbations, we conducted a randomized clinical trial in which patients were assigned to either placebo or treatment arms. The control group follows distribution P which is the 'o' in (a) For the treatment group, we implement soft deterministic interventions on XX specifically do(X 0.5). The resulting distribution P (post- ppostiintervention) visualized by 'e' [b) (n=514). When intervention interacts with selection criteria, we observe that (Y) * PY ) despite the absence of causal relations between X and Y. This discrepancy arises because the consistent application of selection constraints filters out specific samples, altering the distribution of Y. Similarly, under a hard stochastic intervention on Y, do(Y U1-0.., 1],, the distribution Ptv is represented by "o' in (c) (n=386). We observe that ('v(X) * Px) The distribution changes induced by selection after different interventions on both sides, along with the corresponding changes in sample size n, illustrate the symmetry perturbation effect of selection processes and its underlying mechanisms through criterion-based filtering.
With all the notions ready, the joint probability density of intervention over X' is follows:
TTh = p²)(x=
II (1X,TT))
II (X((x(Xaa())) GIXGET})
where pd and pdh denotes the probability density of observational distributions and interventional dis- tributions with intervention target Ti C X, f=(A) indicates the selection constraints on observational sets X, and (ig( indicates the parents of X; in 9 [52]. The joint probability density of observation is TT (X)). T 0. Note that p ("x((XXXXXxxgp (X)))a(())
px(X, X pag pag vX, $ TT.
3 Methodology
Having laid out the necessary preliminaries, now turn interpreting the observed statistical dependencies. We first establish the rationale for leveraging gene perturbation data to detect selection bias, uncover latent confounders, and infer regulatory relations ($ 3.1.. Subsequently, we present computational framework and algorithmic solutions tailored to address these challenges ($) 3.2).
3.1 Differentiating causal relations, selection processes, and latent confounders
We now describe how causal relations, selection processes, and latent confounders can be distin- guished via perturbations, focusing two key patterns: symmetry and perturbation effects.
Perturbation symmetry. Causal relationships (.X + Y) are asymmetrie: perturbing X shifts the distribution of Y, but not vice versa (Figure (a)). ()). Selection processes (X + S + Y) are symmetrie: perturbing either shifts the other, as discussed in Example in Section 2.2 (Figure 3 e)). In the case of latent confounders, unlike selection, perturbations do not propagate through the latent confounder L, i.e., perturbing X or Y does not affect the other (Figure ()). These contrasts provide a basis for differentiating causal, selection, and confounding structures. Note that throughout we focus on acyclic graphs, as our primary goal is to provide a proof of concept for addressing selection bias. Cyclic structures introduce additional complexities not essential for this purpose. Once the acyclic case is fully understood, extending the approach to handle cycles becomes more straightforward.
Perturbation effects. However, perturbation symmetry alone is insufficient when multiple dependen- cies coexist. For instance, both pure regulatory relationships (FigureBaa) and regulation involving latent confounders (Figure d)) are perturbation asymmetric, making them indistinguishable In terms of perturbation effects, fortunately, by modeling observational and perturbation data within the augmented DAG framework (V = {x', L. S, H),, where the perturbation indicators I are incorpo- rated as exogenous variables, we can model how perturbation leads to changes in distribution on other genes via capturing the dependencies between Z and affected ones. To see the power of this new framework, referring back to structures (d) and (a), since conditioning on X constrains the distribution of Y via L, conditional distribution changes between P(Y)X, Ix 1) and P(Y Ix =0) will be

---

Reference Table
A V
A V
A
V IS
Y Ix X,S
V
X X IS
X Hly YY,
(a)
(b)
(d)
(c)
Figure 3: Distinguishing causal, selection, and confounding structures via perturbation effect and indicates the targeted gene pairs of CI test. (a) refers to the direct cause structure between X and Y (represented by 'C'). (b) means there is latent confounder between them (*L"). (c) is the structure of selection process ('S"). (d) stands for causation and latent confounders at the same time ('C L'). (e) stands for causation and selection process the same time ('C S*).
symmetry.
contains CI results. (f) serves as reference table summarizing the CI patterns for each target gene pair: different symbols correspond to different CI relations; black symbols (4, V. indicate the conditional independence, while white symbols (A, V, O. о) indicate the conditional dependence. For example, (a) encodes four CI relations: Y I. Ix (A) and Y 11 Tx X, S the top: X IL Ty S and X V. Ty Y (O) at the bottom.
observed in (d), resulting in Ix Y X. K. In contrast, structure (a) lacks this conditional dependence, making the two structures distinguishable. This complementary use perturbation effects alongside symmetry enables us to distinguish different structures.
Having outlined the distinct perturbation symmetries and perturbation effects associated with the three causal structures, we are now able to categorize the observed dependencies accordingly. To achieve this, we apply conditional independence (CI) tests powerfll statistical tools capable of detecting changes in distributions between the interested variables. Their application relies on two mild assumptions, special case is further discussed in 3.4 as follows:
Assumption 3.1 (Causal Markov assumption). Given DAG G over the variable set every variable M in V is probabilistically independent of its non-descendants given its parents in G.
Assumption 3.2 (Faithfulness assumption). Given a DAG G and distribution P over V, P implies no Cl relations not already entailed by the Markov assumption.
Remark 3.3. These two assumptions common and fundamental for connecting causality with statistical tools, as they ensure that CI tests can correctly capture the underlying causal structure.
Under these mild assumptions, we can now reliably employ the CI test tool to capture distribution changes manifested by perturbation symmetry and effects. The resulting CI patterns between I and X' exhibit distinct behaviors across different structural scenarios, as illustrated in Figure 3 Accordingly, causal relationships, selection processes, and latent confounders are distinguishable. Remark 3.4. When selection interacts with CI testing, between indicator and variables in selection pairs, such as Lx 1 X (V ). and ly 1 X Y (O), as illustrated in the gray-shaded columns (c) and (e) of Figure 3 only when the variable is both intervened and is in the conditional set. This is because testing the dependence between Ix and Y given X is equivalent to detecting differences between the distributions P(Y|X, Ix =) and P(Y|X. Ix = 1). Referring to the pre-intervention case (a) and post-intervention case (c) in Figure 2 it is obvious that P(Y| Ix 0) * P(Y|X, Ix 1), which establishes that Ix V. YIX. Although this special case violates both Markov and faithfulness assumptions, the distinctive characteristic provides opportunity to identify selection. The observed four unique CI patterns are shown Figure i.e., (a), (b), (c), (d).
observe extra conditional dependence
3.2 Algorithm GISL: Handling both selection bias and latent confounding
Building on the different CI patterns across distinct causal structures in Section 3.1 Algorithm I named Gene Regulatory Network Inference in the presence of Selection bias and
we develop
5

---

Algorithm 1: Gene Regulatory Network Inference in the presence of Selection bias and Latent confounders (GISL).
Input:
Observational Do and single gene perturbation ( Xx E T over XIN.- Equivalence class G over X. Confounder pairs L, Sclection pairs S
Output:
Step 1. Initialize: Set 9o as fully undirected graph; L { ). = ]. L' }. Unk = ().
Step Recover the skeleton from observational data Do: g + FGES(g Do).
Step 3. Capture CI patterns from observational data and perturbation data:
foreach (X,,)) g¹ and XX, T do
Compute J = -11,, X, S1, ((,, S1, SJ, [1,; ,
If J = Figure 3(a) then G gg
- Reeulatory Relationship
If J = (b) then L Xz)
LLatent Confounders
If J = (c) then S Xy), G2 f G¹.
* Selection Bias
If J == (d) then L' + X)).
+ +'C&
If J = other patterns not included in Figure 3 then Unk + X,)
Step Correct CI patterns by blocking d-separated paths:
repeat Step 3 (X,, Xy) L', S,, Unk, conditioning on the subsets of their non-endpoints.
If (X,, X2) E L' and correct to (a) then L' := = {(X,,
If (X,, x))= E S and correct to (a) (b) (c) then S == (X,,,,))
until no further edges (X,,X) need correcting
return G LU L.
Latent confounders (GISL), a general nonparametric algorithm, to detect the existence of selection bias and latent confounders, and identify regulatory relationships. We first obtain adjacencies from observational data Do. as it provides the sparest skeleton by statistical criteria such as CI tests or the Bayesian Information Criterion (BIC) [53]. The resulting skeleton, representing conditional dependencies, guides further exploration of underlying causal structures (Step 2).
To capture the differences in symmetry and perturbation effects, examine CI patterns between perturbation indicators II and X' from both observational data Do and single gene perturbation data
{Dx X E T} to identify structures that can inform the skeleton as illustrated in Figure 3 for gene pair (X, (X, The CI results are collected and represented by the set J in Step 3. W While examples offer initial intuition in distinguishing different causal structures, more complex scenarios require further careful analysis for completeness, especially when multiple paths connect node pairs, including d- d-separated paths (representing conditional independence) and inducing paths (representing conditional dependence). defined as follows:
Definition 3.5 (d- de-sparratinn [54] ). Let & be DAG, and let A. B, and C be three disjoint sets of nodes in G. We say that AA and are d-separated by C in G if and only every path between a node in A and node in B is blocked by C.
Definition 3.6 (Inducing path (53). In a DAG with and S, X, Y are any two vertices, and L, S are disjoint sets of vertices not containing X, Y. A path between X and Y is called an inducing path relative to (C, S) if and only if every non-endpoint vertex on p is either in or a collider, and every collider p is an ancestor of either X, Y, or a member of S.
If the inferred structure includes both inducing paths and d-separated paths, further correction is required to remove the dependencies induced by d-separated paths (Step 4), as illustrated in Example 2 in Section 3.2 More complex inducing-patii configurations and their implications for recovering the true causal structure are discussed in Section C Given the adjacencies provided by the skeleton, d d-separated paths can be blocked by conditioning detect the presence of selection bias and latent confounders between each pair of genes as well as regulatory relationships. To represent equivalence classes under latent confounding and selection bias, we adopt the edge semantics from the Partial Ancestral Graph (PAG) framework.
adjacent nodes. After correction, GISL can
Example 2 Consider simple case involving two variables X and Y, connected by two paths: an inducing path + L + Y, and d-separated path + - Y. Without conditioning on Z, the observed CI pattern reflects the influence of both spacifsiiii,i.ll........ both confounding via L and dependency via Z, denoted as "C I.". Only by conditioning on Z. can we block the d-separated path and correctly identify the underlying confounding structure, yielding the CI pattern "L".

---

Table Accuracy % of GISL in identifying selection bias on synthetic data. We report the mean and variance values of accuracy across 10 independent graphs for each configuration.
10
15
20
25
10
15
20
25
Hard intervention
Soft intervention
500
23.7
68.7
72.0::22.2
70.0
60.4:1.4.0
66.6.1:6.6
68.2.1 17.3
67.5.1 15.9
1,000
70.4=20.0
70.0=19.6
74.6 18.6
15.9
75.0 16.7
74.8
80.2
72.1 14.2
1,500
71.0=22.3
72.5=18.6
17.5+
76.9
72.5+ 15.6
80.8+11.3
78.24
77.54 16.3
2,000
73.4=22.6
75,4=19.6
15.7=15.4
13.9
80.0#11.0
73.3
78.24 13.6
75.1 14.7
3.3
Identifiability result of the GISL algorithm
We analyze the identifiability of the proposed GISL framework in detecting selection bias, latent confounders, and regulatory relationships and conclude the following theorem:
Theorem 3.7. (Identifiability of GISL) Let the observational and perturbation data be generated from the DAG model 9 defined in Equation (T). Under Markov 3.7 and faithfulness 3.2 assumptions. when the sample size oc, the causal relationships, selection processes, and latent confounders are identifiable up to the equivalence classes of four types of CI patterns in Figure among variables that are subject interventions. Moreover, the presence of selection processes and latent confounders (existing or not) is identifiable.
Remark 3.8. Identifiability can be established intuitively. Take the causal relation X Y as an example. Structurally, X can only point out with tails, and Y can only be pointed to with arrowheads. To reproduce this structure using other inducing paths, one must substitute the tail with a v-structure (S) or the arrowhead with the hidden common cause (L). If selection is used to replace the tail, then Ix and Y would not be conditionally independent given X (discussed in 3.44) which contradicts the CI pattern of 'C' in Figure a). Alternatively, if a latent confounder replaces the arrowhead, since hidden common causes only contribute to two arrowheads, the requirement of tails and removing the effect of the extra arrowhead necessitates an intermediate node to keep the CI patterns with 'C'. As the intermediate node cannot be blocked, this requires v-structure and an edge point to Y following Definition 3.6, Then, X must point to the intermediate node to form the V-structure with the extra arrowhead. Therefore, X is the ancestor of Y. For more details on proofs of the identifiability of GISL in the presence of latent confounders and selection process, please kindly refer to Appendix D
4 Experiments
We begin by evaluating GISL's ability detect selection bias on synthetic data. Next, we present the benefits of considering selection bias in causal discovery, specifically, benchmark its performance in identifying regulatory relationships tasks, against established baselines. Finally, we apply GISL to real-world gene expression data, using Z-scores as proxy ground truth.
4.1 Identify the selection bias on synthetic data
Nonparametric settings. To better reflect the complexity of gene expression, adopt nonpara- metric structural causal model (SCM) that accommodates latent confounders and selection bias with- assuming specific functional distributional forms: L; Ei,; X; (1 - I,) fo (Xa)(X.)) Ez.] + (pa + Ez.] with selection governed by (X,, XX), C. Es, and E,, are Gaus- sian noise terms with randomly selected means and variances, C' is randomly selected threshold, and for fa are sampled from diverse nonlinear functions: linear, square, sin and log. 4, E [0. indicates the gene is perturbed or not, with corresponding perturbation function f:. Following Section 2.1 we simulate gene perturbations using both hard and soft interventions. In an ideal hard intervention (knockout), one would set do do(X, ~ (a, b)). with the interval (a, b] randomly chosen for each intervention. Soft interventions (knockup and knockdown) model up- or down-regulation by adding uniform noise to the original expression: X; + X, + e, where 8 ~ (/, e) for knockup, and e 21(-d, 0) for knockdown. The magnitudes and d define the strength of the up- down-regulation.
0). But to capture off-target variability, we instead sample
7

---

GES
JCI-SSP
IGSP
UT4GGPP
(Ours)
DAG
DAG Precision
530
s 1300
n- 2038
0.8
0.75
0.4
0.2
0.25
0.8
0.75
Soft
0.4
0.25
0.2
DAG Recall
DAG SHD
n= 500
1998
n= 2000
n- 500
1500
2000
0.8
Hard
25
0.4
100
0.6
Salf
0.4
10
25:
15
16
15
dP Varlables)
Voriablest
Figure 4: Comparison results in identifying regulatory relations under four metries: DAG F1, DAG Precision, DAG Recall, and DAG SHD (Structural Hamming Distance). All values are averaged over 10 runs with different random seeds. Error bars represent the 95% confidence interval.
Synthetic data generation. We first instantiate DAGS using Erdös-Rényí model [56] with the number of edges equal to the number of variables. Then groups of synthetic gene expression data are generated following the SCM and intervention protocols detailed above, sweeping the configurations: the number of observation data n ranged in {500, 1500, 2000}. the number of variables |X|| E {10, 15, 20, 25). For each (n, (x() configuration, and separately for hard and soft intervention, we instantiate 10 independent structures by randomly choosing 1-3 selection pairs and 1-3 confounding pairs, reflecting the complexity of the gene expression scenario.
Comparison results (selection bias, Table I and Table 2). We report the evaluation results of GISL (KCI for CI test with
0.05) in identifying selection bias the synthetic data detailed above. The accuracy score is defined as the percentage of identified selection pairs that are aligned with ground truth in all positive predictions. In T GISL accurately captures selection bias, especially when more data points are collected. Take an example with |X' 20, leveraging more data to precisely characterize distribution changes yields a ~ % accuracy gain. Further- more, to evaluate GISL under increasing selec- tion complexity, we fix | = 15 and n 2000, and sweep the number of selection processes (f/a E (1, 2,2,,344) Table 2 alongside the identifi-
Table 2: Accuracy % under more selection pairs; Upper-bound baseline: Collider identification.
Hard
81.6 8288148
68.1 1 18.1
70.8+ 20.6
Soft
15.5
76.1
70.6+
69.5 14.3
Collider
90.0+5.00
85.0+5.3
83.3+-.66
3.55+399
cation of collider as an approximate upper bound. As more genes are subjected to biased selection, the induced observational constraints make detection increasingly challenging for CI-based methods. Despite this, GISL maintains strong performance across all settings, detecting selection bias in around 70% of affected gene pvirs-even when more than half the variables are biased. To the best of our knowledge, existing methods designed to identify selection bias arising from survival constraints in gene regulatory tasks. As such, we do not report baseline comparisons in this section.
4.2 Identify regulatory relationships on synthetic data
Baselines. To rigorously evaluate GISL's ability to disentangle true regulatory relationships from spurious perrddnncees-pattculatt those arising from selection bias-we conduct experiments on synthetic data, comparing GISL. against robust baseline methods: GIES 40), IGSP [7.. UT-IGSP
[57 and JCI-GSP method used in [57], which is an extension of JCI
[58 with GSP 59].
8

---

Benchmarking results (regulatory relationships). Figure illustrates that GISL achieves superior overall performance in identifying regulatory relationships (DAG F,). Among evaluation metries, precision emerges most critical, as it directly reflects GISL's advantage in disentangling regulatory relationships from spurious dependencies induced by selection bias. While existing approaches often misattribute selection biases causal relationship confounders, GISL effectively distinguishes these phenomena, resulting in markedly improved precision. We emphasize that the objective of this paper is not to propose definitive solution to causal discovery, but rather to highlight that ignoring selection bias can lead to flawed casual interpretation. Further comparisons including with FCI [55], ICD [60], and non-causal baselines are presented in Appendix Moreover, the evaluation in the robustness of GISL across noise levels, structural and functional complexity of the SCM, and larger graphs are discussed in Appendix E
4.3 The presence of selection bias on single-cell gene expression data
Datasets. We use three real-world scRNA-seq perturbation datasets: Dixit 32 and Adamson [6], from K562 leukemia cells (5,012 and 5,060 genes under 19 and 86 single-gene perturbations), and Norman [31, from A549 lung carcinoma cells (5,045 genes under 105 perturbations). Each perturbation targets one gene and is profiled across numerous cells.
Selection bias identification on real-world data (Figure 5). We plot the detection re- sult on the Dixit dataset in Figure [ Re- sults for the other two datasets and detailed analysis are attached in Appendix [F.T] Our method can detect both the underlying se- lection processes and latent confounders at the same time, facilitating the explanation of observed dependence among genes (regula- tory relationship. selection bias, and latent confounders). practice, the explanation of the selection process and latent confounders will guide the biologists that perturbing these genes may not have the expected effect, as the selection process may lead to unafford- able consequences after sample filtering, and latent confounders will have no reaction to perturbation.
Figure 5: Experimental result on 19 perturbed key
genes from perturb-seq [32]. S and L imply the de- tected selection and confounded pairs, blue edges are the regulatory interactions priorly known 33
Evaluation using Z-score. evaluate whether the identified gene pairs are truly subject to selection bias- without relying on ground-trul toound-truth labels-we introduce Z-score: score function that captures the changes in growth rate before and after perturbations (29), de- tailed calculation refers to Appendix G Notably, both
Z-Score Distribution of Genes by Selection
Z-Score
the increase (+) and decrease (-) result in a Z-score absolute value larger than O. We have the follow- ing assertion: true selection pairs must exhibit high absolute Z-score values, while the opposite is not nec- essarily true. To see this, consider X; + S - X,, where X; and X, are genes under the selection pro- cess fa. Perturbing either gene causes differential cell survival and, thus, a higher Z-score. However, high Z-score does not necessarily imply direct selection processes: in X; + S + X, + Xt. perturbing X₁ indirectly alters the joint distribution of (X,,X) under S. yielding high Z-score for X₁ even though X₁ is not under direct selection process. Therefore, only gene pairs for which both exhibit high Z-scores can be selection pairs, and simultaneously, an accurate algorithm for detecting selection pairs should also perform well on detecting pairs with high Z-scores. By applying Z-score threshold of 0.15 0.15 to define high-scoring gene pairs, GISL. achieves accuracies of 77.7% on Dixit (see Figure 6 for comparison in Z-score distribution), 80.2% Adamson, and 74.8% on Norman.
In Selectin
Gene Group
Figure 6:
Comparison of Z-score distribution between genes in selection pairs and others.

---

Conclusion and Discussion
introduce a novel perspective to gene regulatory network inference (GRNI); selection bias, partic- ularly the often-overlooked survival constraints, in shaping dependencies. Building on theoretical insights into the distinct perturbation symmetries and effects associated with regulation, confound- ing, and selection, we propose GISL general nonparametric algorithm capable of disentangling regulatory relationships, latent confounders, and, crucially, dnlection-inducee dependencies.
Empirical analyses on large-scale gene expression datasets reveal pervasive selection bias in real- world scenarios. Extensive benchmarking on synthetic data further demonstrates that (1) failure to account for selection bias undermines the validity of GRNI methods, and (2) explicitly modeling such bias, as done by GISL, yields more accurate and robust causal discovery. Limitations and broader implications of this work are discussed in Appendix H
Acknowledgment
We would like to acknowledge the support from NSF Award No. 2229881, AI Institute for Societal Decision Making (AI-SDM), the National Institutes of Health (NIH) under Contract R01HL159805, and grants from Quris AI, Florin Court Capital, and MBZUAI WIS Joint Program, and the Al Deira Causal Education project.
References
[1] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe'er. Using bayesian networks to analyze expression data. Proceedings of the fourth annual international conference on Computational molecular biology, pages 127-135, 2000.
[2] Eric H Davidson and Douglas H Erwin. Gene regulatory networks and the evolution of animal body plans. Science, 311(5762):796-800 800, 2006.
[3] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.
[4] Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge university press, 2000.
[5] Seongho Kim. ppcor: an package for a fast calculation semi-partial correlation coefficients. Communications for statistical applications and methods, 22(6):665, 2015.
[6] Alicia T Specht and Jun Li. Leap: constructing gene co-expression networks for single-cell ma-sequencing data using pseudotime ordering. Bioinformatics, 33(5):764-766, 2017.
[7] Vân Anh Huynh-Thin Thu, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. Inferring regulatory networks from expression data using tree-based methods. PloS one, 5(9):e12776, 2010.
[8] Anne-Claire Haury, Fantine Mordelet, Paola Vera- Leca-Licona, and Jean-Philippe Vert. Tigress: trustful inference gene regulation using stability selection. BMC systems biology, 6:1-17, 2012.
[9] Thomas Moerman, Sara Aibar Santos, Carmen Bravo González-Blas, Jaak Simm, Yves Moreau, Jan Aerts, and Stein Aerts. Grnboost2 and arboreto: efficient and scalable inference of gene regulatory networks. Bioinformatics, 35(12)2119-2211, 2019.
[10] Patrick Meyer, Kevin Kontos, Frederic Lafitte, and Gianluea Bontempi. Information theoretic inference of large transcriptional regulatory networks. EURASIP journal on bioinformatics and systems biology, 2007:1-9, 2007.
[11] Thalia E Chan, Michael PH Stumpf, and Ann C Babtie. Gene regulatory network inference from single-cell data using multivariate information measures. Cell systems, 5(3):251-22 2017.
10

---

[12] Andrea Ocone, Laleh Haghverdi, Nikola S Mueller, and Fabian Theis. Reconstructing gene regulatory dynamies from high-dimensional single-cell snapshot data. Bioinformaties, 31(12):19--16, 2015.
[13] Hirotaka Matsumoto, Hisanori Kiryu, Chikara Furusawa, Minoru SH Ko, Shigeru BH Ko, Norio Gouda, Tetsutaro Hayashi, and Itoshi Nikaido. Scode: an efficient regulatory net- work inference algorithm from single-cell ma-seq during differentiation. Bioinformatics, 33(15)2231-2321) 2017.
[14] Hirotaka Matsumoto and Hisanori Kiryu. Scoup: a probabilistic model based on the ornstein- uhlenbeck process to analyze single-cell expression data during differentiation. BMC bioinfor- matics, 17:1-16, 2016.
[15] Manuel Sanchez-Castille David Blanco, Isabel M Tienda-Luns MC Carrion, and Yufei Huang. A bayesian framework for the inference of gene regulatory networks from time and pseudo-time series data. Bioinformatics, 34(6))964 970, 2018.
[16] Steven Woodhouse, Nir Piterman, Christoph M Wintersteiger, Berthold Göttgens, and Jasmin Fisher. Scns: graphical tool for reconstructing executable regulatory networks from single- cell genomic data. BMC systems biology, 12:1-7, 2018.
[17] Yuhao Wang. Liam Solus, Karren Yang, and Caroline Uhler. Permutation-bn based causal inference algorithms with interventions. Advances Neural Information Processing Systems, 30, 2017.
[18] Anastasiya Belyaeva, Chandler Squires, and Caroline Uhler. Dei: learning causal differences between gene regulatory networks. Bioinformatics, 37(18):3067-306, 2021.
[19] Mathieu Chevalley, Yusuf Roohani, Arash Mehrjou, Jure Leskovec, and Patrick Schwab. Causalbench: A large-scale benchmark for network inference from single-cell perturbation data. arXiv preprint arXiv:2210. 17283, 2022.
[20] Jiaqi Zhang, Louis Cammarata, Chandler Squires, Themistoklis P Sapsis, and Caroline Uhler Active learning for optimal intervention design in causal models. Nature Machine Intelligence, 1066 2023.
[21] Audrey P Gasch, Paul T Spellman, Camilla M Kao, Orna Careel-Harel, Michael B Eisen, Gisela Storz, David Botstein, and Patrick O Brown. Genomic expression programs in the response of yeast cells to environmental changes. Molecular biology of the cell, 11(12):4241- 4257, 2000.
[22] Luisa Statello, Chun-Jie Guo, Ling-Ling Chen, and Maite Huarte, Gene regulation by long non- coding mas and its biological functions. Nature reviews Molecular cell biology, 22(2):96-118, 2021.
[23] Sergey V Razin and Alexey A Gavrilov. Non-coding mas in chromatin folding and nuclear organization. Cellular and Molecular Life Sciences, 78( (811445449-550 2021.
[24] Oliver Stegle, Leopold Parts, Matias Piipari, John Winn, and Richard Durbin. Using proba- bilistic estimation of expression residuals (peer) to obtain increased power and interpretability of gene expression analyses. Nature protocols, 733:5000-507, 2012.
[25] Lingfei Wang and Tom Michoel. Efficient and accurate causal inference with hidden confounders from genome-transerritooee variation data. PLDS computational biology, 13(8):c1005703, 2017.
[26] Muhammad Ammar Malik and Tom Michoel. Restricted maximun-likelihood method for learning latent variance components in gene expression data with known and unknown con- founders. G3, 12(2)jjkab410, 2022.
[27] Chunlin Li, Xiaotong Shen, and Wei Pan. Nonlinear causal discovery with confounders. Journal of the American Statistical Association, 1295546) 2024.
11

---

[28]
James J. Heckman. Sample selection bias as a specification error. Econometrica, page 153, Dec 1978.
[29]
DepMap. Depmap 23q4 public. figshare+. Journal Name, 2023.
[30] Joshua M Dempster, Isabella Boyle, Francisca Vazquez, David Root, Jesse Boehm, William Hahn, Aviad Tsherniak, and James M McFarland. Chronos: a crispr cell population dynamics model. BioRxiv, pages 2021-02, 2021.
[31] Norman Thomas M., Horlbeck Max A., Replogle Joseph M., Ge Alex Y., Xu Albert, Jost Marco, Gilbert Luke A., and Weissman Jonathan S. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science, Nov 2019.
[32] Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, CharlesP. Fulco, Livnat Jerby-Arnon, Ne- manjaD. Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, Britt Adamson, ThomasM. Norman, EricS. Lander, JonathanS. Weissman, Nir Friedman, and Aviv Regev, Perturb-seq: Dissecting molecular circuits with scalable single-cell rna profiling of pooled genetic screens. PMC, Nov 2016.
[33] Maxim V Kuleshov, Matthew R Jones, Andrew D Rouillard, Nicolas F Fernandez. Qiaonan Duan, Zichen Wang. Simon Koplev. Sherry L Jenkins, Kathleen M Jagodník, Alexander Lachmann, et al. Enrichr: comprehensive gene set enrichment analysis web server 2016 update. Nucleic acids research, 44(W1):99-W99) 2016.
[34] Kaori Sasai, Warapen Treekitkammonnkkll, Kazuharu Kai, Hiroshi Katayama, and Subrata Sen. Functional significance of aurora kinases-p53 protein family interactions in cancer. Frontiers
in oncology, 6:247, 2016.
[35] Paraskevi Briassouli, Florence Chan, Kay Savage, Jorge S Reis-Filho, and Spiros Linardopou- los. Aurora-a regulation of nuclear factor-kb signaling by phosphorylation of ikbo. Cancer research, 67(4):1689-1663 1695, 2007.
[36] Zhaofeng Li, Xueyu Li, Han Yang. Meixiang Huang, Zhu Liu, Zongliang Zhang, Kai Zhao, Xinbao Yin, Guanqun Zhu, Yulian Zhang, al. Torl aipl interacts with p53 enhance cell cycle dysregulation in prostate cancer progression. Molecular and Cellular Biochemistry, pages 1-15, 2025.
[37] Michael Levine and Eric H Davidson. Gene regulatory networks for development. Proceedings of the National Academy of Sciences, 102(14):4936 4942, 2005.
[38] Gregory F Cooper. Causal discovery from data in the presence of selection bias. In Pre- proceedings of the Fifth International Workshop on Artificial Intelligence and Statistics, pages 150-150. PMLR, 1995.
[39] Miguel A Hernán, Sonia Hernández-Diaz, and James M Robins. A structural approach to selection bias. Epidemiology, 15(5):615-625, 2004.
[40] Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs. The Journal of Machine Learning Research, 13(1):2409- 2464, 2012.
[41] Karren Yang, Abigail Katcoff, and Caroline Uhler. Characterizing and learning equivalence classes of causal dags under interventions. In International Conference on Machine Learning, pages 5541-5550. PMLR, 2018.
[42] Sara Magliacane, Tom Claassen, and Joris M Mooij. Ancestral causal inference. Advances in Neural Information Processing Systems, 29, 2016.
[43] Kevin B Korb, Lucas Hope, Ann E Nicholson, and Karl Axnick. Varieties of causal intervention. In Pacific Rim international conference artificial intelligence, pages 322-331. Springer, 2004.
[44] Kevin Hoover. The logic of causal inference: Econometrics and the conditional analysis of causation. Economics & Philosophy, 6(2):207-234, 1990.
12

---

[45] Jin Tian and Judea Pearl. Causal discovery from changes. arXiv preprint 2013.
[46] Sofia Triantafillou and Ioannis Tsamardinos. Constraint-based causal discovery from multiple interventions over overlapping variable sets. The Journal of Machine Learning Research, 16(1):1147-2055, 2015.
[47] Nicolai Meinshausen, Alain Hauser, Joris M Mooij, Jonas Peters, Philip Versteeg, and Peter Bühlmann. Methods for causal inference from gene perturbation experiments and validation. Proceedings of the National Academy of Sciences, 113(27):7361-736) 2016.
[48] AmirEmad Ghassami, Saber Salehkaleybar, Negar Kiyavash, and Kun Zhang. Learning causal structures using regression invariance. Advances in Neural Information Processing Systems, 30, 2017.
[49] Felix Elwert and Christopher Winship. Endogenous selection bias: The problem of condition- ing on collider variable. Annual review of sociology. 40(1):31-53, 2014.
[50] Dylan J Taylor, Surya B Chhetri, Michael Tassia, Arjun Biddanda, Stephanie M Yan, Genevieve L Wojcik, Alexis Battle, and Rajiv C McCoy. Sources of gene expression variation in a globally diverse human cohort. Nature, (32(8023):122-1)) 2024.
[51] Niklas Mähler, Jing Wang. Barbara K Terebieniec, Pär K Ingvarsson, Nathaniel R Street, and Torgeir R Hvidsten. Gene co-expression network connectivity an important determinant of selective constraint. PLoS genefics, 13(4):e1006402) 2017.
[52] Steffen L Lauritzen. Graphical models, volume 17. Clarendon Press, 1996.
[53] Joseph Ramsey, Madelyn Glymour, Ruben
and Clark Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. Inter- national journal of data science and analytics, 3:121-129, 2017.
[54] Judea Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference. Morgan Kaufmann, 1988.
[55] Jiji Zhang. On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence, 172(16-17)1183--1899 2008
[56] Paul Erdös, Alfréd Rényi, et On the evolution of random graphs. Publ, math. inst. hung. acad. sci, 5(1):17-60, 1960.
[57] Chandler Squires, Yuhao Wang, and Caroline Uhler. Sermutation-hasee causal structure leam- ing with unknown intervention targets. In Conference on Uncertainty in Artificial Intelligence, pages 1048- PMLR, 2020.
[58] Joris M Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts. Journal of machine learning research, 21 2020.
[59] Liam Solus, Yuhao Wang, and Caroline Uhler. Consistency guarantees for greedy permutation- based causal inference algorithms. Biometrika, 108(4):795-814 2021.
[60] Raanan Y Rohekar, Shami Nisimov, Yaniv Gurwicz, and Gal Novik. Iterative causal discov- ery in the possible presence of latent confounders and selection bias. Advances in Neural Information Processing Systems, 34:2454-2465. 2021.
[61] Britt Adamson. Thomas M Norman, Marco Jost, Min Cho, James K Nuñez, Yuwen Chen, Jacqueline E Villalta. Luke A Gilbert, Max A Horlbeck, Marco Y Hein, et al. A multiplexed single-cell crispr screening platform enables systematic dissection of the unfolded protein response. Cell, (6777)1867-1188 2016.
[62] Peter Kharchenko, Lev Silberstein, and David Scadden. Bayesian approach to single-cell differential expression analysis. Nature methods, 11(7):740-742, 2014.
13

---

[63] Mingchao Li, Qing Min, Matthew C Banton, and Xinpeng Dun. Single-cell regulatory network inference and clustering identifies cell-type specific expression pattern of transcription factors in mouse sciatic nerve. Frontiers in Cellular Neuroscience, 15:676515, 2021.
[64] Atul Deshpande, Li-Fang Chu, Ron Stewart, and Anthony Gitter. Network inference with granger causality ensembles on single-cell transcriptomics. Cell reports, 38(6), 2022.
[65] Shuo Li, Yan Liu, Long-Cher Shen, He Yan, Jiangning Song, and Dong-Jun Yu. Gmfgrm: a matrix factorization and graph neural network approach for gene regulatory network inference. Briefings in Bioinformatics, 25(2)bbbad52, 2024.
[66] Hung Nguyen, Duc Tran, Bang Tran, Bahadir Pehlivan, and Tin Nguyen. A comprehensive survey of regulatory network inference methods using single cell ma sequencing data. Briefings in bioinformatics, 22(3):bbaa190, 2021.
[67] Jiaqi Zhang, Chandler Squires, and Caroline Uhler. Matching desired causal state via shift interventions. Advances in Neural Information Processing Systems, 34: 1999339-99,,,, 2021.
[68] Loka Li, Ignavier Ng. Gongxu Luo, Biwei Huang, Guangyi Chen, Tongliang Liu, Bin Gu, and Kun Zhang. Federated causal discovery from heterogeneous data. arXiv preprini arXiv:2402. 13244,, 2024.
[69] David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(N0v):507-554, 2002.
[70] Loka Li, Haoyue Dai, Hanin Al Ghothani, Biwei Huang, Jiji Zhang, Shahar Harel, Isaac Bentwich, Guangyi Chen, and Kun Zhang. On causal discovery in the presence of deterministic relations. Advances Neural Information Processing Systems, 37:130920-13099
2024.
[71] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non- gon-gasssian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
[72] Patrik Hoyer, Dominik Janzing, Joris M Mooij. Jonas Peters, and Bernhard Schölkopf. Nonlin- ear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008.
[73] Kun Zhang and Aapo Hyvarinen. On the identifiability of the post-nonlinear causal model. arXiv preprint arXiv: /205 2599, 2012.
[74] Klea Ziu, Slavomír Hanzely, Loka Li, Kun Zhang, Martin Takáč, and Dmitry Kamzolov. vdag Projected stochastic approximation iteration for dag structure learning. arXiv preprint arXiv:2410.23862, 2024.
[75] Peter L. Spirtes, Christopher Meek, and Thomas S Richardson. Causal inference in the presence of latent variables and selection bias. Conference on Uncertainty in Artificial Intelligence, 1995.
[76] Wei Chen, Kun Zhang, Ruichu Cai, Biwei Huang, Joseph Ramsey, Zhifeng Hao, and Clark Glymour. Fritl: A hybrid method for causal discovery the presence of latent confounders. arXiv preprint arXiv: 203. 14238, 2021.
[77] David Kaltenpoth and Jilles Vreeken. Nonlinear causal discovery with latent confounders. In International Conference on Machine Learning, pages 15639-15654. PMLR, 2023.
[78] Ruichu Cai, Zhiyi Huang, Wei Chen, Zhifeng Hao, and Kun Zhang, Causal discovery with latent confounders based higher-order cumulants. In International conference on machine learning, pages 3380-3407. PMLR, 2023.
[79] Angli Xue, Seyhan Yazar, Drew Neavin, and Joseph E Powell. Pitfalls and opportunities for applying latent variables in single-cell equi analyses. Genome Biology. 24(1) (::3,, 2023.
[80] Philip Versteeg, Joris Mooij, and Cheng Zhang. Local constraint-based causal discovery under selection bias. In Conference Causal Learning and Reasoning, pages 840-860. Pmlr, 2022.
14

---

[81] David Kaltenpoth and Jilles Vreeken. Identifying selection bias from observational data. AAAI, 2023.
[82] Kun Zhang. Jiji Zhang, Biwei Huang, Bernhard Schölkopf, and Clark Glymour. On the identifiability and estimation of functional causal models in the presence of outcome-dependent selection. In UAI, 2016.
[83] Elias Bareinboim, Jin Tian, and Judea Pearl. Recovering from selection bias in causal and statistical inference. In Probabilistic and Causal Inference: The Works of Judea Pearl, pages 433 450.. Association for Computing MachineryNew YorkNYUnited States, 2022.
[84] Gregory Cooper and Changwon Yoo. Causal discovery from mixture of experimental and observational data. arXiv preprint arXiv:1301.6686, 1999.
[85] Daniel Eaton and Kevin Murphy. Exact bayesian structure learning from uncertain interven- tions. In Artificial intelligence and statistics, pages 107-114, 14. PMLR, 2007.
[86] Jin Tian and Judea Pearl. Causal discovery from changes. arXiv preprint arXiv: I 301.2312, 2001.
[87] Thomas Richardson and Peter Spirtes. Ancestral graph markov models. The Annals of Statistics, 30(4);922-1030, 2002.
[88] Antti Hyttinen, Patrik O Hoyer, Frederick Eberhardt, and Matti Jarvisalo. Discovering cyclic causal models with latent variables: A general sat-based procedure. arXiv preprint arXiv:1309.60 6886666,,,,,66,, 2013.
[89] Fuyuan Cao. Yunxia Wang, Kui Yu, and Jiye Liang. Causal discovery from unknown inter- ventional datasets over overlapping variable sets. IEEE Transactions on Knowledge and Data Engineering, 2024.
[90] Murat Kocaoglu, Alex Dimakis, and Sriram Vishwanath. Cost-optimal learning of causal graphs. In International Conference on Machine Learning, pages 1875-1884. PMLR, 2017.
[91] Adam Li, Amin Jaber, and Elias Bareinboim. Causal discovery from observational and interventional data across multiple environments. Advances in Neural Information Processing Systems, 16:16942-16956, 2023.
[92] Haoyue Dai, Ignavier Ng, Jianle Sun, Zeyu Tang, Gongxu Luo, Xinshuai Dong, Peter Spirtes, and Kun Zhang. When selection meets intervention: Additional complexities in causal discovery, arXiv preprint 20X3v22503.07302. 2025.
[93]
Saliba, Alexander J. Westermann, Stanislaw A. Gorski, and Jörg Vogel. Single-cell ma-seq: advances and future challenges. Nucleic Acids Research, page 8845-8860, Aug 2014.
[94] Jennifer A. Doudna and Emmanuelle Charpentier. The new frontier of genome engineering with crispr-cas9. Science, Nov 2014.
[95] Matthew H Larson, Luke A Gilbert, Xiaowo Wang, Wendell A Lim, Jonathan S Weissman, and Lei S Qi. Crispr interference (crispri) for sequence-specific control of gene expression. Nature protocols, 8(1 1):2180-2196 2013.
[96] Albert W Cheng, Haoyi Wang, Hui Yang, Linyu Shi, Yarden Katz, Thorold W Theunissen, Sudharshan Rangarajan, Chikdu S Shivalila, Daniel B Dadon, and Rudolf Jaenisch. Multi- plexed activation of endogenous genes by crispr-on, ma-guided transcriptional activator system. Cell research, 1631011163-1111, 2013.
[97] Britt Adamson, ThomasM. Norman, Marco Jost, MinY. Cho, JamesK. Nuñez, Yu-Wen Chen, JacquelineE. Villalta, LukeA. Gilbert, MaxA. Horlbeck, Marco'. Hein, RyanA. Pak, AndrewN. Gray, CarolA. Gross, Oren Parnas, JonathanS. Weissman, Atray Dixit, and Aviv Regev. A multiplexed single-cell crispr screening platform enables systematic dissection of the unfolded protein response. PMC,PMC, Nov 2016.
15

---

[98] Alexandra B Keenan, Denis Tore, Alexander Lachmann, Ariel K Leong, Megan L Woj- ciechowicz. Vivian Utti, Kathleen Jagodník, Eryk Kropiwnicki, Zichen Wang, and Avi Ma' ayan. Chea3: transcription factor enrichment analysis by orthogonal omics integration. Nucleic acids research, 47(W1):W212-W22 2019.
[99] Edward Y Chen, Christopher Tan, Yan Kou, Qiaonan Duan, Zichen Wang, Gabriela Vaz Meirelles, Neil R Clark, and Avi Ma ayan. Enrichr: interactive and collaborative html5 gene list enrichment analysis tool. BMC bioinformatics, 14(1):1-14, 2013.
[100] Natalie de Souza. The encode project. Nature methods, 9111::046-104 2012.
[101] Cecilia Marmai, Rachel E Sutherland, Kevin K Kim, Gregory M Dolganov, Xiaohui Fang, Sophia S Kim, Shuwei Jiang, Jeffery A Golden, Charles W Hoopes, Michael A Matthay, et al. Alveolar epithelial cells express mesenchymal proteins in patients with idiopathic pulmonary fibrosis. American Journal of Physiology- Lung Cellular and Molecular Physiology. 301(1):L71-L77) 2011.
[102] Kazuyoshi Kuwano, Ritsuko Kunitake, Masayuki Kawasaki, Yoshitugu Nomoto, Naoki Hag- imoto, Yoichi Nakanishi, and Nobuyuki Hara. P2Iwafl/cipl/sii and p53 expression in association with dna strand breaks in idiopathic pulmonary fibrosis. American journal of respiratory and critical care medicine, 154())477-483, 1996.
[103] Mark Duncan, Teresa A DeLuca, Hsin- Yu Kuo, Minchang Yi, Milan Mrksich, and William M Miller. Sirtl is critical regulator of k562 cell growth, survival, and differentiation. Experi- mental cell research, 344(1):40-52, 2016.
[104] Jana Zecha, Florian P Bayer, Svenja Wiechmann, Julia Woortman, Nicola Berner, Julian Müller, Annika Schneider, Karl Kramer, Mar Abril-Gil, Thomas Hopf, al. Decrypting drug actions and protein modifications by dose-and time-resolved proteomics. Science, 380(6640):93-14 101, 2023.
[105] Michele Olivieri, Tiffany Cho, Alejandro Álvarez-Quilón, Kejiao Li, Matthew J Schellenberg, Michal Zimmermann, Nicole Hustedt, Silvia Emma Rossi, Salomé Adam, Henrique Melo, et al. A genetic map the response to dna damage in human cells. Cell, 182(2):481-466, 2020.
[106] Kathiresan Selvam, Jun Xu, Hannah E Wilson, Juntaek Oh, Qingrong Li, Dong Wang, and John Wyrick. Elfl promotes transcription-coupled repair in yeast by using its c-terminal domain to bind tfihh. Nature Communications, 15(1):6223, 2024.
16

---

Appendix
A Concepts
Definition A.1 (Marginal independence test). Check whether two variables and Y independent of each other without considering any other variables, Mathematically: X IL Y, if and only if X and Y are independent in the overall data distribution.
Definition A.2 (Conditional independence test), Evaluate whether two variables X and Y are independent given third variable or set of variables Z. Mathematically: X Y|Z, and only X and Y are independent conditioned on Z.
Definition A.3 (d-connection). Every path from a node in X to a node Y is d-connected by Z, and only if X and Y are always conditionally dependent given Z.
Definition A.4 (Partial identifiability). The causal graph is partially identified and only not all causal structures be uniquely determined from the available data and assumptions, and only of plausible structures can be uniquely determined.
Definition A.5 (Identifiabilitty. The causal graph is identified if and only if all causal structures are uniquely determined from the observational data.
B
Related work
In this section, we provide comprehensive review of the literature on causal discovery methods and briefly categorize other computational Gene Regulatory Network Inference (GRNI) methods.
B.1 Gene regulatory network inference
Over the past decades, numerous methods for GRNI have been developed, encompassing computa- tional and causal approaches. Computational models, represented by boolean model, differential equation, gene correlation, and correlation ensemble over pseudo- focus exploring dependen- cies among genes 62
66). These methods focus on dependencies without explicitly considering causal relationships. In contrast, causal models go beyond mere dependence, aiming
13
63
64
65
uncover genuine causal relationships within Gene Regulatory Networks (GRNs) [7 18 67].
B.2
Causal discovery
There are constraint-based causal discovery 3 68 score-based ones 66 70 and methods that
utilize properties of functional forms in the true causal process 77 72 73 There subset of work
considering causal discovery continuous optimization problem (74). These pioneering works provide comprehensive frameworks that follow different principles for causal discovery. However, they have not yet addressed the challenges posed by latent confounders and selection bias.
Causal discovery under latent confounders. When latent confounders are present, confounding effects give rise to spurious dependencies, complicating causal discovery. The FCI (75) 55 algorithm was the first address this challenge, producing equivalence classes represented by Partial Ancestral Graphs (PAGs). Building on FCI's results, (76) leveraged the Triad condition to identify shared latent confounders and infer causal relationships between measured variables under linear constraints. [77 utilized an autoencoder framework to reconstruct nonlinear causal relationships among observed variables while simultancously accounting for the presence of latent confounders. [78] estimated the mixing matrix using higher-order cumulants and introduced the testable One-Latens One-Latent-Compont Compnnnet condition to identify latent variables and establish causal orders. In the context of GRNI, some studies have begun addressing the issue of latent confounders, with focus recovering these confounders 79 and uncovering causal relationships among genes 19
Causal discovery under selection bias. Some fundamental works focused on identifying selection and causal relationships by finding Y-structure [8), identifying selection bias under certain parametric assumptions [8]. studying the identifiability and estimation of functional causal models under the outcome-dependent selection structure [82. and recovering the conditional probability from selection- biased data [83). However, these methods are limited to either parametric assumptions, i.c., lincar
17

---

Gaussian, or putcome-dependent selection structure, which are unsuitable for the general setting (nopparametric. latent confounders, selection bias) and the pairwise selection context of GRNI.
Causal discovery under latent confounders and selection bias When latent confounders and selection bias are present, their hidden nature introduces spurious dependencies between observed variables. These dependencies compromise the core properties used in causal discovery, leading to loss of identifiability and making difficult distinguish true causal relationships. The FCI algorithm was the first attempt to discover ancestral relationships, but it is limited to identifying an equivalence class constrained by the structural information of v-structures. Additionally, significant ambiguities persist regarding the selection structure. Similarly, other approaches have resulted in ancestral equivalence classes that are constrained by graphical properties [60].
75 55]
B.3 Interventional Causal discovery
Early Attempts in Interventional Causal Discovery. The earliest Bayesian approaches [4 85 es- timated the posterior distribution of DAGs using both observational and interventional data. However, these methods did not address key challenges such as identifiability or equivalence class characteri- zation. [86] was the first explore identifiability and Markov equivalence of interventional causal discovery. They focused on single-variable interventions with mechanism changes (soft interventions) and provided graphical criterion for determining when two DAGs are indistinguishable, though no formal representation of the resulting equivalence class was introduced.
When latent confounders are involved. In the pure observational data and nonparametric causal discovery setting, the frameworks of MAG and FCI have been well established [87 55 For interventional causal discovery, various methods have been proposed to address latent variables based
on measuring overlapping variables across different interventions [8] 46 89 and invariance [90] 85
They are either lying under the umbrellas of FCI and the augmented DAG frameworks or using parametric assumptions.
42.
When selection is involved In purely observational and nonparametric causal discovery, selection is typically constrained by structural limitations [82]. However, in the context of interventional causal discovery, various methods have been developed to explicitly address selection bias. These methods leverage interventional data to disentangle the effects of selection mechanisms from genuine causal relations [91 [88. Similarly, these methods are still limited to equivalence classes under the [2] discussed how selection interacts with intervention, and built twin interventional graph model the selection that happens before intervention.
umbrella of FCI. Moreover,
C Understanding inducing paths
X
X
Z
Y
Z
(a)
(b)
(d)
Figure 7: Illstrationn of inducing paths, where X and Y are always d-connected.
A most generalized model might include latent confounders, perturbation indicators, and observed variables involved in the selection process. Nonetheless, such generalized assumptions often render causal relationships indeterminate, making the results less informative.
For example, latent confounders can still make the direct causal relations unidentifiable. Consider the case X + Z + Y with a latent confounder L pointing to both Z and Y shown in Figure (b). Adding direct edge available, as the dependence between X and Y cannot solely be explained by Z. With the selection process, the model +- X <- - ( whether or not to add direct edge X Y. the two scenarios unidentifiable. The causal relationship is identified as limited to ancestor
renders the scenarios equivalent, even if perturbation data /x, /y are
18

---

X
lx
1x
(b)
(c)
Figure 8: Required structure for causal relation, latent confounders, and selection process, where means the nodes on the paths cannot be blocked.
relationships or inducing paths. Similarly, Figure (c) cannot be distinguished from X - S and Figure (d) cannot be distinguished from X < L-Y.
After understanding inducing paths. these paths (which are always d-connected) can mimic any structure reflected in conditional independence (CI) patterns. Consequently, the identifiability of causal relationships, latent confounders, and the selection process remains constrained by the true structure or the inducing path. In the representation of ancestral graphs, they are identifiable up to equivalence classes represented by PAG. The details of the proof can be found in the Appendix D
D Proof
The unique structure involving latent confounders is represented by the collider configuration Ix + X + L + Y +- Iy. It needs both and to be pointed to with an arrowhead shown in Figure 8 b). Any nodes with tails will corrupt the conditional dependence between Ix and Y, and Ty and X. Therefore, causal relations and selection processes cannot be directly involved. Only < is allowed, resulting in d-separated paths between X and Y. Then, the latent confounders can be identified, as the CI pattern remains unaffected when these d-separated paths are blocked. Only when cases with inducing paths, such the scenario depicted in Figure d)), cannot be identified. This is because the d-connected paths between and Y mimic the same unique structures associated with latent confounders. However, this d-connection must have latent confounders provide the arrowhead, leading to the identifiability of the existence of latent confounders.
Latent confounders:
The unique structure of the selection process, which needs both X and Y only points out with tails, which are characterized by the paths Ix - X - S and Iy + - S. As discussed in Remark [3.4] this structure exhibits both dependence and conditional dependence between Ix and Y, and between Ty and X. These conditional dependencies need both tails and arrowheads to interact together. However, without selection processes providing vstructure and d-connection, only causal relations and latent confounders will never meet the requirements of interaction of tails and arrowheads for X, Y, and middle nodes that processes forming the inducing path shown in FigureBic were satisfied. Although there is no direct selection work X and Y. they are biased by the selection process between and X, and between Z and Y. Therefore, the existence of selection bias is identified.
Selection process:
block the path. Only cases involving selection
E Experimental setting and results on synthetic dataset
In this section, we show more results of GISI. compared with baselines and evaluate its robustness. The code GISL available at ahre
The accuracy of identifying latent confounders. Beyond verifying GISL's ability to detect selection bias in Table 3 we evaluate its accuracy in identifying latent confounders on synthetic datasets with n=1500, reported in Table On single-cell gene expression datasets, selection bias is operationalized as the pre-post change in proliferation rate summarized by Z-scores (Section 4.3). By contrast, latent confounders are unobserved and lack a direct interventional response or gold standard, precluding systematic real-world evaluation; accordingly, we assess this component exclusively on synthetic benchmarks.
SID evaluation metric The ability of GISL in identifying regulatory relationships is further evaluated by Structural Intervention Distance (SID) as shown in Table 4
19

---

PAG precision
500
1500
2000
0.9
Hard
Soft
0.6
10
15
20
25
10
15
20
25
10
15
20
25
(# Varlables)
Figure 9: Comparison results in identifying regulatory relations under PAG precision. All values are averaged over 10 runs with different random seeds.
Table 3: Accuracy of GISL in identifying latent confounders on synthetic data. We report the mean and variance values of accuracy across 10 independent graphs for each configuration with n=1500.
X
10
15
20
25
Hard intervention
62.5 + 22.8
63.2 + 20.6
68.8 + 15.2
66.7 + 14.7
Soft intervention
66.8 + 25.3
68.4 1 22.9
70.4 + 18.5
71.2 + 15.6
Table Experimental results of GISL evaluated by Structural Intervention Distance (SID) on the synthetic dataset under hard intervention
10
15
20
25
n= 500
4.5 +
1.5
5.1 + 1.2
5.8 +.9
7.7 + 1.0
==1500
4.3 1
1.2
4.8
1
1.9
5.6
1 1.1
7.3 1.3
n=2000
4.2
4.7
5.8 1.0
7.1 0.9
Table 5: Robustness evaluation of GISL in noise level (the mean Gaussian noise μ)) the complexity of structures (# edges), the number of variables |AY, and the realism function. Baseline setting: I U(0 2), edge 15, | 15, hard intervention, and random nonlinear functions.
Settings
Precision
Recall
SHD
F1
Precision
Recall
SHD
n=1500
n=
Baseline ~
69.621.9
80.82-1.
618819
5.5:6.5
61.6=0.9
86.8
50.1
5.8
72.6+1.8
4.7+4.6
67.-118
Ju (4
89.716
69221.11
4.3.3.8
edge
67.622.6
92.820.9
59.42099
7.243.3
90.1
49.4
8.824.7
|X 25
67.3220.
85.82112
56.11:099
10.821.
64.3=0.6
90.81002
50.61111
11.223.6
I 50
71.2+116
90.7+11.
62.4+2211
14.2+668
68.8-14
87.2++0.5
53.7+0.9
16.3+7.0
Normulised Hill differential equations
67.6.11.7
6..1
59.7.13.2
5.9.15.3
1.3
83.7.12.33
48.9
6.4.1.66
Robustness Table 5 evaluatee the robustness of GISL across noise levels, structural/functional complexity of the Structure Causal Model (SCM), and variable counts. Specifically, increasing the district when sampling the mean of Gaussian noise improves the performance, bigger distribution changes are easier for CI test tools to detect. Greater structural and functional complexity causes modest decline in overall performance, whereas holding other factors fixed, sparser graphs yield better performance.
Overall performance in PAG. To represent causal relationships in the presence of latent con- founders and selection bias, we draw inspiration from ancestral graphs and the notations in (55). We employ six types of edges, >, tx, o-o, o-t, represent equivalence classes of conditional
20

---

Table 6: Experimental results of GISL and computational baselines on synthetic data
Methods
Acc
Recall
F1
SHD
GISL. PIDC M
94.7+0.01
95.1+00.01
94.9+00011
1.0:00.54 19.4+12.6
51.6=8.45
95.240
55.1
PPCOR [5
43.4=10.32
97.310.28
49.719.82
26.8:7.66
independence (CI) patterns. Among these, -, ++, -- and corresponds to the structures illustrated in Figure 3 while undirected edges (-) represent there are both latent confounders and selection bias. To assess the ability of GISL in recovering regulatory relationships, selection pro- cesses, and latent confounders, Figure compares the accuracy of the PAG learned by GISL against canonical baselines that accommodate selection bias and latent confounding.
All in all, although the specific structures of causal relationships, latent confounders, and selection bias are indeterminate, every pair of nodes shares the same d-separation properties, exhibiting identical conditional independence relations. This leads to Markov equivalence, which can be uniquely represented by DAG with L and S. Here, the structures involving L and S do not specifically indicate a true common cause or v-structure; rather, they represent the presence of latent confounders selection bias. However, unlike the output of FCI, the output of GISL is more precise. This is because FCI initializes the causal graph with o o, which can represent all types edges. Subsequently, rules applied based the limited information provided by v-structures to determine the specific edge types. For instance, while FCI only utilizes v-structures to infer additional edges, as shown in (b) of Figure where X o Y is output, GISL. refines this X - Y, providing more accurate representation of the equivalence class CI patterns. With CI patterns, the equivalence classes are more precise.
GISL V.S. computational methods under selection bias We rethink the gene regulatory network inference from causal view and focus identifying the causal relationship. latent confounders, and selection process. The setting and output of GISL differ from those of computational methods, which are unable to address the dependencies caused by latent confounders and selection bias. The experimental results of GISL and computational methods on synthetic data are provided for comparison as shown in Table 6 From the table, is evident that computational methods fail to identify causal relations in the presence of selection bias. This failure occurs because the selection process affects not only the directly targeted variables but also those connected through the same causal pathways. For example, (C) in Figure 7 if perturbing Y, the distribution of both Z and X changes. Even without gene perturbation data, computational methods that rely solely on dependence (correlation or co-occurrence) consider these variables as dependent. Consequently, their output often results in a fully connected graph.
Comparison with CDIS 92 CDIS models selection a one-time process that occurs only before intervention and becomes inactive afterward, using an interventional twin graph. Under this assumption, the distribution satisfies P(Y do(X) = x, S) P(Y S), meaning that after an intervention on X is applied, the selection variable S no longer influences Y. In contrast, our work focuses on biological constraints in genes-a form of inherent selection bias that exists both before and after intervention. These constraints come from basic biological rules, such as essential gene functions conditions needed for cell survive, which do not disappear even when gene is perturbed. As a result, selection continues to affect the system even after X.),, violating the CDIS assumption. In our setting, the correct relation is P(Y do(X) x, S) P(Y 5). This fundamental difference means that CDIS cannot account for such persistent selection effects, leading poor performance in settings where biological constraints are present, shown in Table 7
F
Experimental settings and results of real-world dataset
Data availability. With the advent of next- sett-eeeratiinn sequencing (NGS) techniques, such as single- cell RNA-sequencing (scRNA-seq), the availability single-cell data enables deeper analysis of gene expression in biological systems, offering unprecedented resolution at the level of individual cells [93]. Moreover, thanks to the advancement and maturation of gene sequencing and perturbation tools, including CRISPR-Cas 194). CRISPRi 195), and CRISPRa 196), genes transformed
21

---

Figure 10: Comparison of Z-ocrre distribution between genes in selection pairs and others on Table Comparison with CDIS on synthetic dataset with | n=1500, soft intervention. Experi- mental results are evaluated by FI, Precision, Recall, and SHD in regulatory relationships, and by Precision identifying selection bias.
Method
F1
Precision
Recall
SHD
Precision (sc- lection bias)
GISL
80.0 1.7
83.9 + 3.2
74.1 + 1.6
4.1 + 5.7
80.8 + 11.3
CDIS
63.8 + 4.5
61.4 + 4.7
68.1 + 5.3
8.4 + 6.4
34.1 +: 9.1
Z-Score Distribution of Genes by Selection Status
Z-Score Distribution of Genes by Selection Status
2-Score
-1-
Z-Score
-3-
In Selection Pairs
Other Genes
In Selection Pairs
Other Geres
Gene Group
Gene Group
(a)
(b)
Adamson (a) and Norman (b) datasets.
into viable subjects for causal discovery by providing high-qualin single-gene observational and perturbation (interventional) data through the systematic technique Perturb-seq 997 31]31
In real-world datasets, the causal relationships are evaluated based on Enrichr. However, not all perturbed genes are reported in Enrichr, some genes cannot be perturbed processed by biological tools like ChIP-Seq. To evaluate the selection process, z-score is used to verify its existence. The z-score represents the ratio of the growth rate between perturbed genes and normal genes. Changes in growth rates indicate variations in sample size, which align with the characteristics of the selection process. Thus, the z-score serves as an evaluation tool. Distributions of the z-scores for the reported genes are shown in Figure 6 From the figure, we observe that these genes exhibit differences in growth rates between the perturbed and normal cells, indicating the presence of selection process. In some cell lines, the growth rate does not change, suggesting that the gene is not under selection in those cells. This observation is consistent with our explanation regarding differential gene expression.
F1 Experimental results on real-world datasets
We conduct experiments on three representative datasets of single-cell gene expression in the real world to verify the effectiveness of our proposed method, including data from K562 cells [32 61 and Human Lung Epithelial Cells (HLEC) [II. Firstly, the skeleton discovered among 5012 genes. Then, the structure of the perturbed genes is identified. Since GISL outputs equivalence class, the structures represent CI patterns. For example, + denotes equivalence class corresponding to first CI pattern illustrated in Figure The edge indicates the indeterminate edges usually under selection, characterized by CI patterns where I remains conditionally dependent on observed genes, regardless of the conditioning set. The selection variable S and latent confounder specifically indicate the presence selection bias or latent confounders, respectively.
The experimental result on the Dixit dataset is shown in Figure 5 Where blue edges are verified as correct by Enrichr, where the ChEA 2022 198), TF-gene enrichment [99], and (100) are mainly considered. We consider them to be more reliable because ChIP-Seq, used in ChEA, directly detects the binding sites between transcription factors and genes. Additionally, TF-gene enrichment identifies robust pairs supported by enrichment analysis. Moreover, the selection bias is supported by the high 7-scores observed for the genes such as RACGAPI (-0.903), F2F4 (-0. 158), GABPA (-0.543), and N R2C2 (-0.127), respectively. Experimental results indicate that almost all directed causal
22

---

relationships discovered by GISL are correct, while undirected edges remain undetermined. This is because undirected edges suggest the presence of both selection bias and latent confounders, making the causal relationship unidentifiable. This is because, with selection bias and latent confounders, the CI pattern becomes fully dependent, which obscures the unique characteristics of causal relationships in the CI pattern. The experimental results of PPOCR are shown in Figure T2 The results show that our proposed GISL demonstrates greater accuracy and reliability in matching the results of biological experiments. In addition, the experimental results of GISL on the Norman and Adamson datasets are shown in Figure TT and Figure [3] Correspondingly, the Z-score distribution comparison between genes under selection bias and others is illustrated in Figure 10
Figure 11: Experimental result of GISL on the Norman dataset.
TORTAIPI
NIGCC
IRFI
CENIEE
ELFi
RACGAPI
E2F4
CRESH
YY!
Figure 12: Experimental result of PPCOR the Dixit dataset.
Figure 13: Experimental result GISL on the Adamson dataset.
Moreover, biological analysis for the human lung epithelial and k562 lukemia cell lines is discussed as follows:
Human Lung Epithelial Cells GISL reports that there exist biological constraints between gene
CNNT and CDKNIA. CNNI (Calponin 1) is normally low in healthy lung epithelium but markedly increases during metthelialmeseneehnn transition fibrotic remodeling [00]. CDKNIA (p21) is
23

---

cell-cycle inhibitor that surges under DNA damage stress (e.g., p53 activation), causing growth arrest or senescence; indeed, lung injury models and fibrotic lungs (IPF) show abnormally high p21 in alveolar epithelial cells [02 These observations imply that when lung epithelial cells are pushed toward damaged or transitioning state, both CNN1 and p21 tend rise, which suggests that the cell may need to keep its combined activity in check to maintain normal function. Excessive co-elevation of CNNI (indicating mesenchymal/fibrotic shift) and CDKNIA A (indicating cell-cycle arrest) could drive cells into an unsustainable state (senescence or fibrosis), so a balance in their expression is likely required to preserve viability and tissue homeostasis. This suggests form biological constraint whereby lung epithelial cells limit the concurrent upregulation of CNN1 and CDKNIA to avoid tipping into pathological remodeling or loss of proliferative capacity.
K562 Lukemia Cells GISL reports that there exist biological constraints between gene ELF1 and gene GABPA. ELFI, an ETS-family transcription factor, regulates hematopoictic differentiation and immune genes; in K562 cells, it controls the G1 to S phase transition via CDKNIA (p21), with overexpression inducing apoptosis and underexpression promoting unchecked proliferation []03]. GABPA, nuclear respiratory factor, orchestrates mitochondrial biogenesis and metabolism through the P13K/Ak//mTOR axis; its depletion triggers GO to GI phase arrest, and its binding affinity (pKD) and phosphorylation response (pECsn) underscore requirement for dosage matching downstream signaling intensity [I04]. Both factors converge on stress and apoptotic pathways, GABPA targets (Caspase-9, Bcl-2) overlap with ELF1, regulated DNA-damage genes, and respond coordinately to P13K/Akt perturbations (e.g.. LY294002) and autophagy induction (e.g.. in K562 and ADM cells), suggesting K562 cells impose a constraint ELF1 and GABPA expression to maintain proliferative homeostasis and survival [06)
G Understanding Z-score
Z-score is result of a model proposed by (3), which is designed to describe the dynamic proliferation process. For single-guided RNAS (SgRNA) j targeting gene in cell line c, the number of cells (t) with sgRNA at time t after perturbation is modeled as follows:
(o)(1) = (o) (DejeRest + (1- -l-Pec Paj) Ret)
(2)
where = 0 is the time of perturbation, Pej is the probability that the sgRNA j achieves knockout of its target in cell line c, R. is the unperturbed growth rate the cell line, and R is the new growth rate caused by knockout of the targeted gene in the given cell line. The gene fitness effect is the fractional change in growth rate Teg = R R. - 1. Gene fitness effects are the primary desired output for this type of experiment.
Considering the delayed perturbation effect, the model built in practice to estimate the number of cells are as follows:
L De (t d = =0)((0) (1 + p.P, ('P'
( - -1) (t)
where
e indexes cell line, j indexes sgRNA, g indexes gene, T. indexes library or batch, and is the time elapsed since library transduction
u (t) is the model estimate the normalized readcounts of sgRNA j in cell line e screened in batch or library T. at time f
w (0) is the model estimate of the normalized number of cells initially receiving sgRNA j
_ and pj are the estimated CRISPR knockout efficacies in cell line e with sgRNA j
R is the estimated unperturbed growth rate of the cell line
Teg is the estimated relative change in growth rate for that cell line if gene g targeted by sgRNA j is completely knocked out
dg is the delay between infection and the onset of the growth phenotype
z (L) is normalization equal the sum of the numerator over all sgRNAs j the cell
line for the given library and time point
24

---

(a)
(b)
(c)
Figure 14: The distribution of z-scores for the genes FEV, PTPNI, and SET across all cell lines highlights the variation in sample size.
The Z-score describing the ratio of growth rate before and after perturbation is designed as follows:
(1) (e)(0) (0) PoPj
(e (Reree(t-d))) -1)) 14 2 dg
(4)
Zez (t) Uej (0) vt < dg.
(5)
Using FEC, PTPN1, and SET as representative genes, the calculated Z-scores across a panel of cell lines are shown in Figure 14 Cell lines with a Z-score of zero indicate the absence of selection constraints on the genes: perturbation does not induce measurable cell death, and the expression distribution remains unchanged. In contrast, non-zero Z-scores mark genes under selection con- straints, where perturbation produces pronounced shifts in their expression profiles. Moreover, each gene displays distinct 7-score pattern across the cell-line panel, reflecting gene-specific selective dynamics. These results align with the differential gene-expression patterns expected under selective pressure.
H Discussion
Discussion and Limitations. In cells, we argue that the different intracellular environments, acting as selection mechanisms, restrict the expression of genes. When the environment remains, a selection mechanism is always present. Genes stay in cells with the remaining environment, showing the reasonability of our setting. The superiority of GISL is that it can identify the existence of selection bias and latent confounders, as well as causal relationships. Despite its flexibility in integrating obser- vational and perturbation data within unified causal framework, a fundamental limitation remains: whenever selection bias is present, spurious marginal and conditional dependencies inevitably arise (see Section 3.4). Consequently, reliably detecting selection in the presence such bias remains an outstanding challenge.
More Discussion.
Q1: Why does the paper explicitly look at gene regulatory networks when the method is more generally for causal discovery with interventional data?
A1: GISL is broadly applicable to general causal discovery, but chose the gene regulatory network inference (GRNI) framework because it was inspired by a real-world challenge, the pervasive yet often overlooked selection (survival) bias in biological data, and directly addressing this issue offers immediate value for GRNI and beyond.
Specifically, we focus on GRNI because:
GRN inference is classical causal discovery task with strong biological significance, supported by well-established gene perturbation technologies that provide abundant inter- ventional data.
Despite being a classical problem, modeling gene expression from causal perspective remains challenging due to latent confounding and overlooked selection bias, making it a meaningful testbed for our method.
Q2: Why GISL performs the best regarding the recall? A2: As shown in Figure 4 GISL shows a little bit lower recall in identifying regulatory relationships under selection bias. This is because, theory, is challenging to distinguish (c): selection bias and (e): causal selection bias, as they share the same CI patterns (all conditionally dependent). As result, our GISL categorizes both cases
25

---

(c) and (e) as selection bias and thus shows lower recall when identifying regulatory relationships (causal relations).
On the other hand, the baseline algorithm, e.g. GIES, exhibits higher recall, but this does not necessarily indicate better performance. The key issue is that GIES is unable to account for selection bias, which introduces spurious dependencies into the data. When GIES relies on score function to detect conditional dependencies, it may mistakenly interpret these spurious associations as causal relationships. As a result, the algorithm tends to include additional edges regardless of whether true causal relationships exist, leading inflated recall but low accuracy.
26

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Our work pinpoints selection bias as pervasive, yet previously ignored, cause of spurious gene-gene dependencies in expression data. Leveraging the insight that selection pairs respond symmetrically to gene perturbations (unlike asymmetric regulation untouched confounding), we devise GISL, non parametric causal-discovery algorithm that jointly recovers true regulatory edges, latent confounders, and selection pairs. The proposed GISL algorithm is assessed on both synethic data and real-world gene data, with superiority over classical causal discovery methods and computational GRNI baselines. Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We have a section named Discussion in Appendix where we present the limitation of our method.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "I "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
27

---

3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
The main theorical result is the identifiability (Section 3.3 of the proposed GISL framework. We attach its proof in Appendix ??. As for the two implicit assumptions
Justification:
(Assumption 3.1 and (.2), both of them are common and fundametal assumptions for connecting causality with statistical tools.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proof' can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The most important contribution of this work is the proposed GISL framework, which we present the algorithmic designs using a detailed pesudocode (Algorithm IT. The experimental results have folds: synethic data and real-world data. For synethic data experiments, we detailed the data generation, hyperameters and reported our results with sufficient number runs. For results on real-world data, we detailed the dataset source and the rational and novel evaluation protocol (using Z-score as the proxy of ground-truth labels).
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
28

---

(c) If the contribution is a new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: We clearly stated the dataset (both synethic and real-world data) used this work. As for the code, we presented pesucode for reproduce the main experimental results. While we are not able to make public the code, we will realease it after the paper reviewing.
Guidelines:
The answer NA means that paper does not include experiments requiring code. Please see the NeurIPS code and data submission guidelines (https://nips) cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
ssiisnoliccy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification; We inleuded the generation of synethic data before introducing our experiment results as well as the souree of the real-world datset used in our analysis. Our method didn't t include any hyperameters.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
29

---

Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: For the table results, report the mean value and the variance of evaluation metrics across runs. For the figure comparisons, we plot the mean value and the 95% confidence interval of the evaluation metrisc across runs.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2- 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative eтror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Runing our GISL on synethic dataset, e.g. 2000 samples with 10 variables, takes about 10 minutes on single CPU core.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics Ethics/httss:/ httppp/:/paarrrss..coopapa/aataa
Answer: [Yes]
Justification: Our work focuses on uncovering selection bias and causal regulation in gene expression data. We believe it does not violate any ethical guidelines.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
30

---

If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Our work is foundational in research and not tied to particular applications. While there might be many societal consequences of our work, none which we feel must be specifically highlighted here.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We do not foresee any significant risk of misuse arising from our work.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
31

---

Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We explicitly mentioned the sources of datasets used in our experiment.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, wwwpppppppsataeea thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if has changed) should be provided. If this information is available online, the authors are encouraged to reach out the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
This work did not release any new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
32

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy //np:s:////urraie cc. Cc//oonferencee/22555///// for what should or should not be described.
33


Reviews:
[{"review_id": "StxMq015cX", "reviewer": "Reviewer_TNYS", "summary": "This paper is concerned with inferring latent confounders and selection bias in gene regulatory networks with interventional data. The paper shows that interventional data implies different d-separations in the graphical models for direct cause, existence of latent confounder, and existence of selection bias. Next, an algorithm, GISL, is introduced that uses both observational and interventional data to infer the causal graph as well as the existence of confounders and selection bias. The author's test this method in synthetic settings, where they show that it can identify selection bias. They then test the method on synthetic graphs that have selection bias, and on three gene perturbation datasets.", "questions": "- Why are the experiments only done with selection bias and not with latent confounders as well? Most of the paper talks about identifiability of both. \n- It is hard to interpret the results of 4.1 without a proper baseline. Maybe a baseline where the selection criteria is observed (so it reduces to finding a collider) may provide an upper bound on the performance?\n- There is no mention of the actual conditional independence test used or other experimental details in the Appendix.", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths:**\n- The idea of the paper is simple yet effective.\n\n**Weaknesses:**\n- I'm unsure why the paper explicitly looks at gene regulatory networks when the method is more generally for causal discovery with interventional data. I believe this may reduce the perceived significance of the work.\n- The related work should be in the main paper, and not in the appendix. It should also be emphasised in the related work whether methods require only observational data, or also interventional data (like the current work). A lot of the criticisms of related work, that they require functional assumptions is because they do not require interventional data. The related work could also be more complete. The novelty and differences of this work in comparison to other works is not made very clear.\n- Some of the claims made in the paper could be more accurate. It seems that the algorithm only allows for identifying latent confounders with interventions on all variables (this statement is missing in Theorem 3.7). Even in this case, there are certain cases that are not identifiable (like figure 7(d) as mentioned in Appendix D).\n- The claims about unobserved confounder are not backed by experiments (see question).\n- The SID (structural intervention distance) metric might also be a more natural measure of how good the recovered graph is (as opposed to just SHD).\n\n\n**Comments:**\n- The figures need to be bigger to be legible along with the rest of the text. Figure 2 is not decipherable at all.\n- Figure 3 is quite hard to understand and it took me a while to figure out what the shapes were. I'm also fairly sure 3 (b) top is incorrect.\n- L234: Where is example 3.2? Do you mean example 2?\n- L271: It is claimed that you are using nonparametric SCMs for complexity, but the functions chosen are actually quite simple. Would functions drawn from a GP or an NN not be more appropriate?", "was_revised": true, "final_justification": "Though the presentation of the work can be improved (which is partly addressed by the authors), I do believe taking selection bias into account is an overlooked area.\n\nAlthough I agree with the concerns raised by yzJf, they have been partly addressed in my and their responses. I also think these concerns are generally valid for almost all causal discovery methods.", "rebuttal": "We sincerely thank you for the time dedicated to reviewing our paper and valuable feedback. Please find the response to your comments and questions below.\n\n---\n**Q1:** **Why does the paper explicitly look at gene regulatory networks when the method is more generally for causal discovery with interventional data?**\n\nThank you for this insightful suggestion. Indeed, our method is broadly applicable to general causal discovery, but we chose the gene regulatory network inference (GRNI) framework because it was inspired by a real‐world challenge—the pervasive yet often overlooked selection (survival) bias in biological data—and directly addressing this issue offers immediate value for GRNI and beyond. \n\nSpecifically, we focus on GRNI because:\n1. GRN inference is a classical causal discovery task with strong biological significance, supported by well-established gene perturbation technologies that provide abundant interventional data.\n2.  Despite being a classical problem, modeling gene expression from a causal perspective remains challenging due to latent confounding and overlooked selection bias, making it a meaningful testbed for our method.\n\n---\n**Q2:** **Related works should be in the main paper and could be more complete. The difference between GISL and related works should be clear.**\n\nIn the revised version, **we will move the discussion of related work into the main paper and expand it accordingly.** Due to limited space, we list representative causal discovery methods for comparison in settings and highlight the novelty and difference compared with related works as follows:\n\n| Mehtod | Require Interventional data|Allow latent confounder |Allow selection bias | nonparametric assumption |\n|----------|----------|----------|----------|----------|\n| [1] |✅| ❌  | ❌ | ✅| \n| [2] | ❌ | ✅| ✅|✅| \n| [3] | ❌ | ✅| ❌| ❌|\n| [4] | ✅ | ✅|❌|✅ |\n| [5] | ❌ |❌|  ✅|❌|\n| Ours    |✅|✅|✅|✅|\n\nIn comparison, **GISL is novel in utilizing interventional data to identify relationships, latent confounders, and selection bias in a general setting without relying on parametric assumptions**, and is specifically tailored to handle the **biological selection constraints** that arise in gene regulatory networks. \n\n> [1]. Permutation-Based Causal Structure Learning with Unknown Intervention Targets. UAI 2020.\n\n> [2]. Iterative causal discovery in the possible presence of latent confounders and selection bias. NeurIPS 2021.\n\n> [3]. A Versatile Causal Discovery Framework to Allow Causally-Related Hidden Variables. ICLR 2023.\n\n> [4]. Characterization and Learning of Causal Graphs with Latent Variables from Soft Interventions. NeurIPS 2019.\n\n> [5]. Identifying selection bias from observational data. AAAI 2018.\n\n---\n**Q3:** **The statements about the identifiability in Theorem 3.7 could be more accurate. There are certain cases that are not identifiable (like figure 7(d) as mentioned in Appendix D).**\n\nWe agree with you. We have revised the statement in Theorem 3.7 to make it more accurate as follows: \n\nThe causal relationships, selection processes, and latent confounders are identifiable up to the equivalence classes of four types of CI patterns in Figure 3 **among the variables that are subject to intervention.** The presence of selection and latent confounders (existing or not) is identifiable among these variables.\n\nAdditionally, as shown in Figure 7(d) in Appendix D, **certain structures such as $X \\rightarrow L \\leftarrow Y$ differ in causal structure but yield the same Conditional Independence (CI) patterns, placing them in the same Markov equivalence class under PAG representation**. While the exact structure may not be distinguishable, **the presence of latent confounders between X and Y can still be inferred**. Thus, GISL can reliably identify whether selection bias or latent confounders exist among the intervened variables and identify causal structures up to the equivalence class.\n\n---\n**Q4:** **Lacking the accuracy of identifying latent confounders.**\n\n We report GISL’s ability to identify latent confounders on **synthetic datasets** with 1500 samples as follows:\n\n| # Variables| 10| 15 |20| 25|\n|----------|----------|----------|----------|----------|\n| hard| 62.5 $\\pm$ 22.8  | 63.2 $\\pm$ 20.6  | 68.8 $\\pm$ 15.2  | 66.7 $\\pm$ 14.7 | \n| soft  |   66.8 $\\pm$ 25.3   |   68.4 $\\pm$ 22.9   | 70.4 $\\pm$ 18.5    |    71.2 $\\pm$ 15.6  |\n\nWe focus on selection bias in our experiments for the following reasons:\n\n1.  **One of the main contributions of this work is to highlight the often-overlooked role of selection bias in gene expression data**, and to propose GISL specifically to identify it.     \n2.  **Selection bias can be systematically evaluated** on both synthetic and real-world datasets, as the selection processes (e.g., biological constraints) are either known or can be reasonably assumed. In contrast, **latent confounders can only be evaluated on synthetic datasets**, since there is no ground truth or systematic benchmark for latent confounding in real-world gene expression data.\n\n---\n**Q5:** **The SID (structural intervention distance) metric might also be a more natural measure of how good the recovered graph is.**\n\nThank you for the suggestion. GISL focuses on recovering the **graphical structure**, and **SHD**, which directly measures structural differences globally, is a commonly used and standard evaluation metric in this context. While **SID** provides a meaningful assessment of **local interventional reasoning**, it focuses on **ancestral relationships** and is generally more tolerant of structural discrepancies. As a result, SID may overlook structural inaccuracies that do not impact interventional distributions.\n\nNevertheless, we appreciate your comment and are happy to incorporate it for a more comprehensive analysis. According to your suggestion, we evaluate the SID under varying sample size n for hard intervention as follows:\n| # Variables| 10| 15 |20| 25|\n|----------|----------|----------|----------|----------|\n| n= 500| 4.5 $\\pm$ 1.5  | 5.1 $\\pm$ 1.2  | 5.8 $\\pm$ 0.9  | 7.7 $\\pm$ 1.0 | \n| n=1500  |   4.3 $\\pm$ 1.2   |    4.8 $\\pm$ 1.9   | 5.6 $\\pm$ 1.1     |   7.3 $\\pm$ 1.3  |\n| n=2000  |    4.2 $\\pm$ 1.1  |    4.7 $\\pm$ 1.4   | 5.8 $\\pm$ 1.0     |    7.1 $\\pm$ 0.9  |\n\n---\n**Q6:** **Figure 2 is not decipherable.**\n\nWe will revise and reformat Figure 2 to improve its clarity and readability in the final version of the paper.\n\n---\n**Q7:** **Figure 3 is quite hard to understand, and (b) top is incorrect.**\n\nWe thank you for your patience in interpreting Figure 3. We have added a detailed description for better clarity and corrected the mistake in panel (b). Additionally, we have added a detailed explanation to help readers interpret the figure more easily.\n\nSpecifically, Figure 3 (f) serves as a **reference table** summarizing the CI patterns for each target gene pair. In this table, **black symbols** indicate that the corresponding CI relation **holds** (i.e., the variables are conditionally independent, while **white symbols** indicate that the CI relation **does not hold** (i.e., the variables are conditionally dependent).\n\n**Taking Figure 3(a) as an example,**  the structure exhibits conditional dependence $Y \\not\\perp I_X \\mid S$, which is represented by the white upward triangles (△) in the reference table. The conditional independence $Y \\perp I_X \\mid X, S$ is shown in (a) when $X$ is given and is represented by black downward triangles (▼). These CI patterns for (a) are summarized using the corresponding symbols displayed below the causal graph.\n\n---\n**Q8:**  **Incorrect reference to example 2.**\n\nThanks for pointing this out. We have corrected it. \n\n---\n**Q9:** **L271: The function in nonparametric simulation is simple. Would functions drawn from a GP or an NN not be more appropriate?**\n\nAccording to your suggestion, we have expanded the function pool to include a **neural network (NN)** with a hidden dimension of [5, 5] and ReLU activation. The regulatory function is now randomly sampled from a set that includes linear, quadratic, sin, log, and NN functions. The results of number of variables $\\mathcal{X}$ = 15, number of samples n = 1500  for soft intervention are as follows:\n\n| Method | F1| Precision |Recall| SHD | SID|\n|----------|----------|----------|----------|----------|----------|\n| Without NN | 80.0 $\\pm$ 1.7  | 83.9 $\\pm$ 3.2  | 74.1 $\\pm$ 1.6 | 4.1 $\\pm$ 5.7| 3.6 $\\pm$ 1.2| \n| With NN    |  77.1 $\\pm$ 1.5     |  82.4 $\\pm$ 2.5  |  73.1 $\\pm$ 1.6  |   4.4 $\\pm$ 6.4   | 4.1 $\\pm$ 1.1 |\n\n---\n**Q10:** **Baselines for identifying colliders as upper bound performance.**\n\nFollowing your advice, we conducted experiments on synthetic data comprising 15 variables and 2000 samples under soft interventions to evaluate the identifiability of GISL in detecting colliders, as detailed below:\n\n| # | 1 | 2 | 3 | 4 | \n|----------|----------|----------|----------|----------|\n| Selection process | 78.5 $\\pm$ 15.6  | 76.1 $\\pm$ 14.5  | 70.6 $\\pm$ 17.7 | 69.5 $\\pm$ 14.3| \n| Collider    |  90.0 $\\pm$ 5.0     |  85.0 $\\pm$ 5.3  |  83.3 $\\pm$ 5.6  |   82.5 $\\pm$ 3.9   |\n\n---\n**Q11:** **Lacking the actual conditional independence test used or other experimental details.**\n\n We have added detailed descriptions of the experimental setup as follows:\nSpecifically, following Algorithm 1, the first step—**skeleton discovery**—is performed using **FGES** with the **BIC score** to identify the conditional dependence. In Step 3, we capture CI patterns from both observational and interventional data. For the **nonparametric CI testing**, we use the **Kernel-based Conditional Independence (KCI)** test.\n\n---\nIf any aspect remains unclear, please do not hesitate to let us know; we would be more than happy to discuss it further with you."}, {"review_id": "2whr7QmcKx", "reviewer": "Reviewer_yzJf", "summary": "This paper introduces GISL, a nonparametric causal discovery algorithm tailored for gene regulatory network inference (GRNI) in the presence of both selection bias and latent confounders. The main innovation hinges on the key insight that *selection bias* can create spurious dependencies in gene expression data, particularly in contexts where only cells meeting certain survival thresholds might be observed in lab settings. GISL is a perturbation-inspired framework to distinguish selection-induced dependencies from true regulation or confounding. GISL is operationalized as a flexible causal framework that leverages conditional independence to differentiate between regulatory, confounded, and selection-based relationships, with a well formulated algorithm (Algo 1). The method is validated on both synthetic and real single-cell gene expression datasets (scRNA-seq perturbation datasets), and is compared to competing methodologies.", "questions": "The three main concerns I have, which I believe could improve the paper, are:\n\n**Robustness of Simulation Experiments:**\nThe simulation study could benefit from a more thorough sensitivity analysis. Methods like GISL are often highly sensitive to both noise levels and the dimensionality of the system—particularly in high-gene-count settings typical in gene expression studies. However, the current simulations appear to focus only on relatively small systems (10–25 variables). It's unclear whether GISL's performance generalizes to more realistic, higher-dimensional scenarios, or more complicated functional forms of the SCM. A systematic evaluation of performance across varying noise levels, SCM complexities, and number of genes would help validate the method’s robustness.\n\n**Realism of the Synthetic Data Generator:**\nWhile SCMs offer a clean framework for simulation, it's not entirely clear how well they reflect the complexities of real gene expression dynamics. More biologically grounded simulators—such as the one introduced in Bhuva et al. (“Differential co-expression-based detection of conditional relationships...”)—could provide a more realistic benchmark. If this is too difficult, then at least some context or justification for using SCMs over existing biological simulators would be helpful, especially if the goal is to evaluate performance in settings that mirror actual gene regulation.\n\n**Need for More Real-World Evidence:**\nGiven the strong real-world motivation throughout the paper—for instance, the AURKA–TOR1AIP1 example in the introduction—there’s definitely room to further ground the method in real world biological data. One way to improve this would be to shift some theoretical content to the appendix and expand the real data application. Specifically, the paper could benefit from showcasing concrete examples of relationships that GISL correctly flags as arising from selection bias in leukemia and lung cancer—ideally, in pathways with known biological validation. On the flip side, it would also be valuable to highlight where GISL fails or gives ambiguous results. Both successes and limitations will help biologists understand the tangible benefits of GISL.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths: \n- problem is well motivated and biologically meaningful\n- well formulated (and clearly explained) algorithm\n- interesting simulation experiments\n- strong biological application\n\nWeaknesses (elaborated below): \n- simulation experiments can be made more robust\n- need more real world evidence\n\nNOTE: authors performed extensive additional experiments to address my concerns listed above.", "was_revised": true, "final_justification": "Following the authors' rebuttal and in light of the additional experiments and clarifications, I have revised my evaluation. The authors have addressed my primary concerns by demonstrating the robustness of GISL across multiple dimensions, including noise levels, structural complexity, and gene counts. These results provide stronger evidence for the generalizability of the method. Furthermore, the inclusion of real-world biological evidence enhances the relevance and interpretability of the approach, particularly for biologically focused audiences.", "rebuttal": "We first sincerely appreciate your insightful suggestions and valuable feedback, which have enabled us to strengthen our method through a more comprehensive evaluation. We have carefully considered the concerns you raised and provide our detailed responses below.\n\n---\n**Q1:** **Robustness of Simulation Experiments.** \n\nIn light of your suggestion, we have conducted comprehensive experiments to evaluate the robustness of GISL as follows:\n\n**Noise level**\n\nWe conducted experiments to evaluate the robustness of GISL under varying levels of Gaussian noise variance, where the noise mean was randomly sampled from a uniform distribution U(0, 4). **Experimental results demonstrate that GISL’s performance improves as Gaussian noise variance increases, because larger distributional changes are more easily detected by conditional independence (CI) tests, yielding more accurate inference.**\n\nExperimental setting: hard intervention, 15 variables, 1500 samples.\n| Variance of noise | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| U(0,2) | 69.6 $\\pm$ 1.9  | 80.8 $\\pm$ 4.3  | 61.8 $\\pm$ 1.9 | 5.5 $\\pm$ 6.5|\n| U(0,4)    |  72.6 $\\pm$ 1.8     |  93.5 $\\pm$ 1.7  |  62.1 $\\pm$ 3.2  |   4.7 $\\pm$ 4.6   |\n| U(2,4)    |  77.4 $\\pm$ 0.8     |  89.7 $\\pm$ 1.4  |  69.2 $\\pm$ 1.1  |   4.3 $\\pm$ 3.8   |\n\nExperimental setting: hard intervention, 15 variables, 500 samples.\n\n| Variance of noise | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| U(0,2)    |  61.6 $\\pm$ 0.9  |  86.8 $\\pm$ 0.2  |  50.1 $\\pm$ 1.7  |   5.8 $\\pm$ 1.4   |\n| U(0,4) | 67.4 $\\pm$ 1.8  | 87.2 $\\pm$ 2.9  | 56.1 $\\pm$ 1.7 | 5.7 $\\pm$ 6.6|\n| U(2,4)    |  67.9 $\\pm$ 1.0     |  88.8 $\\pm$ 1.3  |  56.1 $\\pm$ 1.5  |   5.6 $\\pm$ 4.0   |\n\n**SCM complexity**\n\nWe conducted experiments to evaluate the robustness of GISL under varying levels of structural complexity, quantified by the number of edges in the graph. **Experimental results indicate that GISL’s performance experiences a slight decline as graph structural complexity increases.**\n\nExperimental setting: hard intervention, 15 variables, 1500 samples.\n| Number of edges | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| 15 | 69.6 $\\pm$ 1.9  | 80.8 $\\pm$ 4.3  | 61.8 $\\pm$ 1.9 | 5.5 $\\pm$ 6.5|\n| 20    |  67.6 $\\pm$ 0.6     |  92.8 $\\pm$ 0.9  |  59.4 $\\pm$ 0.9  |   7.2 $\\pm$ 3.3   |\n\nExperimental setting: hard intervention, 15 variables, 500 samples.\n\n| Number of edges | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| 15    |  61.6 $\\pm$ 0.9  |  86.8 $\\pm$ 0.2  |  50.1 $\\pm$ 1.7  |   5.8 $\\pm$ 1.4   |\n| 20 | 59.2 $\\pm$ 1.3  | 90.1 $\\pm$ 1.2  | 49.4 $\\pm$ 1.5 | 8.8 $\\pm$ 4.7|\n\n**The number of genes**\n\nWe conducted experiments to evaluate the robustness of GISL under a varying number of genes. **Experimental results indicate that GISL’s performance remains stable and even improves as the gene count increases, owing to the greater sparsity of the corresponding graphs.**\n\nExperimental setting: hard intervention, 1500 samples.\n| Number of genes | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| 15 | 69.6 $\\pm$ 1.9  | 80.8 $\\pm$ 4.3  | 61.8 $\\pm$ 1.9 | 5.5 $\\pm$ 6.5|\n| 25    |  67.3 $\\pm$ 0.3     |  85.8 $\\pm$ 1.2  |  56.1 $\\pm$ 0.9  |   10.8 $\\pm$ 1.9   |\n| 50    |  71.2 $\\pm$ 1.6     |  90.7 $\\pm$ 1.3  |  62.4 $\\pm$ 2.1  |   14.2 $\\pm$ 6.8   |\n\nExperimental setting: hard intervention, 15 variables, 500 samples.\n\n| Number of edges | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| 15    |  61.6 $\\pm$ 0.9  |  86.8 $\\pm$ 0.2  |  50.1 $\\pm$ 1.7  |   5.8 $\\pm$ 1.4   |\n| 25 | 64.3 $\\pm$ 0.6  | 90.8 $\\pm$ 0.2  | 50.6 $\\pm$ 1.1 | 11.2 $\\pm$ 3.6|\n| 50 | 68.8 $\\pm$ 1.4  | 87.2 $\\pm$ 0.5  | 53.7 $\\pm$ 0.9 | 16.3 $\\pm$ 7.6|\n\nFor larger datasets, please see Figures 11 and 13, which present real-world experiments on gene expression data encompassing over **5,000 observed genes** and **more than 100 perturbed genes**.\n\n---\n**Q2:** **Realism of the Synthetic Data Generator.** \n\nIn light of your suggestions, we consider biological simulators, which simulate the regulatory relationship by using the normalised-Hill differential equations and gene expression is a truncated distribution. Experimental results are as follows:\n\nExperimental setting: hard intervention, 15 variables, 1500 samples.\n| Regulatory function | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| Random selected nonlinear | 69.6 $\\pm$ 1.9  | 80.8 $\\pm$ 4.3  | 61.8 $\\pm$ 1.9 | 5.5 $\\pm$ 6.5|\n|  Normalised-Hill differential equations   |  67.6 $\\pm$ 1.7     |  83.1 $\\pm$ 2.6  |  59.7 $\\pm$ 3.2  |   5.9 $\\pm$ 5.3   |\n\n\nExperimental setting: hard intervention, 15 variables, 500 samples.\n\n| Regulatory function | F1| Precision |Recall| SHD |\n|----------|----------|----------|----------|----------|\n| Random selected nonlinear   |  61.6 $\\pm$ 0.9  |  86.8 $\\pm$ 0.2  |  50.1 $\\pm$ 1.7  |   5.8 $\\pm$ 1.4   |\n| Normalised-Hill differential equations | 59.4 $\\pm$ 1.3  | 85.7 $\\pm$ 2.3  | 48.9 $\\pm$ 2.1 | 6.4 $\\pm$ 3.6|\n\nExperimental results demonstrate stable performance on data generated by SCM-driven biological simulators.\n\n---\n**Q3:** **Need for More Real-World Evidence.**\n\n According to your suggestion, we have added more real-world evidence to support the results of GISL as follows:  \n\n**Human Lung Epithelial Cells**\n\nGISL reports that there exist biological constraints between gene CNN1 and CDKN1A. **CNN1** (Calponin 1) is normally low in healthy lung epithelium but markedly increases during epithelial–mesenchymal transition or fibrotic remodeling. [1] **CDKN1A** (_p21_) is a cell-cycle inhibitor that surges under DNA damage or stress (e.g., p53 activation), causing growth arrest or senescence; indeed, lung injury models and fibrotic lungs (IPF) show abnormally high p21 in alveolar epithelial cells [2]. These observations imply that when lung epithelial cells are pushed toward a damaged or transitioning state, both CNN1 and p21 tend to rise, which suggests that the cell may need to keep its combined activity in check to maintain normal function. Excessive co-elevation of CNN1 (indicating mesenchymal/fibrotic shift) and CDKN1A (indicating cell-cycle arrest) could drive cells into an unsustainable state (senescence or fibrosis), so a balance in their expression is likely required to preserve viability and tissue homeostasis. This suggests a form of _biological constraint_ whereby lung epithelial cells limit the concurrent upregulation of CNN1 and CDKN1A to avoid tipping into pathological remodeling or loss of proliferative capacity.\n\n**K562 Lukemia Cells**\n\nGISL reports that there exist biological constraints between gene ELF1 and gene GABPA. **ELF1**, an ETS-family transcription factor, regulates hematopoietic differentiation and immune genes; in K562 cells, it controls the G₁/S transition via CDKN1A (p21), with overexpression inducing apoptosis and underexpression promoting unchecked proliferation [3]. **GABPA**, a nuclear respiratory factor, orchestrates mitochondrial biogenesis and metabolism through the PI3K/Akt/mTOR axis; its depletion triggers G₀/G₁ arrest, and its binding affinity (pKD) and phosphorylation response (pEC₅₀) underscore a requirement for dosage matching downstream signaling intensity [4]. Both factors converge on stress and apoptotic pathways—GABPA targets (Caspase-9, Bcl-2) overlap with ELF1-regulated DNA-damage genes—and respond coordinately to PI3K/Akt perturbations (e.g., LY294002) and autophagy induction (e.g., in K562/ADM cells), suggesting K562 cells impose a constraint on ELF1 and GABPA expression to maintain proliferative homeostasis and survival [5,6].\n\n> [1] Alveolar epithelial cells express mesenchymal proteins in patients with idiopathic pulmonary fibrosis. American Journal of Physiology-Lung Cellular and Molecular Physiology 2011.\n\n> [2] P21Waf1/Cip1/Sdi1 and p53 expression in association with DNA strand breaks in idiopathic pulmonary fibrosis. American journal of respiratory and critical care medicine 1996.\n\n> [3] SIRT1 is a critical regulator of K562 cell growth, survival, and differentiation. Expression Cell Research 2016.\n\n> [4] Decrypting drug actions and protein modifications by dose-and time-resolved proteomics. Science  2023.\n\n> [5] A genetic map of the response to DNA damage in human cells, Cell  2020.\n\n> [6] Elf1 promotes transcription-coupled repair in yeast by using its C-terminal domain to bind TFIIH. Nature Communication 2024."}, {"review_id": "uXEnyznOWj", "reviewer": "Reviewer_7UVK", "summary": "This paper proposes a gene regulatory network inference algorithm to learn the causal regulatory relationship from gene expression data. The gene expression data involves both hard and soft interventions. The paper examines a generalized intervention setting that incorporates the presence of a latent confounder and selection bias. The main idea is to use the asymmetry of intervention with selection bias to indicate the presence of a regulatory relationship.", "questions": "- I don't fully understand the symbols in Figure 3. Could you please explain that in more detail with the reference table?\n- This is related to the previous question. What are the terms like $[I_i,X_j|S]$? Are those indicator variables for the conditional independence?\n- How scalable is the proposed method? Since the proposed method involves conditional independence testing, how sensitive is the proposed method to the size of datasets?\n\n\nSome minor details:\n- The option in Algorithm 1 $\\mathcal{J}==(c)$ and $\\mathcal{J}==(d)$ seems to be flipped and inconsistent with the caption of Figure 3.\n- Do the terms regulatory relations and causal relations mean the same thing in this context? The paper seems to use these terms interchangeably.", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "# Strengths\n\n- The paper studies an interesting and relevant problem.\n- The proposed method is evaluated with both synthetic and real-world data experiments. \n\n# Weaknesses\n\n- Some of the notation and explanation in the paper seems a bit confusing (please see questions below). \n- The paper does not compare with some existing work with a very similar setting (such as Dai et al. 2025).\n- The problem setting of the selection variable is somewhat vague. The symmetry of the selection process requires the selection variable to take different values under intervention. Thus, in a sense, the selection bias is mitigated by introducing the intervention variable in the augmented graph. This is a bit different from the standard interpretation of selection bias in the observational setting. Maybe this deserves a more careful discussion.\n\n\nDai, Haoyue, et al. \"When Selection Meets Intervention: Additional Complexities in Causal Discovery.\" The Thirteenth International Conference on Learning Representations.", "was_revised": true, "final_justification": "In the rebuttal, the authors address my primary concern about the clarity of the paper and comparison with existing works. I have raised my score accordingly. The reason for not giving a higher score is that I think the paper could still be improved with better presentation of the problem setup and the methodology.", "rebuttal": "Thank you for the time dedicated to reviewing our paper. Your insightful comments and valuable feedback have greatly improved its clarity. Following your suggestions, we have carefully revised the manuscript. Please see our point-by-point responses below.\n\n---\n**Q1:** **\"I don't fully understand the symbols in Figure 3. Could you please explain that in more detail with the reference table?\"**\n\nWe have added a detailed explanation of Figure 3 with the reference table to improve clarity.\n\nWe begin by explaining the reference table in Figure 3(f). The **first column** lists key conditional independence (CI) patterns used to distinguish different causal structures. **Each pattern is represented by a specific symbol—such as upward/downward triangles, diamonds, and circles—corresponding to different types of CI relations (e.g., between $I_X$ and $Y$).** \n\nA **green check** indicates that a CI relation **holds**, shown by a **black symbol**. Conversely, a **red X** indicates that the CI relation **does not hold** (i.e., conditional dependence), represented by a **white symbol**. This symbolic representation helps visually summarize the CI patterns associated with each structure in Figure 3.\n\n**Taking Figure 3(a) as an example,**  the structure exhibits conditional dependence $Y \\not\\perp I_X \\mid S$, which is represented by the white upward triangles (△) in the reference table. The conditional independence $Y \\perp I_X \\mid X, S$ is shown in (a) when $X$ is given and is represented by black downward triangles (▼). These CI patterns for (a) are summarized using the corresponding symbols displayed below the causal graph.\n\n---\n\n**Q2:** **\"The paper does not compare with some existing work with a very similar setting (such as Dai et al. 2025).\"**\n\nThank you for pointing out this related work. Following your suggestion, we have added a dedicated discussion in the revised manuscript. In short, while the two methods appear similar at a high level, the key assumptions and problem formulations are in fact quite different. Specifically,\n\n\n**Assumptions of CDIS (Dai et al.):**\n\n**CDIS** models selection as a one-time process that occurs **only before intervention** and becomes inactive afterward, using an interventional twin graph. Under this assumption, the distribution satisfies $P(Y \\mid \\text{do}(X) = x, S) = P(Y \\mid S)$, meaning that after an intervention on $X$ is applied, the selection variable S no longer influences $Y$.\n\n**Assumptions of GISL:**\n\nIn contrast, our work focuses on **biological constraints in genes**—a form of inherent selection bias that exists **both before and after intervention**. These constraints come from basic biological rules, such as essential gene functions or conditions needed for a cell to survive, which do not disappear even when a gene is perturbed. As a result, selection continues to affect the system even after $\\text{do}(X)$, violating the CDIS assumption. In our setting, the correct relation is $P(Y \\mid \\text{do}(X) = x, S) \\neq P(Y \\mid S)$. **This fundamental difference means that CDIS cannot account for such persistent selection effects, leading to poor performance in settings where biological constraints are present.**\n\nMoreover, in light of your suggestion, we have conducted additional experimental comparing GISL and CDIS (Dai et al) on synthetic datasets (number of variables $|\\mathcal{X}|$=15, sample size n = 1500, soft intervention). The results are as below:\n\n| Method | F1| Precision |Recall| SHD | Precision of selection bias|\n|----------|----------|----------|----------|----------|----------|\n| GISL | 80.0 $\\pm$ 1.7  | 83.9 $\\pm$ 3.2  | 74.1 $\\pm$ 1.6 | 4.1 $\\pm$ 5.7| 80.8 $\\pm$ 11.3| \n| CDIS   |  63.8 $\\pm$ 4.5     |  61.4 $\\pm$ 4.7  |  68.1 $\\pm$ 5.3  |   8.4 $\\pm$ 6.4   | 34.1 $\\pm$ 9.1 |\n\nWe demonstrate that, within the context of specific biological constraints, our method achieves superior performance due to the correct identification results.\n\n---\n**Q3:**  **The problem setting of the selection variable is vague. The interpretation of the selection variable and perturbation indicator is unclear.**\n\nThanks for your feedback, which has certainly improved the readability of our paper.\n\n**Problem Setting:**\n\n- **Goal:** Identify regulatory relationships, latent confounders, and selection bias (biological constraints) in gene expression data.\n\n-   **Assumption:** Gene expression is affected by latent confounders and persistent selection bias, as supported by biological evidence (see Lines 29–31 and 33–39).\n\n-   **Solution:** We formulate this as a causal discovery problem under latent confounding and selection bias. By leveraging both observational and gene perturbation (interventional) data, GISL identifies regulatory relationships, latent confounders, and selection bias by capturing differences in **perturbation symmetry** and **intervention effects** (see Section 3.1).\n\n**Clarification on the Selection Variable and Perturbation Indicator:**\n\nAs described in Lines 85–88, **the selection variable S imposes a constraint where only samples with S = 1 are observed; data with S = 0 is unobservable.** This formulation is formally discussed in the \"M Bias\" section on page 5 in [1] and has been widely adopted in subsequent works [2–5]. Thus, our analysis is restricted to the selected subpopulation.\n\nImportantly, **selection bias and intervention are independent**—selection can affect both observational and interventional data. The intervention is designed and classified in Sections 4 and 5 (Pages 5-7) in [6].  The  GISL is designed to **identify** selection bias, not eliminate it. To this end, we introduce a perturbation indicator into an **augmented causal graph**, which enables joint modeling of observational and interventional distributions. This allows CI tests to capture both **perturbation symmetry** and **perturbation effects**, facilitating the identification of regulatory relationships, latent confounders, and selection bias.\n\n> [1] Quantifying Biases in Causal Models: Classical Confounding vs Collider-Stratification Bias. Epidemiology, 2003.\n\n> [2] Causal discovery from data in the presence of selection bias. PMLR, 1995.\n\n> [3] Structural approach to selection bias. Epidemiology, 2004.\n\n> [4] Sample selection bias as a specification error. Econometrica: Journal of the Econometric Society, 1979.\n\n> [5] Endogenous Selection Bias: The Problem of Conditioning on a Collider Variable. Annual Review of Sociology, 2014.\n\n> [6] Variables of causal intervention. UAI, 2004.\n\n---\n**Q4:**  **\"What are the terms like $[I_i, X_j \\mid S]$? Are those indicator variables for the conditional independence?\"**\n\nThanks for pointing this out. $[I_i, X_j \\mid S]$ denotes that we need to capture the CI relation between the intervention indicator $I_i$ and the variable $X_j$, given the selection variable $S$. \n\n---\n**Q5:**  **\"How scalable is the proposed method? How sensitive is the proposed method to the size of datasets?\"**\n\nScalability is a key consideration in the design of GISL. GISL is built to handle **large-scale gene expression datasets**, as demonstrated in our real-world experiments (Section 4.3 and Appendix F.1, Lines 760–786), where GISL was applied to datasets with over **5,000 observed genes (variables)** and more than **100 perturbed genes**. In comparison, the original nonparametric CI test KCI [8] was only evaluated on graphs with **5 variables**, and baseline methods such as [7] consider at most **20 variables**. These comparisons highlight the practical scalability of GISL beyond prior work.\n\nIn terms of sensitivity to dataset size, our experiments on synthetic data (10–25 variables, 500–2,000 samples) and real data show that GISL maintains stable performance across varying sample sizes. This is enabled by two design choices: (1) the use of **skeleton discovery** via FGES, and (2) the application of **pairwise CI tests** based on limited conditioning sets, both of which help ensure computational efficiency as dataset size grows.\n\n> [7] Permutation-Based Causal Structure Learning with Unknown Intervention Targets. UAI, 2020.\n\n> [8] Kernel-based Conditional Independence Test and Application in Causal Discovery. UAI, 2011.\n\n---\n**Q6:** **\"The option in Algorithm 1  $\\mathcal{J}==(c)$  and $\\mathcal{J}==(d)$  seems to be flipped and inconsistent with the caption of Figure 3.\"**\n\nThanks for pointing this out. The options for $\\mathcal{J}==(c)$  and $\\mathcal{J}==(d)$ in Algorithm 1 were inadvertently flipped. We have corrected this error to ensure consistency between the algorithm and the figure.\n\n---\n**Q7:**  **\"Do the terms regulatory relations and causal relations mean the same thing in this context? The paper seems to use these terms interchangeably.\"**\n\nYes, in this context, _regulatory relation_ and _causal relation_ refer to the same concept. We primarily use the term _regulatory relationship_ to align with the gene regulatory network inference (GRNI) context, while _causal relationship_ is used when introducing the formal model to remain consistent with standard causal discovery terminology. This interchangeable use is clarified in Lines 79–81 and is consistent with the GRNI definition in Line 22, where regulatory interactions are modeled as causal effects.\n\n---\n If you have any questions, please do not hesitate to let us know. We would greatly appreciate your feedback and are happy to discuss further."}, {"review_id": "Vq9bqLxLHr", "reviewer": "Reviewer_7xy2", "summary": "This paper focuses on gene regulatory network inference (GRNI), where the goal is to identify causal relationships between gene expression (how one gene promotes or represses the expression of another gene in a cell) based on measurements of gene expression levels observed in cells. The authors highlight that the set of cells observed can represent a significant selection bias; two genes may appear to have high levels of co-expression not because of any specific regulatory relationship but rather because cells with both genes are more likely to be selected for sequencing (“satisfying the survival criteria”).\n\nThe authors point out a lack of methods for GRNI that address this selection bias. They then note that unlike bias arising from hidden confounders, selection bias is symmetric: meaning that perturbing gene 1 influences a behavior in gene 2, and perturbing gene 2 also influences gene 1. They propose a novel method, gene regulatory network inference n the presence of selection bias and latent confounders (GISL) that exploits this symmetry to tolerate selection bias when performing GRNI.\n\nContributions:\n- The authors propose a novel method for GRNI in the presence of selection bias, and claim consistency results under minimal assumptions.\n- The authors draw attention to the issue of selection bias and highlight some tractable properties (symmetry) of selection bias.\n- The authors conduct experiments on both synthetic and real-world data to compare their novel procedure to baselines.", "questions": "Regarding Figure 3: What do the set of symbols (up/down triangles, diamonds, circles in black and white) refer to, and what does it mean for them to be “green check” versus “red X” in the reference table? Are these symbols meant to show that certain conditional independence relationships are/aren’t satisfied in each of the columns? Is it important that columns c and e grey-shaded, and is that meant to communicate something?\n\nAs explained in Weaknesses, I did not understand the discussion of mutilated DAGs and I did not understand how they related to the rest of the paper. Could the authors briefly explain?\n\nIn Figure 4, GISL clearly out-performs the baselines in F1, precision, and SHD, but achieves only middling results in recall. Do the authors have any hypotheses for why GISL would perform so well with regards to the other metrics, but not with respect to recall?", "limitations": "yes", "rating": 3, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 1, "significance": 2, "originality": 4, "strengths_and_weaknesses": "Strengths:\n- The authors make a compelling case that selection bias is a serious obstacle to causal discovery and that it is under-addressed.\n- The paper contains several examples that illustrate concepts at a high-level. Section 2.2 and Figure 2 do an excellent job of communicating spurious dependencies arising from selection bias.\n- The authors claim consistency results for their procedure under minimal (and very typical) assumptions.\n\nWeaknesses:\n\nIn several places, the technical content is extremely dense and terms/quantities appear without definition. This seriously limited my ability to understand central parts of the paper, such as the pseudocode of Algorithm 1, the novel proposed method. \n- Figure 3 looks very detailed but I could not interpret it at all. Unfortunately this is an essential figure to understand, because the novel method proposed (Algorithm 1) is defined by referencing sub-panels of Figure 3. I list some clarifying questions about Figure 3 in the Questions section of the review.\n- The preliminaries are full of undefined technical terminology, e.g. mutilated DAG. If these are standard terms within GRNI and the other reviewers disagree with my assessment then perhaps this is not a significant weakness, but as a reader I found this very confusing. These terms are not defined in the main body nor in Appendix A, which is designated for definitions of core concepts. Some terms are explained indirectly much later in the paper (e.g. knockout vs knockup and knockdown for hard/soft interventions), but mutilated DAGs are never referenced again and I still do not understand what they represent nor why they were presented in the preliminaries.\n- I found several sentences un-parsable. E.g. Line 102, “For hard interventions, consider each $T_k$ as factoring in a mutilated DAG…”. $T_k$ is not defined, and given the general lack of definitions I found it impossible to decipher what $T_k$ referred to based on context. Is it an element of $\\cal{T}$? I thought $\\cal{T}\\subseteq \\cal{X}$ so I expected elements of $\\cal{T}$ to be denoted $X_k$.\n\n\nI caveat this by stating that I am not familiar with definitions relating to interventions, so if other more experienced reviewers find the paper reasonably readable then I defer to their opinions.", "was_revised": true, "final_justification": "I maintain my score (borderline reject) because I do not believe the core contributions of this paper are explained in sufficient clarity in the current state of the manuscript. I explained my reasoning to the authors in my comments, and I reiterate my main concern here: multiple reviewers highlighted confusion about Figure 3. Because the main novel method proposed in this work (Algorithm 1) is defined only through reference to Figure 3, in my opinion, it is necessary that Figure 3 be understandable to readers and completely correct. Personally, I am not confident in my level of understanding of the algorithm and thus cannot be confident in my ability to verify the consistency results. I mark my confidence level at the lowest available (1) to reflect this.", "rebuttal": "Firstly, we would like to sincerely thank you for the considerable time and effort invested in reviewing our paper. Your insightful comments and suggestions have significantly helped us enhance the clarity and readability of our manuscript. In response to your valuable feedback, we have carefully revised and clarified our paper:\n\n---\n\n**Q1: What do the symbols and reference table  mean in Figure 3?**\n\nWe add a detailed explanation for Figure 3 to clarify the notation.\n\nThe set of symbols (up/down triangles, diamonds, and circles) represents the four types of conditional independence (CI) relationships listed in the first column of the reference table in Figure 3(f). Each symbol corresponds to a specific CI pattern. The  **green check** indicates that the CI relationship is **satisfied** and is marked with corresponding **black symbols**,  while the **red X** indicates that the CI relationship is **not satisfied** and is marked with **white symbols**. For example, a black upward triangle indicates that $I_X \\perp Y \\mid S$ holds (i.e., I_X and Y are conditionally independent given S), whereas a white upward triangle indicates that $I_X \\not\\perp Y \\mid S$ (i.e., they are dependent). We hope this clarification makes the figure more accessible and easier to interpret.\n\n**Taking Figure 3(a) as an example,**  the structure exhibits conditional dependence $Y \\not\\perp I_X \\mid S$, which is represented by the white upward triangles (△) in the reference table. The conditional independence $Y \\perp I_X \\mid X, S$ is shown in (a) when $X$ is given and is represented by black downward triangles (▼). These CI patterns for (a) are summarized using the corresponding symbols displayed below the causal graph.\n\n**Q2: What do the grey-shaded columns (c) and (e) mean in Figure 3?**\n\nThe grey shading in columns (c) and (e) indicates that both cases yield identical conditional independence (CI) patterns. Specifically, as discussed in Remark 3.4, the selection process introduces additional conditional dependencies into the observed data, resulting in indistinguishable CI patterns for these two scenarios. Consequently, GISL is unable to differentiate a scenario involving only selection bias (c) from one involving both selection bias and a causal relationship (e). However, this limitation does not impact GISL's capability to identify the presence of selection bias itself, even though the precise causal structure remains ambiguous. This point is formalized in Theorem 3.7. The grey shading thus serves to distinguish scenarios that GISL can uniquely identify (white shading: cases a, b, and d) from those it cannot (grey shading: cases c and e, which are indistinguishable).\n\nPlease let us know if you have any questions about Figure 3. We would be delighted to discuss this further and deeply appreciate your thoughtful feedback.\n\n---\n\n**Q3: Technical terminologies & terms are Undefined: mutilated DAG, knockout, knockup, and knockdown.**\n\nWe sincerely appreciate your careful reading and thoughtful comments. We are grateful for the opportunity to clarify these terms and improve the clarity of our presentation.\n\nTo enhance readability, we have added explicit definitions for commonly used terms introduced in the preliminaries, including mutilated DAG (widely used in interventional causal discovery), perturbation indicator (Definition 1 on page 6 of [1]), intervention targets (third paragraph on page 4 of [2]), and hard/soft interventions (illustrated in the figures on pages 4 and 6 of [3]).\n\n\n**Why introduce Mutilated DAG?**\n_Mutilated DAG_ represents the graphical structure for hard intervention (illustrated in Figure 2(a) and described in Lines 102–104), while _Augmented DAG_ is introduced in Lines 105–108 to represent soft interventions. Since our framework accommodates both types of interventions, it is necessary to present both forms of DAG representations.\n\n\n**Logical progression in introducing terms.**\nWe have aimed to reduce cognitive load for the reader by introducing terminology in a gradual, structured manner. In our preliminaries (Lines 100–102), we **first** introduce gene perturbation methodologies—**gene knockout and transcriptional modulation** (gene knockup or knockdown)—which correspond to hard and soft interventions, respectively, in the causal inference framework. We then describe how these two types of interventions are represented graphically. Specifically, a **mutilated DAG** refers to the graph obtained from a causal DAG by **removing all incoming edges to the intervened variables**, effectively breaking the causal dependence of those nodes on their parents. Similarly, the **augmented DAG** is introduced in Lines 105–108 to represent soft interventions, and is illustrated in Figure 2(b).\n\n> [1] Varieties of causal intervention. PRICAI 2004.\n\n> [2] Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. JMLR 2012.\n\n> [3] Interventions and causal inference. PSA 2006.\n\n---\n\n**Q4: $T_k$ in the sentences is Unparsable.**\n\nEach $T_k$ denotes the target variable set for the k-th intervention. We have added the definition of $\\mathcal{T} = \\{T_1, \\ldots, T_K\\}$ in the revised version. \n\n\n--- \n\n**Q5: Why GISL performs not so well regarding the recall?**\n\nAs you pointed out, GISL shows a little bit lower recall in **identifying regulatory relationships** under selection bias. This is because, as discussed in **Q2** and Figure 3, in theory, it is challenging to distinguish `(c): selection bias` and `(e): causal + selection bias` -- they share the same CI patterns (all conditionally dependent). As a result, our GISL categorizes both cases `(c)` and `(e)` as selection bias and thus shows lower recall when identifying regulatory relationships (causal relations). \n\nOn the other hand, the baseline algorithm, e.g., GIES, exhibits higher recall, but this does not necessarily indicate better performance. The key issue is that GIES is unable to account for selection bias, which introduces spurious dependencies into the data. When GIES relies on a score function to detect conditional dependencies, it may mistakenly interpret these spurious associations as causal relationships. As a result, the algorithm tends to include additional edges regardless of whether true causal relationships exist, leading to inflated recall but low accuracy. \n  \n\n--- \nWe hope our responses have clarified the key points and improved the overall clarity of the paper. If any aspects remain unclear, we would be happy to provide further explanation. Thank you again for your thoughtful review and the considerable effort you invested."}]

Author rebuttal:

Meta-review: The submission develops an algorithm for gene regulatory network inference in the presence of selection bias. Such bias is manifested in the form of symmetric perturbation effects, which cannot be explained by latent confounders assuming a causal DAG. Reviewers overall appreciated the efficacy and strong motivation of the idea. Reviewers raised concerns relating to comparisons with existing work (7UVK, TNYS), technical clarity (7UVK), lack of experiments with both confounders and selection bias (TNYS), lack of sensitivity analysis (yzJf), and lack of real-world case studies (yzJf). In discussions, authors provided additional results with latent confounders, more thorough sensitivity analyses, and  real-world case studies. Reviewers indicated that most of these concerns were addressed to their satisfaction, leading to a consensus towards acceptance.
Decision: Accept (poster)