Year: 2021

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. The loci identified by GrID-Net are also enriched for 10–50x more transcription-factor binding motifs than other methods. 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.

heading: 1 INTRODUCTION
text: Understanding the structure of a multivariate dynamical system often boils down to 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 variables can be very valuable. The framework of Granger causality does exactly that: in a dataset where observations are temporally ordered, a time-dependent variable x (with value xt at time t) is said to “Granger cause” a second variable y if the history of x at time t (i.e., x1, . . . , xt−1) is significantly predictive of y at time t even after accounting for y’s own history (Granger, 1969; Shojaie & Fox, 2021). Originating in econometrics, 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 a total ordering along time.
Often, only a 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.
The key conceptual advance of this work is extending the applicability of Granger causality to systems with partially ordered observations (Figure 1a). On a DAG where each node takes on multivariate values, we define a variable x as Granger causing another variable y if the latter’s value yv at the DAG node v can be significantly predicted by the values of x at the ancestor nodes of v, even after accounting for yv’s own ancestral history. We introduce GrID-Net (Granger Inference on DAGs), a
graph neural network (GNN) framework for Granger causal inference in dynamical systems that can be represented as DAGs (Figure 1b). We modify the standard GNN architecture to enable lagged message-passing on the DAG, enabling us to accumulate past information for a variable according to the DAG’s structure. This, combined with a 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 a 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 a pseudotime analysis (Haghverdi et al., 2016) (Figure 1c). Our work addresses a pressing need to identify temporally causal relationships between noncoding loci and gene expression, as traditional approaches for doing so are expensive and slow. They include population genetics techniques like expression quantitative trait locus (eQTL) studies, in which human genetic variability at a particular genomic locus is associated with expression changes in the gene of interest. Alternatively, perturbation-based approaches directly seek to identify gene expression changes in response to CRISPR alterations of specific noncoding genomic loci (Fulco et al., 2016). The emergence of single-cell multimodal assays that profile ATAC-seq and RNA-seq simultaneously offers a new way of discovering such regulatory relationships, on which we capitalize. Rather than requiring the large population sizes of eQTL studies or being limited to specific noncoding loci as in the perturbation-based 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 observations (“peaks”) in a cell are reported as empty (Singh et al., 2021; Zhang et 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 relationship between chromatin accessibility and gene expression. For a noncoding locus to influence a gene’s expression, the locus typically must first become accessible, upon which regulator proteins act on it to drive a change in gene expression (Figure 1d). 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 this 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 GrID-Net substantially outperforms current methods in identifying known peak–gene links, opening a path towards unlocking critical insights about gene regulation at genome scale. To our knowledge, this work presents both the first framework for Granger causal inference on dynamical systems represented as DAGs and the first single-cell method to explicitly model the temporal asynchrony between chromatin accessibility and gene expression.

heading: 2 METHODS

heading: 2.1 BACKGROUND ON GRANGER CAUSALITY: STANDARD FORMULATIONS
text: In time series analyses, the classical Granger causality framework uses the vector autoregressive (VAR) model, which can be expressed as follows (Lütkepohl, 2005).
yt = L∑ `=1 a (`) 1 yt−` + L∑ `=1 a (`) 2 xt−` + t (1)
Here, yt, xt ∈ R are values at time t for two stationary time series y and x, respectively; a(`)1 and a (`) 2 are coefficients that specify how lag ` affects the future of time series y; and t is a zero-mean noise term. In this model, yt is assumed to be a linear combination of the L most recent values of y and x each, with x said to Granger cause time series y if and only if a(`)2 6= 0 for all `. 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 x to forecast y, while a reduced model excludes the effect of x, only containing terms related to the past values of y.
yt = f (full)(yt−1, ..., yt−L;xt−1, ..., xt−L) + t (2)
yt = f (reduced)(yt−1, ..., yt−L) + t (3)
Here, f (full)(·) and f (reduced)(·) are generalized functions that specify how the value of yt 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 are significantly more accurate.

heading: 2.2 GRANGER CAUSALITY ON A DAG: GRAPH NEURAL NETWORK FORMULATION
text: Let the data be represented by a DAG G = (V, E) with n = |V| nodes (i.e., observations) and directed edges E indicating the information flow or partial order between these observations. For instance, when applied to the case of standard Granger causal inference on time series data, G would be the linear graph corresponding to the time series. Let y,x ∈ Rn correspond to the values of variables y and x on the nodes in V , with x putatively Granger-causing y. To infer Granger causality, we compare the full and reduced models as above:
y = g(h̃(full)y + ch̃ (full) x ) + (4)
y = g(h̃(reduced)y ) + (5)
Here, h̃(full)y , h̃ (reduced) y , h̃ (full) x ∈ Rn represent the historical information of y or x (as denoted by the subscript) that is aggregated by graph neural network (GNN) layers. For example, at each node v ∈ V , h̃[v] is the information accumulated from v’s ancestors in G. While h̃(full)y and h̃(reduced)y are outputs determined from similar architectures, their learned weights and bias terms end up being different since h̃(full)y is influenced by its interaction with h̃x. The coefficient c mediates this interaction and describes the effect of x’s history on y. In addition, ∈ Rn is a zero-mean noise term. We also include g(·) as an 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 are non-negative.
For brevity, we only describe how h̃(full)x is computed, with h̃ (full) y and h̃ (reduced) y being independently computed analogously. Also, sincex occurs only in the full model, for notational convenience we drop its superscript (full) below, writing just h̃x. We express h̃x as the mean of the outputs of L consecutive GNN layers, denoting the layerwise outputs as h(`)x :
h̃x = 1
L L∑ `=1 h(`)x (6)
h(`)x =
{ σ(w
(`) x AT+h (`−1) x + b (`) x ) if ` > 1
σ(w (`) x ATx+ b (`) x ) if ` = 1
(7)
Here, w(`)x , b (`) x ∈ R are the per-layer weight and bias terms, respectively. σ(·) represents the nonlinear activation in each of the GNN layers, chosen here to be the hyperbolic tangent function. A and A+ ∈ Rn×n are matrices defined by the DAG structure, designed to introduce a one-step lag in the first GNN layer that is propagated by subsequent layers. A is the column-normalized adjacency matrix of G, with Aij = 1dj if edge (i, j) ∈ E and 0 otherwise, where dj is the in-degree of node j. We note that a DAG does not have self-loops, so the diagonal terms of A are zero. In A+, we include these diagonal terms, setting them to the same value as others in the column (i.e. dj is incremented by 1 inA+ to preserve normalization).
Together,A andA+ allow us to extend the key intuition of Granger causality to DAGs: predicting a value at 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 v but not the information at v itself. Each subsequent layer introduces the preceding set of parents, with the diagonal term in A+ ensuring that information already stored at v 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 L 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. 1) but with the crucial advantage of leveraging the richness of a DAG’s structure in directing the information flow and aggregating sparse data.

