Granger causal inference on DAGs identifies genomic loci regulating transcription
Authors: Alexander P Wu
Rohit Singh
Bonnie Berger
Venue: ICLR 2022 Poster (2022)
Keywords: Granger causality
causal inference
graph neural networks
gene regulation
single-cell genomics
chromatin accessibility
directed acyclic graphs
single-cell multimodal

Abstract: When a dynamical system can be modeled as a sequence of observations, Granger causality is a powerful approach for detecting predictive interactions between its variables. However, traditional Granger causal inference has limited utility in domains where the dynamics need to be represented as directed acyclic graphs (DAGs) rather than as a linear sequence, such as with cell differentiation trajectories. Here, we present GrID-Net, a framework based on graph neural networks with lagged message passing for Granger causal inference on DAG-structured systems. Our motivating application is the analysis of single-cell multimodal data to identify genomic loci that mediate the regulation of specific genes. To our knowledge, GrID-Net is the first single-cell analysis tool that accounts for the temporal lag between a genomic locus becoming accessible and its downstream effect on a target gene's expression. We applied GrID-Net on multimodal single-cell assays that profile chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) in the same cell and show that it dramatically outperforms existing methods for inferring regulatory locus-gene links, achieving up to 71% greater agreement with independent population genetics-based estimates. By extending Granger causality to DAG-structured dynamical systems, our work unlocks new domains for causal analyses and, more specifically, opens a path towards elucidating gene regulatory interactions relevant to cellular differentiation and complex human diseases at unprecedented scale and resolution.

Published as a conference paper at ICLR 2022
GRANGER CAUSAL INFERENCE ON DAGS IDENTIFIES GENOMIC LOCI REGULATING TRANSCRIPTION
Rohit Singh i.l., Alexander P. Wu & Bonnie Berger 122+
I Computer Science and Artificial Intelligence Laboratory. MIT, Cambridge, MA 02139, USA
2 Department of Mathematics, MIT, Cambridge, MA 02139, USA { [rsingh, a alexwu, bab) @csail.mitt.eu
Authors contributed equally to this work
ABSTRACT
When dynamical system be modeled sequence of observations, Granger causality is a powerful approach for detecting predictive interactions between its variables. However, traditional Granger causal inference has limited utility in domains where the dynamics need to be represented as directed acyclic graphs (DAGs) rather than linear sequence, such with cell differentiation trajec- tories. Here, we present GrID-Net, a framework based on graph neural networks with lagged message passing for Granger causal inference on DAG-structured sys- tems. Our motivating application is the analysis of single-cell multimodal data to identify genomic loci that mediate the regulation of specific genes. To our knowl- edge, GrID-Net is the first single-cell analysis tool that accounts for the temporal lag between a genomic locus becoming accessible and its downstream effect on a target gene's 's expression. We applied GrID-Net on multimodal single-cell assays that profile chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) in the same cell and show that it dramatically outperforms existing methods for inferring regulatory locus-gent links, achieving up to 71% greater agreement with independent population genetics-based estimates. By extending Granger causality to DAG-structured dynamical systems, our work unlocks new domains for causal analyses and, more specifically. opens a path towards elucidating gene regulatory interactions relevant to cellular differentiation and complex human diseases at un- precedented scale and resolution
INTRODUCTION
Understanding the structure of multivariate dynamical system often boils down deciphering the causal relationships between its variables. Since inferring true causality is often infeasible, requiring additional experiments or independent mechanistic insight, statistical analysis of observational data to identify predictive relationships between the system's 's variables can be very valuable. The frame- work of Granger causality does exactly that: in a dataset where observations are temporally ordered, a tine-dependem variable z (with value Ia at time f) is said to "Granger cause" second variable y if the history of x at time (i.e., Ti, -4-1) is significantly predictive of y at time even after accounting for y's own history (Granger] 1969 Shojaie & Fox 2021). Originating in economet- rics, Granger causality has been a powerful tool in many domains including biology, finance, and social sciences (Fujita et al.) 2007: Yao et al. 2015 Benhmad] 2012) Rasheed & Tahir 2012). The prerequisite for applying Granger causality, however, is that there be a clear sequential ordering of observations, i.e., the data must conform to total ordering along time.
Often, only partial ordering of the observations is possible. For instance, cell differentiation trajectories may have branches. In text-mining, the citation graph of publications captures the flow of ideas and knowledge. In such cases, the dynamics of the system are more suitably represented as a directed acyclic graph (DAG) corresponding to the partial ordering, with nodes of the DAG representing dynamical states and its edges indicating the flow of information between states.
"o-correeppondng authors: {rsing bab lasaai wwe.saaal...........mmmmmm 1...m.mmme....
'The code for GrIlD-Net is available at http:s:///ttttie com axxw 6/gridnes dnet

---

Published as a conference paper at ICLR 2022
(xw.c)
(Xwyss
KNN grash of cells
(RNA-seq ATAC the same cell)
Multimodal single-cell data
Ossys)
from integrated
multimodial data
Ossuss
Coll 1
edge direction from
Cell 2
paeudotime
Cell
(x.y6)
Cell N
(6..9)
(G..K)
RNA-soq
ATAC-soq
Standard
DAG
Granger Causality
Granger Causality (this paper)
Closed
Open chromati
d
he2
8=
h
Pronwer
The subset ancestors covered
-=1
h(
GriD-Net Message Passing (L=3)
Figure 1: OVERVIEW (a) We extend the concept of Granger causality, previously applied only to sequentially ordered observations, to partial orderings, i.c., where the information flow be described a directed acyelic
graph (DAG). (b) GrID-Net, our approach for inferring Granger causality on DAG- AAG-structured system, is a graph neural network framework with lagged message-passing (f, I., are defined in Sec. 2.2 (c) To identify genomic loci whose accessibility Granger causes the expression of a specific gene, we apply GrID-Net on multimodal single-cell data where chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) have been profiled in the same cell. We model the lag between them by relating the RNA-seq readouts of cell not to the ATAC-seq readouts of same cell to ATAC-seq readouts cells slightly earlier in the profiled biological process, applying GrID-Net on a DAG constructed from the kNN graph of cells with edges oriented as per pscudotime. (d) GrID-Net leverages the biological intuition that the accessibility of specific genomic locus (ATAC-seq) precedes the binding of regulator proteins the locus, thus Granger causing the expression levels for its targeted gene (RNA-scq) to changed.
The key conceptual advance of this work is extending the applicability of Granger causality to sys- tems with partially ordered observations (Figure Πa). On DAG where each node takes on multi- variate values, we define variable x as Granger causing another variable y if the latter's s value yo at the DAG node v can be significantly predicted by the values of æ the ancestor nodes of u, even after accounting for yo is own ancestral history. We introduce GrID-Net (Granger Inference on DAGs), graph neural network (GNN) framework for Granger causal inference in dynamical systems that can be represented DAGs (Figure[p)) We modify the standard GNN architecture to enable lagged message- eaataage-pssing on the DAG, enabling us to accumulate past information for a variable accord- ing to the DAG's structure. This, combined with mathematical formulation of Granger causality more amenable to neural networks, allows GrID-Net to recover nonlinear long-range Granger causal dependencies between the variables.
While our problem formulation is fully general and GrID-Net is broadly applicable, our motivating application is inferring noncoding genomic loci that influence the expression of specific gene. We apply GrID-Net to single-cell multimodal studies that assay chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) in the same cell. seeking to identify Granger causal relationships between the accessibility of individual chromatin regions (ATAC-seq "peaks") and the expression of specific genes. The system dynamics are represented by a DAG where each node corresponds to a cell. Edges connect cells in similar dynamic states (i.e. similar ATAC-seq and RNA-seq profiles), with edge directions estimated by pseudotime analysis (Haghverdi et al. 2016) (Figure Ie).
Our work addresses pressing need to identify temporally causal relationships between noncoding loci and gene expression, as traditional approaches for doing so are expensive and slow. The latter include population genetics techniques like expression quantitative trait locus (eQTL) studies, in which human genetic variability at particular genomic locus is associated with expression changes in the gene of interest. Alternatively, perturbation- expression changes in response to CRISPR alterations of specific noncoding genomic loci (Fulco) al. 2016). In contrast, we capitalize on the emergence of single-cell multimodal assays that
approaches directly seek to identify gene
2

---

Published as a conference paper at ICLR 2022
profile ATg and RNA-seq simultaneously and offer new way of discovering such regulatory relationships. Rather than requiring the large population sizes of eQTL studies or being limited to specific noncoding loci in the perturbation- approaches, these multimodal experiments enable unbiased high-throughput, genome-wide estimations of peak-gene associations.
However, leveraging this data will require algorithmic innovations. Over 99.5% of ATAC-seq ob- servations ("peaks") in cell are reported as empty (Singh al., 2021; Zhang al. 2019). and this sparsity constrains current approaches. A key limitation of existing methods is also that they are based on correlation estimates and cannot capture the dynamics governing the temporally causal re- lationship between chromatin accessibility and gene expression. For a noncoding locus to influence a gene's 's expression, the locus typically must first become accessible, upon which regulator proteins act on it to drive change in gene expression (Figure Tt). This causal mechanism, therefore, entails temporal asynchrony between the ATAC-seq and RNA-seq modalities (Lara-Astiaso et al. 2014] Ostuni et al. 2013). GrID-Net is specifically designed to recover these dynamics between peaks and genes while being robust to the noise and sparsity in single-cell data.
We apply GrID-Net on three single-cell multimodal datasets that characterize a range of dynamic processes, including cancer drug response and cellular differentiation. We show that GrIDD-ett substantially outperforms current methods in identifying known peak-gene links. To our knowledge, this work presents both the first framework for Granger causal inference dynamical systems represented as DAGs and the first single-cell method to explicitly model the temporal asynchrony between chromatin accessibility and gene expression.
2 METHODS
2.1 BACKGROUND ON GRANGER CAUSALITY: STANDARD FORMULATIONS
In time series analyses, the classical Granger causality framework uses the vector autoregressive (VAR) model, which can be expressed as follows (Lütkepohl 2005).
(1) 02 It- +
Vt
y
Here, ye, E R are values at time & for two stationary time series and x, respectively; aj and (f) coefficients that specify how lag affects the future of time series y; and 6, is zero-mean a2
noise term. In this model, yt is assumed to be a linear combination of the most recent values of and x each, with said to Granger cause time series if and only if ag +f 0 for all l.
A more generalized formulation - the one we follow - for performing Granger causal inference is to consider two related models for forecasting time series y. A full model considers the past values of both y and to forecast y. while reduced model excludes the effect of x, only containing terms related to the past values of y.
V/2 = f f(ull)) (yt-1,..., Li D1 Dt It-L) + ce
(2)
y1 = (reducced) (J:-1, .... Vt L) +
(3)
Here, f(full)))) (-) and f(eeduced) (-) are generalized functions that specify how the value of y depends on past values of y and x. The predictions of the two models are then compared, upon which a Granger causal relationship is declared if the full model's predictions significantly more accurate.
2.2 GRANGER CAUSALITY ON A DAG: GRAPH NEURAL NETWORK FORMULATION
Let the data be represented by a DAG ç = EV,8) with n = | nodes (i.e., observations) and directed edges& E indicating the information flow partial order between these observations. For instance, when applied to the case of standard Granger causal inference on time series data, 9 would be the linear graph corresponding to the time series. Let y. E I² correspond to the values of variables and x on the nodes in V. with x putatively Granger-causing y. To infer Granger causality, compare the full and reduced models as above:
3

