DeepReGraph co-clusters temporal gene expression and cis-regulatory elements through heterogeneous graph representation learning

Data pre-processing

Data pre-processing was done in two stages. Firstly, gene expression and cCRE datasets were normalized across temporal data. The second stage included filtering-out temporal profiles where the correlation across replicates or the variance across time was below some pre-defined thresholds. We explain the pre-processing stages in greater detail below.

We firstly downloaded the gene expression values, called peaks, and bam files of ATAC-seq, H3K27ac, and H3K27me3 from the ENCODE database¹ (mice heart fetal tissue). We used logFPKM values of gene expression throughout this paper. To detect cCRE regions, we re-centered ATAC-seq called peaks across time points and replicates using the bioconductor DiffBind package.²⁹ To determine chromatin accessibility in each cCRE region, we took regions of 500 nucleotides long (250 nucleotides on each side from the center of cCREs). We recorded binding intensities for each replicate of each time point separately by using the bioconductor Rsubread package.³⁰ We counted H3K27ac and H3K27me3 marks for cCREs in the same manner, but used regions 3000 nucleotides long instead (1500 nucleotides each side from the center of cCREs). We took the median value of non-peak regions in each experiment separately to remove background counts from real intensity counts. Then, we subtracted these median values from the intensities of each experiment separately. We used logRPKM values of denoised counts across replicates and time points. Removing the background noise caused the distribution of each experiment to resemble a Gaussian distribution.

At the end of this process, we had two replicas of gene expression time-series profiles for a set of genes and two replicas of cCRE activity time-series for each activity marker over a set of cCREs. We then proceeded to filter out invalid genes and cCREs based on two criteria. The first group of filters we applied to this data were proposed in Ref. 31. For each gene, we measured the correlation of the gene expression profiles of both replicates and discarded every gene whose replicate had a negative correlation value. We did the same with the cCREs for each activity marker profile. If only one activity marker had a negative correlation between replicates, the cCRE was discarded. We also discarded every element (gene or cCRE) where two consecutive time probes had an absolute ratio greater than 2.

We also discarded low-expressed genes and low activity-characterized cCREs.³² We computed the mean gene expression value for each gene and discarded every gene whose mean expression value was under the 0.8 percentile of such mean value distribution. We did the same with cCREs: we computed the mean activity value for each activity marker and discarded every cCRE where the mean value was below a minimum percentile threshold. Such thresholds were 0.8 for ATAC-seq and 0.7 for H3K27ac and H3K27me3. Lastly, for gene expression and each activity marker dataset, we created a 3-rd degree polynomial regression to estimate the variance of each time series as a function of its mean value. We discarded all the elements where the actual variance was below the predicted variance. Our pre-processed and cleaned data consisted of 607 genes and 5239 cCREs from the mouse heart dataset. This dataset is available online on our public repository. Our objective was to shed light on regulatory mechanisms among genes and cCREs. Having cleaned our dataset, modeling it as a heterogeneous graph is the strategy we used to combine all the data and learn a new relation type between elements that gives information about gene expression regulation. We explain the graph modelling in the next paragraph.

Heterogeneous graph model

Many industrial data are represented in a graph, as multiple entities with quantifiable relationships between them. Nowadays, graph modelisation is widely used in bioinformatics.³³ One strategy to extract added-value from graphs is the use of graph RL algorithms. Graph RL has been widely studied by the research community over the past few years years.^25,34–36 Graph RL algorithms seek to map the information of a graph on a reduced-dimension latent space where the geometrical distances between elements in such space resemble the relations in the original graph. In other words, graph RL, also referred to as “graph embedding”, consists in finding a reduced dimensional representation of the nodes of a graph while conserving most of the semantic information of such a graph.^37,38

A graph can be homogeneous if it has a unique type of node and a unique type of relation defined between them, or heterogeneous, if it introduces multiple edge and node types. We argue that regulatory networks during embryonic development can also benefit from being modelled as a heterogeneous graph. Reduced-dimension representations of such a graph could contain multiple types of information that, if well combined, could semantically express regulatory mechanisms of gene expression.