heading: 2.3 COMPARING FULL AND REDUCED MODELS TO INFER GRANGER CAUSALITY
text: Let Θ(full)y ,Θ (reduced) y ,Θ (full) x ∈ R2L denote the set of parameters for the full and reduced models, with Θ(full)x = {w(1)x , ..., w(L)x ; b(1)x , ..., b(L)x } and Θ(full)y ,Θ(reduced)y defined analogously. All the Θ parameters are jointly learned by minimizing the combined loss of the full and reduced models (Montalto et al., 2015):
L(total) = ∑
x,y∈P (∑ v∈V L(ŷ(full)v , yv) + ∑ v∈V L(ŷ(reduced)v , yv) ) (8)
where ŷ(full)v and ŷ (reduced) v 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 (see discussion in Appendix A.8). Also, in a multivariate system, not all pairwise combinations of variables may be relevant: P is the subset of x, y variable pairs whose putative Granger-causal interactions are of interest.
To infer Granger causal interactions between y and x, we compare the set of loss terms associated with the full model L(ŷ(full)v , yv) to those of the reduced model L(ŷ(reduced)v , yv), 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 a one-tailed F-test, comparing the residual sum of squares of the reduced and full models; an alternative assessment using Welch’s t-test yielded very similar results (Appendix A.7). 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.

heading: 2.4 DOMAIN-SPECIFIC MODEL CUSTOMIZATION
text: The loss L 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 number of GNN layers L corresponds to the maximum amount of past information desired. Similarly, to introduce a k-hop lag, the first k GNN layers would use the matrix A while the later layers would useA+. Note that in the above, we described a one-hop lag, which we used for all analyses.

heading: 2.5 TRAINING DETAILS AND HYPERPARAMETERS
text: GrID-Net models were trained using the Adam optimizer with a 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 initialization (Bengio & Glorot, 2010). All GrID-Net models consisted of L = 10 GNN layers; the architectures of the three sub-models (h̃(reduced)y , h̃ (full) y , and h̃ (full) x ) were identical but separate. All models were implemented in PyTorch and trained on a single NVIDIA Tesla V100 GPU.