---

Published as a conference paper at ICLR 2022
- 9 (h (ful)) + ccllfult))
(4)
y - ( (reduced)) +
(5)
Here, ((fult)
bgg ((reduced) [(/ult) RR represent the historical information of y or J. (as denoted
by the subscript) that is aggregated by graph neural network (GNN) layers. For a history h, the value h(v] at node E V is the information accumulated from e's ancestors in g. While hy (fall) and hy (reduced) are outputs determined from similar architectures, their learned weights and bias terms end up being different since ht (full) is influenced by its interaction with h The coefficient mediates this interaction and describes the effect æ's history on y. addition, 6 IR a zero- mean noise term. We also include g(-) as optional link function that maps its input to the support of the distribution of variable y. In our single-cell application, we set g as the exponential function since the target variable (y) represents normalized transcript counts, which
non-negative.
(reduced) being indepen- For brevity, we only describe how (fult) is computed, with hg (Fall) and hy dently computed analogously. Also, since æ occurs only the full model, for notational convenience we drop its superscript (full) below, writing just hg. We express h as the mean the outputs of consecutive GNN layers, denoting the layerwise outputs as h²
1 L
h
h
x
(6)
o (w ATT
+ (())
if(>1
1
(w A )
if
Here, w 1 E are the per-layer weight and bias terms, respectively. a(-) represents the nonlinear activation in each of the GNN layers, chosen here be the hyperbolic tangent function. A and A_ E gnxa are matrices defined by the DAG structure. A is the column-normalized adjacency matrix of G, with A, 2; edge (i.j) E & and 0 otherwise, where d, the in-degree of node j. We note that DAG does not have self-loops, so the diagonal terms of are zero. In A., include these diagonal terms, setting them to the same value others in the column (i.e., d, is incremented by in A+ to preserve normalization).
Together, A and A allow us to extend the key intuition of Granger causality to DAGs: predicting a value node v from the values of its ancestors that are within L steps. The first layer introduces the lag, using information at the parents of u but not the information at v itself. Each subsequent layer introduces the preceding set of parents, with the diagonal term in A4 ensuring that information already stored at is integrated with the new information propagated via v's parents. The sequence of GNN layers therefore reflects the successive aggregation of past information by traversing steps backward along the DAG for each node in the graph. Thus, we are able to conceptually match time series-based formulations (e.g., Eqn. T but with the crucial advantage of leveraging the richness of DAG's s structure in directing the information flow and aggregating sparse data.
2.3
COMPARING FULL AND REDUCED MODELS TO INFER GRANGER CAUSALITY
Let @.ful)) @ ffredwoed) O (fult) E 2t. denote the set of parameters for the full and reduced models, with (full)
(we) (1) (
(L), ( b(4) and (ffull)) ( (reduced) defined analogously. All the
④ parameters are jointly learned by minimizing the combined loss of the full and reduced models (Montalto et al. 2015):
(ttotal)
Σ
Σ C(gSruln, yv) + E(ggfreduccn yu)
(8)
Z.DEP
NEEV
DEV
where (fall) and g/ (reduced) correspond to the predictions for observation v by the full and reduced models, respectively. We note that the full and reduced models have completely separate parameters
4

---

Published as a conference paper at ICLR 2022
Also, in a multivariate system, not all pairwise combinations of variables may be relevant: P is the subset of æ, y variable pairs whose putative Granger-causal interactions are of interest.
(see discussion in Appendix A.8)
To infer Granger causal interactions between y and x; compare the set of loss terms associated with the full model Cgj) . to those of the reduced model L(je ((ffredueed) yy),, noting that the precise functional form of L depends on the domain. If x does not Granger cause y, the full and reduced models will have similar explanatory power and hence similar loss values. Otherwise, the full model will have a lower loss. We assess this by one-tailed F-test, comparing the residual sum of squares of the reduced and full models; an alternative assessment using Welch's 1-test yielded very similar results (Appendix (.77). We rank the set of candidate x-y interactions by the F-statistic, with a higher score corresponding to stronger evidence for a Granger causal interaction.
2.4 DOMAIN-SPECIES MODEL CUSTOMIZATION
The loss & should be chosen as appropriate for the domain. In our single-cell context, y corresponds to gene expression. We preprocessed RNA-seq transcript counts, normalizing and log transforming them so that the mean squared error loss was appropriate; we note that most methods based on Pearson correlation of gene expression also seek to minimize squared loss, explicitly or implicitly.
The GrID-Net model can also be customized to account for different lags and lookbacks. The num- ber of GNN layers I. corresponds to the maximum amount of past information desired. Similarly, to introduce a k-hop lag, the first & GNN layers would use the matrix A while the later layers would use A In the sections above, we described one-hop lag that we have used for all analyses here.
2.5 TRAINING DETAILS AND HYPERPARAMETERS
GrID-Net models were trained using the Adam optimizer with learning rate of 0.001 for 20 epochs or until convergence (defined to be the point at which the relative change in the loss function is less than 0.1/ (P| across consecutive epochs). A minibatch size of 1024 candidate peak- gene pairs was used during training, and trainable parameters in the model were initialized using Glorot ini- tialization (Bengio & Glorot 2010). All GrIID-Net models consisted of I. 10 GNN layers; the architectures of the three sub-models (hy (h (reduced) hy (fall)) and h fffll))) were identical but separate. All models were implemented in PyTorch and trained on a single NVIDIA Tesla V100 GPU.
3 RELATED WORK
Single-cell multimodal methods: The emergence of single-cell multimodal assays has led to the development of tools to leverage the joint profiling of multiple genomic, epigenetic, and functional features in each cell. Many of these tools aim to synthesize the multiple modalities within a dataset to perform downstream analyses that elucidate cell state or gene programs (Singh et al.) 2021 Arge- laguet et al.) 2020 Gayoso et al. 202I). Meanwhile, efforts to detect associations at the level of specific features across modalities have thus far been limited. Current approaches for inferring peak-gent associations predominately rely on simply calculating correlations between these two modalities (Ma et al.) 2020: Zhu al. 2019 Liu et al. (219)). Not only do these sirelltioo-based approaches fail to account for the temporal asynchrony between chromatin accessibility and gene ex- pression, but they are also sensitive to the inherent noisiness of single-cell data (Zhang al. 2019). These limitations point to the need for going beyond analyzing mere correlations to leveraging the joint dynamics of these features in causal inference frameworks, like Granger causality.
Granger causality: Recent approaches for extending Granger causality have primarily focused on enabling the detection of nonlinear or graph graph-based interactions in multivariate time series data, One category of methods involves nonlinear kernel-based regression models (Marinazzo al., 2008] Ren et al. 2020): another seeks to learn graphs that describe causal relationships between variables (Basu et al.) 2015, Shojaie & Michailidis) 2010). Neural network- Granger causality models have also more recently been proposed to account for more generalized nonlinear interactions (Tank) jt al. 2018 Marcinkevičs & Vogi, 2021). Of these methods, efforts have been made to leverage specific neural network architectures more fully take advantage of the sequential ordering of observations in time series data (Khanna & Tan, 2020 Nauta et al.) 2019). Customized Granger causal approaches have also been designed for specific biological applications (Finkle et al.! 2018).
5

---

