A method for constructing a gene regulatory network and related devices

Abstract

The application discloses a gene regulation network construction method and related device, and relates to the technical field of gene regulation network construction. The method comprises the following steps: downloading position data of candidate regulation elements from a SCREEN database, combining a genome annotation file and scATAC-seq data to accurately position promoters and enhancers, introducing a variational autoencoder to accurately determine enhancer-promoter pairs that exist potential interaction, associating enhancers with target genes, processing a co-expression gene network based on potential binding sites of each transcription factor of a to-be-detected object on an enhancer, enhancer-promoter pairs that exist potential interaction, and target genes corresponding to each promoter, and obtaining a gene regulation network. The application can solve the problems of failing to accurately position promoters and enhancers, being difficult to associate enhancers with target genes, being difficult to distinguish direct and indirect relationships, and having too many false positive regulation relationships.

Description

Technical Field

[0001] This application relates to the field of gene regulatory network construction technology, and in particular to a method and related apparatus for constructing an enhancer-driven gene regulatory network based on EnhancerNet. Background Technology

[0002] Gene Regulatory Networks (GRNs) record the relationships between transcription factors (TFs) and their regulated target genes, revealing the mechanisms of gene expression regulation within cells. This helps to understand the molecular basis of cell and tissue behavior, and a deeper understanding of GRNs is extremely important for revealing the molecular mechanisms controlling cell behavior.

[0003] Currently, gene regulatory networks are often inferred using experimental methods such as chromatin immunoprecipitation sequencing and interference experiments with specific transcription factors combined with transcriptome sequencing. However, large-scale inference of gene regulatory networks remains challenging. The advent of single-cell RNA sequencing (scRNA-seq) has provided unprecedented dynamic resolution of the transcriptome, offering new insights for gene regulatory network inference and promoting the development of algorithms based on single-cell sequencing data. Although a series of algorithms for inferring gene regulatory networks based on single-cell sequencing data have been developed, such as TENET, SCENIC, and Dictys, problems such as inaccurate localization of promoters and enhancers, difficulty in associating enhancers with target genes, difficulty in distinguishing direct and indirect relationships, and excessive false positives in regulatory relationships remain serious challenges for these algorithms. Summary of the Invention

[0004] The purpose of this application is to provide a method and related apparatus for constructing a gene regulatory network, which can solve problems such as failure to accurately locate promoters and enhancers, difficulty in associating enhancers with target genes, difficulty in distinguishing direct and indirect relationships, and excessive false positive regulatory relationships.

[0005] To achieve the above objectives, this application provides the following solution:

[0006] In a first aspect, this application provides a method for constructing a gene regulatory network, the method comprising:

[0007] Acquire scATAC-seq data, scRNA-seq data, candidate regulatory element location data, and genome annotation files of the target organism; the target organism is a tissue of the target organism; the candidate regulatory elements include candidate promoters and candidate enhancers, and the candidate regulatory element location data is downloaded from the SCREEN database, including the chromosome number where the candidate regulatory element is located, the start site of the candidate regulatory element on its chromosome, and the termination site of the candidate regulatory element on its chromosome;

[0008] Based on the location data of candidate regulatory elements and the genome annotation file, the candidate promoters are annotated to determine the target genes corresponding to the candidate promoters, thereby obtaining the annotation data of the candidate regulatory elements; the annotation data of the candidate regulatory elements includes the location data of each candidate promoter and the corresponding target gene, as well as the location data of each candidate enhancer.

[0009] Based on scATAC-seq data and annotation data of candidate regulatory elements, regulatory elements are identified, and the readings of the regulatory elements in each cell of the target object are calculated to obtain the annotation data of the regulatory elements. The regulatory elements are open candidate regulatory elements, which include promoters and enhancers. The annotation data of the regulatory elements include the location data of each promoter, the corresponding target gene, and the readings in each cell of the target object, as well as the location data of each enhancer and the readings in each cell of the target object.

[0010] Using the readout matrix as input, a latent variable matrix is ​​calculated using a variational autoencoder, and enhancer-promoter pairs with potential interactions are determined based on the latent variable matrix; the readout matrix includes the readout of each regulatory element in each cell of the target object; the latent variable matrix includes the embedded representation of each regulatory element;

[0011] Based on the motif database, the position data of each enhancer and the readings in each cell of the target object are scanned to determine the potential binding sites of each transcription factor of the target object on the enhancer.

[0012] Based on the potential binding sites of each transcription factor on the enhancer of the target object, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated, and the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain the transcription factor-target gene regulatory network.

[0013] Mutual information between genes was calculated based on scRNA-seq data, and a co-expression gene network was constructed based on the mutual information between genes; the genes included coding genes and target genes of transcription factors.

[0014] The intersection of the transcription factor-target gene regulatory network and the co-expressed gene network is calculated, and the intersection is expanded based on the potential binding sites of each transcription factor of the target object on the enhancer, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, to obtain the gene regulatory network of transcription factor-enhancer-promoter-target gene.

[0015] Secondly, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above-described method for constructing a gene regulatory network.

[0016] Thirdly, this application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for constructing a gene regulatory network.

[0017] Fourthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the above-described method for constructing a gene regulatory network.

[0018] According to the specific embodiments provided in this application, this application has the following technical effects:

[0019] This application provides a method and related apparatus for constructing a gene regulatory network. Positional data of candidate regulatory elements are downloaded from the SCREEN database. Further combining this data with genome annotation files and scATAC-seq data enables precise localization of promoters and enhancers, addressing the problem of inaccurate promoter and enhancer localization. A variational autoencoder is introduced to accurately identify potentially interacting enhancer-promoter pairs. Since the promoter is already associated with the target gene, the enhancer can be further associated with the target gene, solving the problem of difficulty in associating enhancers and target genes. Based on the potential binding sites of each transcription factor on the enhancer of the target gene, the potentially interacting enhancer-promoter pairs, and the target gene corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated. This network is further simplified to obtain the transcription factor-target gene regulatory network. The intersection of the transcription factor-target gene regulatory network and the co-expressed gene network is calculated, further expanding to obtain the transcription factor-enhancer-promoter-target gene gene regulatory network, solving the problems of difficulty in distinguishing direct and indirect relationships and excessive false positive regulatory relationships. This application can solve problems such as failure to accurately locate promoters and enhancers, difficulty in associating enhancers with target genes, difficulty in distinguishing direct and indirect relationships, and excessive false positive regulatory relationships. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a flowchart illustrating a method for constructing a gene regulatory network as provided in Embodiment 1 of this application.