We defined a heterogeneous graph and described it as $\mathcal{G}=\left(\mathcal{V},\mathrm{ℰ},\mathcal{T}\right)$, where $\mathcal{V}$ is the set of nodes, $\mathrm{ℰ}$ is the set of edges and $\mathcal{T}={\mathcal{T}}_{\mathrm{ℰ}}\cup {\mathcal{T}}_{\mathcal{V}}$ is the set containing all the node types and edge types. The set of node types consists of genes and cCREs, and is denoted by ${\mathcal{T}}_{\mathcal{V}}=\left\{\mathit{\text{genes}},\mathit{\text{ccres}}\right\}$. Each gene $g_i,\forall i\in \left|{\mathcal{V}}_{\mathit{\text{genes}}}\right|$ in the dataset is characterized by a gene expression time-series profile that will be encoded in a vector denoted by ${\mathbf{g}}_i$. Note that the length of ${\mathbf{g}}_i$ coincides with the number of time-points considered in the gene expression data. The vectors ${\mathbf{g}}_i,\forall i\in \left|{\mathbf{V}}_{\mathit{\text{genes}}}\right|$ can be concatenated as rows to form the feature matrix represented in red in panel A of Figure 1.

Figure 1. Schematic representation of DeepReGraph.

Our framework combines multiple temporal enhancer activity markers with temporal gene expression data and base-pair distances to build a heterogeneous graph of genes and cis-regulatory elements (CRE). The proposed framework then applies an adapted version of AdaGAE²⁶ to find a low-dimensional, easy-to-cluster embedding representation of data through an iterative optimisation process. In the final embedding produced by our method, the spatial distribution of genes and cCREs resembles candidate gene-expression regulatory Mechanisms.

Every candidate cis-regulatory element, instead, will be denoted as $c_j,\forall j\in \left|{\mathcal{V}}_{\mathit{\text{genes}}}\right|$ and will be characterized by $M$ activity time-series profiles, that will be encoded in corresponding vectors denoted as ${\mathbf{c}}_j^1,{\mathbf{c}}_j^2,\dots ,{\mathbf{c}}_j^M$. The time-points under which gene expression and cCRE activity values are available are the same. In this work, $M=3$. Specifically, we considered the ATAC-seq, H3K27ac, and H3K27me3 activity time-series profiles for each cCRE, but one could include more cCRE activity markers in the set of features. Note that the vectors ${\mathbf{c}}_j^m,\forall j\in \left|{\mathbf{V}}_{\mathit{\text{cCREs}}}\right|$ can be concatenated as rows to form the matrices represented in blue in panel A of Figure 1.

Each node of our graph was associated to one or multiple feature vectors that can be seen as points on a multi-dimensional feature space. The set of feature-spaces defined for the nodes in $\mathcal{G}$ is denoted as $R=\left\{G,C_1,C_2,\dots C_M\right\}$ where $G$ is the gene expression feature space and $C_1,C_2,\dots C_M$ are the feature spaces of the $M$ cCRE activity data. Each one of these feature spaces will be referred to as original feature space, and the vectors ${\mathbf{g}}_{\mathbf{i}},{\mathbf{c}}_{\mathbf{j}}^{\mathbf{1}},{\mathbf{c}}_{\mathbf{j}}^{\mathbf{2}},\dots ,{\mathbf{c}}_{\mathbf{j}}^{\mathbf{M}}$ as original feature vectors from now on. Note that such vectors correspond to the rows of the matrices in panel A of Figure 1. It is well-known that the probability of the existence of a regulatory mechanism between a gene and a cis-regulatory element is inversely proportional to the base-pair distance between them.³⁹ For this reason, in our work, we used another data source, called the Link Matrix, that gives us information about the base-pair distance between genes and cCREs. We represented such a matrix in panel C of Figure 1. Interaction between cCREs and genes is generally possible in the range of 1 Mega base-pair.⁴⁰ Consequently, we considered only distances lower than ${10}^6$ base-pairs. For practical purposes, we scaled the base-pair distance values in the interval $\left[0,1\right]$ to define a scaled base-pair distance function, $D^{\mathit{BP}}$:

Given the scaled distance function, we created a base-pair proximity relationship $S_{\mathit{BP}}$ using a parametric transformation:

The interaction between CREs and the genes they regulate is complex. Some CRE markers are correlated with gene expression and some others are anti-correlated.⁴¹ This correlation or anti-correlation is governed by the mechanism that a specific marker affects the gene expression. For example, chromatin should be opened prior to the gene expression. Therefore, chromatin accessibility is commonly directly correlated with gene expression. H3K27ac is also directly correlated with gene expression, as this epigenetic modification happens in active transcription regions. On the other hand, H3K27me3 shows epigenetic modifications that happen at polycomb repressed regions, and are consequently inversely correlated with the expression of the corresponding regulated genes. We needed to capture the relevance of these and other temporal covariances between gene expression and the correspondent regulatory elements’ activity markers.

To this end, we created a joining score, $J_T$, as a function of the temporal slopes of gene expression and cCRE activity time-series:

We computed a trend-aware score multiplying (1) and (2) for each gene-cCRE pair in our graph:

We represented the computation of Equation (3) in panel D of Figure 1. Having defined the trend-aware score, we could differentiate between any pair of genes, $g_j$ and $g_k$, when these were equally proximal to any cCRE $c_i$ but had gene expression vectors with different temporal trends. In this case, the original base-pair proximity score in (1) assigned the same value to the association $\left(g_j,c_i\right)$ and $\left(g_k,c_i\right)$. Notice that the trend-aware score in (2), instead, distinguishes the most plausible association of $c_i$ with respect to the alternatives $g_j$ and $g_k$ as a function of the temporal slopes of the gene expressions elements. In our experiment, we found the best results setting ${\omega }_{\mathit{\text{ATAC}}}=1,{\omega }_{H3K27\mathit{ac}}=0$ and ${\omega }_{H3k27\mathit{me}3}=0$ in Equation (2).

The set of edge types in the graph is denoted as ${\mathcal{T}}_{\mathrm{ℰ}}=\left\{S_{\mathit{BP}|T},S_G,S_{C_1},S_{C_2},\dots ,S_{C_M}\right\}$. As explained before, $S_{\mathit{BP}|T}$ indicates the trend-aware base-pair proximity relationship explained above; $S_G$ stands for the gene expression profile similarity between genes, and $S_{C_1}$, $S_{C_2}$, …, $S_{C_M}$ are the relationships that express the different cCRE activity time-series similarities. We aimed to combine all edge types in ${\mathcal{T}}_{\mathrm{ℰ}}$ to find a unique multi-semantic edge type that resembled gene expression regulation mechanisms. The process of combining edge types is represented in panel E of Figure 1 and is explained later in this section. To perform the combination of all edge types in ${\mathcal{T}}_{\mathrm{ℰ}}$, we modelled the data as a heterogeneous graph and designed a graph RL algorithm capable of extracting such a dependency between elements using prior biological knowledge and the data available. We explain the graph RL algorithm we used to fulfil our objective in the next paragraph.

Graph representation learning

Graph RL algorithms aim at finding a matrix: $\mathbf{Z}\in {\mathrm{ℝ}}^{\left|\mathcal{V}\right|\times d}$ that contains information about every feature-space in $R$. In our case, we used a graph RL algorithm to find a unique feature space for the nodes of $\mathcal{G}$ that comprises most of the information encoded by the feature spaces contained in $R$. The rows of $\mathbf{Z}$ correspond to a new set of feature vectors for each node in $\mathcal{G}$. Each row of $\mathbf{Z}$ is called an “embedding vector” and is a reduced-dimension representation of the whole set of original feature vectors of a given node.^37,34 The embedding vectors are “multi-semantic” in that they contain information from every feature space in $R$. The dimension of embedding vectors, $d$, is generally lower with respect to the dimension of the original feature spaces. The embedding matrix $\mathbf{Z}$ is represented in panel F of Figure 1.

