MultiGRNFormer: A Transformer-Based Multi-Omics GRN Inference Framework

1. Introduction

Gene Regulatory Networks (GRNs) describe the interactions between transcription factors (TFs) and their target genes that control gene expression levels [1]. Understanding GRNs is crucial for deciphering developmental biology mechanisms, disease pathogenesis, and drug target discovery [2,3].

In recent years, deep learning techniques, based on their powerful nonlinear modeling capabilities, have demonstrated significant potential in the field of GRN inference [4,5]. However, existing deep learning-based approaches face a fundamental limitation: the extreme scarcity of high-quality labeled data [5,6,7]. Most datasets contain only hundreds to thousands of TF-gene regulatory pairs, significantly constraining the learning capacity of deep learning models. More critically, existing methods generally follow a “single-sample single-label” training paradigm, failing to fully exploit the latent information within limited datasets. Additionally, their insensitivity to input feature order makes it difficult to capture the dynamic and sequential associations in multi-omics data. For example, methods such as CNNC [8] and DGRNS[7] convert input data into image representations, which provide intuitive spatial relationships but lose positional information, which however is crucial for fully leveraging data and enhancing model performance.

Data augmentation plays a crucial role in deep learning, especially when the available data is limited. By applying random transformations to the original data—such as rotation, scaling, cropping, flipping, and adding noise—it effectively expands the dataset size, enhances the model’s generalization ability, and reduces the risk of overfitting [9]. Additionally, data augmentation helps models learn more robust features and improves their adaptability to complex scenarios. It has been widely applied in various fields, including computer vision, natural language processing, and bioinformatics, serving as a key technique for boosting the performance of deep learning models [10,11,12]. However, despite its extensive research and application in other domains, data augmentation remains underexplored in the field of GRN inference.

Furthermore, groundbreaking advances in single-cell multi-omics technologies now enable the simultaneous profiling of the transcriptome (scRNA-seq) and epigenome (scATAC-seq) within individual cells [13,14]. This provides a more comprehensive perspective for understanding gene transcription and regulation at the cellular level. Current research has demonstrated that integrating gene expression data with chromatin accessibility data for GRN inference can better distinguish direct from indirect regulatory relationships, reduce false positives and false negatives in inference results, and improve overall accuracy [15,16,17].

This study proposes MultiGRNFormer, a multi-omics inference framework based on the Transformer architecture [18]. This framework introduces two key innovations:

Multi-omics modeling of gene regulatory features: Unlike traditional single-omics inference models that rely only on transcriptomic data, MultiGRNFormer integrates both transcriptomic and chromatin accessibility data. By leveraging a Transformer encoder, our model simultaneously captures regulatory information from gene expression and chromatin accessibility to infer GRNs more accurately.

Data augmentation for GRN inference: We introduce positional encoding to make the model sensitive to input feature order. Additionally, we apply dynamic input reordering to generate multiple training samples from single-label data, effectively augmenting the dataset and improving learning efficiency.

To evaluate the performance of MultiGRNFormer, we conducted experiments on seven different datasets, comparing our approach with existing single-omics inference models. Experimental results across these datasets demonstrate that MultiGRNFormer achieves superior performance in GRN inference. Furthermore, our experiments confirm that the proposed data augmentation method significantly enhances model prediction performance. We also explored the feasibility of cross-dataset joint training. Cross-dataset experiments revealed that training across datasets yields significant advantages on most benchmarks, indicating that increased data diversity helps the model extract universal regulatory patterns across tissues. Finally, ablation studies confirm that integrating multi-omics data provides MultiGRNFormer with richer regulatory information, thereby improving inference accuracy.

2. Materials and Methods

2.1. Datasets

The study used seven publicly available datasets, which consist of scRNA-ATAC-seq data from different tissues of both human and mouse [19,20,21].

In the field of GRN, there has been a persistent lack of real-world networks for model evaluation. As a result, in the GRN inference literature, a common and expedient practice is to evaluate the accuracy of a resulting network by comparing its edges to an appropriate database of TFs and their targets. Therefore, We utilized the non-cell type-specific transcriptional regulatory network ChIP-seq data provided by BEELINE[22] and the functional interactions from the STRING[23] database as ground truth networks. Statistical data are shown in Table 1.

2.2. Data Preprocessing