heading: 3 RELATED WORK
text: 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; Argelaguet et al., 2020; Gayoso et al., 2021). Meanwhile, efforts to detect associations at the level of specific features across modalities have thus far been limited. Current approaches for inferring peak–gene associations predominately rely on simply calculating correlations between these two modalities (Ma et al., 2020; Zhu et al., 2019; Liu et al., 2019). Not only do these correlationsbased approaches fail to account for the temporal asynchrony between chromatin accessibility and gene expression, but they are also sensitive to the inherent noisiness of single-cell data (Zhang et al., 2019). These limitations point to the need for methods that go beyond analyzing mere correlations to leverage 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-based interactions in multivariate time series data. One category of methods involves nonlinear kernel-based regression models (Marinazzo et 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-based Granger causality models have also more recently been proposed to account for more generalized nonlinear interactions (Tank et al., 2018; Marcinkevičs & Vogt, 2021). Of these methods, efforts have been made to leverage specific neural network architectures to 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).
As a 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 cell with an associated multivariate observation, and edges describe the flow of information between cells. In contrast, a causal graph describes a data generation process with nodes corresponding to the dynamical system’s variables. In our single-cell peak–gene inference application, performing Granger causal inference on a DAG allows
us to capture 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 a single-cell dataset into a sequential ordering of cells along pseudotime, a 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 inability to rule out hidden confounders that might mediate the interaction between x and y (Mastakouri et al., 2021). Nonetheless, such analysis remains valuable in our peak–gene inference task, as it enables the identification of peak–gene relationships that are temporally causal, upon which a shortlist of hypotheses can be generated for subsequent perturbational or population genetics studies. 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)) including Bayesian inference (Brodersen et al., 2015; Glymour et al., 2019) or structural equation models (Spirtes & Zhang, 2016) could be applied to address this limitation.

heading: 4 RESULTS
text: Datasets and Preprocessing We analyzed three single-cell multimodal datasets with per-cell ATAC-seq and RNA-seq profiles that characterize a range of dynamical systems, including cell differentiation and drug-treatment responses (Cao et al., 2018; Chen et al., 2019; Ma et al., 2020). The sci-CAR dataset profiled 3,260 human lung adenocarcinoma-derived A549 cells after dexamethasone 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 a 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 perturbation-response 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/10) 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 a 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-distance cutoff because almost all enhancer-gene pairs reported in a recent benchmark study were within 1 Mb (Moore et al., 2019). 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–6 hours and required a 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.
DAG construction For each of the three multimodal datasets, a kNN graph was constructed on representations of cells that unified information from both gene expression and chromatin accessibility via Schema (Appendix A.1, Singh et al. (2021)); other approaches for integrating data modalities (e.g., manifold alignment (Cao et al., 2020)) could also be used here. A pseudotime value was then inferred for each cell, after which edges in the kNN graph are retained if they align with the direction of increasing pseudotime. Since each edge is oriented in the direction of increasing pseudotime, no cycles exist and the resulting graph is a DAG. In our results below, we used k=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 (Appendix A.2) to the choice of k 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 a 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-gene pair that associates naturally occurring genetic variability at a locus to gene expression variations observed 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-C/5C) 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 chromatin accessibility are both manifestations of the underlying gene regulatory process (Figure 1d): a noncoding region becomes accessible (ATAC-seq) for binding by a gene regulator and the chromatin 3D structure is arranged (chromatin interaction) so that the regulator is proximal to the target gene’s promoter, thus enabling transcriptional control (Delaneau et al., 2019). As such, compared to random peak–gene pairs, causal pairs should more likely be in close spatial 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 underlying biology only indirectly and are themselves prone to errors: eQTLs are confounded by correlated mutations elsewhere in the genome and Hi-C/5C data suffers from similar sparsity issues as ATACseq. 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 chromatin 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-based techniques used in previous single-cell multimodal studies (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 use the absolute value of the correlation scores to 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 autoregression (VAR Granger) and a state-of-the-art 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). 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 adenocarcinoma-derived 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 (Gjoneska 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–gene pair in the eQTL 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 predicting the true eQTLs based on the ranking of association scores generated by each method for each peak–gene pair. Across all datasets, GrID-Net outperformed the other methods in predicting eQTLs, achieving the highest AUPRC (area under precision-recall curve) and AUROC (area under receiver operating characteristic) in all cases, indicating that peak–gene pairs prioritized by GrID-Net more closely align with evidence from population genetics (Table 1). The non-linear Granger approach (GVAR) performs slightly better than linear approach (VAR Granger); however, both methods performed 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.
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 a peak–gene pair by identifying the pairs of genomic windows in the corresponding chromatin interaction 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 GrID-Net and alternative methods in predicting these pairs. We found that GrIDNet 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 2).
GrID-Net is more robust to sparse data than existing approaches With the sparsity of singlecell multimodal data expected to remain a concern, we assessed GrID-Net’s robustness to such sparsity. In particular, 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 sparsity. Because of differing true-positive rates between the full set of candidate peak–gene pairs and the sparse-data subset (Appendix A.9, Table 15), 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 outperformance compared to the alternatives was accentuated (Table 3, additional details in Appendix A.9). GrID-Net also noticeably outperforms the other methods in predicting significant chromatin interactions associated with peak–gene pairs consisting of sparse genes or peaks. These results provide further evidence of GrID-Net’s utility in extracting insights from sparse single-cell data.
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 aimed to relate peaks in high scoring peak–gene pairs to the transcription factors (TFs) that might bind there. Accordingly, we selected the highest scoring 1% of peak–gene pairs from each of GrID-Net and its alternatives. We then used Homer (Heinz et al., 2010) to identify TF binding motifs that are co-enriched in sets of peaks linked to a particular gene. The presence of these co-enriched motifs is potentially indica-
tive 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, we observed that the peak–gene pairs detected by GrID-Net were associated with 10–50 times more enriched TF binding motifs than the other methods (Table 4), 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 a valuable role in unveiling mechanistic insights into TF-gene regulatory control.
We further evaluated the functional importance of detected peak–gene pairs by comparing peaks in these peak–gene pairs to TF binding sites supported by TF ChIP-seq experimental data (Oki et al., 2018). We first obtained TF ChIP-seq data for the specific cell type(s) represented in each singlecell multimodal dataset. For each of these datasets, we 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 ChIP-seq-derived 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 4). As TF binding sites are indicative of functional regulatory roles for noncoding loci (Valouev et al., 2008), GrID-Net’s effectiveness in prioritizing peaks bound by TFs provides additional evidence for its ability to reveal important gene regulatory roles for specific loci in the noncoding genome.
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 L (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.11; Tables 16, 17).