Graph RL is also referred to as “graph embedding”. Giving a node $v\in \mathcal{V}$, and a specific edge-type $\tau \in {\mathcal{T}}_{\mathrm{ℰ}}$, the neighborhood of $v$ in $\tau$ is denoted as ${\mathcal{N}}_{\tau }\left(v\right)$ and corresponds to the set of nodes that are connected to $v$ by $\tau$-type edges. We denote with $\mathcal{N}\left(v\right)$ instead the set containing every neighborhood of $v$:

The heterogeneous graph embedding process can be represented by an encoding function that takes in input a node $v\in \mathcal{V}$ and its neighborhood set, $\mathcal{N}\left(v\right)$, and returns a new representation for that node:

One of the most common methodologies for implementing (4) is the encoder-decoder paradigm.²⁵ This paradigm also models an auxiliary decoding function. Considering a homogeneous graph, i.e. a graph with a unique node feature space $X$, the encoder function in (4) takes as input original feature vectors of $X$, and outputs the corresponding embedding vectors:

The decoding function, instead, takes as input a pair of node embeddings produced by the encoder, e.g. ${\mathbf{z}}_i$ and ${\mathbf{z}}_j$, and outputs a real number. Such a number is the prediction of the value of a pre-defined pair-wise relationship between ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$:

The encoding and decoding functions should minimize a reconstruction loss of the form:

Graph auto-encoders

A convenient framework that implements the encoding-decoding paradigm is the deep graph auto-encoder (deep GAE).⁴² In this framework, (4) was implemented through a graph neural network (GNN), an artificial neural network that takes the original feature vectors and the adjacency matrix of $\mathcal{G}$ as input and produces the node embedding vectors. GNNs are non-linear parametric functions, and after a pre-defined initialization schema, the parameters of these functions can be optimized to minimize (6) through gradient descent.

In this work, we proposed the use of the basic GNN model presented in Ref. 25. In a homogeneous graph context, the basic GNN takes as input the node feature matrix $\mathbf{X}\in {\mathrm{ℝ}}^{\left|\mathcal{V}\right|\times m}$ and a pre-defined graph adjacency matrix denoted by $\mathbf{A}\in {\mathrm{ℝ}}^{\left|\mathcal{V}\right|\times \left|\mathcal{V}\right|}$, and performs non-linear parametric operations with this information to produce $\mathbf{Z}$. In matrix notation, the operations made by the basic GNN model can be represented as follows:

We defined a parametric “architecture” for the encoding function using (7), but still needed to define $\mathbf{A}$ and $\mathbf{X}$, to drive the optimization of the parameters of $f$ to minimize (6). The initialization of these matrices can follow multiple strategies. For example, in the neural message passing framework,⁴³ given a homogeneous graph context, one generally constructs $\mathbf{X}$ stacking the original feature vectors as the rows of such a matrix. However, it is not the only valid strategy. In our heterogeneous graph context, we decided to initialize $\mathbf{X}$ as the identity matrix because we had multiple original feature matrices. We were confident to use the identity matrix to initialize $\mathbf{X}$ because we carefully initialized the adjacency matrix $\mathbf{A}$ combining the information of every relationship type in our graph, as we will explain in the next paragraph.

Finally, the similarity measure that the decoder is meant to predict needs to be specified to implement (5). In this paper, we proposed the use of the decoding similarity function proposed by Ref. 26 to converge into a clustering-friendly embedding, i.e., an embedding with dense and easily identifiable clusters. In the next paragraph, we give a brief explanation of the algorithm in,²⁶ which performs graph RL on a single-feature space. For a more detailed exposition of the work of Xuelong et al., the reader is referred to the original paper.²⁶

Adaptive-sparsity graph auto-encoder

In this work, we proposed the use of AdaGAE,²⁶ which is a state-of-the-art method for RL and clustering through graph modelling of a dataset. The whole AdaGAE algorithm is wrapped in a gray rectangle in Figure 1. This algorithm is defined in a homogeneous data context. In other words, given a unique original feature matrix $\mathbf{X}\in {\mathrm{ℝ}}^{\left|\mathcal{V}\right|\times m}$, the objective was to produce a low-dimensional embedding matrix $\mathbf{Z}\in {\mathrm{ℝ}}^{\left|\mathcal{V}\right|\times d}$ where $d\ll m$. Authors in²⁶ used a deep GAE to embed high-dimensional datasets in a graph-like, low-dimensional format. Xuelong et al. generated a sparse graph as a function of the Euclidean distances between the original feature vectors using a k-nearest neighbors⁴⁴ (k-NN) model.

