Construction method of gene regulation and control network and related device

Abstract

The invention discloses a gene regulatory network construction method and a related device, and relates to the technical field of gene regulatory network construction.The method comprises the steps that position data of candidate regulatory elements are downloaded from an SCREEN database, a genome annotation file and scATAC-seq data are combined, a promoter and an enhancer are accurately positioned, a variational auto-encoder is introduced, and a gene regulatory element is constructed. The method comprises the following steps: accurately determining an enhancer-promoter pair with potential interaction, associating an enhancer with a target gene, and processing a co-expression gene network based on a potential binding site of each transcription factor of a to-be-detected object on the enhancer, the enhancer-promoter pair with potential interaction and the target gene corresponding to each promoter. And obtaining the gene regulatory network. The method can solve the problems that the promoter and the enhancer cannot be accurately positioned, the enhancer and the target gene are difficult to associate, direct and indirect relationships are difficult to distinguish, false positive regulation relationships are excessive and the like.

Description

Technical Field

[0001] The present application relates to the technical field of gene regulatory network construction, and in particular to a method and related device for constructing an enhancer-driven gene regulatory network based on EnhancerNet. Background Art

[0002] The Gene Regulatory Network (GRN) records the relationship between transcription factors (TFs) and their regulated target genes, reveals the gene expression regulation mechanism inside cells, and can help understand the molecular basis of cell and tissue behavior. In-depth understanding of the gene regulatory network is extremely important for revealing the molecular mechanism that controls cell behavior.

[0003] At present, gene regulatory networks are often inferred through experimental methods such as chromatin immunoprecipitation sequencing technology, interference experiments of specific transcription factors combined with transcriptome sequencing, but it is still difficult to infer gene regulatory networks on a large scale. The emergence of single-cell RNA sequencing (scRNA-seq) provides unprecedented dynamic resolution of the transcriptome, which also provides new ideas for the inference of gene regulatory networks and promotes the generation and development of algorithm models for inferring gene regulatory networks based on single-cell sequencing data. Although a series of algorithm models for inferring gene regulatory networks based on single-cell sequencing data have been developed, such as TENET, SCENIC, Dictys, etc., the failure to accurately locate promoters and enhancers, the difficulty in associating enhancers with target genes, the difficulty in distinguishing direct and indirect relationships, and the excessive number of false positive regulatory relationships are still severe challenges faced by algorithm models for inferring gene regulatory networks. Summary of the invention

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

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

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

[0007] Obtaining scATAC-seq data, scRNA-seq data, location data of candidate regulatory elements, and a genome annotation file of an object to be detected; the object to be detected is a tissue of an organism to be detected; the candidate regulatory elements include candidate promoters and candidate enhancers, and the location data of the candidate regulatory elements are downloaded from the SCREEN database, and the location data of the candidate regulatory elements include the number of the chromosome where the candidate regulatory element is located, the starting site of the candidate regulatory element on the chromosome where the candidate regulatory element is located, and the ending site of the candidate regulatory element on the chromosome where the candidate regulatory element is located;

[0008] Based on the position data of the candidate regulatory element and the genome annotation file, the candidate promoter is annotated, the target gene corresponding to the candidate promoter is determined, and the annotation data of the candidate regulatory element is obtained; the annotation data of the candidate regulatory element includes the position data of each candidate promoter and the corresponding target gene and the position data of each candidate enhancer;

[0009] Based on the scATAC-seq data and the annotation data of the candidate regulatory elements, the regulatory elements are determined, and the readings of the regulatory elements in each cell of the object to be detected are calculated to obtain the annotation data of the regulatory elements; the regulatory elements are open candidate regulatory elements, the regulatory elements include promoters and enhancers, and the annotation data of the regulatory elements include the position data of each promoter, the corresponding target gene and the readings in each cell of the object to be detected, and the position data of each enhancer and the readings in each cell of the object to be detected;

[0010] Taking the read matrix as input, using the variational autoencoder to calculate the latent variable matrix, and determining the enhancer-promoter pairs with potential interactions based on the latent variable matrix; the read matrix includes the reads of each regulatory element in each cell of the object to be detected; 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 object to be detected are subjected to motif scanning to determine the potential binding site of each transcription factor of the object to be detected on the enhancer;

[0012] Based on the potential binding sites of each transcription factor of the object to be detected 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 a transcription factor-target gene regulatory network;

[0013] The mutual information between genes is calculated based on the scRNA-seq data, and a co-expression gene network is constructed based on the mutual information between genes; the genes include 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 object to be detected 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] In a second aspect, the present application provides a computer device, comprising: 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-mentioned method for constructing a gene regulatory network.

[0016] In a third aspect, the present application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-mentioned method for constructing a gene regulatory network.

[0017] In a fourth aspect, the present application provides a computer program product, including a computer program, which implements the above-mentioned method for constructing a gene regulatory network when executed by a processor.

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