[0022] Figure 2 This is a schematic diagram of the enhancer-driven gene regulatory network with PBX1 as the core in the artificial blood system provided in Embodiment 1 of this application.

[0023] Figure 3 This is a schematic diagram of the structure of a computer device provided in Embodiment 2 of this application. Detailed Implementation

[0024] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0025] Example 1

[0026] This embodiment provides a method for constructing a gene regulatory network, which includes the following steps.

[0027] Step S1: Obtain scATAC-seq data, scRNA-seq data, candidate regulatory element location data, and genome annotation file of the object to be tested; the object to be tested is the tissue of the organism to be tested; the candidate regulatory elements include candidate promoters and candidate enhancers, and the location data of the candidate regulatory elements is downloaded from the SCREEN database. The location data of the candidate regulatory elements includes the chromosome number where the candidate regulatory element is located, the start site of the candidate regulatory element on its chromosome, and the termination site of the candidate regulatory element on its chromosome.

[0028] Step S2: Based on the location data of the candidate regulatory elements and the genome annotation file, the candidate promoters are annotated to determine the target genes corresponding to the candidate promoters, thereby obtaining the annotation data of the candidate regulatory elements. The annotation data of the candidate regulatory elements includes the location data of each candidate promoter and its corresponding target gene, as well as the location data of each candidate enhancer.

[0029] Step S3: Based on scATAC-seq data and annotation data of candidate regulatory elements, determine the regulatory elements and calculate the readings of the regulatory elements in each cell of the target object to obtain the annotation data of the regulatory elements; the regulatory elements are open candidate regulatory elements, and the regulatory elements include promoters and enhancers. The annotation data of the regulatory elements include the location data of each promoter, the corresponding target gene and the readings in each cell of the target object, as well as the location data of each enhancer and the readings in each cell of the target object.

[0030] Step S4: Using the readout matrix as input, a latent variable matrix is ​​calculated using a variational autoencoder, and enhancer-promoter pairs with potential interactions are determined based on the latent variable matrix; the readout matrix includes the readout of each regulatory element in each cell of the object to be tested; the latent variable matrix includes the embedded representation of each regulatory element.

[0031] Step S5: Based on the motif database, perform motif scanning on the position data of each enhancer and the readings in each cell of the target object to determine the potential binding sites of each transcription factor of the target object on the enhancer.

[0032] Step S6: Based on the potential binding sites of each transcription factor of the target object on the enhancer, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated, and the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain the transcription factor-target gene regulatory network.

[0033] Step S7: Calculate the mutual information between genes based on scRNA-seq data, and construct a co-expression gene network based on the mutual information between genes; the genes include coding genes and target genes of transcription factors.

[0034] Step S8: Calculate the intersection of the transcription factor-target gene regulatory network and the co-expressed gene network, and expand the intersection based on the potential binding sites of each transcription factor of the target object on the enhancer, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, to obtain the gene regulatory network of transcription factor-enhancer-promoter-target gene.

[0035] By implementing steps S1 to S8 above, this embodiment downloads the location data of candidate regulatory elements from the SCREEN database. Further combining this with genome annotation files and scATAC-seq data enables precise localization of promoters and enhancers, resolving the problem of inaccurate promoter and enhancer localization. Introducing a variational autoencoder accurately identifies potentially interacting enhancer-promoter pairs. Since the promoter is already associated with the target gene, the enhancer can be further associated with the target gene, resolving the difficulty in associating enhancers and target genes. Based on the potential binding sites of each transcription factor on the enhancer of the target object, the potentially interacting enhancer-promoter pairs, and the target gene corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated. This network is further simplified to obtain the transcription factor-target gene regulatory network. The intersection of the transcription factor-target gene regulatory network and the co-expressed gene network is calculated, further expanding to obtain the transcription factor-enhancer-promoter-target gene gene regulatory network, resolving the problems of difficulty in distinguishing direct and indirect relationships and excessive false-positive regulatory relationships.

[0036] To address the problems of existing algorithms used for inferring gene regulatory networks, such as failure to accurately locate promoters and enhancers, difficulty in associating enhancers with target genes, difficulty in distinguishing direct and indirect relationships, and excessive false positive regulatory relationships, this embodiment develops the algorithm model EnhancerNet (i.e., the gene regulatory network construction method used in this embodiment). It is used to infer large-scale gene regulatory networks of transcription factors-enhancer-promoter-target genes based on scATAC-seq data (Single-cell Assay for Transposase-Accessible Chromatinusing sequencing) and scRNA-seq data.

[0037] Below, we will briefly introduce the solutions to the four problems addressed in this embodiment:

[0038] (1) Solutions for failure to accurately locate promoters and enhancers

[0039] Promoters and enhancers are the two most important cis-regulatory elements of gene expression. Promoters are usually located near the transcription start site and initiate gene transcription by binding to RNA polymerase and transcription factors. Enhancers are long-range regulatory elements that can function upstream, downstream, and even within introns of genes, enhancing the expression of specific genes by altering chromatin conformation or forming DNA loops. The challenge in accurately locating cis-regulatory elements lies in the complexity and diversity of these two key regulatory regions: typical promoter regions contain some conserved sequence elements, such as TATA boxes, CAAT boxes, and GC boxes, which facilitate the recognition and binding of transcription factors. However, promoter sequences can vary significantly among different genes, increasing the difficulty of accurate localization. Even within the same species, the length, sequence composition, and position of functional elements of promoter sequences can differ. In particular, some gene promoters may lack clearly conserved sequence elements, making accurate identification difficult through sequence alignment or conventional annotation methods. Unlike promoters, enhancers do not have strict location restrictions or specific sequence features, which makes their identification more complex. Enhancer regions are usually associated with specific chromatin markers such as H3K27ac and H3K4me1. The presence of these chromatin markers helps in the identification of enhancers. However, enhancers have a high degree of sequence diversity, and enhancers of different genes may have significant differences in sequence composition and regulatory characteristics, which further increases the difficulty of their precise location.

