Non-paired single cell multi-omics data gene regulatory network inference method

Abstract

The invention discloses a non-paired single cell multi-omics data gene regulatory network inference method, which comprises the following steps: collecting single cell sequencing non-paired data, and preprocessing the data; inputting the pre-processed single cell sequencing non-paired data into a multi-layer variational auto-encoder structure to generate advanced potential variables and reconstructed single cell sequencing non-paired data; training a layered GAN by adopting a layered adversarial alignment mechanism according to the advanced potential variables and the reconstructed single cell sequencing non-paired data; according to the advanced potential variables, adopting a mutual nearest neighbor strategy to train a mutual nearest neighbor module; fusing the advanced potential variables of different modes by adopting a gating fusion mechanism, constructing a unified advanced potential variable, and completing adaptive feature integration of the multi-mode advanced potential variables; according to prior information, an initial adjacency matrix is established, a gene regulation network is constructed in combination with unified advanced potential variables, and the gene regulation network with biological rationality is directly deduced from non-pairing input.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and in particular relates to a method for inferring gene regulatory networks from unpaired single-cell multi-omics data. Background Art

[0002] Normal cellular function relies on the precise regulation of gene transcription. This process enables cells to sense and respond to changes in the external environment and coordinates key biological processes such as cell proliferation, differentiation, and apoptosis. Transcription factors (TFs) are core regulators of this regulatory system. They recognize and bind to transcription factor binding sites (TFBSs) with specific sequence patterns in the genome, thereby regulating the promoters or enhancer regions of target genes, influencing their expression. Furthermore, TFs often form complexes with cofactors and chromatin remodelers, participating in the regulation of chromatin structure and transcription initiation. To systematically characterize these regulatory mechanisms, researchers often use graph-based approaches to construct directed networks of TFs and their regulated targets, namely gene regulatory networks (GRNs). GRNs are hierarchical, sparse, and dynamic in structure, making them essential tools for understanding cell fate determination, lineage differentiation, and disease states.

[0003] In recent years, with the rapid development of next-generation sequencing (NGS) and single-cell multi-omics platforms, researchers have been able to obtain multi-level information at single-cell resolution, including transcription, chromatin, and protein expression. This has enabled high-dimensional, mechanism-driven modeling of gene regulatory networks. Numerous methods have been proposed for inferring gene regulatory networks, primarily encompassing two approaches. First, co-expression-based statistical modeling methods, such as WGCNA, ARACNe, and GENIE3, calculate regulatory relationships between genes based on correlations or mutual information within expression matrices. While computationally efficient, these methods generally fail to capture true causal regulatory relationships or integrate spatial chromatin state information. Second, multi-omics integration approaches, incorporating prior knowledge of mechanisms, construct regulatory networks by integrating information such as transcription factor binding motifs, chromatin accessibility, and expression profiles. Methods such as SCENIC, scMTNI, and Dictys incorporate motif enrichment analysis, transcription factor activity scores, Bayesian modeling, or meta-learning mechanisms, enhancing the biological interpretability of the networks. However, existing technologies generally have problems such as poor scalability, insufficient generalization ability, and high requirements for modal pairing, making them difficult to adapt to real-world scenarios with large-scale, unpaired data.

[0004] In actual single-cell sequencing experiments, due to sample heterogeneity, sequencing efficiency differences and experimental costs and other factors, paired modal data of the same cell cannot be obtained, resulting in a lack of sample one-to-one correspondence between scRNA-seq and scATAC-seq. This limitation makes it impossible to directly apply a large number of existing methods, thereby giving rise to two types of non-paired multi-omics integration strategies: one is the "cut-and-match" strategy, which maps different modalities to the same latent space through joint dimensionality reduction, and uses anchor cells or shared features to achieve modality alignment, such as the joint PCA method proposed by Buenrostro et al. to reduce noise to some extent, but often at the expense of modality-specific information; the other is the integration method based on deep generative model, which learns the latent representation from non-paired data by introducing variational autoencoder (VAE), generative adversarial network (GAN) and other mechanisms, and realizes modality alignment and missing modality completion. Existing methods such as scCross combine VAE and GAN framework to realize embedding alignment between non-paired single-cell modalities, but their inference range is limited to the cell-level latent space, and they cannot analyze the gene-level regulatory edges or form clear gene regulatory network structures.

[0005] At the same time, existing methods generally have low computational efficiency, insufficient integration ability of biological priors, and lack of interpretability of output network structure when dealing with high-dimensional multi-omics data. Most frameworks rely on fixed motif prior graph or nearest neighbor strategy based on similarity, making it difficult to adapt to different data distributions, cell type diversity or dynamic changes in developmental trajectories. Especially in the face of real data without pairing information, there is still a lack of a unified, universal and scalable computational framework that can effectively integrate non-paired modalities while outputting sparse, semantically clear and biologically interpretable gene regulatory networks. Therefore, how to construct a regulatory network inference method for non-paired single-cell multi-omics data, which integrates omics features, introduces prior knowledge and has structural interpretability without relying on paired samples, is still a key problem that needs to be solved in current technology. SUMMARY

[0006] In view of the above problems in the prior art, the non-paired single-cell multi-omics data gene regulatory network inference method provided by the present application solves the problems of insufficient integration of non-paired data, lack of structural interpretability of regulatory edge inference and difficulty in effectively utilizing prior information in the prior art.