[0019] The present application provides a method for constructing a gene regulatory network and a related device, downloading the position data of candidate regulatory elements from the SCREEN database, further combining the genome annotation file and scATAC-seq data, and being able to accurately locate promoters and enhancers, solving the problem of failing to accurately locate promoters and enhancers. The introduction of a variational autoencoder can accurately determine the existence of a potential interacting enhancer-promoter pair, and since the promoter is associated with the target gene, the enhancer can be further associated with the target gene, solving the problem of difficulty in associating enhancers with target genes. Based on the potential binding sites of each transcription factor of the object to be detected on the enhancer, the existence of a potential interacting enhancer-promoter pair and the target gene corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated, and the transcription factor-target gene regulatory network is further simplified, and the intersection of the transcription factor-target gene regulatory network and the co-expression gene network is calculated, and the gene regulatory network of the transcription factor-enhancer-promoter-target gene is further expanded to solve the problem of difficulty in distinguishing direct and indirect relationships and too many false positive regulatory relationships. The present application can solve the problems of failure to accurately locate promoters and enhancers, difficulty in associating enhancers with target genes, difficulty in distinguishing direct and indirect relationships, and too many false positive regulatory relationships. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0021] Figure 1 A schematic diagram of a process for constructing a gene regulatory network provided in Example 1 of the present application.

[0022] Figure 2 This is a schematic diagram of the gene regulatory network driven by the enhancer with PBX1 as the core in the human hematopoietic system provided in Example 1 of the present application.

[0023] Figure 3 A schematic diagram of the structure of a computer device provided in Example 2 of the present application. DETAILED DESCRIPTION

[0024] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0025] Example 1

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

[0027] Step S1, obtaining scATAC-seq data, scRNA-seq data, location data of candidate regulatory elements and genome annotation files of the object to be detected; the object to be detected is the tissue of the organism to be detected; the candidate regulatory elements include candidate promoters and candidate enhancers, the location data of the candidate regulatory elements are downloaded from the SCREEN database, and the location data of the candidate regulatory elements include the number of the chromosome where the candidate regulatory element is located, the starting site of the candidate regulatory element on the chromosome where it is located, and the ending site of the candidate regulatory element on the chromosome where it is located.

[0028] Step S2, based on the position data of the candidate regulatory element and the genome annotation file, annotate the candidate promoter, determine the target gene corresponding to the candidate promoter, and obtain the annotation data of the candidate regulatory element; the annotation data of the candidate regulatory element includes the position data of each candidate promoter and the corresponding target gene and the position data of each candidate enhancer.

[0029] Step S3, based on the scATAC-seq data and the annotation data of the candidate regulatory elements, determine the regulatory elements, and calculate the readings of the regulatory elements in each cell of the object to be detected 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 position data of each promoter, the corresponding target gene and the readings in each cell of the object to be detected, and the position data of each enhancer and the readings in each cell of the object to be detected.

[0030] Step S4, using the reading matrix as input, using the variational autoencoder to calculate the latent variable matrix, and determining the enhancer-promoter pairs with potential interactions based on the latent variable matrix; the reading matrix includes the readings of each regulatory element in each cell of the object to be detected; 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 object to be detected, and determine the potential binding site of each transcription factor on the enhancer of the object to be detected.

[0032] Step S6, based on the potential binding sites of each transcription factor of the object to be detected on the enhancer, the enhancer-promoter pairs with potential interactions and the target genes corresponding to each promoter, generate a transcription factor-enhancer-promoter-target gene regulatory network, and simplify the transcription factor-enhancer-promoter-target gene regulatory network to obtain a transcription factor-target gene regulatory network.

[0033] Step S7, calculating the mutual information between genes based on the scRNA-seq data, and constructing 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 object to be detected 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] Implement the above-mentioned steps S1 to S8, the present embodiment downloads the position data of the candidate regulatory elements from the SCREEN database, and further combines the genome annotation file and scATAC-seq data, so that the promoter and enhancer can be accurately located, solving the problem of failing to accurately locate the promoter and enhancer. Introducing a variational autoencoder, it is possible to accurately determine the existence of a potential interacting enhancer-promoter pair, and since the promoter is associated with the target gene, the enhancer can be further associated with the target gene, solving the problem of being difficult to associate the enhancer with the target gene. Based on the potential binding site of each transcription factor of the object to be detected on the enhancer, there are potential interacting enhancer-promoter pairs and the target gene corresponding to each promoter, a transcription factor-enhancer-promoter-target gene regulatory network is generated, and the transcription factor-target gene regulatory network is further simplified, and the intersection of the transcription factor-target gene regulatory network and the co-expression gene network is calculated, and the gene regulatory network of the transcription factor-enhancer-promoter-target gene is further expanded to obtain, solving the problem that direct and indirect relationships are difficult to distinguish and false positive regulatory relationships are too many.

[0036] In view of the problems existing in the algorithm models for inferring gene regulatory networks, such as failure to accurately locate promoters and enhancers, difficulty in associating enhancers and target genes, difficulty in distinguishing direct and indirect relationships, and too many false positive regulatory relationships, this embodiment develops an algorithm model EnhancerNet (i.e., a method for constructing a gene regulatory network used in this embodiment), which is used to infer the gene regulatory network of transcription factor-enhancer-promoter-target gene on a large scale based on scATAC-seq data (Single-cell Assay for Transposase-Accessible Chromatinusing sequencing) and scRNA-seq data.

[0037] The following is a brief introduction to the solutions to the four problems in this embodiment:

[0038] (1) Solutions that fail to accurately locate promoters and enhancers

[0039] Promoters and enhancers are the two main 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 act upstream, downstream, or even in introns of genes to enhance the expression of specific genes by changing chromatin conformation or forming DNA loops. The challenge of accurately locating cis-regulatory elements lies in the complexity and diversity of these two key regulatory regions: promoters and enhancers. Typical promoter regions contain some conserved sequence elements, such as TATA boxes, CAAT boxes, and GC boxes, which help to identify and bind transcription factors. However, the promoter sequences of different genes may differ significantly, which increases the difficulty of accurate positioning. Even in the same species, the length, sequence composition, and location of functional elements of promoter sequences may vary. In particular, the promoters of some genes may lack clear conserved sequence elements, making it difficult to accurately identify them through sequence alignment or conventional annotation methods. Unlike promoters, enhancers do not have strict positional restrictions and specific sequence characteristics, which makes their identification more complicated. Enhancer regions are usually associated with specific chromatin marks such as H3K27ac and H3K4me1. The presence of these chromatin marks helps to identify enhancers. However, enhancers have large sequence diversity, and enhancers between different genes may have significant differences in sequence composition and regulatory properties, further increasing the difficulty of their precise positioning.

[0040] The SCREEN database in the Encyclopedia of DNA Elements (ENCODE) project provides accurate positioning information of promoters and enhancers: it integrates and analyzes a large amount of multi-omics data, defines candidate promoters as DNA sequences with high DNase-seq signals, high H3K4me3 signals, and whose centers are located within 200bp of the transcription start site, and defines candidate enhancers as DNA sequences with high DNase-seq signals and high H3K27ac signals and not overlapping with the region of the candidate promoter. At present, no algorithm model uses this database. This embodiment introduces the SCREEN database, extracts the position data of candidate regulatory elements (i.e., candidate promoters and candidate enhancers) in the SCREEN database, identifies candidate regulatory elements near the transcription start site as candidate promoters according to the genome annotation file, identifies other candidate regulatory elements as candidate enhancers, and combines scATAC-seq data to obtain open candidate regulatory elements, and uses the open candidate regulatory elements as regulatory elements. The regulatory elements include promoters and enhancers, thereby solving the problem that the current relevant algorithm models fail to accurately locate promoters and enhancers.

[0041] (2) Solutions to the difficulty of associating enhancers with target genes

[0042] At present, it is still a huge challenge to link enhancers to their target genes. Although Hi-C technology can capture the three-dimensional structure of chromatin and reveal the interactions between remote regions in the genome, it performs poorly in the precise positioning of remote interactions. In addition, Hi-C technology is difficult to achieve single-cell resolution, resulting in mixed data of different cell types, making it difficult to identify specific chromatin structures and gene regulation patterns in different cell types. This embodiment can decipher the long-range interactions between enhancers and promoters at the single-cell level by algorithmically modeling scATAC-seq data. First, the original high-dimensional sparse data determined based on scATAC-seq data is input into the encoder of the variational autoencoder (VAE), which compresses the high-dimensional sparse data features of each regulatory element (also called chromatin region) into a low-dimensional latent space representation, thereby retaining the main features of the input data, and at the same time, through regularization constraints, ensuring that the generated latent variable matrix conforms to the predetermined latent distribution. Then, the decoder of the variational autoencoder reconstructs the chromatin accessibility data according to 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 embedded representation of a regulatory element in the latent space, and each column is a dimension of the latent space. Through the processing of the above encoder and decoder, the variational autoencoder can deal with the sparsity and noise problems of the data while capturing the intrinsic structure of the chromatin accessibility data. Then, the extracted latent variable matrix is ​​used to calculate the Wasserstein distance between the regulatory elements to obtain the distance matrix. The Wasserstein distance is sensitive to the difference in probability distribution and can better reflect the difference in the distribution of different regulatory elements in a complex latent space. Unlike the Euclidean distance, the Wasserstein distance is more suitable for processing asymmetric probability distributions, which is common in the latent space of the variational autoencoder. Then, the distance matrix is ​​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 interactions are more common and efficient. Therefore, a distance penalty term is added to the calculated affinity matrix, so that the affinity between regulatory elements that are physically close is relatively enhanced, while regulatory elements that are farther away are subject to additional penalties. The penalized affinity (i.e., adjusted affinity) is used to measure the potential interaction between regulatory elements. Finally, by setting an affinity threshold, enhancer and promoter pairs on the same DNA chain with an adjusted affinity greater than the affinity threshold are selected to indicate the existence of potential long-range interactions between the two, and enhancers and target genes are further associated based on the association between promoters and target genes.

[0043] (3) Solutions for the difficulty in distinguishing direct and indirect relationships and excessive false positive regulatory relationships