[0040] The SCREEN database within the ENCODE (Encyclopedia of DNA Elements) project provides accurate location information for promoters and enhancers. It integrates and analyzes a large amount of multi-omics data, defining candidate promoters as DNA sequences with high DNase-seq and H3K4me3 signals, centered within 200 bp of the transcription start site. Candidate enhancers are defined as DNA sequences with high DNase-seq and H3K27ac signals that do not overlap with the regions of candidate promoters. Currently, no algorithmic models use this database. This embodiment introduces the SCREEN database, extracting the location data of candidate regulatory elements (i.e., candidate promoters and candidate enhancers) from the SCREEN database. Based on genome annotation files, candidate regulatory elements near the transcription start site are identified as candidate promoters, and other candidate regulatory elements are identified as candidate enhancers. Combined with scATAC-seq data, open candidate regulatory elements are obtained. These open candidate regulatory elements are used as regulatory elements, including promoters and enhancers, thus solving the problem that current related algorithmic models fail to accurately locate promoters and enhancers.

[0041] (2) Solutions to the difficulty in associating enhancers and target genes

[0042] Currently, linking enhancers to their target genes remains a significant challenge. While Hi-C technology can capture the three-dimensional structure of chromatin and reveal interactions between long-range regions in the genome, it performs poorly in the precise localization of long-range interactions. Furthermore, Hi-C technology struggles to achieve single-cell resolution, resulting in mixed data from different cell types, making it difficult to identify specific chromatin structures and gene regulatory patterns in different cell types. This embodiment uses algorithmic modeling on scATAC-seq data to decipher long-range interactions between enhancers and promoters at the single-cell level. First, the original high-dimensional sparse data determined based on scATAC-seq data is input into the encoder of a variational autoencoder (VAE). The encoder compresses the high-dimensional sparse data features of each regulatory element (also known as a chromatin region) into a low-dimensional latent space representation. By compressing the high-dimensional sparse data features into a low-dimensional latent space representation, the main features of the input data are preserved. Simultaneously, regularization constraints ensure that the generated latent variable matrix conforms to a predetermined latent distribution. Then, the decoder of the variational autoencoder reconstructs the chromatin accessibility data based on the latent variable matrix, ensuring that the extracted latent variable matrix can effectively reconstruct the original high-dimensional sparse data structure. Each row of the output latent variable matrix represents the embedding representation of a regulatory element in the latent space, and each column represents one dimension of the latent space. Through the above encoder and decoder processing, the variational autoencoder can capture the intrinsic structure of the chromatin accessibility data while addressing the sparsity and noise issues of the data. The extracted latent variable matrix is ​​then used to calculate the Wasserstein distance between regulatory elements, resulting in a distance matrix. The Wasserstein distance is sensitive to differences in probability distributions and can better reflect the differences in the distribution of different regulatory elements in the complex latent space. Unlike Euclidean distance, the Wasserstein distance is more suitable for handling asymmetric probability distributions, which are common in the latent space of the variational autoencoder. This distance matrix is ​​then further converted into an affinity matrix using a Gaussian kernel function to represent the similarity between regulatory elements in the latent space. Then, although enhancers can act on relatively distant genes, close-range interactions are more common and efficient. Therefore, a distance penalty term is added to the calculated affinity matrix, which relatively enhances the affinity between physically close regulatory elements, while imposing an additional penalty on distant regulatory elements. The penalized affinity (i.e., adjusted affinity) is used to measure the potential interaction between regulatory elements. Finally, by setting an affinity threshold, enhancer-promoter pairs on the same DNA strand with adjusted affinity greater than the threshold are selected to indicate the existence of potential long-range interactions. Based on the association between the promoter and the target gene, enhancers and target genes are further associated.

[0043] (3) Solutions to the difficulty in distinguishing direct and indirect relationships and the excessive number of false-positive regulatory relationships

[0044] Currently, most algorithms for inferring gene regulatory networks based on single-cell sequencing data focus on measuring co-expression relationships between genes. However, these algorithms lack information on the interactions between transcription factors and target genes, resulting in inherent drawbacks such as difficulty in distinguishing direct and indirect relationships and an excessive number of false positives. The EnhancerNet developed in this embodiment can construct a co-expression gene network based on scRNA-seq data. It measures co-expression relationships by calculating partial correlation coefficients (i.e., mutual information) between pairs of genes in large batches. Compared to other methods for calculating co-expression relationships, this method has the advantage of assessing the association between the expression levels of two genes while eliminating the influence of other genes. However, the co-expression gene network generated at this time still contains false positive regulatory relationships. In this embodiment, the gimmeotifs algorithm is further used to identify the binding site of each transcription factor in the precisely located enhancer based on the transcription factor binding motifs in the CisBPv.2 database, and a transcription factor-target gene regulatory network is obtained. The intersection of this transcription factor-target gene regulatory network and the co-expression gene network can be used to obtain the gene regulatory network of the biological process under study. Furthermore, the target genes directly regulated by the transcription factor can be directly identified by the identification of the transcription factor binding motif, thereby avoiding the generation of too many indirect regulatory relationships and removing most of the false positive regulatory relationships.

[0045] The following, combined with Figure 1 The algorithm model of this embodiment is described in detail, mainly condensed into six steps: precise localization of open promoters and enhancers, use of variational autoencoders for feature extraction, calculation of the similarity of open patterns of regulatory elements, inference of co-expressed gene networks, identification of transcription factor binding sites, and construction of enhancer-driven gene regulatory networks.

[0046] (1) Precise location of open promoters and enhancers, corresponding to steps S1-S3.

[0047] For the biological processes being studied, i.e., the object of detection is the tissue of the organism to be tested, specifically tissues from human and mouse species, or tissues from other species, single-cell transposase-accessible chromatin sequencing is required to obtain scATAC-seq data, or relevant scATAC-seq data can be downloaded from public databases. The raw scATAC-seq data is preprocessed through the following steps:

[0048] 1) Use the `fastq-dump` command in the `sratoolkit` to convert the SRA format data (also known as the format file) to FASTQ format. If the FASTQ format data contains sequencing data from multiple cells, a custom script is needed to extract the FASTQ format data of individual cells based on the barcode information.

[0049] 2) Use Trimmomatic to remove headers from fastq format data.

[0050] 3) Use bowtie2 to align the fastq format data to the reference genome corresponding to the object to be tested, and obtain the sam format data.

[0051] 4) Use samtools to convert SAM format data to BAM format data.

[0052] The scATAC-seq data in this embodiment is in BAM format.