[0007] In order to achieve the above-mentioned purpose of the application, the technical scheme adopted by the present application is as follows: a non-paired single-cell multi-omics data gene regulatory network inference method, comprising the following steps: S1, collect single-cell sequencing non-paired data of human hematopoietic cell differentiation, preprocess the single-cell sequencing non-paired data to generate preprocessed single-cell sequencing non-paired data; S2, input the preprocessed single-cell sequencing non-paired data into a multi-layer variational autoencoder structure to generate high-level latent variables and reconstructed single-cell sequencing non-paired data; S3, train a hierarchical GAN using a hierarchical adversarial alignment mechanism according to the high-level latent variables and the reconstructed single-cell sequencing non-paired data, and the trained hierarchical GAN is used to guide the alignment of high-level latent variables of different modalities in an abstract space; S4, train a nearest neighbor module using a nearest neighbor strategy according to the high-level latent variables, and the trained nearest neighbor module is used to constrain the consistency of local structures between modalities in an embedding space; S5, fuse the high-level latent variables of different modalities using a gating fusion mechanism to construct unified high-level latent variables and complete adaptive feature integration of multi-modal high-level latent variables; S6, establish an initial adjacency matrix according to prior information ， Construct a gene regulatory network based on the unified high-level latent variables to complete gene regulatory network inference.

[0008] Further, in S1, the single-cell sequencing non-paired data includes scRNA-seq data and scATAC-seq data. The preprocessing method includes row normalization, principal component analysis, and high-variation feature screening processing.

[0009] Further, in S2, the multi-layer variational autoencoder structure includes a first-level encoder, a second-level encoder, and a zero-inflated decoder; S2 includes the following steps: S21, input the count matrix of the single-cell sequencing non-paired data into the first-level encoder to generate omics-specific latent variables; S22, input the omics-specific latent variables into the second-level encoder to generate high-level latent variables; S23, input the omics-specific latent variables into the zero-inflated decoder to generate reconstructed single-cell sequencing non-paired data.

[0010] Further, in S21, the expression of the generated omics-specific latent variables is specifically: In the formula, k is the omics, RNA represents the RNA modality, and ATAC represents the ATAC modality, is a softplus activation function, is a mean value learned by the first-level encoder, a scale vector learned by the first-level encoder, a Hadamard product, a first Gaussian noise for differentiable sampling; wherein, a first-level encoder, a count matrix of the single-cell sequencing unpaired data, denotes a multi-dimensional standard normal distribution with mean 0 and variance 1 in each dimension and independent between dimensions, a regularization, a batch normalization, a LeakyReLU activation function, a linear layer; In S22, the expression of the generated high-level latent variable is specifically: wherein, a mean learned by the second-level encoder, a scale vector learned by the second-level encoder, a second Gaussian noise for differentiable sampling; wherein, a second-level encoder, a tanh activation function; In S23, the expression of the generated count matrix of the reconstructed single-cell sequencing unpaired data is specifically: wherein, a zero-inflated decoder, a zero-inflated log-normal distribution model, a linear transformation layer for predicting logits corresponding to the zero-inflated probability, used to determine whether an element is structural 0 in the zero-inflated model, a linear transformation layer for predicting the mean parameter μ of the log-normal main distribution in the ZILN distribution, used to determine the center position of the expression value, a linear layer for predicting the scale parameter δ of each gene in the log-normal main distribution, thereby modeling the variability of the expression data,b is an auxiliary variable used to estimate the zero-inflated probability.

[0011] Furthermore, in S3, the hierarchical GAN ​​includes a first discriminator and a second discriminator, and S3 includes the following sub-steps: S31, concatenate high-level latent variables to generate the first modal label, and train the first discriminator through binary cross entropy loss; S32. Connect the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data to generate a second modality label, and train the second discriminator through binary cross entropy loss.

[0012] Further: In S31, the first loss function of the first discriminator is trained by binary cross entropy loss The specific expression is: Where, For the i The modal original labels of samples, N is the sample size, For the i The first modal label of samples, calculate the first modal label The specific expression is: Where, is the result of the connection of high-level latent variables; Where, For the connection operation, is the high-level latent variable of the cell in the RNA modality, It is the high-level latent variable of cells in ATAC mode; In S32, the second loss function of the second discriminator is trained by binary cross entropy loss The specific expression is: Where, For the i The original sequencing data expression level of each sample, For the i The second modal label of samples, calculate the second modal label The specific expression is: Where, is the connection result of the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data, .

[0013] Further, in S4, the method for training the mutual nearest neighbor module is specifically: According to the advanced latent variable, pairs of cells that are mutual nearest neighbors of each other are screened, and an alignment loss is calculated according to all pairs of cells that are mutual nearest neighbors of each other, and the mutual nearest neighbor module is trained; wherein the expression of the alignment loss is specifically: In the formula, is the i-th cell in the RNA modality, i is the j-th cell in the RNA modality, is the i-th cell in the ATAC modality, j is the j-th cell in the ATAC modality, MNN is a pair of cells that are mutual nearest neighbors of each other, represents and is a pair of cells that are mutual nearest neighbors of each other.

[0014] Further, in S5, a unified advanced latent variable is constructed, and the expression of the unified advanced latent variable is specifically: In the formula, is a dynamic gating, is the number of omics types; In the formula, is a sigmoid activation function.