We first preprocessed the raw count matrices for scRNA-seq and scATAC-seq of the dataset. For each count matrix, we denoted rows as features (genes or peak regions) and columns as cells throughout the paper below. Each data matrix was removed if a row or column contained less than 0.1% non-zero values. Data quality control was performed by Seurat V4, including but not limited to total read counts, mitochondrial gene ratios, and blacklist ratios [24]. After quality control, we selected the top 1,000 highly variable genes for subsequent analyzes based on gene variance, and the final gene expression matrix and chromatin accessibility matrix were obtained. The gene expression matrix was denoted as ${\mathit{X}}^R=\left\{x_{ij}^R\right|i=\mathrm{1,2},\dots ,I;j=\mathrm{1,2},\dots ,J\}$ , with totally $I$ genes and $J$ cells. The chromatin accessibility matrix was denoted as ${\mathit{X}}^A=\left\{x_{kj}^A|k=\mathrm{1,2},\dots ,K;j=\mathrm{1,2},\dots ,J\right\}$, which has $K$ peak regions in $J$ cells.

To apply the chromatin accessibility data to GRNs inference, we used the method described in MAESTRO[25] to calculate the peak region in $X^A$ as the regulatory potential of the corresponding gene. Specifically, based on the distance between peak $k$ and gene $i$ in the genome, the regulatory potential weight of peak $k$ for gene $i$ is denoted as $w_{ik}$ and calculated as follows:

$$w_{ik}=\left\{\begin{array}{l}0, {\mathrm{d}}_{\mathrm{i}\mathrm{k}}>150\mathrm{kb} \mathrm{or} \mathrm{peak} k \mathrm{located} \mathrm{in} \mathrm{any} \mathrm{nearby} \mathrm{genes}\\ {\displaystyle \frac{1}{\mathrm{L}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}\left(\mathrm{e}\mathrm{x}\mathrm{o}\mathrm{n}\right)}}, \mathrm{peak} k \mathrm{located} \mathrm{at} \mathrm{the}\mathrm{exon} \mathrm{regions} \mathrm{of} \mathrm{the} \mathrm{gene} i\\ 2^{-{\displaystyle \frac{d_{ik}}{d_0}}},\mathrm{otherwise}\end{array}\right.$$

(1)

where $d_{ik}$ denotes the distance from the center of peak $k$ to the transcriptional starting site of gene $i$, and $d_0$ is the half-decay of that distance (set to 10kb). The regulatory potential $w_{ik}$ of peak $k$ for gene $i$ is usually calculated as $2^{-{\displaystyle \frac{d_{ik}}{d_0}}}$. If $d_{ik}>150kb$, $w_{ik}$ should be less than 0.0005, and in order to save computational time, we set it to 0. In MAESTRO, to better fit the gene expression model, if peak $k$ is located in the exon region of the gene, $w_{ik}$ should be 1 according to the formula. But since reads tend to be located in longer exons than shorter exons, to normalize the likelihood of background reads, the total exon reads were normalized by the total exon length of each gene exon. The peak $k$ regulatory potential of gene $i$ in cell $j$ can then be calculated as:

$$r_{ik|j}=w_{ik}\times x_{kj}^A$$

(2)

Finally, the scATAC-seq matrix ${\mathit{X}}^A$ is transformed into the gene regulatory potential matrix ${\mathit{X}}^P$ by summing the regulatory potential scores of peaks regulating the same gene:

$$x_{ij}^P=\sum _k{r}_{ik|j}$$

(3)

2.3. The MultiGRNFormer framework

MultiGRNFormer is inspired by the design concept of STGRNS[4], utilizing a Transformer architecture to capture intercellular gene expression variations. However, unlike STGRNS, which aggregates multiple cells into motifs, we designed a novel input structure to enhance the model’s sensitivity to intercellular gene expression changes. Additionally, our input design enables MultiGRNFormer to capture the regulatory state of genes within each cell by integrating multimodal data. Specifically, in addition to transcriptomic data traditionally used by existing models, we incorporate chromatin accessibility data to provide an alternative perspective on the gene regulatory process. MultiGRNFormer consists of four key components: an Input Constructor , a Positional Encoder, a Multiple Transformer Encoders, and a Regulation Predictor , as illustrated in Figure 1c.

Input Data: MultiGRNFormer takes two modalities of input data: the gene expression matrix $X^R=\{x_{ij}^R\left|i=\mathrm{1,2},\dots ,I;j=\mathrm{1,2},\dots ,J\right\}$ and the gene regulatory potential matrix $X^P=\{x_{ij}^P\left|i=\mathrm{1,2},\dots ,I;j=\mathrm{1,2},\dots ,J\right\}$, where $I$ denotes the number of genes, and $J$ represents the number of cells.