[0053] Single-cell sequencing was performed on the subjects to be tested to obtain scRNA-seq data.

[0054] Download the genome annotation file in GTF format from the GENCODE database.

[0055] Location data for two types of candidate regulatory elements—candidate promoters and candidate enhancers—were downloaded from the publicly available ENCODE project's SCREEN database. The location data for candidate regulatory elements included the chromosome number on which the candidate regulatory element resides, the start site of the candidate regulatory element on its chromosome, and the termination site of the candidate regulatory element on its chromosome. In this SCREEN database, candidate promoters are described as promoter-like signatures (PLS), and candidate enhancers are described as enhancer-like signatures (ELS).

[0056] Based on this, scATAC-seq data, scRNA-seq data, candidate regulatory element location data, and genome annotation files of the object to be tested can be obtained.

[0057] This embodiment further annotates candidate promoters based on the location data of candidate regulatory elements and genome annotation files, determines the target genes corresponding to the candidate promoters, and obtains annotation data of candidate regulatory elements. The annotation data of candidate regulatory elements includes the location data of each candidate promoter and the corresponding target gene, as well as the location data of each candidate enhancer.

[0058] Specifically, the `annotatePeak` function in the R package `ChIPseeker` can directly associate PLS with its target gene. In the `annotatePeak` function's parameter settings, `tssRegion = c(-1500, 500)` indicates that the region 1500 bp upstream to 500 bp downstream of the gene transcription start site is designated as a candidate promoter region. In other words, PLS falling within this candidate promoter region are directly assigned as candidate promoters for that gene. Further processing outputs a BED format file containing candidate regulatory elements (candidate promoters and candidate enhancers). This BED file contains five columns: the chromosome number of the candidate regulatory element, the start site of the candidate regulatory element on its chromosome, the termination site of the candidate regulatory element on its chromosome, the type of the candidate regulatory element (PLS or ELS), and the corresponding target gene (only the target gene corresponding to PLS).

[0059] In this embodiment, candidate promoters are annotated based on the location data of candidate regulatory elements and the genome annotation file to determine the target genes corresponding to the candidate promoters. Specifically, this includes using the location data of candidate regulatory elements and the genome annotation file as input, and using the annotatePeak function to annotate the candidate promoters to determine the target genes corresponding to the candidate promoters.

[0060] This embodiment further identifies regulatory elements based on scATAC-seq data and annotation data of candidate regulatory elements, and calculates the readings of the regulatory elements in each cell of the subject to be tested to obtain the annotation data of the regulatory elements. The regulatory elements are open candidate regulatory elements, which include promoters and enhancers. The annotation data of the regulatory elements includes the location data of each promoter, the corresponding target gene and the readings in each cell of the subject to be tested, as well as the location data of each enhancer and the readings in each cell of the subject to be tested.

[0061] Specifically, the scATAC-seq data in BAM format and the BED format containing candidate regulatory elements are used as input. The BED tools are used to detect open candidate regulatory elements in this biological process and calculate the read counts of open candidate regulatory elements in a single cell. This generates multiple BED format files containing regulatory elements, each representing a cell. Each BED format file contains 7 columns of data: the regulatory element's ID (composed of a string of chromosome ID, start site, and end site), the chromosome ID of the regulatory element, the start site of the regulatory element on the chromosome, the end site of the regulatory element on the chromosome, the type of regulatory element (PLS or ELS), the corresponding target gene (only the target gene corresponding to PLS), and the read count of the regulatory element in the cell. After merging these BED format files into a single BED format file, a column of cell numbering information is added. The resulting BED format file, named SingleCell_Merged.bed, can be converted into matrix form using Pandas' pivot_table() function, yielding a readout matrix M. Here, the regulatory element number serves as the row index, the cell number as the column index, and the readout as the matrix value. Therefore, this step effectively obtains the precise location information of promoters and enhancers of a specific biological process on the chromosome, as well as the relationship between the promoter and its target gene.

[0062] In this embodiment, based on scATAC-seq data and annotation data of candidate regulatory elements, the regulatory elements are determined, and the readings of the regulatory elements in each cell of the target object are calculated. Specifically, this includes: using scATAC-seq data and annotation data of candidate regulatory elements as input, the bedtools tool is used to determine the regulatory elements, and the readings of the regulatory elements in each cell of the target object are calculated.

[0063] (2) Identify transcription factor binding sites, corresponding to step S5.

[0064] The identified enhancers are individually extracted and output as BED format files. The `run_cistarget` function of the `pycisTarget` package in Python is used to perform motif scanning on the enhancer BED files, thereby searching for transcription factor binding motifs in a motif database. The reference genome version for humans is hg38, and for mice it is mm10. The results retain four targeting relationships: "Direct_annot" (direct evidence that the motif is associated with a specific transcription factor TF), "Motif_similarity_annot" (based on tom-tom motif similarity), "Orthology_annot" (based on orthogonality of transcription factor TFs directly associated with the motif), or "Motif_similarity_and_Orthology_annot" (based on motifs with similarity thresholds and orthologous annotations with homology thresholds). After scanning and processing, a BED format file containing the identification results is obtained. This BED format file describes the potential binding sites of each TF on the input enhancer. Therefore, this step effectively obtains information on the binding of transcription factors to enhancers in the studied biological process.

[0065] In this embodiment, based on a motif database, a motif scan is performed on the position data of each enhancer and the readings in each cell of the target organism to determine the potential binding sites of each transcription factor of the target organism on the enhancers. Specifically, the motif scan is performed using the run-cistarget function based on the motif database to determine the potential binding sites of each transcription factor of the target organism on the enhancers.

[0066] (3) The variational autoencoder is used for feature extraction, corresponding to step S4.

[0067] The readout matrix M is further subjected to dimensionality reduction using a variational autoencoder to capture nonlinear relationships in the data and extract features, thereby obtaining a more informative representation than the original features. The model construction and training of the variational autoencoder are accomplished using PyTorch. The readout matrix M first undergoes preprocessing, which involves performing a logarithmic transformation on the matrix elements, followed by standardization of the transformed data to achieve numerical stability and improve model performance, as follows:

[0068] M' = log(M+1);

[0069]

[0070] Where M' is the logarithmically transformed matrix; M” is the preprocessed matrix; μ is the mean of the logarithmically transformed matrix M'; and σ is the standard deviation of the logarithmically transformed matrix M'.