The main novelty of AdaGAE is the adaptiveness in the sparsity parameter $k$. The graph auto-encoder is, in fact, iteratively optimized with respect to various objective functions like (6), and each objective function is constructed using a different sparsity parameter. With a custom decoding function, such a dynamic model leads to a sparse graph. In other words, they produce a clustering-friendly embedding of the original data.

In other words, at each iteration, $l\in \left\{\mathrm{0,1,2}\dots ,L-1\right\}$, a sparse and weighted graph is generated. In this graph, the nodes correspond to the samples of the original dataset. Instead, the weights of the edges are inversely proportional to the Euclidean distances between the original feature vectors. Specifically, at iteration $l$, a graph is generated in which the weight of the link between nodes $v_i$ and $v_j$ is denoted as $p_{i,j}^l$ and is a function of a sparsity parameter $k^l$:

Notice that, in (8), the parameter $k^l$ might have a different value at each iteration. Such a parameter induces the sparsity of the generated graph: for each node $v_i$, only up to $k^l$ elements will have a non-zero weighted link with that node. Thus, at each iteration $l$, AdaGAE generates a $k^l$-sparse graph. Notice also that the function (8) is a discrete probability density function (PDF):

The embedding matrix $\mathbf{Z}$ produced by AdaGAE was iteratively optimized. We denoted with ${\mathbf{Z}}^l$ the embedding matrix after iteration $l$. The embedding vector of node $v_i$ after iteration $l$, instead, coincides with the $i$-th row of ${\mathbf{Z}}^l$ and is denoted as ${\mathbf{z}}_{\mathbf{i}}^{\mathbf{l}}$. At each iteration, AdaGAE proposes to implement the decoding function in (5) as a function of the embeddings of the previous iteration:

At each iteration $l$, AdaGAE proposes to use the Kullback-Leibler divergence (KL-divergence)⁴⁵ of ${\mathbf{P}}^{\mathbf{l}}$ with respect to ${\mathbf{Q}}^{\mathbf{l}}$ to implement (6). Explicitly, the reconstruction loss that AdaGAE seeks to minimize is the following:

Xuelong et al. relied on the graph convolutional network (GCN) model⁴⁶ to implement the encoder of their deep GAE. However, in this paper, we used a simpler GNN model, as explained in the previous paragraph. We represented the iterative computation of (8) and (9) and the iterative optimization of (10) in panel E of Figure 1.

Notice that the main dynamic criterion that causes the objective functions to change is the sparsity parameter $k^l$. However, the sparsity is not the only thing that changes from iteration to iteration: In the first iteration, the distribution family ${\mathbf{P}}^0$ is generated through (8) using the Euclidean distances between the original feature vectors. After the first iteration, instead, the distribution family ${\mathbf{P}}^{\mathbf{l}},\forall 1<l<L-1$ is generated using (8) as a function of the Euclidean distances between the embedding vectors of the ${\mathbf{Z}}^{l-1}$ embedding matrix. By doing so, AdaGAE considers high-order neighborhoods of each node to construct the final embedding matrix. In other words, AdaGAE is said to “exploit the high-level information” present in the data.

Different embeddings for various iterations are plotted in panel G of Figure 1. In those plots, one can see the embedding of a small set of genes and cCREs extracted from the original mouse heart dataset. Large points represent genes, while cCREs are the small points. Point colors, instead, correspond to single-modality or homogeneous cluster assignments, i.e., cluster assignments of genes and cCREs, taking into account their corresponding feature spaces separately. Note that the embedding turns more cluster-friendly at each iteration. In other words, the density and separation of clusters are increased at each iteration, and clusters can be caught visually in the final embedding in panel H of Figure 1. In AdaGAE, the sparsity parameter initialization, $k_0$, the increment of this parameter from iteration to iteration, $\mathit{\delta k}$, and the number of iterations, $L$, regulate the number of clusters in the final embedding. In our experiment, we tested various parameter configurations and found eight significant GERM clusters with $L=11$, $k_0=350$, $\mathit{\delta k}=25$.