Published as a conference paper at ICLR 2022
As clarification, we distinguish the DAGs referred to in this paper from the concept of "causal graphs". Here, each node in the DAG corresponds to a state with an associated multivariate obser- vation, and edges describe the ordering between states. In contrast, a causal graph describes data generation process with nodes corresponding to the dynamical system's 's variables. In our single-cell peak-gene inference application, performing Granger causal inference on DAG allows us to cap- ture the diversity of cell states as well as their dynamics. We note that using the standard Granger causality framework would require forcing the rich heterogeneity of single-cell dataset into a se- quential ordering of cells along pseudotime, suboptimal choice due to the possibility of multiple differentiation trajectories, substantial biological and technical variation in single-cell data (Zhang) et al. 2019), and multiple cells representing very similar cell states (Baran et al.) 2019).
We also note that a limitation of most Granger causal analyses - including ours is the inabil- ity to rule out hidden confounders that might mediate the interaction between x and y [Mastakouri et al. 202I). Nonetheless, such analysis remains valuable in our peak-gene inference task, as cnables the identification of peak-gene relationships that are temporally causal, upon which short- list of hypotheses can be generated for subscquent perturbational More importantly, even indirect peak-gene interactions are biologically meaningful. In future work, other approaches for causal inference (e.g.. |Entner& & Hoyer (2010); Pfister et al. (2019)) includ- ing Bayesian inference (Brodersen et al. 2015]
population geneties studies.
| Glymour et al. 2019) or structural equation models
(Spirtes& & Zhang 2016) could be applied to address this limitation.
4 RESULTS
Datasets and Preprocessing We analyzed three single-cell multimodal datasets with per-cell ATAC-seq and RNA-seq profles profiles that characterize a range of dynamical systems, including cell dif- ferentiation and drug-treatment responses (Cao et al. 2018 Chen al. 2019 Ma al.) (2020). The sci-CAR dataset profiled 3,260 human lung adeeoocacinoma--erivea A549 cells after dexametha- sone treatment. The SNARE-seq dataset evaluated 5,081 cells from the mouse neonatal cerebral cortex. Lastly, the SHARE-seq dataset contained 34,774 cells from differentiating skin cells in mouse. We applied geometric sketching (Hie et al., 2019) to identify representative subset of 5,000 cells for the SHARE-seq dataset, which we use for all downstream analyses. These datasets are high-dimensional (about 200, 000 peaks and 30, 000 genes per cell) and sparse (over 99.5% of peaks in a cell have zero counts). For each dataset, GrID-Net produces results that are specific to the cell type/state composition of the study. While we chose these studies to be tissue-specific, there remains some diversity in cell type and state due to cell-differentiation and cntuuurbatio-resppnae variability. To flexibly hone in on specific cell types within a study, users can adapt GrID-Net by simply filtering the dataset to preserve only the cells of interest.
Following recommended practice (Hie et al., 2020), we preprocessed each dataset by applying log(1 + CPM and log(1 + +CPM/100) transformations to the raw RNA-seq and ATAC-seq count data, respectively; here CPM indicates counts per million. The transformed RNA-seq data for each gene was also divided by the maximum value for that gene. For each single-cell multimodal dataset, we then compiled set of candidate peak-gene pairs (i.e.. P in Eqn. 8 by selecting ATAC-seq peaks within 1 Mb from each gene (Appendix A.3). We chose this genomic disttnce cutoff because almost all enhancer-gent gene pairs reported a recent benchmark study were within Mb (Moore et al.) [211)) After this filtering, we evaluated 507,408, 916,939, and 872.039 candidate peak-gene pairs for the sci-CAR, SHARE-seq, and SNARE-seq datasets, respectively. Training GrID-Net on these sets of candidate peak-gene pairs took roughly 2 2-6 hours and required maximum of 40 GB of RAM. We further discuss scalability, runtime and memory usage details along with recommendations for using geometric sketching (Hie et al. 2019) for large datasets in Appendix A.12
For each of the three multimodal datasets, a kNN graph was constructed on cell representations that unified information from both gene expression and chromatin accessibility via Schema (Appendix A. Singh et al. (202I)); other approaches for integrating data modalities (e.g., manifold alignment (Cao et al. 2020 could also be used here. We next inferred a pseudotime value for each cell. Edges in the KNN graph were then retained only if they aligned with the direc- tion of increasing pseudotime. thus ensuring that cycles exist and the resulting graph is DAG. In our results below, we used 1=15 to build the kNN graph and the diffusion pseudotime algorithm (Haghverdi et al., 2016) to infer pseudotime stamps, although we found GrID-Net to be robust (Ap-
DAG construction

---

Published as a conference paper at ICLR 2022
pendix A.2)) to the choice of ls and other approaches for inferring pseudotime, specifically, Palantir Setty et al. 2019) and Monocle (Qiu et al. 2017).
Evaluating causality by using independent information A key contribution of this work is framework for systematically evaluating putative Granger causal peak-gene relationships; to our knowledge, no such framework currently exists. Since ground truth genome-scale perturbational data is unavailable, we propose that the prediction of peak-gene associations be compared against independent, cell type-matched information sources. One information source we use are eQTLs from cell types matched to those profiled in the single-cell studies. An eQTL is a locus-gems pair that associates naturally occurring genetic variability at a locus to gene expression variations ob- served in large populations, thus serving as a proxy for direct perturbation of the locus (Nica & Dermitzakis, 2013). We also evaluate against cell type-matched chromatin interaction (Hi-C75C) data that report pairs of genomic loci that are spatially close in the 3D structure of the genome (Schoenfelder & Fraser 2019). Intuitively, chromatin interaction and accessibility are both mani- festations of the underlying gene regulatory process so that, compared to random yeak-gens pairs, causal pairs should more likely be in close spatial (i.e., Hi-C) proximity.
We note that, in absolute terms, a high overlap of Granger causal peak-gene pairs with these datasets should not be expected. Like ATAC-seq, both eQTLs and chromatin interaction capture the underly- ing biology only indirectly and are themselves prone to errors: eQTL. s are confounded by correlated mutations elsewhere in the genome and Hi-C/5C data suffers from similar sparsity issues as ATAC- seq. Therefore, we focus here on a relative evaluation, comparing GrID-Net's overlap with these datasets against that of alternative approaches. Also, while good overlap with both eQTL and chro- matin interaction datasets is desirable, we believe emphasis should be given to eQTLs since they are the better proxy for an actual perturbation.
Benchmarking We evaluated two correlation- techniques used in previous single-cell multi- modal studics (Cao et al. 2018 Chen et al. 2019; Ma et al. 2020): a) Pearson correlation (Pearson Corr) of peak counts and gene expression, with each cell considered an independent observation; and b) Pearson correlation computed after averaging the RNA-seq and ATAC-seq profiles of each cell and its 50 neighbors in the kNN graph of cells (Pseudocell Corr). The pseudocell approach, currently the state-of-the-art seeks robustness by replacing point estimates in the ATAC-seq (or RNA-seq) space with local averages. We used the absolute value of the correlation scores rank peak-gene pairs in our analyses below. To assess the advantage of operating on a DAG for partially ordered data, we also computed a traditional, linear Granger causal estimate using vector autoregres- sion (VAR Granger) and a state-of-the-ant generalized vector autoregression method for nonlinear Granger causal inference (GVAR). For both these methods, data was forced into a total ordering along pseudotime by partitioning the pseudotime range into 100 bins of equal size and averaging the ATAC-seq and RNA-seq profiles across cells assigned to each bin (Appendix (A.4). (A.44). We note that the Pseudocell Corr method uses information from just the kNN graph structure, while both the VAR Granger and GVAR methods use only pseudotime information. By making use of both the KNN graph and pseudotime information, GrID-Net is able to combine their strengths.
Predicting eQTLs We obtained eQTL data from the GTEx project for the human cell types that most closely correspond to the cells from the single-cell multimodal assays (GTEx Consortium, 2013). For the lung adenocarcinomm--dervell A549 cells in the sci-CAR dataset, we used eQTL data for human lungs. The SHARE-seq and SNARE-seq studies profiled mouse tissues, for which matching eQTL data is limited. Since human and mouse epigenomes show strong conservation (Goneekka et al. 2015) Xiao et al.) 2012), we used eQTL data for human skin and the cerebral cortex, respectively (see Appendix A.5). and mapped this data to the mouse genome. For each dataset, we retained peak-gene pairs with a matching locus- locus-gene gene pair in the eQT'L data.
To assess the effectiveness of GrID-Net in predicting eQTLs,, we labeled peak- gene pairs associated with eQTLs as true (eQTL, p < 10-10) or false (eQTL p > 0.9), discarding pairs that were not in either category. GrID-Net and the alternative methods were evaluated on their accuracy in predict- ing the true eQTLs based the ranking of association scores generated by each method for each peak-gen gene pair. Across all datasets, GrID-Net outperformed the other methods in predicting eQTLs, achieving the highest AUPRC (area under precision-recee curve) and AUROC (area under receiver operating characteristic) in all cases, indicating that peak- -enn pairs prioritized by GrID-Net more closely align with evidence from population genetics (Table T]. The non-linear Granger approach
7

---