Input Constructor: The format of data input plays a critical role in maximizing the model’s feature extraction capabilities. Thus, we carefully designed an Input Constructor to structure the input data effectively. For a given transcription factor (TF) $a$ and target gene $b$, their gene expression values in the $j$th cell are $x_{aj}^R\in X^R$ and $x_{bj}^R\in X^R$, respectively, while their regulatory potential values are $x_{aj}^P\in X^P$ and $x_{bj}^P\in X^P$. To ensure the model accurately captures the regulatory state of a gene within a single cell, we concatenate all relevant data from the same cell into a single feature vector: $X_j^{ab}=\left\{x_{aj}^R,x_{aj}^P,x_{bj}^R,x_{bj}^P\right\}$. We retain only the cells where both TF $a$ and gene $b$ are expressed to minimize the impact of missing values. Then, we randomly sample $N$ cells from this set (with replacement if the available number of cells is less than $N$). Finally, the model input is structured as: $X^{ab}=\left\{X_j^{ab}|j=\mathrm{1,2},\dots ,N\right\}$.

Positional Encoding: Since the Transformer encoder lacks an inherent understanding of positional information in the input data, we apply positional encoding to $X^{ab}$ in order to obtain $X_{pos}^{ab}$, computed as follows:

$$X_{pos}^{ab}=X^{ab}+\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\left(X^{ab}\right),$$

(1)

where $\mathrm{P}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}(·)$ is the positional encoding function. We adopt a sinusoidal positional encoding scheme, where sine functions are applied to odd-indexed elements and cosine functions are applied to even-indexed elements:

$$PE\left(m,2n\right)=\sin \left(m/{10000}^{2n/4}\right),$$

(2)

$$PE\left(m,2n+1\right)=\cos \left(m/{10000}^{2n+1/4}\right),$$

(3)

where $m$ denotes the position of a cell in $X^{ab}$, and $n$ represents the feature index, even indices corresponding to $2n$ and odd indices to $2n+1$.

Multiple Transformer Encoders: We employ a multi-layer Transformer encoder [18] with residual connections to extract multi-omics features that define gene regulatory interactions. The self-attention mechanism in the Transformer encoder allows the model to focus more on cells that provide strong regulatory evidence. The hierarchical multi-layer structure enhances the model’s ability to capture higher-order interactions, while residual connections mitigate the issues of gradient vanishing and explosion commonly observed in deep networks.

Specifically, our encoder consists of multiple identical Transformer layers, each comprising two sublayers: multi-head self-attention and a position-wise feed-forward network, as illustrated in Figure 1b. The transformation between layers is formally defined as:

$$H^{\left(l\right)}=\mathrm{L}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}\mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\left(H^{\left(l-1\right)}+\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{d}\left(H^{\left(l-1\right)}\right)\right)$$

(4)

$$H^{\left(l\right)}=\mathrm{L}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}\mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\left(H^{\left(l\right)}+\mathrm{F}\mathrm{F}\mathrm{N}\left(H^{\left(l\right)}\right)\right)$$

(5)

where $H^{\left(l\right)}\in {\mathbb{R}}^{{N\times d}_{model}}$ denotes the feature representation at the $l$th encoder layer. In our model, $d_{model}=4$, since each sample includes two gene expression values and two regulatory potential values.

$$hea{d}_i=\mathrm{A}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\left(Q{W}_i^Q,K{W}_i^K,V{W}_i^V\right),$$

(6)

where $W_i^Q,W_i^K,W_i^V\in {\mathbb{R}}^{d_{model}\times d_k}$ are the query, key, and value projection matrices, respectively, with $d_k=d_{model}/h$, in which $h$ is the number of attention heads. The scaled dot-product attention is computed as:

$$\mathrm{A}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\left(Q,K,V\right)=\mathrm{s}\mathrm{o}\mathrm{t}\mathrm{f}\mathrm{m}\mathrm{a}\mathrm{x}\left({\displaystyle \frac{Q{K}^T}{\sqrt{d_k}}}\right)V$$

(7)

The outputs from all attention heads are concatenated and linearly projected to form the final representation:

$$\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{d}\left(Q,K,V\right)=\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{a}\mathrm{t}\left(head1,\dots ,hea{d}_h\right)W$$

(8)

The position-wise feed-forward network (FFN) consists of two linear transformations with a $\mathrm{G}\mathrm{E}\mathrm{L}\mathrm{U}$ activation function:

$$\mathrm{F}\mathrm{F}\mathrm{N}\left(x\right)=W_2\left(\mathrm{G}\mathrm{E}\mathrm{L}\mathrm{U}\left(W_{1}x+b_1\right)\right)+b_2,$$

(9)

where $W_1\in {\mathbb{R}}^{d_{model}\times h{d}_{model}}$ expands the dimension by a factor of 4, and $W_2\in {\mathbb{R}}^{{h{d}_{model}\times d}_{model}}$ projects it back to the original size.

Regulation Predictor: We employ a multi-layer perceptron (MLP) as the regulation predictor. First, the regulatory features obtained from the Transformer encoder, $H\in {\mathbb{R}}^{N\times 4}$, are flattened into a one-dimensional vector:

$$H^{\prime }=\mathrm{F}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{n}\left(H\right)$$

(10)

Then, $H^{\prime }$ is passed through the MLP to predict the regulatory probability $p_{ab}$ between TF $a$ and gene $b$:

$$p_{ab}=\sigma \left(W_2·\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left(W_1{H}^{\prime }+b_1\right)+b_2\right),$$

(11)

where $\sigma$ denotes the $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{d}$ activation function.

Parameter Optimization: To train the model, we minimize the Binary Cross Entropy (BCE) loss between the predicted regulatory interactions and ground truth labels, enabling the model to distinguish fine-grained regulatory relationships:

$$BCE=-{\displaystyle \frac{1}{N}}\sum _{i=1}^N y_i\cdot \mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}\left(x\right)\right)+\left(1-y_i\right)\cdot \mathrm{l}\mathrm{o}\mathrm{g}\left(1-\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}\left(x\right)\right)$$

(10)

where $i$ denotes the $i$th TF-gene pair and $y_i$ denotes the label of the $i$th TF-gene pair.

2.4. Experimental Setting and Hyperparameter Optimization

For each real network, we use 4/5 of the set of positive labels as the training set and 1/5 of the positive labels as the test set. Due to the high imbalance of the dataset (as in Table 1), we randomly sampled negative individuals in the unobserved links. The training set was constructed with an equal number of positive and negative individuals. For the test set, we employed the following sampling strategy to simulate real-world sparse regulatory networks:

$${\displaystyle \frac{N_{positive}}{{N_{positive}+N}_{negative}}}={\displaystyle \frac{\left|E\right|}{\frac{1}{2}{N}_{TF}\times N_{gene}}},$$

(11)

where $\left|E\right|$ represents the number of edges in the regulatory network of the dataset, $N_{TF}$ and $N_{gene}$ denote the number of transcription factors (TFs) and genes, respectively. $N_{positive}$ and $N_{negative}$ represent the number of positive and negative individuals in the test set. This widely adopted approach aims to evaluate model performance under biologically realistic conditions where interaction networks are extremely sparse.

The baseline methods included in this paper are:

• CNNC [8]: predicts GRNs using deep convolutional neural networks.
• STGRNS [4]: a supervised learning method based on Transformer architecture.
• GENIE3 [26]: an unsupervised learning method based on random forests that constructs GRNs using regression coefficient weights.
• GRNBoost2 [27]: an unsupervised learning method for GRN inference using random gradient boosting regression and early stopping regularization.

For a fair comparison of the models, all models were trained on the same training sets and model performance was compared on the same test sets. Due to the high imbalance of government nutrients in the test set, we use the AUPRC metric, in addition to the AUROC metric, to measure the ability of the models to identify positive individuals.

3. Results

3.1. Parameter analysis

MultiGRNFormer is influenced by several key parameters: the number of sampled cells $N$, the number of encoder layers, and the hidden layer size of the MLP. To determine an optimal set of parameters, we conducted a parameter analysis on two datasets, Human-bone and Human-breast. We selected these datasets of different sizes and performed five-fold cross-validation on the training sets of each dataset using various parameter combinations. Specifically, the training set was equally divided into five subsets, with one subset iteratively held out as the validation subset while the others were used for training. Final performance metrics were averaged across all validation rounds.