heading: 5 DISCUSSION
text: We extend the applicability of Granger causal analysis to dynamical systems represented as DAGs, introducing GrID-Net to perform such inference. GrID-Net takes advantage of the expressive power of graph neural networks, enabling the detection of long-range nonlinear interactions between variables. We focused on applying it to study multimodal single-cell gene regulatory dynamics, given the prevalence of graph-based representations of cellular landscapes in the field. GrID-Net demonstrated substantial improvements in accuracy relative to 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, positions 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 applicability of GrID-Net. The ability of GrID-Net to capture long-range, nonlinear interactions in DAGs opens the door to novel analyses on the many existing datasets that are characterized by these graph-based dynamics as well.

heading: ACKNOWLEDGEMENTS
text: We thank the anonymous reviewers of ICLR 2022 for their valuable feedback and suggestions. Figure 1 was created using biorender.io.

heading: REPRODUCIBILITY STATEMENT
text: In order to ensure that the full set of results presented in this paper are reproducible, we have included source code for our method, which can be found in the supplementary materials. In addition, all analyses described in the paper were performed on publicly available 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 detailed pre-processing 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/5C, and TF binding motif data in the “Appendix” section. We also made sure to provide details on 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”.

heading: A APPENDIX

heading: A.1 DAG CONSTRUCTION DETAILS
text: For the three single-cell multimodal datasets, we applied a standard log(1+CPM/10) transformation to the raw RNA-seq counts data and a log(1+CPM/100) 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 et al., 2015) were retained. In the original SHARE-seq and SNARE-seq studies, topic modelingbased 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 tf-idf, a procedure derived from text-based topic modeling that has been applied extensively to analyze single-cell ATAC-seq data (Cusanovich et al., b;c;a).
The normalized RNA-seq and ATAC-seq data were then projected onto their top 100 principal components (excluding components that have Spearman correlation ρ > 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 a secondary modality. For the sci-CAR dataset, the time of data collection post-drug treatment was used as a tertiary modality. Meanwhile, ATAC-seq served as the primary modality for the SHARE-seq dataset, as the original SHARE-seq study identified strong cell cycle effects represented 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 (k = 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 a DAG.

heading: A.2 DAG ROBUSTNESS ANALYSIS
text: In addition to the diffusion pseudotime algorithm, we also evaluated Palantir and Monocle for inferring pseudotime. The performance of GrID-Net was robust to the choice of algorithm (Tables 5, 6,7,8).

heading: A.3 SINGLE-CELL DATA PREPROCESSING AND CANDIDATE PEAK–GENE PAIR SELECTION
text: 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 are 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.1% 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.

heading: A.4 BENCHMARKING
text: For the VAR Granger model, peak–gene pairs were scored on the − log10(p)-value associated with the statistical significance of Granger causality, as calculated using the grangercausalitytests function with default settings in the statsmodels package (Seabold & Perktold, 2010).

heading: A.5 MATCHING EQTLS WITH SINGLE-CELL ATAC-SEQ & RNA-SEQ DATA
text: The procedure to match a 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 statistically 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 et al., 2021) with default
settings to convert human genomic coordinates represented in the GTEx datasets to mouse coordinates. 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 it is organized around different strains of mice rather than tissues (Doss et al., 2005; Hofstetter et 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 (Gjoneska et al., 2015; Xiao et al., 2012).

heading: A.6 CHROMATIN INTERACTION DATA PREPROCESSING
text: 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% of interactions as the set to be predicted.

heading: A.7 COMPARISON OF F-TEST AND WELCH’S T-TEST TO RANK PUTATIVE PEAK–GENE INTERACTIONS
text: For each candidate peak–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-squared-errors of the reduced and full models as per an F distribution with (2L + 1, n − 4L − 1) degrees of freedom. Here, n is the number of cells, and for each peak-gene pair, each of the L GNN layers includes two parameters per peak/gene, and c is an additional parameter in the full model). Another approach would be to compare the full and reduced models’ losses with the Welch’s t-test, where the alternative hypothesis is that the mean of the loss terms for the full model is 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.
We compared the performance of the two tests on the eQTL and chromatin interaction prediction tasks, finding that they yielded very similar results (Tables 9,10). This similarity is likely due to our use of the test statistics to rank peak–gene pairs, upon which the top-ranking pairs are selected. Moreover, the degrees of freedom are identical across all candidate pairs, so their ranking depends only on the relative losses of the full and reduced models.

heading: A.8 COMPARISON OF COMBINED AND SEPARATED LOSS FUNCTIONS FOR FULL AND REDUCED MODELS
text: 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 as intended.

heading: A.9 ROBUSTNESS OF GRID-NET TO SPARSITY
text: Tables 13 and 14 expand upon the results in Table 3 by also showing standard deviations (across three runs) of GVAR’s and GrID-Net’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.

heading: A.10 TF BINDING MOTIF ENRICHMENT ANALYSIS
text: For each gene that was associated with more than 5 peaks in the highest scoring 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 BenjaminiHochberg-corrected 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.

heading: A.11 GRID-NET SUGGESTS POSSIBLE LINK BETWEEN GENOMIC DISTANCE AND TEMPORAL REGULATORY CONTROL
text: 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 (≥ 100 kb). We then evaluated different architectures of GrID-Net that vary in the number of GNN layers L, hypothesizing 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% of peak–gene pairs detected by different architectures, we found this to indeed be the case: across all datasets, architectures with higher L detected larger ratios of distal-to-proximal peak–gene pairs (Tables 16, 17), suggesting that distal peak–gene interactions may be marked by longer pseudotime difference between peak and genes. This finding suggests the possibility of a previously unreported link between genomic distance and the temporal lag between a regulatory element’s accessibility and its target gene’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 between distal peaks and genes may be mediated by other regulatory mechanisms rather than being a direct influence. One hypothesis is that the complex coordination of regulatory factors required to enable distal regulatory interactions is associated with greater lags between distal regulatory elements’ accessibility and the gene expression changes they affect (Panigrahi & O’Malley, 2021).

heading: A.12 RUNTIME AND MEMORY USAGE
text: GrID-Net’s runtime and memory usage are dependent on a variety of factors, most notably the number of nodes in the DAG describing the system of interest. In our model, the DAG is represented by the N × N adjacency matrices A and A+, where n 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 singlecell datasets, we recommend the use of sketching techniques (Hie et al., 2019; DeMeo & Berger, 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, 35 GB, and 39 GB for these respective datasets.

Review summaries:
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.   Also, 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.   

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.

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.   

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.