Published as a conference paper at ICLR 2022
(GVAR) performs slightly better than linear approach (VAR Granger); however, both methods per- formed substantially worse than GrID-Net, pointing to the advantage of the DAG-based inference approach for such systems over standard Granger-causal inference on forcibly-ordered data.
Table 1: eQTL prediction accuracy on three multimodal datasets
scl-CAR
SHARE-seq
SNARE-seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Pearson Corr
0,0594
0.514
0.0306
0,515
0.0096
0.619
Pseudocell Corr
0.0591
0.512
0.0385
0.541
0,0084
0.573
VAR Granger
0.0542
0.493
0.0321
0.506
0.0037
0.439
GVAR (nonlinear)
0,0556 + 0.0016
0.497 + 0.016
0.0330 + 0.0039
0,536 + 0.049
0,0040
0.0004
01,485 + 0.070
GrID-Net
0.0645 : 0.0004
0.530 + 0.002
0.0604 = 0.0041
0.606 0.019
0.0198 0.0088
0.681 0.006
Predicting chromatin interactions We obtained chromatin interaction data for the specific cell types and experimental settings profiled in the three single-cell multimodal datasets (Davis et al. 2018 Dixon et al. 2012 Poterlowicz et al., 2017). Chromatin interactions were determined for peak-gene gene pair by identifying the pairs of genomic windows in the corresponding chromatin inter- action dataset that overlapped the peak and the transcription start site (TSS) of the gene. We then identified the peak- gene pairs with significant chromatin interactions (Appendix [A.6) and compared the accuracy of GrIDD-Net and alternative methods in predicting these pairs. We found that GrIDD Net achieved the highest AUPRC and AUROC values in this prediction task across all datasets, demonstrating its ability to better identify these key functional features in gene regulation (Table
Table 2: Chromatin interaction prediction accuracy on three multimodal datasets
sci-CAR
SHARE-seq
SN SNARE-Seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Pearson Corr
0.0105
0.502
0.0093
0,356
0,0106
0,513
Pscudocell Corr
0.0103
0.496
0.0113
0.506
0,0108
0,517
VAR Granger
0.0101
0.491
0.0207
0.474
0,0102
0.504
GVAR (nonlinear)
0.0103 - 0.0001
0.509
+ 0.005
0.0108
: 0.0014
0.440 1: 0.124
0.0104 = 0.0003
0.013
GrID--Net
0.0113 + 0.0000)
0.550 + 0.001
0.0402
0.0023
0.753 0.010
0.0115 0.0002
0.540 0.006
GrID-Net is more robust to sparse data than existing approaches We hypothesized that GrID- Net's message passing procedure, which aggregates sparse and noisy measurements across cells using the DAG's structure, would offer it robustness to the sparsity of single-cell multimodal data. Because of differing true-positive rates between the full set of candidate peak-gend pairs and the sparse data subset (Appendix [A.9 Table [5). the relative performance of methods are compared within rather than across settings. We focused on the subset of genes or peaks for which data is especially sparse (fewer than 5% or 2% of cells for genes or peaks, respectively). Evaluating the various methods on the eQTL prediction task on this subset, we found that GrID-Net's outperfor- mance compared to the alternatives was accentuated (Table additional details in Appendix (A.9). GrID-Net also noticeably outperforms the other methods in predicting significant chromatin inter- actions associated with peak-gene pairs consisting of sparse genes or peaks. These results provide further evidence of GrID-Nee's utility in extracting insights from sparse single-cell data.
Table 3: eQTL and chromatin interaction prediction accuracy for sparse genes and peak (AUPRC)
eQTL
Chromatin interactions
sci-CAR
SHARE-seq
SNARE-seq
sci-CAR
SHARE-seq
SNARE-seq
Method
Sparse Genes
Sparse
Sparse
Sparse Peaks
Spurse Genes
Spurse Peaks
Sparse Genes
Sparse
Sparse
Sparse
Sparse
Sparse
Peaks
Genes
Peaks
Genes
Peaks
Genes
Peaks
Peurson Corr
0.0607
0.0566
0.0472
0.0195
0.0136
0.0078
0.0116
0.0106
0.0140
0.0184
0.0102
0.0088
Pseudocell Corr
0.0589
0.0562
0.0567
0.0167
0.0122
0.0058
0.0114
0.0104
0.0116
0.0154
0.0101
0.0089
VAR Granger
0.0494
0.0520
0.0445
0.0168
0.0045
0.0029
0.0115
0.0101
0.0049
0.0664
0.0097
0.0087
(nonlineur]
0.0512
0.0533
0.0496
0.0261
0.0047
0.0027
0.0113
0.0104
0.0068
0.0117
0.0098
0.0091
GrID-Net 0.0668
0.0628
0.0818
0.1450
0.0218
0.0169
0.0121
0.0113
0.0455
0.0596
0.0116
0.0098
Peak-gene pairs prioritized by GrID-Net are supported by TF motif and ChIP-seq data To assess the functional relevance of the detected peak-gene pairs, we sought to relate high-scoring peaks to the transcription factors (TFs) that might bind there. We selected the highest scoring 1% of peak-gena pairs from each of GrID-Net and its alternatives. We next applied Homer (Heizz et al. 2010) to identify TF binding motifs that are co-enriched in sets of peaks linked to particular gene.
8

---

Published as a conference paper at ICLR 2022
The presence of these co-enriched motifs is potentially indicative of coordinated regulatory control, whereby a gene is modulated by multiple TF binding events featuring the same or related sets of TFs (Pott & Lieb] 2015 Carleton et al., 2017). Consequently, the co-occurrence of these motifs serves as a proxy for the functional significance of the proposed peak-gene pairs. Across all datasets, observed that the peak-gene pairs detected by GrID-Net were associated with 10-50 times more enriched TF binding motif's than the other methods (Tablef4), with GrID-Net's outperformance being statistically significant (p < 10-10 one-sided binomial test with Bonferroni correction, details in Appendix [A.10]. These results suggest that GrID-Net may serve valuable role in unveiling mechanistic insights into TF-gene regulatory control.
We further evaluated the functional importance of detected peak-gene gene pairs by comparing peaks in these peak-gene pairs to TF binding sites supported by TF ChIP-seq experimental data (Oki et al.) each single-
2018). We first obtained TF ChIP-seq data for the specific cell type(s) represented cell multimodal dataset. For each of these datasets, then identified the 0.1 % of peaks involved in the highest scoring peak-gene pairs and determined the proportion of such peaks that overlapped with SolP-seq-deriver TF binding sites. We observed that peaks prioritized by GrID-Net more consistently overlapped TF binding sites relative to other methods across all three datasets (Table As TF binding sites are indicative of functional regulatory roles for noncoding loci (Valouev al.) 2008). GrID-Net's effectiveness in prioritizing peaks bound by TT' provides additional evidence for its ability reveal important gene regulatory roles for specific loci in the noncoding genome.
Table 4: Functional relevance of peak-gene pairs. Left: Number of putative TF-target gene relations detected
from top I % of peak-gene pairs. Right: Proportion of prioritized peaks associated with TF binding sites.
Number of putative TF- Targeg gene relations
Proportion of peaks overlapping TF binding sites
Method
sci-CAR
SHARE-seq
SNARE-seq
sci-CAR
SHARE-seq
SNARE-seq
Pearson Corr
2
95
0.630
0.521
0.232
Pseudocell Corr
16
68
23
0.565
0.573
0.261
VAR Granger
45
56
0.652
0.479
0.268
GVAR
2.9
11.9
19,9
+1.1
10.2
113
0,703 20.070
0.528
10.049
0.245 10.033
GrID-Net
=
71
226
- 31
1275
0.862 0.013
0.774 0.012
0.573 0.004
Investigating a link between genomic distance and regulatory control We also explored if larger peak-gene genomic distances correspond to greater temporal lags between peak accessibility and gene expression. By varying (the number of GNN layers), we allowed for greater pseudotime difference between peaks and genes, finding that architectures with larger L did produce a larger proportion of distal (vs. proximal) peak-gene pairs in their top hits (Appendix A.IT Tables I6
DISCUSSION
We extend the applicability of Granger causal analysis to dynamical systems represented as DAGs, introducing GrID-Net perform such inference. GrID-Net takes advantage of the expressive power of graph neural networks, enabling the detection of long-range nonlinear interactions between vari- ables. We focused on applying it to study multimodal single-cell gene regulatory dynamics, given the prevalence of graph-based representations of cellular landscapes the field. GrID-Net demon- strated substantial improvements in accuracy relative current methods in predicting peak-gene links from independent, cell type-matched information sources. Capitalizing on the high resolution of ATAC peaks (~ 1 kb), regulatory links detected by GrID-Net from single-cell multimodal assays may be used to precisely hone in on functional regions within noncoding loci or to complement existing chromatin interaction data towards this goal (Fulco et al. 2019 Nasser et al. 2021). This increased resolution, combined with the genome-wide scale of single-cell multimodal assays, po- sitions GrID-Net to serve as a vital tool for furthering our understanding of key aspects of gene regulatory dynamics relevant to fundamental biology and disease.
Additionally, the ever growing quantities of network data across numerous non-biological domains, including social media networks and financial transaction networks, also points to the broader ap- plicability of GreD-Net. The ability of GrID-Net to capture long-range range, nonlinear interactions in DAGs opens the door to novel analyses on the many existing datasets that are characterized by these graph-based dynamics well.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
We thank the anonymous reviewers of ICLR 2022 for their valuable feedback and suggestions. Fig- ure I was created using biorender.o.. All authors were supported by the NIH grant R35GM141861.
REPRODUCIBILITY STATEMENT
In order to ensure that the full set of results presented in this paper are reproducible, we have in- cluded source code for our method, which can found at pntttps://githhh gridnet In addition, all analyses described in the paper datasets, for which we have provided references throughout the main text. Similarly, the software tools used for biology-specific analyses (e.g., Homer) are publicly available, and we have described the specific settings that were used for this study in the "Appendix" section. We have also de- tailed pre-proue pr-pccessssng steps for the various single-cell datasets in the main text and included both pre-processing and analysis procedures related to the eQTL, Hi-C/SC, and TF binding motif data in the "Appendix" section. We also made
com lexw 6/ performed on publicly available
provide details the specific training details and hyperparameter choices related to training our neural network model in a dedicated subsection in the "Methods" section entitled "Training Details and Hyperarameters".
REFERENCES
Richard Argelaguet, Damien Arnol, Danila Bredikhin, Yonatan Deloro, John C. Marioni, and Oliver Stegle. Mofa+: a statistical framework for comprehensive integration of multi-men modal single-cell data. Genome Biology, 21(111), 2020.
Yael Baran, Akhiad Bercovich, Arnau Sebe-Pedros, Yaniv Lubling, Amir Giladi, Elad Chomsky, Zohar Meir, Michael Hoichman, Aviezer Lifshitz, and Amos Tanay. Metacell: analysis of single- cell ma-seq data using k-nn graph partitions. Genome Biology, 20(206), 2019.
Sumanta Basu, Ali Shojaie, and George Michailidis. Network granger causality with inherent group- ing structure. The Journal of Machine Learning Research, 16(1):417-453, 2015.
Yoshua Bengio and Xavier Glorot. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of AISTATS 2010, 256, 2010.
François Benhmad. Modeling nonlinear Granger causality between the oil price and US dollar: A wavelet based approach. Economic modelling, 2944) 2012.
Kay H Brodersen, Fabian Gallusser, Jim Koehler, Nicolas Remy, and Steven I. Scott. Inferring causal impact using bayesian structural time-series models. The Annals of Applied Statistics, (1):247-274, 2015.
Junyue Cao, Darren A. Cusanovich, Vijay Ramani, Delasa Aghamirzaie, Hannah A. Pliner, An- drew J. Hill, Riza M. Daza, Jose L.
Jonathan S. Packer, Lena Christiansen, Frank J. Steemers, Andrew C. Adey, Cole Trapnell, and Jay Shendure. Joint profiling of chromatin accessibility and gene expression thousands of single cells. Science, (640999)) 380-1385, 2018.
K. Cao, X. Bai, Hong, and L. Wan. Unsupervised topological alignment for single-cell multi- omics integration. Bioinformatics, 36(Suppl1) i48 -156, 072020.
Julia B. Carleton, Kristofer C. Berrett, and Jason Gertz. Multiplex enhancer interference reveals col- laborative control of gene regulation by estrogen receptor a-bound enhancers. Cell Systems, 5(4): 333-344, 2017.
Song Chen, Blue Lake, and Kun Zhang. High-throughput sequencing of the transcriptome and chromatin accessibility in the same cell. Nature Biotechnology, 37:1452-1457, 2019.
Darren A. Cusanovich, Riza Daza, Andrew Adey, Hannah A. Pliner, Lena Christiansen, Kevin L. Gunderson, Frank J. Steemers, Cole Trapnell, and Jay Shendure. Multiplex single cell profiling of chromatin accessibility by combinatorial cellular indexing. Science, 348(6237).
10

---

Published as a conference paper at ICLR 2022
Darren A. Cusanovich, Andrew J. Hill, Delasa Aghamirzaie, Riza M. Daza, Hannah A. Pliner, Joel B. Berletch, Galina N. Filippova, Xingfan Huang, Lena Christiansen, William S. DeWitt, Choli Lee, Samuel G. Regalado, David F. Read, Frank J. Steemers, Christine M. Disteche, Cole Trapnell, and Jay Shendure. A single-cell atlas of in vivo mammalian chromatin accessibility. Cell, 174(5), b.
Darren A. Cusanovich, James P. Reddington, David A. Garfield, Riza M. Daza, Delasa Aghamirzaie, Raquel Marco-Ferrerge Hannah A. Pliner, Lena Christiansen, Xiaojie Qiu, Frank J. Steemers, Cole Trapnell, Jay Shendure, and Eileen E. M. Furlong. The cis-regulatory dynamics of embryonic de- velopment at single-cell resolution. Nature, 555(7697),
Carric A. Davis, Benjamin C. Hitz, Cricket A. Sloan, Esther T. Chan, Jean M. Davidson, Idan Gab- dank, Jason A. Hilton, Kriti Jain, Ulugbek K. Baymuradov, Aditi K. Narayanan, Kathrina C. Onate, Keenan Graham, Stuart R. Miyasato, Timothy R. Dreszer, J. Seth Strattan, Otto Jolanki, Forrest Y. Tanaka, and J. Michael Cherry. The encyclopedia of dna elements (encode): data portal update. Nucleic Acids Research, 46(Database issue):D794-8801, 2018.
B. DeMeo and B. Berger. Hopper: a mathematically optimal algorithm for sketching biological data. Bioinformatics, 36(Suppl)) : 1236 - -i241, 072020.
Jesse R. Dixon, Siddarth Selvaraj. Feng Yue, Audrey Kim, Yan Li, Yin Shen, Ming Hu, Jun S. Liu, and Bing Ren. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature, 485(7398) 2012.
S. Doss, F. E. Schadt, T. A. Drake, and J. Lusis. Cis-acting expression quantitative trait loci in mice. Genome Res, 15(5):681 1000)68116999 1-99,, May 2005.
Doris Entner and Patrik O Hoyer. On causal discovery from time series data using FCI. Probabilistic graphical models, pp. 121-128, 2010.
J. Finkle, J. J. Wu, and Bagheri. Windowed Granger causal inference strategy improves discovery of gene regulatory networks. Proc Natl Acad Sci U A, (1509)22522-2557) 02 2018.
Andr Fujita, Joo R. Sato, Miguel Garay, Rui Yamaguchi, Satoru Miyano, Mari Sogayar, and Car- los E. Ferreira. Modeling gene expression regulatory networks with the sparse vector autoregressive model. BMC Systems Biology, 1(39), 2007.
Charles P. Fulco, Mathias Munschauer, Rockwell Anyoha, Glen Munson, Sharon R. Grossman, Eliz- abeth M. Perez, Michael Kane, Brian Cleary, Eric S. Lander, and Jesse M. Engreitz. Systematic mapping of functional enhancer-promoter connections with crispr interference. Science, 354(6313): 769-773, 2016.
Charles P. Fulco. Joseph Nasser, Thouis R. Jones, Glen Munson, Drew T. Bergman, Vidya Subrama- nian, Sharon R. Grossman, Rockwell Anyoha, Benjamin R. Doughty, Tejal A. Patwardhan, Tung H. Nguyen, Michael Kane, Elizabeth M. Perez, Neva C. Durand, Caleb A. Lareau, Flena K. Stamen- ova, Frez Lieberman Aiden, Eric S. Lander, and Jesse M. Engreitz. Activity-by--conacc model of unhancer-promote regulation from thousands of crispr perturbations. Nature Genetics, 51:1664- 1669, 2019.
Adam Gayoso, Zoö Steier, Romain Lopez, Jeffrey Regier, Kristopher L. Nazor, Aaron Streets, and Nir Yosef. Joint probabilistic modeling of single-cell multi-omic data with totalvi. Nature Methods, 18: 272-282, 2021.
E. Gjoneska, A. R. Pfenning, H. Mathys, G. Quon, A. Kundaje, L. H. Tsai, and M. Kellis. Conserved epigenomic signals in mice and humans reveal immune basis of Alzheimer's disease. Nature, 518 (7539):365-369, Feb 2015.
C. Glymour, Zhang, and Spirtes. Review of Causal Discovery Methods Based on Graphical Models. Front Genet, 10:524, 2019.
Jairo Navarro Gonzalez, Ann S. Zweig, Matthew L. Speir, Daniel Schmelter, Kate R. Rosenbloom, Brian J. Raney, Conner C. Powell, Luis R. Nassar, Nathan D. Maulding, Christopher M. Lee, Brian T. Lee, Angie S. Hinrichs, Alastair C. Fyfe, Jason D. Fernandes, Mark Diekhans, Hiram Clawson, Jonathan Casper, Anna Benet-Pagès. Galt P. Barber, David Haussler, Robert M. Kuhn, Maximilian Haeussler, and W. James Kent. The ucse genome browser database: 2021 update. Nu- cleic Acids Research, 49:D1046-D1057 2021.
11