[0071] Next, we construct the encoder, which is responsible for compressing the input data into the latent space. Essentially, the structure of the encoder can be represented by the following equation:

[0072] h = ReLU(W1M″+b1);

[0073] μ(M″) = W2h + b2;

[0074] logσ 2 (M”)=W3h+b3;

[0075] Where h is the hidden layer output obtained through the ReLU activation function; W1, W2, W3 are the encoder weight matrices; b1, b2, b3 are the encoder bias terms; μ(M”) is the mean vector of the latent space; logσ 2 (M”) represents the variable log-variance of the latent space. The weight matrix and bias terms are randomly generated during initialization and dynamically adjusted during training using an optimization algorithm to minimize the loss function.

[0076] The latent variable matrix Z is sampled from the latent distribution using the reparameterization technique, as shown in the following formula:

[0077] Z = μ(M”) + ∈·σ(M”);

[0078] Where ∈ represents the standard normal distribution. Noise in the mid-sample; denoted as the standard deviation of the potential space.

[0079] Next, a decoder is constructed. The goal of the decoder is to map the latent variable matrix Z back to the data space to reconstruct the input matrix M". The output of the decoder can be expressed by the following formula:

[0080]

[0081] in, To reconstruct the matrix; W d and b d These represent the weights and biases of the decoder, respectively; the sigmoid function is used to constrain the output value between 0 and 1. d and b d It is dynamically learned during the training process, with the initial value being randomly distributed, and finally updated to the value that best adapts to the data through an optimization algorithm.

[0082] The role of the decoder is not only to generate a reconstruction matrix similar to the input data. It is also an important mechanism for effectively constraining the latent space during training. Through comparison... The difference between "M" and "M" is that the decoder provides optimization signals to the model, thereby guiding the encoder to compress high-dimensional data into a low-dimensional potential space while minimizing information loss, achieving the dual goals of dimensionality reduction and preservation of important data information.

[0083] During training, the loss function of the variational autoencoder consists of two parts: reconstruction loss and KL divergence loss. Its purpose is to constrain the distribution of the latent space to a standard normal distribution while maintaining the quality of the input data reconstruction. The overall form of the loss function L is:

[0084]

[0085] Reconstruction loss Used to measure the reconstruction matrix The difference between the preprocessed matrix M and the preprocessed matrix M' is used to ensure that the data generated by the decoder is as close as possible to the real data. This is expressed using the binary cross-entropy form as follows:

[0086]

[0087] Where, n rows and n cols Reconstructed matrices The number of rows and columns of the preprocessed matrix M.

[0088] KL divergence D KL (q(Z|M”)||p(Z)) is used to measure the difference between the latent coding distribution q(Z|M”) and the standard normal distribution p(Z). The specific calculation method is as follows:

[0089]

[0090] Where k is the dimension index of the latent variable matrix Z; The square of the mean of the values ​​in the k-th dimension of the latent variable matrix Z; Let be the variance of the values ​​in the k-th dimension of the latent variable matrix Z.

[0091] When the loss function L is below the threshold δ = 10 -2 When the training ends, the encoder outputs the latent variable matrix Z. The latent variable matrix Z is in matrix form, where each row represents the embedding representation of a control element in the latent space, and each column is a dimension of the latent space.

[0092] (4) Calculate the similarity of the open modes of the control element, corresponding to step S4.

[0093] The latent variable matrix Z output by the variational autoencoder is used to calculate the similarity of the open modes of the control elements. First, each row of the latent variable matrix Z is normalized so that the sum of each row is 1. Assuming Z is an n×d matrix, where n is the number of control elements and d is the dimension of the features, the normalization calculation formula is as follows:

[0094]

[0095] Where Z'[i][k] is the element in the i-th row and k-th position of the standardized matrix Z'; Z[i][k] is the element in the i-th row and k-th position of the latent variable matrix Z; S i Let Z be the sum of the elements in the i-th row of the latent variable matrix Z.

[0096] Each row of the standardized matrix Z' is a valid probability distribution, suitable for calculating the Wasserstein distance between control elements. To calculate the Wasserstein distance between each pair of control elements, i.e., Z'[i] and Z'[j], a matrix D needs to be constructed, where D... ij The Euclidean distance between row i and row j in the standardized matrix Z' is defined as:

[0097]

[0098] The Wasserstein distance can then be calculated using the following optimization problem:

[0099]

[0100] Where Π(p,q) is the set of all possible transportation plans, and p and q are vectors representing the probability distributions of the i-th and j-th rows of the standardized matrix Z', respectively. i =Z'[i][k], q j =Z'[j][k];T ij Let T be the amount of transport to move the quantity from row i to row j; the symbol inf denotes the "minimum lower bound", which means that among all possible transport plans T, the transport cost is the minimum, i.e., the optimal transport plan. The optimal transport plan is solved by using the POT package in Python or the optimize module in the SciPy package.

[0101] Calculating the Wasserstein distance between any two rows of the standardized matrix Z' yields an n×n distance matrix W, where the row and column indices are control elements, and the values ​​are the Wasserstein distances between control elements. Next, a local adaptive Gaussian kernel is used to transform the distance matrix W into an affinity matrix A. Specifically, the kernel bandwidth of each control element is defined by its kNN Wasserstein distance, and the affinity between control element i and control element j is defined as:

[0102]

[0103] Among them, A ij To regulate the affinity between control element i and control element j; W ij σ is the Wasserstein distance between control element i and control element j; i and σ j These are the local kernel bandwidths of control element i and control element j, respectively. The local kernel bandwidth is determined based on the Wasserstein distance of the k nearest neighbor, with k = 5 by default.

[0104] Next, a distance penalty term is introduced to adjust affinity:

[0105]

[0106] Among them, A' ij The adjusted affinity between control element i and control element j; d ij τ represents the distance (in kb) between regulatory elements i and j on the genome; τ is the distance scale parameter used to control the strength of the distance penalty, which defaults to 100.

[0107] The calculated adjusted affinity matrix A' contains information about the similarity of the open modes of regulatory elements. By setting an affinity threshold (default 0.3), pairs of regulatory elements with potential interactions in this biological process can be selected. Based on the promoter or enhancer to which the regulatory element pair belongs, enhancer-promoter pairs with potential interactions can be selected.