[0015] Further, S6 includes the following steps: S61, constructing an initial adjacency matrix based on scRNA-seq data A ; In the formula, is an adjacency matrix representing a prior regulatory network after mapping processing, is a randomly generated adjacency matrix, is a learnable coefficient for balancing the influence of prior knowledge and data-driven initialization; In the formula, is a binary matrix for representing prior TF-to-gene regulatory links, and T is the number of TFs; S62, inputting the initial adjacency matrix, the unified advanced latent variable, and the count matrix of the scRNA-seq data into a multi-layer perceptron, training the multi-layer perceptron using a reconstruction loss, and constructing a gene regulatory network; wherein the expression of the reconstruction loss is specifically: wherein, is the reconstruction evaluation by sampling the latent variable from the output distribution of the encoder, is the output of the decoder, generating the original data with the fused latent variable, is a hyper-parameter for adjusting the weight of the KL divergence term, for balancing the reconstruction loss and the regularization term, is the approximate posterior distribution of the fused latent variable given the input data X, is the prior distribution of the fused latent variable , is a hyper-parameter for regulating the weight of the sparse regularization term to control the sparsity of the model parameters, is the weight of the sparse regularization term, is the L1 norm.

[0016] The beneficial effects of the present application are: (1) The present application provides a non-paired single-cell multi-omics data gene regulatory network inference method, which extracts non-paired modal features through a multi-level variational autoencoder, extracts shared and exclusive latent representation features between multiple modalities, aligns the modal features, and retains the semantic features within the modalities; at the same time, an adversarial learning strategy is combined to guide the alignment of the latent distributions of different modalities in the abstract space. By introducing a gating fusion mechanism to weight and integrate the information of different modalities, a nearest neighbor alignment mechanism is introduced to correct the local structure of non-paired cells, and prior knowledge is injected into the adjacency matrix to integrate the prior knowledge and the nonlinear model, balance the prior domain knowledge and data-driven, effectively improve the structural sparsity and biological interpretability of the regulatory network, and thus accurately infer the transcription factor-target gene regulatory relationship in non-paired single-cell data, solve the problems of insufficient integration ability of non-paired data, lack of structural explanation of regulatory edge inference and difficulty in effectively using prior information in the prior art, and realize direct inference of a gene regulatory network with biological rationality from non-paired input.

[0017] (2) The present application ingeniously handles the heterogeneity between multiple omics data sets while maintaining the internal consistency of a single omics data, providing better performance compared to current methods.

[0018] (3) The present application proves its effectiveness in the gene regulatory network inference task of non-paired single-cell multi-omics data through explainability analysis, demonstrating its effectiveness. Based on the comparative experiments and case analysis designed in the technical solution, the framework proposed in the present application achieves leading results. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1A flow chart of a non-paired single-cell multi-omics data gene regulatory network inference method of the present application.

[0020] Figure 2 Results of cross-modal alignment of scRNA-seq and scATAC-seq by the present application; Figure 2 a is a UMAP plot of scRNA-seq cells, Figure 2 b is a UMAP plot of scATAC-seq cells, Figure 2 c is a joint UMAP embedding plot colored by cell type, Figure 2 d is a joint embedding plot colored by modality.

[0021] Figure 3 A clustering result plot; Figure 3 a is a UMAP plot of original scRNA-seq profiles, Figure 3 b is a UMAP plot of RNA expression profiles predicted from scATAC-seq input by the present application.

[0022] Figure 4 A comparison plot of GLUE-inferred regulatory interaction subnetworks.

[0023] Figure 5 A comparison plot of regulatory interaction subnetworks inferred by the present application.

[0024] Figure 6 A node centrality indicator highlighting key transcription factors in the topology and functional analysis of the regulatory network inferred by the present application Figure 1 .

[0025] Figure 7 A node centrality indicator highlighting key transcription factors in the topology and functional analysis of the regulatory network inferred by the present application Figure 2 .

[0026] Figure 8 A plot of the correlation between TF centrality and GO enrichment breadth in the topology and functional analysis of the regulatory network inferred by the present application.

[0027] Figure 9 A functional enrichment plot of TF targets inferred by the present application. DETAILED DESCRIPTION

[0028] The specific embodiments of the present application are described below to enable those skilled in the art to understand the present application, but it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, any changes that are obvious within the spirit and scope of the present application as defined and determined by the appended claims are obvious, and all applications utilizing the concept of the present application are within the scope of protection.

[0029] The basic idea of ​​the present invention is to address the problem that single-cell sequencing unpaired data cannot be directly aligned, and to design a gene regulatory network reconstruction method that integrates generative latent representation learning and graph structure inference mechanisms. Taking into account that most of the scRNA-seq and scATAC-seq data commonly seen in practical applications are unpaired, and there are data structure heterogeneity and information redundancy problems between modalities, traditional integration methods that rely on paired data are difficult to apply directly. To this end, the present invention introduces a mutual nearest neighbor alignment mechanism to perform local structural correction on unpaired cells; at the same time, it constructs a trainable regulatory edge adjacency matrix, introduces the TF-target prior graph obtained by motif inference as a soft constraint, and realizes interpretable learning of the regulatory network. Finally, through an end-to-end deep neural network structure, the three processes of latent representation alignment, modality fusion and structure generation are combined to complete cross-modal feature integration and gene regulatory network reconstruction without the need for sample-level pairing.