[0044] At present, most of the algorithm models for inferring gene regulatory networks based on single-cell sequencing data are based on the measurement of co-expression relationships between genes. The gene regulatory networks inferred by such algorithm models lack interaction information between transcription factors and target genes, and have inherent defects such as difficulty in distinguishing direct and indirect relationships and too many false positive regulatory relationships. The EnhancerNet developed in this embodiment can construct a co-expression gene network based on scRNA-seq data, and measure the co-expression relationship by calculating the partial correlation coefficient (i.e., mutual information) between paired genes in large quantities. Compared with other methods for calculating co-expression relationships, the method of calculating the partial correlation coefficient has the advantage of being able to evaluate the association between the expression levels of two genes under the condition of eliminating the influence of other genes. However, the co-expression gene network generated at this time still contains false positive regulatory relationships. This embodiment further uses the gimmemotifs algorithm to identify the binding sites of each transcription factor in the precisely located enhancer based on the transcription factor binding motifs in the CisBPv.2 database, and obtains a transcription factor-target gene regulatory network. The gene regulatory network in the biological process under study can be obtained by intersecting the transcription factor-target gene regulatory network with the co-expression gene network. In addition, the target gene directly targeted and regulated by the transcription factor can be directly identified through the identification of the transcription factor binding motifs, thereby avoiding the generation of too many indirect regulatory relationships and removing most of the false positive regulatory relationships.

[0045] The following, combined Figure 1 The algorithm model of this embodiment is introduced in detail, which is mainly summarized into six steps: precise positioning of open promoters and enhancers, use of variational autoencoders for feature extraction, calculation of 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) Accurate positioning of open promoters and enhancers, corresponding to steps S1 to S3.

[0047] For the biological process being studied, that is, the object to be detected is the tissue of the organism to be detected, which can be tissues of two species, human and mouse, or tissues of other species, single-cell transposase-accessible chromatin sequencing is required to obtain scATAC-seq data, or relevant scATAC-seq data is downloaded from a public database. The raw data of scATAC-seq data is preprocessed by the following steps:

[0048] 1) Use the fastq-dump command of sratoolkit to convert the sra format data (also called format file) into fastq format data. If the fastq format data contains sequencing data of multiple cells, a customized script is required to separate the fastq format data of a single cell from it according to the barcode information.

[0049] 2) Trimmomatic was used to remove adapters from fastq format data.

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

[0051] 4) Use samtools to convert the sam format data into bam format data.

[0052] The scATAC-seq data in this example is data in bam format.

[0053] Single-cell sequencing is performed on the object to be detected to obtain scRNA-seq data.

[0054] The genome annotation file in GTF format was downloaded from the GENCODE database.

[0055] The positional data of two types of candidate regulatory elements, candidate promoters and candidate enhancers, were downloaded from the SCREEN database of the public ENCODE project. The positional data of the candidate regulatory elements included the chromosome number of the candidate regulatory element, the start site of the candidate regulatory element on the chromosome where it was located, and the end site of the candidate regulatory element on the chromosome where it was located. In the SCREEN database, candidate promoters were described as promoter-like signatures (PLS) and candidate enhancers were described as enhancer-like signatures (ELS).

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

[0057] This embodiment further annotates the candidate promoter based on the position data of the candidate regulatory element and the genome annotation file, determines the target gene corresponding to the candidate promoter, and obtains the annotation data of the candidate regulatory element. The annotation data of the candidate regulatory element includes the position data of each candidate promoter and the corresponding target gene and the position data of each candidate enhancer.

[0058] Specifically, the annotatePeak function of the R package ChIPseeker can be used to directly associate PLS with its target gene. In the parameter setting of the annotatePeak function, tssRegion=c(-1500,500), it means that the region from 1500bp upstream to 500bp downstream of the gene transcription start site is designated as the candidate promoter region, that is, the PLS falling in the candidate promoter region is directly assigned as the candidate promoter of the gene. Through further arrangement, a bed format file containing candidate regulatory elements (candidate promoters and candidate enhancers) can be output. The bed format file contains 5 columns of data, which are the chromosome number of the candidate regulatory element, the starting site of the candidate regulatory element on the chromosome where it is located, the termination site of the candidate regulatory element on the chromosome where it is located, the type of the candidate regulatory element (PLS or ELS), and the corresponding target gene (only the target gene corresponding to the PLS).

[0059] At this time, in this embodiment, based on the position data of the candidate regulatory element and the genome annotation file, the candidate promoter is annotated to determine the target gene corresponding to the candidate promoter, specifically including: taking the position data of the candidate regulatory element and the genome annotation file as input, using the annotatePeak function to annotate the candidate promoter to determine the target gene corresponding to the candidate promoter.

[0060] This embodiment further determines the regulatory elements based on the scATAC-seq data and the annotation data of the candidate regulatory elements, and calculates the readings of the regulatory elements in each cell of the object to be detected 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 position data of each promoter, the corresponding target gene and the readings in each cell of the object to be detected, and the position data of each enhancer and the readings in each cell of the object to be detected.