Notice that the original work of Xuelong et al. defined (8) as a function of the Euclidean distances over a unique feature space. We extended the AdaGAE framework to deal with multiple node feature spaces. In particular, at each iteration $l$, we combined the Euclidean distances of multiple node feature spaces - the gene expression and the cCRE activity time-series - to generate a unique distance function that is given to (8) to generate ${\mathbf{P}}^l$. In the next paragraph, we explain how we combined the feature spaces to adapt AdaGAE to heterogeneous graph RL.

Extending AdaGAE to heterogeneous graphs

As explained previously, we modelled the whole set of gene expression and cCRE-activity data as a heterogeneous graph. We proposed to extend the work in Ref. 26 to heterogeneous graph RL. We ran the same process described in the previous paragraph, with some differences. The first difference was the construction of the reference distance function to feed to (8). At each iteration $l$, we combined the Euclidean distances defined in each one of the original feature spaces to form a unique distance function ${\mathbf{D}}^l$ that generated the sparse distribution family ${\mathbf{P}}^l$.

Recall that $R=\left\{G,C_1,C_2,\dots ,C_M\right\}$ is the set of feature spaces defined, where $G$ is the gene-expression feature space and $C_1,C_2,\dots ,C_M$ are the cCRE activity feature spaces, for example ATAC-seq, H3k27ac, aH3k27me3, etc. At each iteration $l$, for each feature space $r\in R$, we computed the Euclidean distances between the original feature vectors in $r$, scaled these distances into the interval $\left[0,1\right]$, and denoted them with $d_{i,j}^r,\forall i,j\in \left|\mathcal{V}\right|$. Notice that the Euclidean distances between the original feature vectors were only defined between same-class elements: The Euclidean distances in the gene-expression feature space $G$, were defined between genes, and the distances in $C_1,C_2,\dots ,C_M$, only between cCREs. Consequently, we set the distance between different-class elements to the maximum value, i.e., 1, for each feature space. Finally, we create a custom linear combination of these distance functions to form a unique distance function ${\widehat{\mathbf{D}}}^l$:

At the beginning of the $l$-th iteration, the Euclidean distances between the embedding vectors of ${\mathbf{Z}}^{\mathbf{l}-\mathbf{1}}$ are denoted by ${\mathbf{D}}_Z^{l-1}$. We combined such information with ${\widehat{\mathbf{D}}}^l$ defined in (11), to generate the unified distance function, ${\mathbf{D}}^l$ as follows:

In (13), ${\alpha }_Z^l$ is the fixed parametric importance weight of the Euclidean distances in ${\mathbf{D}}_Z^{l-1}$. Notice that, in (12), ${\widehat{\mathbf{D}}}_Z^{l-1}$ allows to exploit the high-level information in the data and ${\widehat{\mathbf{D}}}^l$ helps to reduce the loss of information with respect to the original feature spaces at each iteration.

At each iteration $l,\forall l\in \left[0,L-1\right]$, we computed ${\mathbf{D}}^l$ using (12). This computation is represented in panel E of Figure 1. We then used ${\mathbf{D}}^l$ and the corresponding sparsity parameter $k^l$ to compute a connectivity distribution family ${\widehat{\mathbf{P}}}^l$ using (8). Notice that this distribution family was generated using only the same-class relation types in ${\mathcal{T}}_{\mathrm{ℰ}}-\left\{S_{\mathit{BP}|T}\right\}$. In other words, we only used the node feature spaces in $R$ to create ${\widehat{\mathbf{P}}}^l$. We needed to add the trend-aware base-pair proximity relationship $S_{\mathit{BP}|T}$ defined in (3) to compute the definitive distribution family ${\mathbf{P}}^l$ that takes into account all information we have:

After computing ${\widehat{\mathbf{P}}}^l$, we computed $Q^l$ using (9), and optimized the embedding matrix $\mathbf{Z}$, minimizing a reconstruction loss of the form (6). AdaGAE proposed the use of (10) to implement the reconstruction loss. In this paper, however, we propose to extend such a loss function to better separate elements that are distant in the original feature spaces, as explained in the next paragraph.

Regularization of the loss function

The minimization of (10) has the same effect of an attractive force that tends to collapse points in the embedding space, reducing the distances between $z_j$ and $z_j$ proportionally to $p_{i,j}^l$. When considering a single node $v_i$, such an objective function could lead to errors in the manipulation of the position of the most distant neighbors of such node.⁴⁷ Inspired by Ref. 48, in this work, we added a regularization term that acts as a repulsive force. Such a repulsive force tends to increase ${\widehat{d}}_{i,j}^l$ proportionally to the inverse of $p_{i,j}^l$. In other words, we pushed elements away from each other proportionally to their distance ${\mathbf{D}}^l$ in (12). To create this force, we added the KL-divergence of $1-{\mathbf{P}}_i^l$ with respect to $1-{\mathbf{Q}}_i^l$ to (10).

Notice that the union of the attractive and repulsive force is equivalent to the minimization of the binary cross-entropy of ${\mathbf{P}}^l$ with respect to ${\mathbf{Q}}^l$ as defined in Ref. 48. Lastly, we added parametric importance weights to the attractive and repulsive force components and constructed our final objective function:

We iteratively optimized the parameters of our GNN-based encoder. At each iteration, we minimized (15) with a different set of parameters. After running the optimization for a predefined number of iterations, we processed to a clustering-friendly embedding where the clusters reflect plausible gene regulatory networks. We then ran k-means clustering on the embedding to find plausible gene expression regulatory mechanisms. This clusterization is represented in panel H of Figure 1.

DeepReGraph generated high-quality co-clusters

We employed a principle component analysis (PCA) to visualize the cCRE clusters on dimension reduced plots as shown in panel A of Figure 2. This figure highlights a clear separation of clusters on $\mathit{PC}0$, $\mathit{PC}1$, and $\mathit{PC}2$ which explains the 57% variability in the cCRE datasets. We also performed k-means clustering using $k=8$. Panel B of Figure 2 shows a confusion matrix that compares cCRE clusters from DeepReGraph with the uni-modal clustering produced with k-means clustering. It is clear from this confusion matrix that the result of the multi-modal clustering of DeepReGraph was highly similar to the uni-modal clustering one. This similarity indicates that DeepReGraph clustering does take into account the same-class pattern similarities when defining clusters. We also visualized the PCA plot of gene expression data in panel C of Figure 2. PC0 by itself explains 76% of the variability in gene expression. Basically, gene expression has two distinct patterns throughout mouse fetal heart tissue development: genes with increasing expression and genes with decreasing expression. Interestingly, these two major patterns have been divided into sub-patterns based on the cCRE clusters that plausibly drive their regulation.

Figure 2. Intersections between single-modality k-means clusters and co-clusters induced by DeepReGraph.

A) Principal component analysis (PCA) reduced dimension of candidate cis-regulatory elements (cCRE) profiles colored by the correspondent DeepReGraph cluster. B) Intersection between only-cCRE agglomerative clustering and cCRE extracted from DeepReGraph heterogeneous clusters. C) PCA reduced dimension of gene expression profiles colored by the corresponding DeepReGraph cluster.

The clustered patterns produced by DeepReGraph are presented in Figure 3. In this figure, the y-axis contains the mean-centered values for gene expression in the left-most column, while mean-centered CRE features are in the last three columns. The mean value reduction process mentioned was done as follows: given a time-series vector $\mathbf{x}=\left[x_0,x_1,\dots ,x_t\right]$, with a mean value $\widehat{x}=\frac{{\sum }_{i=0}^t{x}_i}{\left|\mathbf{x}\right|}$, the mean-reduced vector is $\tilde{\mathbf{x}}=\left[x_0-\widehat{x},x_1-\widehat{x},\dots ,x_t-\widehat{x}\right]$. The x-axes of the plots instead correspond to the considered time-points of mouse fetal heart development. Notice that each cluster contains a set of genes and cCREs. The area between the first and 0.75 quantile is colored for each trend plot, to help visualize the trend of each cluster.