[0030] like Figure 1 As shown, in one embodiment of the present invention, a method for inferring gene regulatory networks from unpaired single-cell multi-omics data includes the following steps: S1. Collecting unpaired single-cell sequencing data of human hematopoietic cell differentiation and preprocessing the data to generate preprocessed unpaired single-cell sequencing data. S2. Input the preprocessed single-cell sequencing unpaired data into a multi-layer variational autoencoder structure to generate high-level latent variables and reconstructed single-cell sequencing unpaired data; S3. A hierarchical GAN ​​is trained using a hierarchical adversarial alignment mechanism based on high-level latent variables and reconstructed single-cell sequencing unpaired data. The trained hierarchical GAN ​​is used to guide the alignment of high-level latent variables of different modalities in an abstract space. S4. Based on the high-level latent variables, a mutual nearest neighbor strategy is used to train the mutual nearest neighbor module. The trained mutual nearest neighbor module is used to constrain the consistency of the local structure between modalities in the embedding space, thereby enhancing the semantic alignment capability of high-level latent variables across modalities. S5. Use the gated fusion mechanism to fuse high-level latent variables of different modalities, construct a unified high-level latent variable, and complete the adaptive feature integration of multi-modal high-level latent variables; S6. Establish the initial adjacency matrix based on prior information ， Combine unified high-level latent variables to construct gene regulatory networks and complete gene regulatory network inference.

[0031] In S1, single-cell sequencing unpaired data include scRNA-seq data (single-cell RNA sequencing data) and scATAC-seq data (single-cell chromatin accessibility data); Preprocessing methods include row normalization, principal component analysis, and high-variance feature screening.

[0032] In the embodiment, the application collects scRNA-seq data and scATAC-seq data, and performs gene regulation network inference from multiple data dimensions, which helps better characterize complex gene regulation processes. The collected single-cell sequencing non-paired data is from Integrated single-cell analysis maps the continuous regulatory landscape of human hematopoietic differentiation.

[0033] In S2, the multi-layer variational autoencoder structure includes a first-level encoder, a second-level encoder and a zero-inflated decoder; S2 includes the following sub-steps: S21, inputting the count matrix of the single-cell sequencing non-paired data into the first-level encoder to generate a omics-specific latent variable; S22, inputting the omics-specific latent variable into the second-level encoder to generate a high-level latent variable; S23, inputting the omics-specific latent variable into the zero-inflated decoder to generate reconstructed single-cell sequencing non-paired data.

[0034] In the embodiment, in order to effectively integrate unpaired single-cell multi-omics data, the application introduces a multi-layer variational autoencoder (VAE) structure, combined with a generative adversarial learning strategy, which gradually extracts and aligns latent representations across omics layers, and maps from omics-specific latent space to shared latent space for downstream tasks.

[0035] In S21, the omics-specific latent variable is generated The expression is specifically: In the formula, k is the omics, RNA represents the RNA mode, and ATAC represents the ATAC mode, is a softplus activation function, is the mean learned by the first-level encoder, is the scale vector learned by the first-level encoder, is a Hadamard product, is the first Gaussian noise for differentiable sampling; In the formula, is the first-level encoder, is a count matrix of single-cell sequencing non-paired data, wherein, represents expression of cells on genes, represents openness of cells on chromatin regions, represents a multi-dimensional standard normal distribution, each dimension of which has a mean of 0, a variance of 1, and is independent of each other, is a regularization, is batch normalization, is a LeakyReLU activation function, is a linear layer;

[0036] In this embodiment, the Dropout rate is set to 0.2 to alleviate overfitting, and all layer orders are applied to achieve robust feature learning. In S22, the expression of the high-level latent variable is specifically: In the formula, is the mean learned by the second-level encoder, is the scale vector learned by the second-level encoder, is the second Gaussian noise for differentiable sampling; wherein the high-level latent variable captures abstract and generalized semantics across omics, enhancing cross-modal alignment; In the formula, is the second-level encoder, is a tanh activation function; In S23, the expression of the count matrix of the reconstructed single-cell sequencing non-paired data is specifically: In the formula, is a zero-inflated decoder, is a zero-inflated lognormal distribution model, is a linear transformation layer for predicting logits corresponding to zero-inflation probability, used to determine whether an element is structural 0 in a zero-inflation model (such as ZINB or ZILN), is a linear transformation layer for predicting the mean parameter of the log-normal main distribution in the ZILN distributionμ a linear transformation layer for determining the center position of the expression value, for predicting the scale parameter of each gene in the log-normal main distribution δ a linear layer to model the variability of the expression data, b auxiliary variables for estimating the zero inflation probability.

[0037] In this embodiment, in order to reconstruct The present application uses a decoder based on Zero-Inflated Log-Normal (ZILN) to jointly capture the count and dropout events through a Zero-Inflated Log-Normal distribution model, thereby improving the reconstruction quality.

[0038] In S3, the hierarchical GAN includes a first discriminator and a second discriminator, and S3 includes the following steps: S31, connect the high-level latent variable to generate the first modality label, and train the first discriminator through binary cross-entropy loss; S32, connect the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data to generate the second modality label, and train the second discriminator through binary cross-entropy loss.

[0039] In this embodiment, since the extracted high-level latent variable has different modality semantic information, the present application introduces a hierarchical GAN, which has two discriminators operating at different semantic depths. This hierarchical adversarial setting enhances the latent alignment and generation fidelity, so that the present application can perform high-quality modality conversion and gene regulation inference.