[0061] Specifically, the bam format file of scATAC-seq data and the bed format file containing the candidate regulatory elements are used as input, the bedtools tool is used to detect the candidate regulatory elements that are open in the biological process, and the readings (read counts) of the open candidate regulatory elements in a single cell are calculated to generate multiple bed format files containing regulatory elements. Each bed format file represents a cell, and each bed format file contains 7 columns of data, namely the number of the regulatory element (composed of a string of chromosome number, start site, and end site), the number of the chromosome where the regulatory element is located, the start site of the regulatory element on the chromosome where it is located, the end site of the regulatory element on the chromosome where it is located, the type of regulatory element (PLS or ELS), the corresponding target gene (only the target gene corresponding to PLS), and the reading of the regulatory element in the cell. After integrating these bed format files into one bed format file, a column of information about the cell number will be added. The bed format file generated by merging these bed format files representing single cells is called SingleCell_Merged.bed. The SingleCell_Merged.bed file can be converted into a matrix form through the pivot_table() function of Pandas to obtain the reading matrix M, in which the number of the regulatory element is used as the row index, the number of the cell is used as the column index, and the reading is used as the value of the matrix. Therefore, this step actually obtains the precise positioning information of the promoter and enhancer of a specific biological process on the chromosome and the relationship between the promoter and the target gene.

[0062] At this point, in this embodiment, based on the scATAC-seq data and the annotation data of the candidate regulatory elements, the regulatory elements are determined, and the readings of the regulatory elements in each cell of the object to be detected are calculated, specifically including: using the scATAC-seq data and the annotation data of the candidate regulatory elements as input, using the bedtools tool to determine the regulatory elements, and calculating the readings of the regulatory elements in each cell of the object to be detected.

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

[0064] The identified enhancers are picked out separately and output as bed format files. The run_cistarget function of the pycisTarget package in Python language is used to perform motif scanning on the bed format files of the enhancers, thereby searching for transcription factor binding motifs in the motif database. The reference genome version of humans in the motif database is hg38, and the reference genome version of mice is mm10. The results retain four targeting relationships, namely "Direct_annot" (direct evidence that the motif is associated with that transcription factor TF), "Motif_similarity_annot" (based on tomtommotif similarity), "Orthology_annot" (based on the orthogonality of the transcription factor TF directly associated with the motif) or "Motif_similarity_and_Orthology_annot" (based on the motif containing the similarity threshold and the orthologous annotation containing the homology threshold). After scanning and sorting, a bed format file containing the identification results will be obtained, which describes the potential binding sites of each TF on the input enhancer. Therefore, this step actually obtains the binding information of transcription factors to enhancers in the biological process under study.

[0065] In this embodiment, based on the motif database, the position data of each enhancer and the readings in each cell of the object to be detected are subjected to motif scanning to determine the potential binding site of each transcription factor of the object to be detected on the enhancer. Wherein, the position data of each enhancer and the readings in each cell of the object to be detected are subjected to motif scanning to determine the potential binding site of each transcription factor of the object to be detected on the enhancer, specifically comprising: based on the motif database, using the run-cistarget function to perform motif scanning on the position data of each enhancer and the readings in each cell of the object to be detected, and determine the potential binding site of each transcription factor of the object to be detected on the enhancer.

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

[0067] The reading matrix M is further reduced in dimension 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 completed using PyTorch. The reading matrix M is first preprocessed by performing a logarithmic transformation on the matrix elements and then normalizing the transformed data to achieve numerical stability and improve the performance of the model, as follows:

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

[0069]

[0070] Wherein, M' is the matrix after logarithmic transformation; M" is the matrix after preprocessing; μ is the mean of the matrix M' after logarithmic transformation; σ is the standard deviation of the matrix M' after logarithmic transformation.