The optimal parameter combinations differed significantly between the two datasets: for Human-bone, the optimal combination was that $N=1000$, encoder layers is 2, and MLP hidden size is 1024; while for Human-breast, it was that $N=500$, encoder layers is 1, and MLP hidden size is 512. We hypothesize that this discrepancy is related to dataset size. To select an appropriate set of parameters, we averaged the results across different parameter settings, as shown in Figure 2. The results indicate that the model is not highly sensitive to the number of encoder layers and the hidden layer size of the MLP, as different values for these parameters had minimal impact on performance. However, the number of sampled cells $N$ had a significant impact. The optimal value for $N$ was 800 for Human-bone and 500 for Human-breast, likely due to the larger dataset benefiting from a higher $N$. The number of encoder layers exhibited similar trends across both datasets, with one layer performing the worst and two to three layers yielding competing results. The hidden layer size of the MLP showed consistent performance across all tested values. Given the fact that most datasets are similar in size to or smaller than Human-breast, we adopted the parameter combination of $N=500$, encoder layers is 2, and MLP hidden size is 128 for all experiments unless otherwise specified.

3.2. Performance of MultiGRNFormer in Gene Regulatory Network Inference

To evaluate the performance of our model, we conducted a comprehensive benchmark test against four established methods across seven datasets, with results shown in Figure 3. MultiGRNFormer achieved the best performance on most datasets, particularly excelling in the AUPRC metric, where it outperformed all other models across all datasets. However, on some datasets, the AUROC metric approached random comparison outcome. This phenomenon may be attributed to the highly imbalanced test sets, where the overwhelming presence of negative samples allows AUROC to be dominated by the true negative rate (TN). As a result, slight model performance improvements may not be reflected in the AUROC scores. In the TFS+1000 tests, MultiGRNFormer outperformed others on both AUROC and AUPRC. Specifically, our model achieved higher AUROC scores than in the TFS+500 setting, whereas AUPRC scores exhibited a slight decline. This suggests that an increase in training data improved model performance in AUROC, but the more imbalanced nature of the TFS+1000 test set counteracted the benefits by reducing recall.

3.3. Enhancing MultiGRNFormer Performance through Data Augmentation

A primary challenge in deep learning-based gene regulatory network inference is the limited availability of labeled data, which constrains the scale of deep learning models. Deep learning models are highly sensitive to training data size, with larger training sets generally leading to better performance. As shown in Figure 3, all deep learning models achieved the best results on the largest dataset (human bone). MultiGRNFormer incorporates a positional encoding module, enabling the model to be sensitive to input order. Leveraging this feature, we designed a novel data augmentation strategy: without altering the gene regulatory relationships, we generate multiple training samples by randomly shuffling the input order of a TF-gene label. This approach modifies the relative positions of cells in the input data while preserving the regulatory relationships between genes, thereby expanding the size of the training dataset.

We expanded the training sets from three datasets by different augmentation factors and evaluated performance using their respective test sets. Applying this approach in the TFS+1000 setting, as shown in Figure 4, we observed performance improvements as the augmentation factor increased. For the largest dataset, Human-bone, applying a fivefold data augmentation led to a 5.6% increase in AUROC and a 2.2% increase in AUPRC compared to no augmentation. For Human-breast, the respective improvements were 3.7% and 2%. On the smallest dataset, Human-jejunum, augmentation effects were negligible for factors of 2–4 but became significant at a factor of 5. We further applied fivefold data augmentation across all seven datasets (Figure 5) and found consistent performance improvements, suggesting that our model captures regulatory relationships by leveraging gene expression variations across different cells.

3.4. Cross-Dataset Learning to Improve MultiGRNFormer Performance

Based on the effectiveness of data augmentation, we investigated whether MultiGRNFormer could leverage features learned from other datasets to enhance performance on a target dataset. We designed an experiment where the model was pre-trained on multiple datasets (25 epochs) before being fine-tuned on the target dataset (10 epochs) and evaluated on the target test set. As shown in Figure 6a, this approach yielded limited improvements, with most datasets showing no significant advantage. We hypothesize that the current training data scale does not sufficiently mitigate batch effects across datasets. To address this, we repeated the experiment using fivefold data augmentation for all training datasets. The results in Figure 6b indicate that the combined approach (cross-dataset & Augmented 5) significantly improved model performance compared to augmentation alone. This suggests that integrating data augmentation and cross-dataset training not only mitigates batch effects but also enables the model to learn more generalizable gene regulatory relationships.

On most datasets, the cross-dataset & Augmented 5 approach substantially improved model performance, with AUROC increasing by an average of 4.7% and AUPRC by 2.7%. However, performance gains varied across species. In the five human datasets, AUROC improved by an average of 6.4% and AUPRC by 3.7%. In contrast, the two mouse datasets showed minimal improvement (AUROC +0.5%, AUPRC +0.1%), with a performance decline observed in the mouse-brain dataset. We attribute this to two factors: (1) the mouse datasets were significantly smaller than the human datasets, containing only a few hundred regulatory interactions, which remained insufficient even after augmentation, and (2) interspecies differences likely exceeded batch effects, which our model struggled to account for.