[0040] In S31, the first loss function of the first discriminator in S31 is trained through binary cross-entropy loss The expression of the first loss function is specifically: In the formula, is the modality original label of the i th sample, i is the number of samples, which is the number of cells in this embodiment, N is the first modality label of the i th sample, and the expression of the first modality label The expression of the first modality label is specifically: i In the formula, is the connection result of the high-level latent variable; In the formula, is a connection operation, ​​is the high-level latent variable of the cell in the RNA modality, It is the high-level latent variable of cells in ATAC mode; In S32, the second loss function of the second discriminator is trained by binary cross entropy loss The specific expression is: Where, For the i The original sequencing data expression level of each sample, For the i The second modal label of samples, calculate the second modal label The specific expression is: Where, is the connection result of the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data, .

[0041] In S4, the method for training the nearest neighbor modules is as follows: According to the high-level latent variables, the cell pairs that are the nearest neighbors of each other (Mutual Nearest Neighbors, MNN) are selected, the alignment loss is calculated based on all the cell pairs that are the nearest neighbors of each other, and the mutual nearest neighbor module is trained; among them, the alignment loss The specific expression is: Where, The first in RNA mode i cell, The first j Cells, MNN is a pair of cells that are each other’s nearest neighbors, express and are pairs of cells that are each other's nearest neighbors.

[0042] In this embodiment, in order to improve the consistency of local structure, a mutual nearest neighbor module based on the mutual nearest neighbor strategy is introduced. The mutual nearest neighbor module effectively enhances the cross-modal semantic alignment capability by constraining the consistency of local structures between modalities in the embedding space, and improves the robustness of the model to local perturbations between heterogeneous modalities.

[0043] In S5, a unified high-level latent variable is constructed The specific expression is: Where, For dynamic gating, High-level latent variables reflecting the specificity of each omics in the fusion process The relative importance of , and is calculated independently for each omics type to avoid mutual interference, The number of omics types, with a small constant added to the denominator This is to ensure numerical stability during training and avoid division by zero; Where, is the sigmoid activation function, which is used to map values ​​to scope.

[0044] In this embodiment, to further enhance the integration of multimodal high-level latent variables, the present invention introduces an adaptive gated fusion mechanism. This mechanism dynamically adjusts the contribution of each omics-specific latent representation during the fusion process, addressing the sparsity problem caused by direct concatenation and improving the performance of downstream tasks. The adaptive gated fusion mechanism regulates the contribution of each omics: noisy or redundant features are given a near-zero weight and suppressed, while informative features are emphasized, thereby improving the expressiveness and robustness of the fused representation. The final fused representation retains cross-omics consistency and can be interpreted through gating weights, revealing the contribution of each omics modality to cell state modeling.

[0045] In steps S1~S5, a unified cross-modal latent representation is effectively constructed While semantic alignment of unpaired data is achieved, translating these high-dimensional abstract features into biologically interpretable regulatory relationships remains a significant challenge. To address this issue, step S6 of the present invention introduces a unified end-to-end framework that integrates prior biological knowledge and data-driven learning for GRN inference from multi-omics embeddings. By combining variational autoencoders with graph-based structure learning, this framework leverages existing regulatory priors while simultaneously discovering new interactions from heterogeneous omics data, thereby improving the accuracy and interpretability of GRN inference.

[0046] S6 includes the following sub-steps: S61. Constructing an initial adjacency matrix based on scRNA-seq data A , ,in G represents the total number of genes, RRepresentation matrix; Given that previous TF-target gene networks usually cannot cover all experimentally observed genes, this paper adopts a conversion strategy to map TF-gene regulatory relationships into gene-gene prior connections. These priors are integrated through weighted fusion and random initialization: Where, is a binary matrix used to represent the regulatory links from prior TFs to genes, and T is the number of TFs; Where, is the adjacency matrix representing the prior regulatory network after mapping, is a randomly generated adjacency matrix, is a learnable coefficient used to balance the influence of prior knowledge and data-driven initialization; S62. Input the initial adjacency matrix, unified high-level latent variables, and count matrix of scRNA-seq data into a multilayer perceptron, train the multilayer perceptron using reconstruction loss, and construct a gene regulatory network. In this embodiment, the unified high-level latent variables are combined with the initial adjacency matrix in the decoder. A The overall training objective combines the reconstruction loss, KL divergence regularization and Sparsity constraint, reconstruction loss The specific expression is: Where, To sample latent variables from the encoder's output distribution for reconstruction evaluation, For the output of the decoder, use the fusion latent variables to generate the original data, is a hyperparameter used to adjust the weight of the KL divergence term, used to balance the reconstruction loss and regularization term, Fusion of latent variables for given input data X The approximate posterior distribution of is the fusion latent variable The prior distribution of To control the weight of the sparse regularization term hyperparameters to control the sparsity of model parameters, is the sparse regularization term weight, is the L1 norm.

[0047] where the first term evaluates the reconstruction error of gene expression; the second term enforces latent variables to follow standard Gaussian priors; the third term imposes regularization on the weight matrix W associated with A regularization, encouraging sparsity and aligning with prior regulatory knowledge.

[0048] During training, the adjacency matrix A is dynamically updated, preserving known regulatory links while allowing the model to discover new regulatory edges. This enables joint optimization under the guidance of both knowledge priors and multi-omic data.