[0071] Next, we build the encoder, which is responsible for compressing the input data into the latent space. In essence, the structure of the encoder can be expressed 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 by the ReLU activation function; W1, W2, W3 are the weight matrices of the encoder; b1, b2, b3 are the bias terms of the encoder; μ(M") is the mean vector of the latent space; logσ 2 (M”) is the variable logarithmic variance of the latent space. The weight matrix and bias terms are randomly generated during initialization and dynamically adjusted during training through the optimization algorithm to minimize the loss function.

[0076] The latent variable matrix Z is sampled from the latent distribution using the reparameterization trick, as follows:

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

[0078] Among them, ∈ is from the standard normal distribution Noise sampled in is the standard deviation of the latent space.

[0079] Next, we construct a decoder. 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 as follows:

[0080]

[0081] in, is the reconstruction matrix; W d and b d are the weight and bias of the decoder respectively; the sigmoid function is used to limit the output value between 0 and 1. d and b d It is learned dynamically during the training process. The initial value is randomly distributed and is eventually updated to the value that best adapts to the data through the optimization algorithm.

[0082] The decoder is not only used to generate a reconstruction matrix similar to the input data It is also an important mechanism to effectively constrain the latent space during training. The decoder provides an optimization signal for the model, thereby guiding the encoder to compress high-dimensional data into a low-dimensional latent space while minimizing information loss, achieving the dual goals of dimensionality reduction and retaining important information of the data.

[0083] During training, the loss function of the variational autoencoder consists of two parts: reconstruction loss and KL divergence loss, which aims to constrain the distribution of the latent space to a standard normal distribution while maintaining the quality of the reconstruction of the input data. 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" ensures that the data generated by the decoder is as close to the real data as possible. It is expressed in binary cross entropy form as follows:

[0086]

[0087] Among them, n rows and n cols Reconstructed matrices are with 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 potential coding distribution q(Z|M”) and the standard normal distribution p(Z), and is calculated as follows:

[0089]

[0090] Among them, k is the dimension index of the latent variable matrix Z; is the square of the mean of the kth dimension of the latent variable matrix Z; is the variance of the values ​​of the kth dimension of the latent variable matrix Z.

[0091] When the loss function L is below the threshold δ = 10 -2 When , the training is ended, and the encoder is used to output the latent variable matrix Z. The latent variable matrix Z is in matrix form, each row represents the embedded representation of a regulatory element in the latent space, and each column is a dimension of the latent space.

[0092] (4) Calculate the similarity of the open patterns of regulatory elements, corresponding to step S4.

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

[0094]

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

[0096] Each row of the normalized matrix Z' is a valid probability distribution, which is suitable for calculating the Wasserstein distance between regulatory elements. In order to calculate the Wasserstein distance between each pair of regulatory elements, namely Z'[i] and Z'[j], it is necessary to construct a matrix D, where D ij represents the Euclidean distance between row i and row j in the normalized matrix Z', defined as:

[0097]

[0098] Then, the Wasserstein distance can be calculated by the following optimization problem:

[0099]

[0100] Among them, Π(p,q) is the set of all possible transportation plans, p and q are vectors, representing the probability distribution of the i-th and j-th rows of the standardized matrix Z', that is, p i = Z'[i][k],q j = Z'[j][k]; T ij is the transportation volume from row i to row j; the symbol inf represents the “minimum lower bound”, which here means the minimum transportation cost among all possible transportation plans T, that is, the optimal transportation plan. The optimal transportation plan is solved by using the POT package in Python or the optimize module in the SciPy package.

[0101] The Wasserstein distance between any two rows of the normalized matrix Z' is calculated to obtain an n×n distance matrix W, whose row and column indices are regulatory elements and whose values ​​are Wasserstein distances between regulatory elements. Next, the local adaptive Gaussian kernel is used to convert the distance matrix W into the affinity matrix A. Specifically, the kernel bandwidth of each regulatory element is defined by the Wasserstein distance of its kNN, and the affinity between regulatory element i and regulatory element j is defined as:

[0102]

[0103] Among them, A ij is the affinity between regulatory element i and regulatory element j; W ij is the Wasserstein distance between regulatory element i and regulatory element j; σ i and σ j are the local kernel bandwidths of regulatory element i and regulatory element j, respectively. The local kernel bandwidths are determined based on the Wasserstein distance of the k nearest neighbors, and the default value is k=5.

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

[0105]

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

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

[0108] This embodiment uses a reading matrix as input, uses a variational autoencoder to calculate a latent variable matrix, and determines enhancer-promoter pairs with potential interactions based on the latent variable matrix. The reading matrix includes the readings of each regulatory element in each cell of the object to be detected, and the latent variable matrix includes an embedded representation of each regulatory element.

[0109] Among them, the reading matrix is ​​used as input, and the latent variable matrix is ​​calculated using the variational autoencoder, which specifically includes:

[0110] 1) Preprocessing the reading matrix to obtain a preprocessed matrix, wherein the preprocessing includes logarithmic transformation and standardization.

[0111] 2) Take the preprocessed matrix as input and process it using the encoder of the variational autoencoder to obtain the latent variable matrix.

[0112] Among them, the enhancer-promoter pairs with potential interactions are determined based on the latent variable matrix, including:

[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 a distance matrix, where the number of rows and columns of the distance matrix is ​​equal to the total number of regulatory elements, and the element in the i-th row and j-th column of the distance matrix is ​​the Wasserstein distance between the i-th regulatory element and the j-th regulatory element.

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

[0116] 4) applying a distance penalty to the affinity matrix to obtain an adjusted affinity matrix, wherein 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 i-th row and j-th column of the adjusted affinity matrix is ​​the adjusted affinity between the i-th regulatory element and the j-th regulatory element.

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

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

[0119] EnhancerNet constructs a gene co-expression network by calculating the mutual information between genes. The data used in this step is scRNA-seq data. The scRNA-seq data and scATAC-seq data belong to the same biological process. The scRNA-seq data is a matrix X, where X ijrepresents the expression level of the jth gene in the ith cell. The goal is to calculate the mutual information between the expression levels of two genes to construct a co-expression network between genes. The formula for calculating the mutual information is:

[0120]

[0121] Among them, I(G1; G2) is the mutual information between gene G1 and gene G2; and are the value sets of gene G1 and gene G2 expression levels 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, EnhancerNet can use kernel density estimation to obtain its corresponding probability density, and its estimated form is:

[0123]

[0124] Where n is the number of cells; h is the bandwidth parameter, the default h = 1, which controls the degree of smoothing; K is the Gaussian kernel function, and its formula is:

[0125]

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

[0127]

[0128] The formula can be further simplified as:

[0129]

[0130] Among them, x i is the expression level of gene G1 in the ith cell.

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

[0132]

[0133] Among them, h x and h y It is also the bandwidth parameter, the default value is h x =h y =1;y i is the expression level of gene G2 in the ith cell.

[0134] After calculating the mutual information between any two genes, EnhancerNet selects a percentile of the mutual information as the mutual information threshold (the default is 80 percentile). If the mutual information between genes is greater than the mutual information threshold, it is considered that there is a potential co-expression relationship and is used to construct a gene co-expression network, in which the 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 the mutual information between genes based on scRNA-seq data, and constructs a co-expression gene network based on the mutual information between genes, and 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, so the association between the coding genes of transcription factors and target genes can be determined, which is also the association between transcription factors and target genes.