---

Published as a conference paper at ICLR 2022
Clive W. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: Journa of the Econometric Society. pp. 424 438, 1969.
GTEX Consortium. The genotype-tissue expression (gtex) project. Nature Genetics, 45(6), 2013.
Laleh Haghverdi, Maren Büttner, F. Alexander Wolf, Florian Buettner, and Fabian J. Theis. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods, 13:845-848, 2016.
Sven Heinz, Christopher Benner, Nathanael Spann, Eric Bertolino, Yin C. Lin, Peter Laslo, Jason X. Cheng, Cornelis Murre, Harinder Singh, and Christopher K. Glass. Simple combinations of lineage- determining transcription factors prime cis-regulatory elements required for macrophage and b cell identities. Molecular Cell, 38:576-589, 2010.
Brian Hie, Hyunghoon Cho, Benjamin DeMeo, Bryan Bryson, and Bonnie Berger. Geometric sketch- ing compactly summarizes the single-cell transcriptomic landscape. Cell Systems, 14483-493, 2019.
Brian Hie, Joshua Peters, Sarah K Nyquist, Alex Shalek, Bonnie Berger, and Bryan D Bryson. Computational methods for single-cell ma sequencing. Annual Review of Biomedical Data Science, 3:339-364, 2020.
J. Hofstetter, D. A. Svihla-Jones, and A. R. Mayeda. A QTL on mouse chromosome 12 for the ge- netic variance in free-running circadian period between inbred strains of mice. Circadian Rhythms, 5:7, Oct 2007.
Saurabh Khanna and Vincent Y.F. Tan. Economy statistical recurrent units for inferring nonlinear granger causality. In International Conference on Learning Representations, 2020.
David Lara-Astiaso, Assaf Weiner, Erika Lorenzo- Verenzo-Vivas, Irina Zaretsky. Diego Adhemar Jaitin, Eyal David, Hadas Sereu,SShal,, Alexander Mildner, Deborah Winter, Steffen Jung. Nir Friedman, and Ido Amit. Chromatin state dynamics during blood formation. Science. 345(6199):943-949, 2014. doi: 10.1 126/science. 1256271.
Longgi Liu, Chuanyu Liu, Andrés Quintero, Liang Wu, Yue Yuan, Mingyue Wang, Mengnan Cheng, Lizhi Leng, Liqin Xu, Guoyi Dong, Rui Li, Yang Liu, Xiaoyu Wei, Jiangshan Xu, Xiaowei Chen, Haorong Lu, Dongsheng Chen, Quanlei Wang. Qing Zhou, Xinxin Lin, Guibo Li, Shiping Liu, Qi Wang, Hongru Wang, Lynn Fink, Zhengliang Gao, Xin Liu, Yong Hou, Shida Zhu, Huan- ming Yang, Yunming Ye, Ge Lin, Fang Chen, Carl Herraan, Roland Bils, Zhouchun Shang, and Xun Xu. Deconvolution of single-cell multi-omics layers reveals regulatory heterogeneity. Nature Communications, 10:470, 2019.
Helmut Lüttepphl. New Introduction to Multiple Time Series Analysis. MIT Press, 2005.
Sai Ma, Bing Zhang, Lindsay M. LaFave, Andrew S. Earl, Zachary Chiang, Yan Hu, Jiarui Ding, Alison Brack, Vinay K. Kartha, Law Travis Tay, Tristan, Caleb Lareau, Ya-Chieh Hsu, Aviv Regev, and Jason D. Buenrostro. Chromatin potential identified by shared single-cell profiling of ma and chromatin. Cell, 103 2020.
Ričards Marcinkeviis and Julia Vogt. Interpretable models for granger causality using self- explaining neural networks. In 9th International Conference on Learning Representations, ICLR 2021, 2021.
Daniele Marinazzo, Mario Pellicoro, and Sebastiano Stramaglia. Kernel method for nonlinear granger causality. Physical Review Letters, 100(14103), 2008
Atalanti Mastakouri, Bernhard Schölkopf, and Dominik Janzing. Necessary and sufficient conditions for causal feature selection in time series with latent common causes. In International Conference on Machine Learning, pp. 7502-7511. PMLR, 2021.
Alessandro Montalto, Sebastiano Stramaglia, Luca Faes, Giovanni Tessitore, Roberto Preveteg, and Daniele Marinazzoa. Neural networks with non-uniform embedding and explicit validation phase assess granger causality. Neural Networks, 71:159-171, 2015.
12

---

Published as a conference paper at ICLR 2022
Jill E. Moore, Henry E. Pratt, Michael J. Purcaro, and Zhiping Weng. A curated benchmark of unhancer-gene interactions for evaluating enhancer- target gene prediction methods. Genome Bi- ology, 20119,2019
Joseph Nasser, Drew Bergman, Charles P. Fulco, Philine Guckelberger, Benjamin R. Doughty, Te- jal A. Patwardhan, Thouis R. Jones, Tung H. Nguyen, Jacob C. Ulirsch, Fritz Lekschas, Kristy Mualim, Heini M. Natri, File M. Weeks, Glen Munson, Michael Kane, Helen Y. Kang, Ang Cui, John P. Ray. Thomas M. Fisenhaure, Ryan I.. Collins, Kushal Dey, Hanspeter Pfister, Alkes I.. Price, Charles B. Epstein, Anshul Kundaje, Ramnik J. Xavier, Mark J. Daly, Hailiang Huang, Hilary K. Finucane, Nir Hacohen, Eric S. Lander, and Jesse M. Engreitz. Genome-wide enhancer maps link risk variants to disease genes. Nature, 593:238 243, 2021.
Meike Nauta, Doina Bucur, and Christin Seifert. Causal discovery with attention-based convolutional neural networks. Machine Learning and Knowledge Extraction, 1:312-3444 2019.
Alexandra C Nica and Emmanouil T. Dermitzakis. Expression quantitative trait loci: present and future. Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 368(162))2012036,, 2013.
Shinya Oki, Tazro Ohta, Go Shioi, Hideki Hatanaka, Osamu Ogasawara, Yoshihiro Okuda, Hideya Kawaji, Ryo Nakaki, Jun Sese, and Chikara Meno. Chip-atlas: data-mining suite powered by full integration of public chip-seq data. EMBO Reports, 19:e46255, 2018.
Renato Ostuni, Viviana Piccolo, Iros Barozzi, Sara Polletti, Alberto Termanini, Silvia Bonifacio, Alessia Curina, Elena Prosperini, Serena Ghisletti, and Gioacchino Natoli. Latent enhancers ac- tivated by stimulation differentiated cells. Cell, 1-2:157-171, 2013.
Anil Panigrahi and Bert W. O' O'Malley. Mechanisms of enhancer action: the known and the unknown. Genome Biology, 22(108), 2021.
Niklas Pfister, Peter Bühlmann, and Jonas Peters.
Invariant causal prediction for sequential data. Journal of the American Statistical Association, 1 114(52)):1244-12766 2019.
Krzysztof Poterlowicz, Joanne L. Yarker, Igor Malashchuk, Brian R. Lajoie, Andrei N. Mardaryev, Michal R. Gdula, Andrey A. Sharov, Terumi Kohwi-Shigematsu, Vladimir A. Botchkarev, and Michael Y. Fessing. 5c analysis of the epidermal differentiation complex locus reveals distinct chromatin interaction networks between gene-rich and gene-poor tads skin epithelial cells. PLoS Genetics, 13(9); 1006966, 2017.
Sebastian Pott and Jason D. Lieb. What are super-enhancers? Nature Genetics, 47:8-12, 12,8-1,,00 2015.
X. Qiu, Q. Mao, Y. Tang, L. Wang, R. Chawla, H. A. Pliner, and C. Trapnell. Reversed graph embed- ding resolves complex single-cell trajectories. Nat Methods, 14(10):979-98) 982, Oct 2017.
Hafsa Rasheed and Muhammad Tahir. Fdi and terrorism: Co-integration & granger causality. Interna- tional Affairs and Global Strategy, 4(0):155, 2012.
Weijie Ren, Baisong Li, and Min Han. A novel granger causality method based on hsic-lasso for revealing nonlinear relationship between multivariate time series. Physica A: Statistical Mechanics and its Applications, 541 2020.
Rahul Satija, Jeffrey A. Farrell, David Gennert, Alexander F. Schier, and Aviv Regev. Spatial recon- struction of single-cell gene expression data. Nature Biotechnology. 33:495-502, 2015.
2015.
Stefan Schoenfelder and Peter Fraser. Long-range control. Nature Reviews Genetics, 20:437-455. 2019.
promoter contacts in gene expression
Skipper Seabold and Josef Perktold. statsmodels:
Econometric and statistical modeling with python. Proceedings of the 9th Python in Science Conference, 2010.
M. Setty, V. Kiseliovas, J. Levine, A. Gayoso, L. Mazutis, and D. Pe' er. Characterization of cell fate probabilities in single-cell data with Palantir. Nat Biotechnol, 37(4):451 -460, 04 2019.
13

---