[0108] In this embodiment, the reading matrix is ​​used as input, the latent variable matrix is ​​calculated using a variational autoencoder, and the enhancer-promoter pairs with potential interactions are determined based on the latent variable matrix. The reading matrix includes the readings of each regulatory element in each cell of the object to be tested, and the latent variable matrix includes the embedded representation of each regulatory element.

[0109] The latent variable matrix is ​​calculated using a variational autoencoder, with the reading matrix as input. Specifically, this includes:

[0110] 1) Preprocess the reading matrix to obtain the preprocessed matrix. The preprocessing includes logarithmic transformation and standardization.

[0111] 2) Using the preprocessed matrix as input, the encoder of the variational autoencoder is used to process it to obtain the latent variable matrix.

[0112] Specifically, the identification of enhancer-promoter pairs with potential interactions based on the latent variable matrix includes:

[0113] 1) Standardize each row of the latent variable matrix so that the sum of the elements in each row is 1, and obtain the standardized matrix.

[0114] 2) Calculate the Wasserstein distance between any two rows in the standardized matrix to obtain the distance matrix. The number of rows and columns of the distance matrix is ​​equal to the total number of control elements. The element in the i-th row and j-th column of the distance matrix is ​​the Wasserstein distance between the i-th control element and the j-th control element.

[0115] 3) The distance matrix is ​​processed using the Gaussian kernel function to obtain the affinity matrix. The number of rows and columns of the affinity matrix is ​​equal to the total number of control elements. The element in the i-th row and j-th column of the affinity matrix is ​​the affinity between the i-th control element and the j-th control element.

[0116] 4) Apply a distance penalty to the affinity matrix to obtain the adjusted affinity matrix. The number of rows and columns of the adjusted affinity matrix is ​​equal to the total number of control elements. The element in the i-th row and j-th column of the adjusted affinity matrix is ​​the adjusted affinity between the i-th control element and the j-th control element.

[0117] 5) For each element in the adjusted affinity matrix, determine whether the element is greater than the preset affinity threshold. If so, the two regulatory elements corresponding to the element are regarded as a pair of regulatory elements with potential interaction. If the pair of regulatory elements with potential interaction includes a promoter and an enhancer, then the pair of regulatory elements with potential interaction is regarded as an enhancer-promoter pair with potential interaction.

[0118] (5) Infer the gene co-expression network, corresponding to step S7.

[0119] EnhancerNet constructs a gene co-expression network by calculating mutual information between genes. This step uses scRNA-seq data, which, like scATAC-seq data, pertains to the same biological process. The scRNA-seq data is represented by a matrix X, where X... ijLet represent the expression level of the j-th gene in the i-th cell. The goal is to calculate the mutual information between the expression levels of two genes to construct a co-expression network. The formula for calculating mutual information is:

[0120]

[0121] Wherein, I(G1;G2) represents the mutual information between gene G1 and gene G2; and Let x and y be the sets of expression levels of genes G1 and G2 in all cells, respectively; p(x,y) is the joint probability density function of gene G1 expression level x and gene G2 expression level y; p(x) and p(y) are the marginal probability density functions of gene G1 expression level x and gene G2 expression level y, respectively.

[0122] Using kernel density estimation, EnhancerNet can estimate these probability density functions. For each gene expression level, EnhancerNet can obtain its corresponding probability density using kernel density estimation, with the following form:

[0123]

[0124] Where n is the number of cells; h is the bandwidth parameter, with a default value of h=1, which controls the smoothness; and K is the Gaussian kernel function, whose formula is:

[0125]

[0126] here Substituting the Gaussian kernel function into the kernel density estimation formula, we get:

[0127]

[0128] This formula can be further simplified to:

[0129]

[0130] Where, x i denoted as the expression level of gene G1 in the i-th cell.

[0131] The joint probability density p(x,y) can be estimated similarly:

[0132]

[0133] Among them, h x and h y Also a bandwidth parameter, the default is h. x =h y =1; y i denoted as the expression level of gene G2 in the i-th cell.

[0134] After calculating the mutual information between any two genes, EnhancerNet selects a percentile of mutual information as the mutual information threshold (default 80th percentile). If the mutual information between genes is greater than this threshold, they are considered to have a potential co-expression relationship and are used to construct a gene co-expression network. In this network, nodes represent genes, and the weights of the edges are given by the mutual information values. The higher the value, the stronger the co-expression relationship between genes.

[0135] This embodiment calculates gene-gene mutual information based on scRNA-seq data and constructs a co-expression gene network based on this mutual information. The genes include the coding genes and target genes of transcription factors. It should be noted that transcription factors are genes before transcription and translation; therefore, the association between the coding genes and target genes of transcription factors can be determined, which is also the association between transcription factors and target genes.

[0136] The scRNA-seq data includes the expression level of each gene in each cell of the target organism. Based on the scRNA-seq data, the mutual information between genes is calculated, and a co-expression gene network is constructed based on this mutual information. Specifically, this includes:

[0137] 1) For each gene, determine the marginal probability density function of the gene based on the expression level of the gene in each cell of the subject to be tested.

[0138] 2) For any two gene pairs, the joint probability density function of the gene pair is determined based on the expression levels of the two genes in each cell of the target organism. The mutual information of the gene pair is calculated based on the joint probability density function and the marginal probability density functions of the two genes in the gene pair. The mutual information of all gene pairs constitutes the mutual information between genes.

[0139] 3) If the mutual information of a gene pair is greater than a preset mutual information threshold, the two genes in the gene pair are connected to obtain the connection relationship between any two genes.

[0140] 4) Using each gene as a node, connect the nodes based on the connection relationship between any two genes to obtain a co-expressed gene network.

[0141] (6) Construct an enhancer-driven gene regulatory network, corresponding to steps S6 and S8.

[0142] By integrating information obtained from scATAC-seq data regarding transcription factors to enhancers, enhancers to promoters, and promoters to target genes, EnhancerNet can construct a basic gene regulatory network. This network describes the regulatory pathway of transcription factors-enhancer-promoters-target genes and can be further simplified to a regulatory network of transcription factors regulating target genes. Based on the gene co-expression network obtained from scRNA-seq data, EnhancerNet can remove redundancy from the regulatory network of transcription factors regulating target genes using an intersection method, the formula of which is:

[0143] R′ ij =R ij ∩C ij ;