[0136] Among them, the scRNA-seq data includes the expression level of each gene in each cell of the object to be detected, and the mutual information between genes is calculated based on the scRNA-seq data, and the co-expression gene network is constructed based on the mutual information between genes, which specifically includes:

[0137] 1) For each gene, the marginal probability density function of the gene is determined based on the expression level of the gene in each cell of the object to be detected.

[0138] 2) For a gene pair consisting of any two genes, the joint probability density function of the gene pair is determined based on the expression levels of the two genes included in the gene pair in each cell of the object to be detected, and the mutual information of the gene pair is calculated based on the joint probability density function of the gene pair and the marginal probability density function of the two genes included in the gene pair. The mutual information of all gene pairs constitutes the mutual information between genes.

[0139] 3) 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 to obtain a connection relationship between any two genes.

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

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

[0142] By integrating the information from transcription factors to enhancers, enhancers to promoters, and promoters to target genes obtained based on scATAC-seq data, EnhancerNet can construct a basic gene regulatory network that describes the regulatory pathway of transcription factors-enhancers-promoters-target genes, and can be further simplified into the regulatory network of transcription factors regulating target genes. Based on the gene co-expression network obtained based on scRNA-seq data, EnhancerNet can remove redundancy in the regulatory network of transcription factor-regulated target genes by finding the intersection method, and the formula is:

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

[0144] Among them, R' ij is the intersection; R ij is the regulatory relationship between transcription factor i and target gene j inferred based on scATAC-seq data, i.e., the transcription factor-target gene regulatory network; C ij is the co-expression relationship between transcription factor i and target gene j inferred based on scRNA-seq data, i.e., the gene co-expression network. Finally, the regulatory network obtained by EnhancerNet is the "enhancer-driven gene regulatory network", which describes the regulatory network pathway by which specific transcription factors regulate target genes by binding to enhancers in the biological process 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 of the object to be detected on the enhancer, the enhancer-promoter pairs with potential interactions, and the target genes corresponding to each promoter, and simplifies the transcription factor-enhancer-promoter-target gene regulatory network to obtain the transcription factor-target gene regulatory network, calculates the intersection of the transcription factor-target gene regulatory network and the co-expression gene network, and expands the intersection based on the potential binding sites of each transcription factor of the object to be detected 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 the transcription factor-enhancer-promoter-target gene.

[0146] Among them, based on the potential binding sites of each transcription factor of the object to be detected 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 a transcription factor-target gene regulatory network, which specifically includes:

[0147] 1) Based on the potential binding sites of each transcription factor of the object to be detected on the enhancer, the connection relationship between the transcription factor and the enhancer is determined.

[0148] 2) Based on the potential interacting enhancer-promoter pairs, determine the connection relationship 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 relationship between transcription factors and enhancers, the connection relationship between enhancers and promoters, and the connection relationship between promoters and target genes, the 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 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 to determine the connection relationship between the transcription factor and the target gene, 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 the transcription factor-target gene regulatory network.

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

[0153] Example 2

[0154] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Figure 3 As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, referred to as I / O) and a communication interface. Among them, the processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store data. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, a method for constructing a gene regulatory network is implemented.

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

[0156] In an exemplary embodiment, a computer device is provided, including a memory and a processor, wherein a computer program is stored in the memory, and when the processor executes the computer program, the method for constructing a gene regulatory network in Example 1 is implemented.

[0157] Example 3

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

[0159] Example 4

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