Published as a conference paper at ICLR 2022
A. Shojaie and G. Michailidis. Discovering graphical Granger causality using the truncating lasso penalty. Bioinformatics, 26(18)i517-522 Sep 2010.
Ali Shojaie and Emily Fox. Granger causality: A review and recent advances. arXiv preprint Sxxiv21005006655 2021.
Rohit Singh, Brian Hie, Ashwin Narayan, and Bonnie Berger. Schema: metric learning enables inter-
pretable synthesis of heterogeneous single-cell modalities. Genome Biology, 22(131), 2021.
P. Spirtes and K. Zhang. Causal discovery and inference: concepts and recent methodological ad- vances. Appl Inform (Berl), 3:3, 2016.
Alex Tank, Ian Covert, Nick Foti. Ali Shojaie, and Emily B. Fox. Neural granger causality for nonlinear time series. arXiv preprint arXiv:/802.05842 2018.
Anton. Valouev, David S. Johnson, Andreas Sundquist, Catherine Medina, Elizabeth Anton, Scrafim Batzoglou, Richard M. Myers, and Arend Sidow. Genome-wide analysis of transcription factor binding sites based on chip-seq data. Nature Methods, 5:829-834, 2008.
S. Xiao, D. Xie, X. Cao, P. Yu, X. Xing, C. C. Chen, M. Musselman, M. Xie, F. D. West, H. A. Lewin, T. Wang, and S. Zhong. Comparative epigenomic annotation of regulatory DNA. Cell, 149(6): 1381-1392 Jun 2012.
Shun Yao, Shinjae Yoo, and Dantong Yu. Prior knowledge driven granger causality analysis on gene regulatory network discovery. BMC Bioinformatics, 16(273), 2015.
Xiuwei Zhang, Chenling Xu, and Nir Yosef. Simulating multiple faceted variability in single cell ma sequencing. Nature Communications, 10(2611), 2019.
Chenxu Zhu, Miao Yu, Hui Huang, Ivan Juric, Armen Abnousi, Rong Hu, Jacinta Lucero, M. Margarita Behrens, Ming Hu, and Bing Ren. An ultra eightthroughpt method for single-cell joint analysis of open chromatin and transcriptome. Nature Structural Molecular Biology, 36:1063-1077 2019.
A APPENDIX
A.1 DAG CONSTRUCTION DETAILS
For the three single-cell multimodal datasets, we applied a standard log(I+CPM/10) transformation the raw RNA-seq counts data and log( 1/00) transformation to the raw ATAC-seq counts data, where CPM indicates counts per million. Genes and accessible regions detected in fewer than 0.1% of cells in each dataset were excluded. For each dataset, the expression of each gene was scaled to unit variance and zero mean, after which the top 2000 most highly variable genes (Satija al. 2015) were retained. In the original SHARE-seq and SNARE-seq studies, topic modeling- based approaches were used to transform ATAC-seq counts data for downstream trajectory analyses, prompting us to transform the ATAC-seq counts in these datasets using if-idf, a procedure derived from text-based topic modeling that has been applied extensively to analyze single-cell ATAC-seq data (Cusanovich et al.
The normalized RNAA-sq and ATAC-seq data were then projected onto their top 100 principal com- ponents (excluding components that have Spearman correlation p 0.9 with the total counts per cell). Schema (Singh et al., 2021) was then used to generate unified representations for each cell that synthesize information from both the RNA-seq and ATAC-seq representations. For the sci-CAR and SNARE-Seq datasets, RNA-seq was used as the primary modality in Schema with ATAC-seq acting as secondary modality. For the sci-CAR dataset, the time of data collection post-drug treat- ment was used as tertiary modality. Meanwhile, ATAC-seq served as the primary modality for the SHARE-sey dataset, as the original SHARE-seq study identified strong cell cycle effects repre- sented in the RNA-seq data for certain subpopulations of cells and consequently used ATAC-seq as the modality for analyzing lineage dynamics for the dataset (Ma et al. 2020). We set the minimum correlation between the primary and the remaining modalities to be 0.9 in Schema. For each dataset, a KNN graph (A: = 15) based on distances between the unified representations generated by Schema was constructed. Pseudotime stamps were then inferred for each cell based on this kNN graph using diffusion pseudotime (Haghverdi et al., 2016). Finally, edges in the kNN graph that cohere with the direction of increasing pseudotime are retained to generate DAG.
14

---

Published as a conference paper at ICLR 2022
A.2 DAG ROBUSTNESS ANALYSIS
In addition to the diffusion pseudotime algorithm, we also evaluated Palantir and Monocle for in- ferring pseudotime. The performance of GrID-Net was robust to the choice of algorithm (Tables
Table 5: eQTL prediction accuracy by pscudotime algorithm
sci-CAR
SHARE--eq
SNARE-Seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Diffusion
0.065 + 0.0004
0.530 +: 0.0019
0.060 + 0041
0.606 + 0.0194
0.020 + 0.0088
0,681 +: 0.0055
Palantir
0.064 + 0.0004
0.530 + 0.0018
0.060 = 0.0028
0.611 + 0.0083
0.019 + 0.0068
0.682 + 0.0106
Monocle
0.065 1: 0.0004
0.531 1: 0.0019
0.061 4 0.0037
0.605 1 0.0266
0.018 1 0.0068
0.675 1 0.0099
Table Chromatin intcraction prediction accuracy by pscudotime algorithm
sci-CAR
SHARE-sQq
SNARE-Seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Diffusion
0.011 + 0,0000
0.550 +: 0.0008
0.040 + 0.0023
0.753 + 0.0105
0.011 0.0002
0.540 0.0060
Palantir
0.011 + 0.0000
0.550 + 0.0009
0.057 = 0.0193
0.783 + 0.0373
0.011 0.0001
0.543 0.0039
Monocle
0.011 + 0.0000
0.548 + 0.0010
0.057 + 0.0209
0.793 + 0.0525
0.011 + 0.0002
0.544 + 0.0052
Table 7: cQTL predictio accuracy by k in kNN -gaph construction
sci-CAR
SHARE-seq
SNARE-SQ
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
km30
0.0645 + 0.0004
0.5303 0.0019
0.0630 0.0051
0.6151 + 0.0162
0.0196
0.0087
0.6788 + 0.0101
k=50
0.0645 + 0.0004
0.5302 0.0017
0.0629 + 0.0069
0.6159 + 0.0160
0.0198
0.0088
0.6778 + 0.0108
k=15
0.0645 + 0.0004
0.5299 = 0.0019
0.0604 1: 0.0041
0.6060 + 0.0194
0.0198 0.0088
0.6805 + 0.0055
A.3 SINGLE-CELL DATA PREPROCESSING AND CANDIDATE PEAK-GENE PAIR SELECTION
We applied the same CPM and log transformations as above to the raw RNA-seq and ATAC-seq counts data for all three single-cell multimodal datasets. We then scaled the RNA-seq data for each gene by the maximum value for that gene, so that expression values for all genes in each dataset in the same [0,1] range. The peak gene pairs considered for each dataset were determined by selecting genes and peaks detected in greater than 0. % of cells and then retaining all peak-gene pair combinations in which the peak is within 1 Mb of the gene. For the SHARE-seq dataset, we only considered peaks detected in greater than 1% of cells so as to have a comparable set of candidate peak-gene pairs across the three datsets.
A.4 BENCHMARKING
For the VAR Granger model, peak-gene pairs were scored on the - 10) as- sociated with the statistical significance of Granger causality, as calculated using the waangercaassllttytess function with default settings in the stat smodels package (Seabold & Perktold 2010).
A.5 MATCHING EQTLS WITH SINGLE-CELL ATAC-SEQ RNA-SEQ DATA
The procedure to match peak-gene pair with a eQTL variant-gene pair involved identifying all catalogued variants positioned in the peak of interest and selecting the variant with the most statis- tically significant association with the gene in the peak-gene pair. To perform this procedure for the SNARE-seq and SHARE-seq datasets, we used liftOver (Gonzalez al. 2021) with default settings to convert human genomic coordinates represented in the GTEx datasets to mouse coordi- nates. We chose to use human eQTL data here because it is important to obtain tissue-specific eQTL data. Unfortunately, there is relatively little eQTL data available in mouse and, to our knowledge, when available is organized around different strains of mice rather than tissues (Doss et al. 2005 Hofstetter al.) 2007). On the other hand, the transfer of eQTL knowledge from human to mouse is well-motivated because the epigenome is largely conserved across mouse and human GGoonskkl al. 2015 Xiao al. 2012).
15

---