[0144] Among them, R' ij For intersection; R ij This refers to the regulatory relationship between transcription factor i and target gene j inferred from scATAC-seq data, i.e., the transcription factor-target gene regulatory network; C ij This represents the co-expression relationship between transcription factor i and target gene j, inferred from scRNA-seq data, i.e., a gene co-expression network. Ultimately, the regulatory network obtained by EnhancerNet is the "enhancer-driven gene regulatory network," which describes the regulatory network pathways by which specific transcription factors exert their effects on target genes by binding enhancers in the biological processes under study.

[0145] This embodiment generates a transcription factor-enhancer-promoter-target gene regulatory network based on the potential binding sites of each transcription factor on the enhancer of the target object, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter. The transcription factor-enhancer-promoter-target gene regulatory network is then simplified to obtain a transcription factor-target gene regulatory network. The intersection of the transcription factor-target gene regulatory network and the co-expressed gene network is calculated. This intersection is then expanded based on the potential binding sites of each transcription factor on the enhancer of the target object, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, resulting in a transcription factor-enhancer-promoter-target gene regulatory network.

[0146] Specifically, based on the potential binding sites of each transcription factor on the enhancer of the target object, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated. This network is then simplified to obtain the transcription factor-target gene regulatory network, which specifically includes:

[0147] 1) Based on the potential binding sites of each transcription factor on the enhancer of the target object, determine the connection relationship between the transcription factor and the enhancer.

[0148] 2) Based on the existence of potentially interacting enhancer-promoter pairs, determine the connectivity between enhancers and promoters.

[0149] 3) Based on the target gene corresponding to each promoter, determine the connection relationship between the promoter and the target gene.

[0150] 4) Based on the connection relationships between transcription factors and enhancers, between enhancers and promoters, and between promoters and target genes, transcription factors, enhancers, promoters, and target genes are connected to generate a transcription factor-enhancer-promoter-target gene regulatory network.

[0151] 5) If the transcription factors in the transcription factor-enhancer-promoter-target gene regulatory network can be linked to the target gene through the enhancer and the promoter, then the transcription factors and the target gene are linked to determine the connection relationship between the transcription factors and the target gene. Based on the connection relationship between the transcription factors and the target gene, the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain the transcription factor-target gene regulatory network.

[0152] This embodiment uses the artificial hematopoietic process as an example. The scRNA-seq data used comes from the Human Cell Atlas: https: / / data.humancellatlas.org / explore / projects / 091cf39b-01bc-42e5-9437-f419a66c8a45, and the scATAC-seq data comes from the GEO database, accession number GSE96769. EnhancerNet identified an enhancer-driven gene regulatory network with PBX1, one of the main regulators of hematopoietic stem and progenitor cells (HSPCs), at its core, such as... Figure 2 As shown, the downstream target genes of PBX1 (highlighted in blue) are all regulated by PBX1 binding to enhancers, and the transcriptional activity of downstream target genes is regulated through the connection between enhancers and their target genes. Taking MEIS1, a downstream target gene of PBX1 and also an HSPC transcription regulator, as an example, the transcription factor PBX1 binds to the enhancer region of stain 2 (located at chr2: 66362408-66362565), and further regulates it remotely by interacting with the promoter region of MEIS1 (located at chr2: 66433420-66433767).

[0153] Example 2

[0154] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 3 As shown, this computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides the environment for the operating system and computer programs stored in the non-volatile storage media. The database stores data. The I / O interfaces are used for exchanging information between the processor and external devices. The communication interface is used for communicating with external terminals via a network. When executed by the processor, the computer program implements a method for constructing a gene regulatory network.

[0155] Those skilled in the art will understand that Figure 3 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0156] In one exemplary embodiment, a computer device is provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the method for constructing a gene regulatory network in Embodiment 1.

[0157] Example 3

[0158] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the method for constructing a gene regulatory network in Embodiment 1.

[0159] Example 4

[0160] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the method for constructing a gene regulatory network in Embodiment 1.

[0161] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with relevant regulations.

[0162] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0163] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A method for constructing a gene regulatory network, characterized in that, The method for constructing the gene regulatory network comprises the following steps: obtaining scATAC-seq data, scRNA-seq data, position data of candidate regulatory elements and a genome annotation file of a to-be-detected object; the to-be-detected object is a tissue of a to-be-detected organism; the candidate regulatory elements comprise candidate promoters and candidate enhancers; the position data of the candidate regulatory elements are obtained by downloading from a SCREEN database; the position data of the candidate regulatory elements comprise a number of a chromosome where the candidate regulatory elements are located, a start site of the candidate regulatory elements on the chromosome and an end site of the candidate regulatory elements on the chromosome; annotating the candidate promoters based on the position data of the candidate regulatory elements and the genome annotation file, determining target genes corresponding to the candidate promoters, and obtaining annotation data of the candidate regulatory elements; the annotation data of the candidate regulatory elements comprise position data of each candidate promoter, corresponding target genes and position data of each candidate enhancer; determining regulatory elements based on the scATAC-seq data and the annotation data of the candidate regulatory elements, and calculating reads of the regulatory elements in each cell of the to-be-detected object, to obtain annotation data of the regulatory elements; the regulatory elements are open candidate regulatory elements; the regulatory elements comprise promoters and enhancers; the annotation data of the regulatory elements comprise position data of each promoter, corresponding target genes and reads in each cell of the to-be-detected object, and position data of each enhancer and reads in each cell of the to-be-detected object; using a variational autoencoder to calculate a latent variable matrix based on a read matrix as input, and determining enhancer-promoter pairs with potential interaction based on the latent variable matrix; the read matrix comprises reads of each regulatory element in each cell of the to-be-detected object; the latent variable matrix comprises embedded representations of each regulatory element; performing motif scanning on position data of each enhancer and reads in each cell of the to-be-detected object based on a motif database, to determine potential binding sites of each transcription factor of the to-be-detected object on the enhancer; generating a transcription factor-enhancer-promoter-target gene regulatory network based on the potential binding sites of each transcription factor of the to-be-detected object on the enhancer, the enhancer-promoter pairs with potential interaction and target genes corresponding to each promoter, and simplifying the transcription factor-enhancer-promoter-target gene regulatory network to obtain a transcription factor-target gene regulatory network; calculating mutual information between genes based on the scRNA-seq data, and constructing a co-expression gene network based on the mutual information between the genes; the genes comprise coding genes of transcription factors and target genes. An intersection of the transcription factor-target gene regulatory network and the co-expression gene network is calculated, and the intersection is expanded based on potential binding sites of each transcription factor of the subject on enhancers, enhancer-promoter pairs with potential interactions, and target genes corresponding to each promoter, to obtain a transcription factor-enhancer-promoter-target gene gene regulatory network.