[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, stored data, displayed data, 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 relevant data must comply with relevant regulations.

[0162] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, 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 article uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, according to the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A method for constructing a gene regulatory network, characterized in that: The method for constructing the gene regulatory network comprises: Obtaining scATAC-seq data, scRNA-seq data, location data of candidate regulatory elements, and a genome annotation file of an object to be detected; the object to be detected is a tissue of an organism to be detected; the candidate regulatory elements include candidate promoters and candidate enhancers, and the location data of the candidate regulatory elements are downloaded from the SCREEN database, and the location data of the candidate regulatory elements include the number of the chromosome where the candidate regulatory element is located, the starting site of the candidate regulatory element on the chromosome where the candidate regulatory element is located, and the ending site of the candidate regulatory element on the chromosome where the candidate regulatory element is located; Based on the position data of the candidate regulatory element and the genome annotation file, the candidate promoter is annotated, the target gene corresponding to the candidate promoter is determined, and the annotation data of the candidate regulatory element is obtained; the annotation data of the candidate regulatory element includes the position data of each candidate promoter and the corresponding target gene and the position data of each candidate enhancer; Based on the scATAC-seq data and the annotation data of the candidate regulatory elements, the regulatory elements are determined, and the readings of the regulatory elements in each cell of the object to be detected are calculated to obtain the annotation data of the regulatory elements; the regulatory elements are open candidate regulatory elements, the regulatory elements include promoters and enhancers, and the annotation data of the regulatory elements include the position data of each promoter, the corresponding target gene and the readings in each cell of the object to be detected, and the position data of each enhancer and the readings in each cell of the object to be detected; Taking the read matrix as input, using the variational autoencoder to calculate the latent variable matrix, and determining the enhancer-promoter pairs with potential interactions based on the latent variable matrix; the read matrix includes the reads of each regulatory element in each cell of the object to be detected; the latent variable matrix includes the embedded representation of each regulatory element; Based on the motif database, the position data of each enhancer and the readings in each cell of the object to be detected are subjected to motif scanning to determine the potential binding site of each transcription factor of the object to be detected on the enhancer; Based on the potential binding sites of each transcription factor of the object to be detected 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 a transcription factor-target gene regulatory network; The mutual information between genes is calculated based on the scRNA-seq data, and a co-expression gene network is constructed based on the mutual information between genes; the genes include coding genes and target genes of transcription factors; 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 object to be detected 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.

2. The method for constructing a gene regulatory network according to claim 1, characterized in that: Based on the position data of the candidate regulatory element and the genome annotation file, annotating the candidate promoter and determining the target gene corresponding to the candidate promoter specifically includes: Taking the position data of the candidate regulatory element and the genome annotation file as input, annotating the candidate promoter using the annotatePeak function to determine the target gene corresponding to the candidate promoter; Based on the scATAC-seq data and the annotation data of the candidate regulatory elements, the regulatory elements are determined, and the reading of the regulatory elements in each cell of the object to be detected is calculated, specifically including: Taking scATAC-seq data and annotation data of candidate regulatory elements as input, bedtools tool was used to identify regulatory elements and calculate the read count of the regulatory elements in each cell of the subject to be detected.

3. The method for constructing a gene regulatory network according to claim 1, characterized in that: Taking the reading matrix as input, the variational autoencoder is used to calculate the latent variable matrix, which includes: Preprocessing the reading matrix to obtain a preprocessed matrix; the preprocessing includes logarithmic transformation and standardization; The preprocessed matrix is ​​used as input and processed by the encoder of the variational autoencoder to obtain a latent variable matrix.

4. The method for constructing a gene regulatory network according to claim 1, characterized in that: Determining enhancer-promoter pairs with potential interactions based on the latent variable matrix specifically includes: Standardize each row of the latent variable matrix so that the sum of the elements in each row is 1, and 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 i-th row and j-th column of the distance matrix is ​​the Wasserstein distance between the i-th regulatory element and the j-th 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 i-th row and j-th column of the affinity matrix is ​​the affinity between the i-th regulatory element and the j-th regulatory element; Applying a distance penalty 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 i-th row and j-th column of the adjusted affinity matrix is ​​the adjusted affinity between the i-th regulatory element and the j-th regulatory element; For each element in the adjusted affinity matrix, determine whether the element is greater than a preset affinity threshold; if so, the two regulatory elements corresponding to the element are regarded as a pair of regulatory elements with potential interaction, and if the pair of regulatory elements with potential interaction includes a promoter and an enhancer, the pair of regulatory elements with potential interaction is regarded as an enhancer-promoter pair with potential interaction.

5. The method for constructing a gene regulatory network according to claim 1, characterized in that: Based on the motif database, the position data of each enhancer and the readings in each cell of the object to be detected are scanned for motifs to determine the potential binding sites of each transcription factor on the enhancer of the object to be detected, specifically including: Based on the motif database, the run-cistarget function is used to perform motif scanning on the position data of each enhancer and the readings in each cell of the object to be detected to determine the potential binding sites of each transcription factor of the object to be detected on the enhancer.

6. The method for constructing a gene regulatory network according to claim 1, characterized in that: Based on the potential binding sites of each transcription factor of the object to be detected 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 a transcription factor-target gene regulatory network, which specifically includes: Based on the potential binding site of each transcription factor of the object to be detected on the enhancer, determining the connection relationship between the transcription factor and the enhancer; Based on the enhancer-promoter pairs with potential interactions, determine the connection relationship between enhancers and promoters; Based on the target gene corresponding to each promoter, determining the connection relationship between the promoter and the target gene; Based on the connection relationship between transcription factors and enhancers, the connection relationship between enhancers and promoters, and the connection relationship between promoters and target genes, the transcription factors, enhancers, promoters and target genes 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 to determine the connection relationship between the transcription factor and the target gene, 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 for constructing a gene regulatory network according to claim 1, characterized in that: scRNA-seq data includes the expression level of each gene in each cell of the object to be detected. Based on the scRNA-seq data, the mutual information between genes is calculated, and a co-expression gene network is constructed based on the mutual information between genes, which includes: For each gene, determining a marginal probability density function of the gene based on the expression level of the gene in each cell of the subject to be detected; For a gene pair consisting of any two genes, based on the expression levels of the two genes included in the gene pair in each cell of the object to be detected, the joint probability density function of the gene pair is determined, and based on the joint probability density function of the gene pair and the marginal probability density functions of the two genes included in the gene pair, the mutual information of the gene pair is calculated; the mutual information of all the gene pairs constitutes the 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 to obtain a connection relationship between any two genes; Each gene is taken as a node, and the nodes are connected based on the connection relationship between any two genes to obtain a co-expressed gene network.

8. A computer device comprising: 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 method for constructing a gene regulatory network according to any one of claims 1 to 7.

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

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

Images