We conducted additional experiments to explore the feasibility of training a universal model capable of inferring gene regulatory networks across datasets. In these experiments, we trained models on the six remaining datasets (with fivefold data augmentation) and directly evaluated them on the target dataset’s test set, without fine-tuning on its training set. The resulting models performed only slightly better than random predictions. This suggests that training on the target dataset remains crucial for learning dataset-specific batch effects, indicating that our current model cannot fully eliminate these effects.

3.5. Leveraging Multi-Omics Data to Enhance MultiGRNFormer Performance

In all previous experiments, we utilized multi-omics datasets for model training. To assess whether our model benefits from multi-omics data, we conducted experiments using only transcriptomic data, with results shown in Figure 7. In most cases, multi-omics data outperformed transcriptomics alone, particularly in the data augmentation (Figure 7b) and cross-dataset training (Figure 7c) experiments, where AUROC increased by 1.4% and 2.5%, respectively, and AUPRC increased by 1.1% and 1.9%. However, in standard training (Figure 7a), multi-omics data did not provide a significant advantage. This suggests that without data augmentation or cross-dataset training, the traditional “single-sample, single-label” method is insufficient to capture dynamic correlations across omics layers. By contrast, data augmentation and cross-dataset training provided sufficient samples, allowing MultiGRNFormer to leverage multi-omics data more effectively, leading to superior performance in these experimental settings.

4. Discussion

This study systematically validates the superior performance of MultiGRNFormer newly proposed in gene regulatory network (GRN) inference driven by multi-omics data. Experimental results demonstrate that by integrating transcriptomic (scRNA-seq) and chromatin accessibility (scATAC-seq) data, a well-trained model can effectively capture dynamic changes in gene expression and the open states of regulatory elements, thereby reducingfalse positives in unimodal inference. In data augmentation and cross-dataset joint training scenarios, the multi-omics model achieved average AUPRC improvements of 1.1% and 1.9%, respectively, compared to using transcriptomic data alone. To address the challenge of limited labeled data, we propose an input sequence randomization-based augmentation strategy, which enhances data utilization efficiency by generating multiple perspectives from a single sample. This technique increases the effective use of limited labeled data by fivefold, leading to AUROC and AUPRC improvements of 5.6% and 2.2%, respectively, on the largest dataset (human bone). Meanwhile, for smaller datasets such as human-jejunum, the enhancement effect follows a nonlinear accumulation pattern, highlighting the role of data augmentation in improving the generalization ability of deep learning models in GRN inference.

Further investigation reveals that the synergy between cross-dataset training and data augmentation can overcome the limitations of single-dataset training. When combined with data augmentation, cross-dataset training improves AUROC by an average of 6.4% and AUPRC by 3.7% across five human datasets, indicating that the model can capture regulatory patterns across different tissues. However, cross-species training (e.g., human-to-mouse) does not perform well and may even lead to performance degradation, underscoring the challenge of improving model performance through cross-species training.

Although fully cross-dataset “zero-shot” inference has not yet been achieved, the combination of data augmentation and cross-dataset training has demonstrated superior performance over single-dataset training in most cases, suggesting that MultiGRNFormer has the potential to learn generalizable regulatory patterns. The model performs particularly well in multi-omics-driven scenarios. While the advantage of multi-omics is limited under conventional training, its performance surpasses unimodal approaches significantly when combined with data augmentation and cross-dataset training, confirming that leveraging complementary multimodal features requires sufficient training samples. These findings provide new insights for future GRN inference model development.

This study demonstrates the superior performance of MultiGRNFormer in multi-omics-driven gene regulatory network (GRN) inference; however, several promising research directions remain open for future exploration. Subsequent work may proceed along the following avenues: First, to construct a more comprehensive and accurate representation of the regulatory landscape, future versions of MultiGRNFormer could be extended to incorporate additional modalities—such as DNA methylation [28], histone modifications [29], and 3D chromatin conformation [30]. A major challenge will be to develop effective methods for aligning and integrating these heterogeneous data types within a unified deep learning framework. Second, future research could investigate more advanced attention mechanisms or graph neural network architectures to better capture the complex, non-linear interactions between transcriptomic and epigenomic data.