2. The method of constructing a gene regulatory network according to claim 1, wherein, The candidate promoters are annotated based on the position data of the candidate regulatory elements and the genomic annotation file, and target genes corresponding to the candidate promoters are determined, specifically including: The candidate promoters are annotated using the annotatePeak function with the position data of the candidate regulatory elements and the genomic annotation file as input, and target genes corresponding to the candidate promoters are determined; Based on the scATAC-seq data and the annotation data of the candidate regulatory elements, regulatory elements are determined, and the reads of the regulatory elements in each cell of the subject are calculated, specifically including: The regulatory elements are determined using the bedtools tool with the scATAC-seq data and the annotation data of the candidate regulatory elements as input, and the reads of the regulatory elements in each cell of the subject are calculated.

3. The method of constructing a gene regulatory network according to claim 1, wherein, The latent variable matrix is calculated using the variational autoencoder with the read matrix as input, specifically including: The read matrix is preprocessed to obtain a preprocessed matrix; the preprocessing includes logarithmic transformation and standardization processing; The latent variable matrix is obtained by processing the preprocessed matrix using the encoder of the variational autoencoder.

4. The method of constructing a gene regulatory network according to claim 1, wherein, The enhancer-promoter pairs with potential interactions are determined based on the latent variable matrix, specifically including: Each row of the latent variable matrix is standardized so that the sum of the elements of each row is 1, to obtain a standardized matrix; The Wasserstein distance between any two rows in the standardized matrix is calculated to obtain a distance matrix; the number of rows and columns of the distance matrix is equal to the total number of regulatory elements, and the element in the ith row and jth column of the distance matrix is the Wasserstein distance between the ith regulatory element and the jth regulatory element; The distance matrix is processed using a Gaussian kernel function to obtain an affinity matrix; the number of rows and columns of the affinity matrix is equal to the total number of regulatory elements, and the element in the ith row and jth column of the affinity matrix is the affinity between the ith regulatory element and the jth regulatory element; Distance penalty is applied to the affinity matrix to obtain an adjusted affinity matrix; the number of rows and columns of the adjusted affinity matrix is equal to the total number of regulatory elements, and the element in the ith row and jth column of the adjusted affinity matrix is the adjusted affinity between the ith regulatory element and the jth regulatory element; For each element in the adjusted affinity matrix, it is judged whether the element is greater than a preset affinity threshold; if yes, the two regulatory elements corresponding to the element are taken as a regulatory element pair with potential interaction, and if the regulatory element pair with potential interaction includes a promoter and an enhancer, the regulatory element pair with potential interaction is taken as an enhancer-promoter pair with potential interaction.

5. The method of constructing a gene regulatory network according to claim 1, wherein, Based on the motif database, motif scanning is performed on the position data of each enhancer and the read count in each cell of the to-be-detected object to determine the potential binding site of each transcription factor of the to-be-detected object on the enhancer, specifically including: Based on the motif database, motif scanning is performed on the position data of each enhancer and the read count in each cell of the to-be-detected object to determine the potential binding site of each transcription factor of the to-be-detected object on the enhancer, specifically including:

6. The method of constructing a gene regulatory network according to claim 1, wherein, Based on the potential binding site of each transcription factor of the to-be-detected object on the enhancer, the enhancer-promoter pair with potential interaction and the target gene corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated, and the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain a transcription factor-target gene regulatory network, specifically including: Based on the potential binding site of each transcription factor of the to-be-detected object on the enhancer, the connection relationship between the transcription factor and the enhancer is determined; Based on the enhancer-promoter pair with potential interaction, the connection relationship between the enhancer and the promoter is determined; Based on the target gene corresponding to each promoter, the connection relationship between the promoter and the target gene is determined; Based on the connection relationship between the transcription factor and the enhancer, the connection relationship between the enhancer and the promoter and the connection relationship between the promoter and the target gene, the transcription factor, the enhancer, the promoter and the target gene are connected to generate a transcription factor-enhancer-promoter-target gene regulatory network; If the transcription factor in the transcription factor-enhancer-promoter-target gene regulatory network can be connected to the target gene through the enhancer and the promoter, the transcription factor and the target gene are connected, the connection relationship between the transcription factor and the target gene is determined, and based on the connection relationship between the transcription factor and the target gene, the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain a transcription factor-target gene regulatory network.

7. The method of constructing a gene regulatory network according to claim 1, wherein, If the transcription factor in the transcription factor-enhancer-promoter-target gene regulatory network can be connected to the target gene through the enhancer and the promoter, the transcription factor and the target gene are connected, the connection relationship between the transcription factor and the target gene is determined, and based on the connection relationship between the transcription factor and the target gene, the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain a transcription factor-target gene regulatory network. If the transcription factor in the transcription factor-enhancer-promoter-target gene regulatory network can be connected to the target gene through the enhancer and the promoter, the transcription factor and the target gene are connected, the connection relationship between the transcription factor and the target gene is determined, and based on the connection relationship between the transcription factor and the target gene, the transcription factor-enhancer-promoter-target gene regulatory network is simplified to obtain a transcription factor-target gene regulatory network. For each gene, based on the expression amount of the gene in each cell of the to-be-detected object, the marginal probability density function of the gene is determined; For each gene pair consisting of any two genes, a joint probability density function of the gene pair is determined based on expression amounts of the two genes included in the gene pair in each cell of the subject to be detected, and mutual information of the gene pair is calculated based on the joint probability density function of the gene pair and marginal probability density functions of the two genes included in the gene pair; mutual information of all the gene pairs constitutes mutual information between genes; If the mutual information of the gene pair is greater than a preset mutual information threshold, the two genes included in the gene pair are connected, and a connection relationship between any two genes is obtained; Each gene is taken as a node, and the nodes are connected based on the connection relationship between any two genes, and a co-expression gene network is obtained.

8. A computer device comprising: A memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that the processor executes the computer program to implement the method for constructing a gene regulatory network according to any one of claims 1-7.

9. A computer readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the method for constructing a gene regulatory network according to any one of claims 1-7.

10. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the method for constructing a gene regulatory network according to any one of claims 1-7.

Images