[0049] To verify the effectiveness of the method of the present application, the following experimental data is provided in this embodiment: I. Comparative experiment Inference of gene regulatory networks (GRNs) from single-cell multi-omic data is crucial for decoding cellular regulation. However, a major challenge lies in the ubiquitous unpaired datasets, where different omic modalities (e.g., scRNA-seq and scATAC-seq) are measured in different cells. This mismatched structure hinders traditional GRN inference frameworks that rely on tightly aligned modalities. The present application addresses this challenge by aligning modality-specific latent spaces through variational autoencoders and adversarial training, and then inferring TF-target interactions through gated fusion and graph-guided regulatory priors integration.

[0050] The present application was benchmarked against seven representative GRN inference methods: scMTNI, DeepSEM, scGLUE, LINGER, SCENIC, scGPT, and UnpairReg. Evaluation was conducted on two widely used benchmarks: (1) the Buenrostro dataset, which is unpaired and represents hematopoietic lineage differentiation; (2) a paired PBMC dataset from 10x Genomics. Model performance was evaluated using AUC, AUPR, and F1-score metrics.

[0051] Table 1. Comparison of GRN inference performance (AUC, AUPR, F1-score) As shown in Table 1, on the unpaired dataset, the present application achieved the highest AUPR (0.688) and F1-score (0.773), outperforming the suboptimal method (scGLUE) by 5.9% AUPR and 4.2% F1. Compared with LINGER, which is designed specifically for paired data, the present application improved AUPR by 24.9% and F1-score by 21.9%, verifying the ability to achieve high-quality regulatory edge prediction under unpaired structure.

[0052] Notably, to enable LINGER to be evaluated on the unpaired Buenrostro dataset, we applied cell-level matching and down sampling, reducing the scRNA-seq sample size from 14,432 to 2,034 cells. This pre-processing likely compromised the complete regulatory landscape and led to LINGER’s underperformance in this setting (AUPR: 0.439, F1: 0.554). The necessity of this modification further highlights the importance of developing methods like the present invention that natively support unpaired data without relying on aggressive heuristic alignment or data reduction. Compared to UnpairReg, which also supports unpaired multi-omic data, the present invention consistently achieves superior performance, underscoring the advantage of generative models in revealing the latent regulatory relationships inherent to unpaired data.

[0053] On the paired PBMC dataset, the present invention maintains top performance, achieving the highest AUPR (0.672) and F1 score (0.770), while LINGER achieves the highest AUC (0.708). This indicates that the present invention effectively balances global ranking and local edge recovery, while robustly generalizing across paired and unpaired settings.

[0054] Furthermore, DeepSEM, scGPT, and SCENIC are designed for single-modality scRNA-seq data. Their relatively low performance on both benchmarks—especially in terms of AUPR and F1 score—highlight the limitations of RNA-only approaches in reconstructing transcriptional regulation. For example, on the PBMC dataset, the present invention’s AUPR is significantly higher compared to SCENIC, with nearly a doubling in F1 score, indicating improved performance in regulatory network inference. These gains underscore the advantage of leveraging multi-omic signals, including chromatin accessibility, to reveal context-specific TF-target relationships.

[0055] To ensure biologically grounded evaluation, we use ChIP-seq derived TF-target relationships from the Cistrome project as ground truth. This high-confidence resource aggregates experimentally validated transcription factor binding sites across different human cell types and tissues, providing a reliable benchmark for regulatory interaction evaluation. The present invention’s consistent advantage over this stringent ground truth further reinforces its practical utility in GRN reconstruction.

[0056] II. Leadership in single-cell multi-omic integration capability Accurate integration of single-cell multi-omic data is crucial for generating unified cell representations that support downstream tasks such as cell type classification and gene regulatory network (GRN) inference. The present invention addresses this challenge through a multi-level variational autoencoder (VAE) framework augmented with a hierarchical generative adversarial network (GAN) that facilitates robust alignment of unpaired modalities.

[0057] We evaluated the present invention on the human hematopoietic differentiation dataset from Buenrostro et al., which includes 14,432 scRNA-seq profiles spanning 7 cell types, and 2,034 scATAC-seq profiles spanning 10 cell types. Both modalities exhibit significant heterogeneity in sample size and cell type distribution. The scATAC-seq data covers approximately 150,000 chromatin accessible peaks, while the scRNA-seq data contains 12,558 genes, including 1,212 highly variable genes for downstream analysis.

[0058] To assess integration quality, we performed UMAP projection on the latent space generated by the present invention. As shown in Figure 2 the present invention effectively aligns scRNA-seq and scATAC-seq profiles while preserving biologically meaningful cell type structure, achieving coherent modality alignment while preserving biologically meaningful clusters.

[0059] In addition to integration, the present invention supports cross-modality translation, allowing for the inference of missing omic layers. In this task, we used the present invention to predict scRNA-seq expression profiles from scATAC-seq input and mapped the predicted profiles into the transcriptomic latent space. This translation process also involved aligning 10 ATAC-defined cell types with their transcriptomic counterparts.

[0060] Figure 3 The clustered results of the translation were visualized and confirmed that the present invention accurately reproduced RNA-like clustering structure from chromatin accessibility data. These results demonstrate the ability of the model to perform biologically coherent cross-modality generation under unpaired conditions. Figure 3 (a) is a UMAP plot of the original scRNA-seq profiles. Figure 3(b) UMAP plot of predicted RNA expression profiles from scATAC-seq input of the present application. The predicted embedding preserves clusters consistent with real RNA profiles, demonstrating the ability of the present application to infer transcriptome states from epigenetic signals.