We compared the performance of the two tests on the eQTL and chromatin interaction prediction tasks, finding that they yielded very similar results (Tables I0). This similarity is likely due to our use of the test statistics rank peak-gent gene pairs, upon which the top-ranking pairs are selected. Published as a conference paper at ICLR 2022
Table 8: Chromatin interaction prediction accuracy by ls in kNN-graph construction
sci-CAR
SHARE-seq
SNARE-seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
k=30
0.0113 + 0,0000
0.5499
0.0009
0,0362 + 0.0036
0.7233 + 0.0264
0.0114 0.0002
0.5408 + 0.0053
km50
0.0113 + 0.0000
0.5503
0.0009
0.0351 0.0046
0.7152 : 0.0272
0.0113 0.0002
0.5373 + 0.0057
k=15
0.0113 + 0.0000
0.5497 = 0.0008
0.0402 + 0.0023
0.7527 + 0.0105
0.0115 0.0002
0,5400 + 0.0060
A.6
CHROMATIN INTERACTION DATA PREPROCESSING
For the sci-CAR, SHARE-seq, and SNARE-Seq datasets, the corresponding chromatin interaction datasets contained scores indicating the frequency of interactions between genomic windows. We retained peak-gene pairs that were matched with the genomic windows presented in the chromatin interaction datasets and chose the most frequent 1% interactions as the set to be predicted.
A.7 COMPARISON OF F-TEST AND WELCH'S S T-TEST TO RANK PUTATIVE PEAK-GENT INTERACTIONS
For each candidate peak-gene gene pair, we compare the full and reduced models using an F-test, as is common in Granger causal inference. Specifically, we evaluate the sum--of-quaree-errro of the reduced and full models as per an F distribution with (2L + 1, 4L - 1) degrees of freedom. Here, is the number of cells, each of the L GNN layers includes two parameters (the full model has two such architectures while the reduced model has one) and is an additional parameter in the full model. Another approach would be to compare the full and reduced models' losses with the Welch's 's t-test, where the alternative hypothesis is that the mean of the loss terms for the full model less than the mean of the loss terms for the reduced model. In both approaches, candidate interactions are then ranked by the test statistic with higher scores corresponding to a stronger likelihood of Granger causal interaction.
Since the degrees of freedom are identical across all candidate pairs, their ranking depends only on the relative losses of the full and reduced models.
Table 9: Comparison of F-test and t-test on eQTI. prediction accuracy
sci-CAR
SHARE-seq
SNARE-Seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Hitest
0.0645 + 0.0004 0.5300 + 0.0020
0.0604
0.0041 0.6060 + 0.0190
0.0198 + 0.0088
0.6810 + 0.0060
(-test
0.0697 + 0.0006 0,5400 + 0.0020
0.0573 0.0047
0.6100 + 0.0180
0.0175 + 0.0062
0.6770 :: 0.0080
Table 10: Comparison of F-test and (-test chromatin interaction prediction accuracy
Sci-CAR
SHARE-seq
SNARE-Seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
F-test
0.0113 # 0.0000
0.5500 + 0.0010
0.0402 = 0.0023
0.7530 + 0.0100
0.0115 + 0.0002 0,5400 1: 0.0060
f-test
0.0109 = 0.0000
0,5460 + 0.0010
0.0375 0.0066
0.7620 + 0.0110
0.0115 0.0003
0.5370 +: 0.0070
A.8 COMPARISON OF COMBINED AND SEPARATED LOSS FUNCTIONS FOR FULL AND REDUCED MODELS
During the GNN's training, the sum of the losses of the full and reduced models is optimized. Since the two models do not share any parameters, optimizing the sum of their losses corresponds to optimizing each loss separately. We chose this joint-optimization strategy for efficiency and to avoid stochastic variations if the models were trained separately. Nonetheless, to confirm that optimizing the sum of losses does not introduce any unintended dependencies, we also evaluated GrID-Net's performance on the eQTL and chromatin interaction prediction tasks when the models were optimized in completely separate training runs. The results for separate and joint training were essentially identical (Tables 11 12. suggesting that our joint training works intended.
16

---

Published as a conference paper at ICLR 2022
Table 11: Comparison of combined and separated loss functions on eQTI. predictiin accuracy
sci-CAR
SHA
SNARE-seg
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Separuted
0.0645 + 0.0004
0.530 = 0.002
0.0603 + 0.0043
0.605 + 0.020
0.0201 + 0.0085
0.684 + 0.006 0.681 = 0.006
Combined
0.0645 + 0.0004
0.530 = 0.002
0.0604 + 0.0041
0.606 :: 0.019
0.0198 + 0.0088
Table 12: Comparison of combined and separated loss functions on chromatin interaction prediction accuracy
sci-CAR
SHARE--eq
SNARE-seq
Method
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
Separated
0.0113 + 0.0000
0.550 = 0.001
0.0404 + 0.0019
0.752 + 0.010
0.0115 + 0.0002
0.539 + 0.006
Combined
0.0113 + 0.0000
0.550 = 0.001
0.0402 0.0023
0.753 0.010
0.0115 + 0.0002
0.540 + 0.006
A.9 ROBUSTNESS OF GRID-NET TO SPARSITY
Tables 13 and 14 expand upon the results in Table also showing standard deviations (across three runs) of GVAR's and GrID-Net's s AUPRC results; the other methods are fully deterministic and their results do not change across runs. Table 13 shows the baseline proportions of true positives across the different comparison conditions to provide context for differences in the range of AUPRC values across these conditions.
Table 13: eQTL prediction accuracy for sparse genes and peak (AUPRC). Accompanies Table 3
sei-CAR
SHARE-sqq
SNARE-seq
Method
Sparse Genes
Sparse Posks
Sparse Genes
Sparse Peaks
Sparse Genes
Sparse Poaks
Pearson Cor
0,0607
0,0566
0,0472
0.0195
0.0136
0,0078
Pseudocell Corr
0.0589
0.0562
0.0567
0.0167
0.0122
0.0058
VAR Granger
0.0494
0,0520
0,0445
0.0168
0.0045
0,0029
GVAR (nonlinear)
0.0512
0.0014
0.0533
0.0015
0.0496
0.0049
0.0261
0.0040
0.0047
0.0004
0.0027
0.0002
GrID-Net
0.0668
0.0008
0.0628 0.0004
0,0818 0.0051
0,1450
0.0310
0,0218 0.0096
0.0169
0.0082
A. TF BINDING MOTIF ENRICHMENT ANALYSTS
For each gene that was associated with more than 5 peaks in the highest scoring 1 1% peak gene pairs, we applied HOMER with default parameters on the genomic windows defined by the set of peaks to test for enrichment of TF binding motifs (Heinz et al. 2010). We used a Benjamini- Hochberg-correcte
p-value cutoff of 0.01 to call enriched TF binding motifs. Against each baseline, the statistical significance of the difference between GrID-Net and the baseline in the number of enriched motifs was assessed via a one-sided binomial test, with Bonferroni correction for multiple hypothesis testing. All comparisons were significant at the p 10 10 level.
A.11 GRID-NET SUGGESTS POSSIBLE LINK BETWEEN GENOMIC DISTANCE AND TEMPORAL REGULATORY CONTROL
GrID-Net's explicit consideration of the temporal asynchrony between chromatin accessibility and gene expression enables the genome-wide study of lags in gene regulation for the very first time. We classified peak-gene pairs as proximal (genomic distance < 10kb) or distal (2 100 kb). We then evaluated different architectures of GrID-Net that vary in the number of GNN layers hy- pothesizing that architectures with higher L (i.e., allowing greater pseudotime difference between peaks and genes) may detect more distal peak-gene pairs. Inspecting the top 1 1% of greak-gene pairs detected by different architectures, we found this to indeed be the case: across all datasets. archi- tectures with higher detected larger ratios of distal-to- -to--ooxiaaal peak- gene pairs (Tables 16 17). suggesting that distal peak-gene interactions may be marked by longer pseudotime difference be- tween peak and genes. This finding suggests the possibility of previously unreported link between genomic distance and the temporal lag between a regulatory element's accessibility and its target gene's 's expression. However, like other Granger causal analyses, our approach also can not account for hidden confounders; this is especially pertinent here since the greater pseudotemporal lag be- tween distal peaks and genes may be mediated by other regulatory mechanisms rather than being direct influence. One hypothesis is that the complex coordination of regulatory factors required
17

---

Published as a conference paper at ICLR 2022
Table 14: Chromatin interaction prediction accuracy for sparse genes and peak (AUPRC). Accompanies Table
sci-CAR
SHARE
SNARE-seg
Method
Sparse Genes
Sparse
Sparse Gienes
Sparse
Sparse Gieee
Sparse
Pearson Corr
0.0116
0.0106
0.0140
0.0102
Pscadocell Corr
0.0114
0.0104
0.0116
0.0154
0.0101
0.0089
Granger
0.0049
GVAR (nonlincar)
0.0113 + 0.0000
0,0104 + 0.0001
0,0068
0.0010
0.0117 + 0.0034
0.0098 + 0.0003
0.0091 + 0.0003 0.0098 0.0002
GrID-Net
0.0121 0.0001
0.01 0.0000
0.0455
0.0144
0.0596 0.0177
0.01 + 0.0002
Table 15: Baseline proportions of true positives in full, sparse gene, and sparse peak comparison conditions.
sci-CAR
SHARE-seq
SNARE-seq
Category
<QTL
Chrom Int
eQTL
Chrom Int
eQTL
Chrom Int
Full
0.0558
0.0101
0.0302
0.0111
0.00401
0,0100
Sparse Genes
0.0515
0.01
0.0452
0.0071
0.00488
0.0095
Sparse Peaks
0.0535
0.0102
0.0192
0.0111
0.00284
0.0085
enable distal regulatory interactions is associated with greater lags between distal regulatory ele- ments' accessibility and the gene expression changes they affect (Panigrahi & O O'Malley 2021).
Table 16: Percentages of peak gene pairs within 10kb
sci-CAR SHARE-seq SNARE-seq
2 GNN Layers 14.5%
6.04%
2.53%
GNN Layers
14.1%
5.20%
2.35%
10 GNN Layers
13.7%
5.33%
2.25%
A.12 RUNTIME AND MEMORY USAGE
GrID-Net's runtime and memory usage are dependent number of nodes the DAG describing the system interest, In our model, the DAG is represented by the N x N adjacency matrices A and A+, where is the number of nodes in the DAG (i.e. number of cells in a multimodal single-cell dataset). These matrices serve as the primary bottlenecks for both runtime and memory usage. Memory usage scales quadratically with n, while calculations involving A or A+ can be parallelized on a GPU to reduce the runtime. For very large single- cell datasets, 2020) so that the matrices will fit into GPU memory. In addition, training time scales linearly with the number of candidate peak-gene pairs. In our results above, we evaluated 507,408 916,939, and 872,039 candidate peak gene pairs for the sci-CAR, SHARE-seq, and SNARE-seq datasets, respectively. For these same datasets, we considered systems of 3,260, 5,000, and 5,081 cells. The average total runtime for GrID-Net was 1.8 hours, 5.5 hours, and 5.3 hours, while peak memory usage was approximately 31 GB, GB, and 39 GB for these respective datasets.
variety of factors, most notably the
recommend the of sketching techniques (Hie et al.] 2019 DeMeo & Berger)
18

---

Published as a conference paper at ICLR 2022
Table 17: Ratio of distal (> 100kb) to proximal 10kb) peak-gene pairs in the top hits
sci-CAR
SNARE-seq SHARE-seq
2 GNN Layers
6.10 + 0.062
12.53
0.128
35.97 + 1.322
5 GNN Layers
6.23 + 0.119
15.33 + 0.522
37.45 + 1.774
10 GNN Layers
6.05 + 0.055
15.77 0.150
41.16 + 1.228
19