Figure 3. Co-clusters identified by DeepReGraph and the relative profile trends.

The first column of plots contains gene expression time profiles, and the rest of the columns contain enhancer activity time profiles. The space between the first and third quartile for each plot was colored to better show the trend.

DeepReGraph revealed the regulatory signature of mouse fetal heart development

We investigated the enriched function of gene expression clusters and enriched transcription factor binding site motifs of cCRE clusters to decipher the general signature of mouse fetal heart development. Figure 4 summarizes these enrichment analysis. It clearly shows that all gene expression clusters have clear functional annotation and all cCRE clusters entail clear transcription factor binding site motifs.

Figure 4. A) Enriched function of gene expression extracted from DeepReGraph clusters. B) Enriched transcription factor binding site motifs of cCRE for the same clusters.

In general, expression of genes related to cell proliferation functions decrease during mouse fetal development, while expression of genes related to heart functions increase.⁴ However, if we look at the enriched terms for the detected gene expression clusters in panel A of Figure 3, and the pattern of gene expression the cCREs that drive them in Figure 3, we can observe that the story is not so simple. Two of the largest gene expression clusters are cluster2 (192 genes) which represents the heart functional genes (enriched for heart contraction function), and cluster7 (249 genes) which represents the cell proliferation genes. The smaller gene expression clusters have more specific enriched functions. Considering the other gene expression clusters down-regulated during development, cluster0 was more enriched for DNA replication, while cluster3 was enriched for non-heart developmental processes. Similarly, the smaller gene expression clusters up-regulated during development gained specific functional enrichment: cluster1 was enriched in metabolic processes to generate energy for the heart to function. Cluster4 was enriched for ventricular cardiac muscle cell membrane repolarization, and cluster6 was enriched for heart contraction. Cluster5 contained genes enriched for regulation of muscle system process. Here, co-clustering of gene expression and cCREs enabled us to derive smaller and more specific clusters. Otherwise, as it is clear from Figure 2 panel C, the smaller clusters of gene expression with more focused functions could not be deduced from gene expression alone.

Multi-modal clustering can also describe how cCREs can drive gene expression during development. First, smaller gene expression clusters gained more cCREs per genes, as Figure 3 shows. Secondly, cCREs with different epigenetic patterns have been linked to similar pattern of gene expression, as can be seen also in Figure 3. For example, in the major cluster with up-regulated genes during development, i.e. cluster2, the trends of ATAC-seq and H3K27ac increased as expected. Similarly, for cluster7 which entails a large set of down-regulated genes throughout development, the trend of ATAC-seq and H3K27ac decreased. However, the H3K27me3 pattern for both cluster2 and cluster7 showed no changes in Figure 3. We can observe polycomb removal events in cluster1, cluster3, and cluster4 as H3K27me3 levels decreased in these clusters.

Multi-modal clustering can further clarify how gene expression changes through development, as cCRE clusters exhibited clearly enriched motifs in panel B of Figure 4. These enriched motifs can prioritize the candidate transcription factors which drive the development. For example, the MEF2 motif was enriched in the clusters related to up-regulated genes. The only exception was cluster1, for which the AP1 motif was enriched. We can assume that a MEF2 binding transcription factor, and more probably MEF2C, was the major transcription factor which caused the increase in gene expression values. Other binding site motifs were been enriched for the clusters entailing temporally up-regulated genes in Figure 4 panel B. Interestingly, these motifs were highly cluster-specific. A similar dynamic was seen for clusters with temporally decreased gene expression. Although zinc finger motifs were enriched in these clusters, the GATA motif was only enriched in cluster0 and the basic helix-loop-helix (BHLH), nuclear receptor motifs are only enriched in cluster3.

Data and software availability

Source code, Data and Interactive Notebook available at: https://github.com/QwertyJacob/DeepReGraph. Archived source code, data and notebook at time of publication: https://doi.org/10.5281/zenodo.6416055

Data are available under the terms of the Apache License, Version 2.0