[0061] III. Case Study To evaluate the ability of the present application to reveal gene regulatory mechanisms, we constructed a directed, weighted transcriptional regulatory network from its output and performed systematic analysis of network topology and functional annotation. We used the human hematopoietic differentiation single-cell multi-omics dataset from Buenrostro et al. as a case study, applying the present application and GLUE for regulatory inference. For each method, the top 50 regulatory edges ordered by EdgeWeight were extracted to construct subnetworks as shown in Figure 4 to Figure 5 To ensure the reliability of interactions, ChIP-seq data were integrated, and only experimentally supported regulatory edges were kept as high-confidence links.

[0062] The regulatory network predicted by GLUE exhibits a simple radial topology with SPI1 as the only central node. In contrast, the present application constructed a more complex network that not only identified SPI1 but also multiple highly interconnected key transcription factors (TFs), including FOS, JUN, JUNB, KLF6, SPI1, and YBX1, forming a multi-centered regulatory module. This suggests that the present application more effectively captures complex regulatory structures and potential co-regulatory mechanisms, significantly enhancing the interpretability and biological relevance of the network.

[0063] We further used classical centrality metrics, PageRank, Degree, and Betweenness, to characterize the network inferred by the present application as shown in Figure 6 and Figure 7 Key TFs such as FOS, JUN, YBX1, and SPI1 occupy central positions, indicating their regulatory importance. To assess their functional regulatory capacity, we integrated GO enrichment analysis and plotted the relationship between PageRank scores and the breadth of functional enrichment as shown in Figure 8JUN exhibited high centrality and extensive GO enrichment, representing a typical "high centrality - high function" TF. FOS and JUNB, as core members of the AP-1 family, also located in the central position and enriched in multiple biological processes, implying their synergistic roles in hematopoietic differentiation. Although YBX1 had relatively low centrality, its targets were significantly enriched in stress response and metabolic pathways, indicating functional diversity and context-dependent regulation. SPI1 showed high PageRank but limited GO enrichment, implying its key hub role in the structure but narrow functional range. KLF6 exhibited moderate centrality and enrichment, suggesting its auxiliary role in specific regulatory modules.

[0064] To further elucidate functional specialization, GO enrichment of target genes of high centrality TFs was visualized by heatmaps. Figure 9 showed unique but complementary enrichment patterns: SPI1 and KLF6 were enriched in neutrophil-mediated immunity, implying immune regulation; JUNB was associated with negative regulation of multiple cellular processes, indicating the modulation of hematopoietic microenvironment, Figure 9 In the middle, a heatmap of -log10 (adjusted p-value) of enriched GO Biological Process terms associated with the top six TFs' target genes. The strong enrichment on immune regulation, apoptosis, and cell cycle pathways demonstrated the biological relevance of the inferred network. Multiple TFs (FOS, JUN, JUNB, KLF6, SPI1) were co-enriched in cell cycle regulation, reflecting the coordinated control of proliferation and differentiation. Additional enrichment in programmed cell death and differentiation further highlighted the diverse regulatory roles in hematopoiesis.

[0065] Finally, based on the top three GO terms of each TF, we reconstructed TF-Target-Function subgraphs, systematically delineating the functional regulatory architecture, highlighting candidate genes and modules for mechanistic exploration.

[0066] By capturing diverse TF roles, the present invention revealed biologically meaningful, context-specific regulatory modules. This enabled a deeper understanding of gene regulation in hematopoiesis and related processes.

[0067] In the description of the application, it needs to be understood that the terms "center", "thickness", "upper", "lower", "horizontal", "top", "bottom", "inner", "outer", "radial" and the like indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the application and simplifying the description, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as a limitation on the application. In addition, the terms "first", "second", "third" are only for the purpose of description, and cannot be understood as indicating or implying relative importance or implying the number of technical features indicated. Therefore, the features defined by "first", "second", "third" can explicitly or implicitly include one or more of the features.

Claims

1. A method for inferring gene regulatory networks from unpaired single-cell multi-omics data, characterized by: The following steps are involved: S1. Collecting unpaired single-cell sequencing data of human hematopoietic cell differentiation and preprocessing the data to generate preprocessed unpaired single-cell sequencing data. S2. Input the preprocessed single-cell sequencing unpaired data into a multi-layer variational autoencoder structure to generate high-level latent variables and reconstructed single-cell sequencing unpaired data; S3. A hierarchical GAN ​​is trained using a hierarchical adversarial alignment mechanism based on high-level latent variables and reconstructed single-cell sequencing unpaired data. The trained hierarchical GAN ​​is used to guide the alignment of high-level latent variables of different modalities in an abstract space. S4. Based on the high-level latent variables, a mutual nearest neighbor strategy is adopted to train the mutual nearest neighbor module. The trained mutual nearest neighbor module is used to constrain the consistency of the local structure between modalities in the embedding space. S5. Use the gated fusion mechanism to fuse high-level latent variables of different modalities, construct a unified high-level latent variable, and complete the adaptive feature integration of multi-modal high-level latent variables; S6. Establish the initial adjacency matrix based on prior information ， Combine unified high-level latent variables to construct gene regulatory networks and complete gene regulatory network inference.

2. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 1, characterized in that: In S1, single-cell sequencing unpaired data include scRNA-seq data and scATAC-seq data; Preprocessing methods include row normalization, principal component analysis, and high-variance feature screening.

3. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 1, characterized in that: In S2, the multi-layer variational autoencoder structure includes a first-stage encoder, a second-stage encoder, and a zero-expansion decoder; S2 includes the following sub-steps: S21. Input the count matrix of single-cell sequencing unpaired data into the first-level encoder to generate omics-specific latent variables; S22, input the omics-specific latent variables into the second-level encoder to generate high-level latent variables; S23. Input the omics-specific latent variables into the zero-inflated decoder to generate reconstructed single-cell sequencing unpaired data.

4. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 3, characterized in that: In S21, generate omics-specific latent variables The specific expression is: Where, k For omics, , RNA represents RNA modality, ATAC represents ATAC modality, is the softplus activation function, is the mean learned by the first-level encoder, is the scale vector learned by the first-level encoder, For Hadamard, is the first Gaussian noise used for differentiable sampling; Where, is the first-level encoder, is the count matrix of single-cell sequencing unpaired data, Represents a multidimensional standard normal distribution, where each dimension has a mean of 0 and a variance of 1, and the dimensions are independent of each other. is regularization, is batch normalization, is the LeakyReLU activation function, is a linear layer; In S22, high-level latent variables are generated The specific expression is: Where, is the mean learned by the second-level encoder, is the scale vector learned by the second-level encoder, is the second Gaussian noise used for differentiable sampling; Where, is the second-stage encoder, is the tanh activation function; In S23, the count matrix of reconstructed single-cell sequencing unpaired data is generated The specific expression is: Where, is a zero-expansion decoder, is the zero-inflated lognormal distribution model, It is a linear transformation layer used to predict the logits corresponding to the zero-inflated probability, and is used to determine whether an element is structurally 0 in the zero-inflated model. is the mean parameter used to predict the log-normal main distribution in the ZILN distribution μ The linear transformation layer is used to determine the center position of the expression value. is the scale parameter used to predict each gene in the log-normal main distribution δ Linear layers, thereby modeling the variability of expression data, b is an auxiliary variable used to estimate the zero-inflated probability.

5. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 4, characterized in that: In S3, the hierarchical GAN ​​includes a first discriminator and a second discriminator. S3 includes the following sub-steps: S31, concatenate high-level latent variables to generate the first modal label, and train the first discriminator through binary cross entropy loss; S32. Connect the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data to generate a second modality label, and train the second discriminator through binary cross entropy loss.

6. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 5, characterized in that: In S31, the first loss function of the first discriminator is trained by binary cross entropy loss The specific expression is: Where, For the i The modal original labels of samples, N is the sample size, For the i The first modal label of samples, calculate the first modal label The specific expression is: Where, is the result of the connection of high-level latent variables; Where, For the connection operation, is the high-level latent variable of the cell in the RNA modality, It is the high-level latent variable of cells in ATAC mode; In S32, the second loss function of the second discriminator is trained by binary cross entropy loss The specific expression is: Where, For the i The original sequencing data expression level of each sample, For the i The second modal label of samples, calculate the second modal label The specific expression is: Where, is the connection result of the reconstructed single-cell sequencing unpaired data and the single-cell sequencing unpaired data, .

7. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 6, characterized in that: In S4, the method for training the nearest neighbor modules is as follows: According to the high-level latent variables, the cell pairs that are each other's nearest neighbors are selected, and the alignment loss is calculated based on all the cell pairs that are each other's nearest neighbors, and the nearest neighbor module is trained; among them, the alignment loss The specific expression is: Where, The first in RNA mode i cell, The first j Cells, MNN is a pair of cells that are each other’s nearest neighbors, express and are pairs of cells that are each other's nearest neighbors.

8. The gene regulatory network inference method based on unpaired single-cell multi-omics data according to claim 7, characterized in that: In S5, a unified high-level latent variable is constructed The specific expression is: Where, For dynamic gating, is the number of omics types; Where, is the sigmoid activation function.

9. The method for inferring gene regulatory networks from unpaired single-cell multi-omics data according to claim 8, characterized in that: S6 includes the following sub-steps: S61. Constructing an initial adjacency matrix based on scRNA-seq data A ; Where, is the adjacency matrix representing the prior regulatory network after mapping, is a randomly generated adjacency matrix, is a learnable coefficient used to balance the influence of prior knowledge and data-driven initialization; Where, is a binary matrix used to represent the regulatory links from prior TFs to genes, and T is the number of TFs; S62. Input the initial adjacency matrix, unified high-level latent variables, and count matrix of scRNA-seq data into a multilayer perceptron, train the multilayer perceptron using reconstruction loss, and construct a gene regulatory network. Among them, the reconstruction loss The specific expression is: Where, To sample latent variables from the encoder's output distribution for reconstruction evaluation, For the output of the decoder, use the fusion latent variables to generate the original data, is a hyperparameter used to adjust the weight of the KL divergence term, used to balance the reconstruction loss and regularization term, Fusion of latent variables for given input data X The approximate posterior distribution of is the fusion latent variable The prior distribution of To control the weight of the sparse regularization term hyperparameters to control the sparsity of model parameters, is the sparse regularization term weight, is the L1 norm.

Images