Reviews:
[{"review_id": "rpJ6OzuH7OQ", "reviewer": "Reviewer_dibr", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Update:** After reading the rebuttal and other reviews, I have increased my original score. Thank you to the authors for addressing the concerns raised in great detail.\n\n**Strengths:**\n- The problem tackled by the authors is, in my opinion, an exceedingly difficult one.\n- The experimental section is extensive and the authors have done their best to analyze different facets of the problem. \n- The paper is very carefully written. \n\n**Weaknesses:**\n- The usefulness of the method heavily relies on how good the DAG construction is. In this paper, the authors construct a kNN graph to define the DAG, which seems a reasonable choice, but then how sensitive is the overall method to this choice?\n- Without confidence bands, many of the results presented are rather speculative. For example, based on the results reported in Table 5, the authors speculate that distal peak-gene interactions may be characterized by longer lags between peaks and genes. While this conclusion seems reasonable, it is rather hard to gauge whether there is a singificant decrease in proximal peak-gene pairs when going from two to ten GNN layers.\n- The authors make a decent attempt at putting their work into context, but if they want to focus on the causal aspect, there is a whole literature on causal Bayesian network inference and structural causal models that seems to be missing in the reference section. Apart from that, some other papers that seem relevant:\n   - https://www.jmlr.org/papers/volume16/basu15a/basu15a.pdf\n   - https://doi.org/10.1093/bioinformatics/btq377\n   - https://www.pnas.org/content/115/9/2252.short\n\nOther comments:\n- page 8, Table 3: How do you explain the huge increase in performance for GrID-Net (eQTL & Sparse Peaks + Chromatin interactions & Sparse Genes + Chromatin Interactions & Sparse Peaks) and VAR Granger (Chromatin interactions & Sparse Peaks) over the correlation methods for the SHARE-seq data set? Based on Tables 1 and 2 it seems to be related to a sudden increase in AUPRC? Could you comment on that?\n- page 11, title of reference by Heinz et al. is incorrect \n- page 13, DAG construction: I would be careful when normalizing the data when used for causal inference because the normalization might destroy those relationships. Here I understand that the DAG is determined post-normalization, so that should be okay, since these relationships are inferred.", "summary_of_the_paper": "The authors propose a method for modeling single-cell multimodal data using Granger causality. They introduce a graph neural network called GrID-Net (Granger Inference on DAGs), where they assume that the network variables (the cells) follow a DAG structure. They focus on applying their method to the task of inferring non-coding genomic loci that causally influence the expression of a specific gene.", "summary_of_the_review": "Overall the method is sound and the authors have put in a lot of work to showcase the merits of their approach. I wouldn't say that the proposed approach \"outperforms existing methods for inferring regulatory locus-gene links\", but it seems to be an incremental step forward. Keeping in mind that the problem tackled is an exceedingly difficult one, I am inclined to give this paper the benefit of a doubt.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "mxDO4WU2AMM", "reviewer": "Reviewer_QYmT", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Detailed Comments:\n\n1. There are some issues with the model or claims about it that we note:\n1.a. \nThere are many instances in the manuscript where the relation between genetic loci and expression is referred to as “casual”. It needs to be made clear that this is “Granger causal” or “temporally causal” as this does not represent true causality. An obvious example would be a true causal factor that affects ATAC peak at time X and later, after some unrelated process takes effect, expression of a gene at time X+t. Clearly the definitions used here would conclude a casual relation that does not exist.\n1.b.\nUsing a t-test to check for differences in means of the full and reduced model seems problematic. A reduced model will always have higher MSE than a full model on the sample data. The metric needs to be adjusted for the number of free parameters/model complexity. An F test would be more appropriate. Alternatively, the authors could evaluate performance on a held out test set. As it stands, this issue makes us doubt the reported results. \n1.c.\nThere is a somewhat misuse of terminology with regards to the L variable defined in the manuscript as this does not represent the same lag as that is defined by Granger causality. The manuscript interprets L correctly but overloads the use of the term “lag”. It would be helpful to show a scenario (i.e reduction) where GrID-Net is equivalent to Granger causality inference using a linear model. For example, using a linear graph (Fig. 1a left) as input with L=1.\n1.d.\nIt was somewhat unclear to us if model implementation and target function in Eq 10 means the parameters of the reduced and full are shared. It seems to be the case. If so, isn’t it possible that the reduced model suffers from the full model optimization problem? I.e. if the reduced model was optimized on its own without shared params it would converge to a better predictive solution? \n\n2. There are several issues with the analysis that we note:\nWe understand that a technical contribution is of prime importance for ICLR. However, the following suggestions for improving the data analysis could significantly improve the practical utility of the manuscript. \n2.a. Data\n2.a.i\nUsing liftover to map human eQTLs to mouse genome is problematic.\n2.a.ii\nThe direction of the eQTL is not clear since they are chosen purely based on a p-value cutoff. Analyzing the relation between variant creation, change in expression and change in ATAC peaks would help clarify what the model is learning (enhancer/repressor creation/disruption).  \n\n2.b Performance and Evaluation \n2.b.i\nThe choice of baselines is not optimal. The Var Granger approach will perform poorly since it is just a simplification of GrID-Net. The correlation methods could have been applied to the graph directly to make the comparison more fair given that the graph itself is not part of the proposed method and is created using existing software. \n2.b.ii \nThe model is almost making random predictions (AUROC hovering around 0.5). It also does not improve significantly compared to the baselines. We understand that the relation to eQTL based calls are problematic but with such poor correspondence it brings doubt to using this as a metric. \n2.b.iii\nThe sparsity analysis is ancillary and does not contribute much to the manuscript. First, it is strange that the performance is marginally better on the sparse data (correlation of 0s inflating the result?). Second, without explicitly modeling single cell sparsity, mentioning robustness to sparsity is not productive. \n\n2.c. Mechanistic Analysis\n2.c.i\nHow is a TF associated with a gene? You need to show that the KO of a TF is actually associated with a gene expression increase/decrease. Presence of a TF motif does not always mean the expression is altered. \n2.c.ii\nThere is no enrichment analysis of TF binding. Does the GrID-Net detect more peak-gene pairs overall compared to the baselines? If so, a raw number of detected pairs doesn’t tell us anything about enrichment. A simple fisher test would suffice.\n2.c.iii\nThe interpretation behind the proximal vs distal analysis is not supported biologically. This is an extremely bold claim without experimental evidence. The authors should rephrase this section to be more speculative.\n2.c.iv\nConsider using enrichment tests here instead of a percentage. It is not clear how different the shown percentage is from the null. For example, is the ratio of distal to proximal labels close to the observed percentage after changing L? If so, the percentage increase is too small to have any biological meaning. \n\n2.d. Lack of clarity\n2.d.i \nHow are the p-value cutoff chosen for classification in AUPRC? Even if a standard 0.05 is used, it should still be mentioned.\n2.d.ii\nHow is Pearson correlation being computed with missing values in the sparsity analysis? If using 0s, this would inflate the correlation. \n2.d.iii\nHow does the model handle multiple terminal lineages (i.e. multiple nodes without any outgoing edges). What is y in this case?\n\n2.e Relation to previous work:\n2.e.i\nThe authors frame the entire work as a novel generalization of Granger causality. While we understand the need to frame the work in terms of ML novelty in a ML conference, we found this to misrepresent the actual work. In practice, Granger causality is not defined exactly as it is here, and the “generalization of Granger causality for DAGs” boils down to standard GNN operations over graphs with a memory of length L, using two competing architectures. The solution seems fine for the problem at hand (given a good answer to the possible issues raised above), but the remote relation to Granger seems to justify not much more than a general comment in the discussion.\n\n\nUPDATE\n==============\nWe read all of the authors detailed answers and reviewed the updated manuscript. The authors did a very good job addressing many of the concerns and comments raised by us and the other reviewers. While some concerns may remain we acknowledge the good work and effort made by the authors to improve the analysis, clarify details and make more precise claims. We have updated our score accordingly.\n", "summary_of_the_paper": "In this work, the authors propose a model framed as an extension to Granger causality for DAGs, and apply their model to infer “causal” relations between changes in gene expression and chromatin accessibility over time. Rather than relying on linear models and ordered data, the authors use a graph neural network to account for partially ordered data (i.e. branching structure). Similar to Granger inference, they compare a full model using paired data to a reduced model which excludes the lagging predictor. The authors then focus on analyzing temporally ordered scRNASeq and scATACSeq in 3 data sets from human and mouse. They test the method on predicting “causality” between the data modalities by using eQTLs and and chromatin interactions as proxy labels and compare against several baseline methods, noting that the method is robust to sparsity in data. Finally, they find enrichment of TF binding motifs in peak-gene pair regions and conclude that the model detects a lag between distal ATACSeq peaks and gene expression.  \n", "summary_of_the_review": "Overall, the paper introduces a model that, while a simple GNN, which appears to be a very reasonable approach for solving the problem of detecting the effects of chromatin accessibility on expression. The problem itself is interesting and important. In addition, the manuscript is generally clear and well written though the authors should be more careful about how claims are interpreted and described. However, the work does not offer many novel technical innovations since the GNN architecture is standard and the domain specific modeling/insights is lacking. Also, aspects of the evaluation of the model are concerning since the authors do not use out of sample data and should consider other baseline.  \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "irtTonFf-F1", "reviewer": "Reviewer_A4M6", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\nOne thing that i would like the authors to mention in their paper is that the presented work assumes causal sufficiency; i.e. that there are no hidden confounders in the data. This is a very strict assumption which is also often violated in real datasets. It is important to point out here that any Granger related method suffers from this issue, and for that reason the interpretation of the findings as 'causal' should be careful, as such methods are not able to differentiate between X \\rightarrow Y, and X \\leftarrow H \\rightarrow Y, where H is a potential hidden node acting as a confounder. \n\nI would like to ask some questions the authors:\n\n1. What is the intuition of the authors about the robustness of their method in sparse settings? I see the data driven conclusion,\nbut i would be interested in the theoretical explanation or some intuition about it from the authors.\n\n2. I appreciate a lot the presentation of the related work and the comparison experiments. Nevertheless, I would like to also see comparison with a non-linear implementation of Granger, as it is not straightforward to understand here whether the poorer performance of VAR Granger comes from the linearity of the method, or from the fact that Granger itselelf does not take into account partial ordering of the data.\n\n3. What is the variance of the results presented in Table 1? Is the difference between GrID-Net and the other methods significant?\n\n4. In the paragraph about methods related to Granger causality (p.5), the authors make an effort to list the related work, however they miss some important causal discovery methods in time series (and lagged panel data in general) methods that outperform Granger and don't require causal sufficiency.\ni.e.\n - Mastakouri et al.2020 Necessary and Sufficient Conditions for Causal feature selection in time series with latent common causes.\n- D. Entner, and P. O. Hoyer (2010). On causal discovery from time series data using FCI (PGM2010)\n- N. Pfister, P. Bühlmann, J. Peters: Invariant Causal Prediction for Sequential Data. Journal of the American Statistical Association, 114(527):1264-1276.", "summary_of_the_paper": "The paper proposes a modification of Granger causality so that it is applicable in partially ordered observations. Authors achieve that by adapting the architecture of GNNs to enable lagged-passed information in the DAG. Finally they apply it on a real data problem with the goal to infer non-coding genomic loci that causally influence the expression of a specific gene. In comparison with other three methods, the proposed method GrID net outperforms them on an open dataset.\n", "summary_of_the_review": "Overall I find the paper to be very well written, with a new method clearly motivated from a real-world problem the field of genomics. For that reason i am inclined for acceptance of the paper. I would like the authors to mention the strict assumption of causal sufficiency that their work assumes and make it clear that their results may suffer from it. I would also like to ask the authors to answer my questions above.\n\n--------------------- update ---------------------------\n\nSince the authors did not respond to any of my comments and concerns above, I decided to reduce my score. Not clearly stating that causal sufficiency is assumed is a major point, as it leads to wrong \"causal conclusions\". Moreover, the lack of explanation from the authors regarding their intuition in some key points of the paper contributes to lack of clarity. In addition to the above, the lack of comparison with linear Granger or at least discussion about it, make me decrease my score significantly. It is a pity because this could have been a good paper, if the authors have taken the time to answer and fix these points by answering to the reviewers. Without their responses/edits the paper contains strong overstatements.\n\n\n------------------- update after authors posted their rebuttal -------------------\nI am reverting my score to its original value, given that the authors answered all my comments and updated their paper accordingly.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "aHgo_EQeMr1", "reviewer": "Reviewer_xWaf", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths\n•\tIntegrative and generalized modeling of DNN and granger causal inference over DAG\n\n•\tThe authors provided details and reasonable justifications on three datasets, data processing, training & hyperparameters and pseudotiming analysis for building DAGs\n\n•\tThe benchmark is comprehensive and convincing, providing different metrics to evaluate GrID vs. other state-of-arts\n\n•\tValidation of gene-region interactions by independent eQTLs and Hi-C data\n\nWeaknesses\n•\tThe authors probably need to elaborate on the GNN architecture to represent historical information, ideally via visualization on Fig 1b. Also, do Hy and Hx use the same architecture (e.g., # of layers)?\n\n•\tIs it possible to compare DAGs by alignment-based methods such as manifold alignment methods? For instance, UnionCom, MMD-MA, etc use manifolds to align scRNA-seq & scATAC-seq and uncover pseudotiming trajectories. I suggest that the authors also benchmark pseudotiming approaches.\n\n•\tHow cell-type specific are those causal gene-region interactions? Can they better classify cell types than gene expression/chromatin accessibility only? Also, the authors still used the bulk tissue eQTLs data in GTEx that can be confounded by different cell types. \n\n•\tThe authors should discuss scalability, esp. on single-cell data applications. \n", "summary_of_the_paper": "This paper developed an exciting method, GrID-net to infer granger causality between multi-modality on directed acyclic graph (DAG) in the framework of deep neural network (GNN), aiming to improve identifying causal interactions across modalities. Primarily, the authors first find the DAG among samples and represent the historical information of modal features via aggregating GNN layers over DAG. They then apply granger causal inference to test the causal interactions of any two features if adding such GNN-based historical information of one feature significantly helps reduce the prediction loss of another feature. \n\nAlso, it is nice that the authors applied GrID-net to multiple recent single-cell multi-modal datasets (scRNA-seq and scATAC-seq) and identified casual chromatin regions for gene expression at the cell type level (grander causal gene-region interactions), providing novel gene regulatory mechanistic insights. Besides, the authors validated gene-region interactions by eQTLs, Hi-C, and TFBS enrichment.  \n", "summary_of_the_review": "In general, the manuscript was well written and organized logically. The GrID-net modeling was also presented clearly. I believe the new model has broad and timely impact on single cell data integration and understanding temproal gene regulatory mechanisms. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The AC and reviewers all agree that the paper proposes a very interesting framework to extend Granger Causality to DAG structured dynamical systems with important applications. 

The submission was the object of extensive discussion, and the AC and reviewers all agree that the author feedback satisfactorily addresses the vast majority of their concerns. We strongly urge the authors to incorporate all the points and revisions mentioned in their feedback. 

We certainly hope that the author will pursue this line of work and consider scaling their approach to tackle larger applications such as those related to social networks.
Decision: Accept (Poster)