Single-cell Hi-C data enhancement method and system and storage medium

Abstract

The invention relates to a bioinformatics technology, in particular to a single-cell Hi-C data enhancement method and system and a storage medium. The method comprises the following steps: preprocessing multi-resolution Hi-C data, extracting Hi-C contact matrixes under different resolutions, and constructing a graph structure to obtain graph structure data under different resolutions; constructing a variational graph auto-encoder VGAE model as a fusion unit model, inputting the graph structure data into the fusion unit model, carrying out feature extraction and matrix reconstruction, realizing cross-scale fusion of the graph structure data under different resolutions, obtaining fusion potential spatial features under different resolutions, and carrying out fusion on the fusion potential spatial features under different resolutions. And reconstructing the contact matrix as an enhanced Hi-C contact matrix. According to the method, information of low-resolution data and information of high-resolution data are effectively fused, and the problem of sparsity of single-cell Hi-C data is solved.

Description

Technical Field

[0001] The present invention relates to bioinformatics technology, and in particular to a single-cell Hi-C data enhancement method, system and storage medium. Background Art

[0002] Genomic DNA is not a simple linear state in the cell nucleus, but exists in the form of chromatin that is highly folded and condensed with a specific high-level spatial conformation and is stored in the cell nucleus as a carrier of genetic information. Genomics is the collective characterization and quantification of all genes of an organism, studying the structure, function, evolution, localization, and other issues of the genome, and analyzing their interrelationships and their impact on organisms. A large number of studies have shown that the spatial structure of chromosomes is closely related to most biological processes in cells, including genome regulation, stable genome expression, DNA replication, and chromosome translocation; some studies have also shown that the 3D structure of chromatin is also related to disease, and there are also studies showing that the insulating properties of the boundaries of TADs (Topologically Associating Domains) have a significant impact on the epigenetic landscape of single cells.

[0003] Since most existing 3D genome modeling methods focus on constructing chromosome / genomic structures of cell populations derived from bulk Hi-C (High-throughput Chromosome Conformation Capture) data (usually millions of cells). The reconstructed structure is the average structure of a cell population or subpopulation, which cannot capture the variability of individual cells well. Although single cells of the same type usually have some common structural properties, different single cells differ in the interactions between topologically related domains TADs, chromosome loops within TADs, and the distribution of chromosome compartments. These structures (such as TADs, chromosome loops, chromosome compartments) play an important role in gene transcription, histone methylation, cell cycle, and cell development.

[0004] The development of single-cell technologies for probing 3D genome conformation, especially the development of single-cell chromosome conformation capture technology, has enabled us to understand genome function better than before; however, due to the extreme sparsity and high noise of single-cell Hi-C data (such as graph data), it is still difficult to study genome structure and function using Hi-C data of single cells. It is very important to use single-cell Hi-C data to analyze and model 3D genome / chromosome structure to study cell variability. Therefore, it becomes very necessary to improve the quality of single-cell Hi-C data, and it is crucial to develop a computational method to enhance contact maps.

[0005] In recent years, some studies have applied deep learning technology to 3D genomics, including TADs recognition, Loops (ring-shaped interaction structures of DNA or chromatin in three-dimensional space, i.e., chromatin loops) recognition, and contact matrix reconstruction. Researchers can reveal the connection between structure and function at the single cell level that is ignored by Bulk Hi-C through the dynamic chromosome conformation at the single cell level.

[0006] Based on single-cell Hi-C data, it is necessary to reveal the structural characteristics of single-cell 3D genomes through effective computational methods. For single-cell Hi-C data, the classic analysis process is to first read the cell barcodes through sequencing, compare them with the reference genome, and merge them into the single-cell Hi-C contact matrix; secondly, pre-process the data, and then use effective computational methods to perform dimensionality reduction, imputation, embedding, and clustering; finally, perform downstream analysis on the processed data to reveal the heterogeneity and dynamics of the three-dimensional genome. Existing technologies can embed single-cell Hi-C data into low-dimensional space, thereby capturing biologically significant changes in 3D chromatin structure, but it is difficult to cluster cells.

[0007] Although existing methods have made some progress in improving the resolution and analysis performance of Hi-C data, they still have the following limitations: 1. Single-resolution enhancement: Most current methods focus on the enhancement of single-resolution data, and it is difficult to simultaneously utilize multi-resolution feature information, resulting in limited enhancement effects at high resolutions (such as 10kb).

[0008] 2. Insufficient adaptability: Some methods, such as HiCPlus and HiCNN, were originally designed to process bulk Hi-C data and have poor adaptability to sparse single-cell Hi-C data.

[0009] 3. Insufficient feature extraction: Existing feature extraction methods such as HiCRep / MDS make it difficult to perform cell clustering and multidimensional analysis at the single-cell level.

[0010] 4. Data noise and sparsity challenges: Single-cell Hi-C data is highly sparse and noisy. Most methods are not ideal when dealing with data interpolation, especially when detecting key structures such as TADs and Loops.

[0011] In short, the contact matrix in the single-cell Hi-C data of existing technologies shows significant sparsity, which limits the in-depth understanding of the three-dimensional genome structure. These deficiencies restrict the application of existing methods in high-resolution single-cell Hi-C data enhancement and downstream analysis, so it is very necessary to develop multi-scale graph neural network models to enhance high-resolution single-cell Hi-C data. Summary of the invention

[0012] In order to solve the sparsity problem of existing single-cell Hi-C data, the present invention proposes a single-cell Hi-C data enhancement method, system and storage medium. By designing a multi-scale variational graph autoencoder (VGAE) model, the information of low-resolution and high-resolution data is effectively fused, thereby solving the sparsity problem of single-cell Hi-C data and enhancing the integrity and quality of the single-cell Hi-C contact matrix.

[0013] In one aspect, the present invention provides a single cell Hi-C data enhancement method, comprising the following steps: S1. Preprocess the multi-resolution Hi-C data, extract the Hi-C contact matrix of each cell at different resolutions, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions; S2. Construct a variational graph autoencoder VGAE model as the fusion unit model, input graph structure data at different resolutions into the fusion unit model, perform feature extraction and matrix reconstruction, realize cross-scale fusion of graph structure data at different resolutions, obtain fusion latent space features at different resolutions, and reconstruct the contact matrix as the enhanced Hi-C contact matrix.

[0014] Preferably, the above enhancement method further comprises the steps of: S3. Identify domain boundaries in Hi-C data using the enhanced Hi-C contact matrix.

[0015] In another aspect, the present invention provides a single cell Hi-C data enhancement system for implementing the above enhancement method, wherein the enhancement system comprises the following modules: A preprocessing module is used to preprocess multi-resolution Hi-C data, extract the Hi-C contact matrix at different resolutions for each cell, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions; The fusion unit model construction module constructs a variational graph autoencoder VGAE model as the fusion unit model, inputs graph structure data at different resolutions into the fusion unit model, performs feature extraction and matrix reconstruction, realizes cross-scale fusion of graph structure data at different resolutions, obtains fusion latent space features at different resolutions, and reconstructs the contact matrix as the enhanced Hi-C contact matrix.

[0016] In yet another aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed on a computer, the computer is enabled to execute the above-mentioned enhancement method.

[0017] Compared with the prior art, the technical effects achieved by the present invention include: 1. Multi-resolution fusion strategy: This invention combines single-cell Hi-C data of different resolutions (e.g., 100kb, 50kb, 10kb) for the first time, and effectively fuses the information of low-resolution and high-resolution data by designing a multi-scale variational graph autoencoder VGAE model, thereby solving the sparsity problem of single-cell Hi-C data. This multi-resolution fusion method not only enhances the integrity and quality of the single-cell Hi-C contact matrix, but also provides richer genomic structural information at a higher resolution, providing a more reliable data basis for the subsequent identification of TADs and Loops domains.

[0018] 2. Enhanced data quality: The model is designed based on graph neural network, and the variational graph autoencoder VGAE is used as the core architecture. The graph structure characteristics of Hi-C data are fully utilized, the contact matrix is ​​represented as an adjacency matrix of the graph structure, and a fusion unit model is created for feature extraction and matrix reconstruction. The fusion unit model can use data features at different scales to perform cross-scale feature fusion of low-resolution and high-resolution graph structure data, enhance the overall quality of the Hi-C contact matrix, realize information transmission and sharing between features of different resolutions, and improve the prediction performance of the model and data enhancement effect; multi-resolution data fusion can effectively fill in the missing parts of sparse data, reduce noise, and improve the ability of downstream analysis to parse genome structure.

[0019] 3. Efficient identification of TADs and Loops structures: Based on the cross-domain analysis results of multi-scale features and through the enhanced contact matrix, the present invention can simultaneously identify three-dimensional genome features (such as TADs and Loops) at multiple resolutions at the single-cell level. Especially at the single-cell level, by accurately locating the boundaries of TADs and Loops structures, it can help researchers understand the spatial distribution of genomic functional regions and their changes in different cellular states, which is of great significance for studying the spatial organization of chromosomes and gene regulation.

[0020] 4. Scalability and adaptability of the algorithm: The algorithm framework designed in the present invention is suitable for sparse and highly heterogeneous single-cell Hi-C data, and can also be extended to other similar high-dimensional graph structure data, with good adaptability and versatility. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 is a flow chart of a single-cell Hi-C data enhancement method based on multi-resolution and multi-scale in an embodiment of the present invention; Figure 2 4 is a block diagram of a single-cell Hi-C data enhancement system based on multi-resolution and multi-scale in an embodiment of the present invention. DETAILED DESCRIPTION

[0022] The technical solution of the present invention is further described in detail below in conjunction with embodiments and drawings, but the implementation manner of the present invention is not limited thereto.

[0023] Example This embodiment is based on a multi-resolution and multi-scale single-cell Hi-C data enhancement method. The core idea is to enhance single-cell Hi-C data at multiple resolutions, especially to achieve effective feature information fusion between low-resolution and high-resolution images.

[0024] In order to ensure the efficient operation of the model, the data input part adopts an optimized data processing solution, using the straw package to efficiently extract single-cell Hi-C data and convert it into a format suitable for graph neural network processing, so that data from multiple single cells can be quickly processed and efficient training and analysis can be performed.

[0025] In order to effectively utilize single-cell Hi-C data at multiple resolutions (such as 100kb, 50kb, and 10kb), the data must first be preprocessed and then enhanced through a graph neural network (GNN). Through preprocessing, the data is converted into a format suitable for input into a fusion model, and then the fusion model is used for learning and reconstruction to further enhance the data quality.

[0026] like Figure 1 As shown, the single cell Hi-C data enhancement method of this embodiment specifically includes the following steps: S1. Preprocess the multi-resolution Hi-C data, extract the Hi-C contact matrix of each cell at different resolutions, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions.

[0027] This step specifically includes: S11. Preprocess the single-cell Hi-C data at different resolutions and extract the Hi-C contact matrix of each cell at different resolutions.

[0028] Preprocessing mainly includes normalization and sparse matrix filling. Normalization ensures that data of different resolutions are numerically comparable, while sparse matrix filling aims to supplement the missing parts in the contact matrix to avoid the negative impact of sparsity on model training.

[0029] In this embodiment, multi-resolution Hi-C data are first converted into single-cell Hi-C contact matrices of different resolutions, and then the single-cell Hi-C contact matrices of different resolutions are normalized. The Hi-C contact matrix at each resolution is mainly composed of chromosome fragments (i.e., segments of genomic regions), and each matrix element in the contact matrix represents the contact frequency between corresponding chromosome fragments; after normalization, the single-cell Hi-C contact matrix is ​​stored as a sparse matrix for subsequent graph neural network processing.

[0030] The purpose of normalization is to make Hi-C contact matrices at different resolutions comparable. Assume that a contact matrix M is given, whose matrix elements represent the contact frequencies between chromosome fragments. The normalization operation can be achieved by the following formula: ; in The contact matrix i Row, No. j The contact frequency of the column, is the mean of all matrix elements in the contact matrix M, is the standard deviation of all elements in the contact matrix M, Represents the normalization of the contact matrix.

[0031] Through normalization, the values ​​of the contact matrix will become a distribution with a mean of 0 and a standard deviation of 1, thereby improving the comparability between data of different resolutions.

[0032] In this embodiment, the data source of the Hi-C contact matrix is ​​usually a public database, such as a mouse data set in the GEO database. This data set contains single-cell Hi-C data of various tissue types. The data enhancement method of this embodiment can be applied to other related single-cell Hi-C data sets to test its applicability and generalization ability in different cell types. In the research process of this embodiment, 1954 files with the suffix .hic were downloaded, which contained Hi-C data from different single cells; these .hic files were generated into a list using the straw package, and code was written to read the Hi-C contact matrix element values ​​of the 19 chromosomes of each cell, and finally 37126 contact matrices were obtained. The data set acquisition process can be summarized as follows: download Hi-C data from a public database; use the straw package to read the .hic file and extract the Hi-C contact matrix of each cell; use the extracted contact matrix data for subsequent model training and analysis.

[0033] S12. Construct the Hi-C contact matrix at each resolution into a graph structure.

[0034] The preprocessed Hi-C contact matrix is ​​constructed into graph structures of different scales, where the nodes of the graph structure represent the chromosome segments of the genome (such as the segmentation of the genomic region), and the edges of the graph structure represent the contact frequency between different chromosome segments. The node features use the features of topologically associated domains (TADs) and chromatin compartments (A / B compartment) at a resolution of 1Mb.

[0035] The Hi-C contact matrix can be converted into a graph structure, where the rows and columns of the matrix correspond to the nodes of the graph structure, and the matrix elements correspond to the edge weights (contact frequencies) between the nodes. For the normalized contact matrix M at each resolution, a graph structure G can be constructed, where the node set of the graph structure is the chromosome segment set V of the genome, and the edge set of the graph structure is the contact frequency between chromosome segments. The adjacency matrix A of the graph structure is determined by the elements in the contact matrix: A(i,j)= M(i,j) ; A(i,j) Represents a node in a graph structure i and nodes j So for each contact matrix M, the index with a value is the edge set of the adjacency matrix A. , the value is the edge weight corresponding to the adjacency matrix A .

[0036] The Hi-C contact matrix at each resolution is considered as a graph structure. The rows and columns in the contact matrix represent different regions of the genome, and the value of each matrix element represents the contact frequency between chromosome segments.

[0037] Through preprocessing, graph structure data of different resolutions such as 10kb, 50kb, and 100kb can be obtained, thereby converting the contact matrix into graph structure data so that the subsequent model can better capture the complex spatial relationship between genomes.

[0038] S2. Construct a variational graph autoencoder VGAE model as the fusion unit model, input graph structure data at different resolutions into the fusion unit model, perform feature extraction and matrix reconstruction, realize cross-scale fusion of graph structure data at different resolutions, obtain fusion latent space features at different resolutions, and reconstruct the contact matrix as the enhanced Hi-C contact matrix.

[0039] This step builds a variational graph autoencoder VGAE model based on graph neural network, inputs graph structure data at different resolutions into the constructed variational graph autoencoder VGAE model, and obtains low-dimensional space representations of different resolutions. The nodes of the graph neural network represent genomic regions, and the edges represent the interactions between these genomic regions (i.e., contact frequency). Through graph convolution operations, graph neural networks can effectively capture the complex relationships between nodes in high-dimensional space and better capture the topological features of three-dimensional genomic data, which is particularly important for processing multi-level spatial structures in Hi-C data.

[0040] The VGAE model adopts an encoder-decoder structure. The encoder is used to extract node features from graph structures of different resolutions, while the decoder reconstructs the contact matrix through the inner product operation of the matrix to learn the genome structure information at different scales.

[0041] In this step, the VGAE model can capture the structural features of the genome at different scales by processing graph structure data of different resolutions. Specifically, in the encoding stage, the graph structure data of multiple resolutions is input simultaneously, and the VGAE model can learn the node features of the graph structure data at each scale, and extract the potential structural information through graph convolution operations; then the graph structure node features at different resolutions are transformed in dimension, weighted fused with appropriate weights, and the graph structure data of multiple resolutions are fused. The fusion of data at multiple resolutions can not only effectively improve the overall quality of the contact matrix, but also enhance the biological significance of the data. In the decoding stage, the VGAE model reconstructs the contact matrix through the inner product operation of the matrix to generate a more accurate Hi-C contact matrix, thereby better revealing the three-dimensional structure and spatial characteristics of the genome.

[0042] This step specifically includes: S21. Design the variational graph autoencoder VGAE model as the fusion unit model.

[0043] Taking advantage of the variational graph autoencoder (VGAE) in processing sparse matrix graph structure data, especially its outstanding performance in link prediction and feature learning, the variational graph autoencoder (VGAE) based on graph neural network was selected to build the fusion unit model.

[0044] S22. Set the encoder of the VGAE model, and input multi-resolution graph structure data at the same time so that the VGAE model can learn the graph structure node features at each resolution, and extract the potential structural information through graph convolution operations; then transform the dimensions of the graph structure node features at different resolutions, enhance the fusion of the graph structure node features at multiple resolutions, and output the fused latent space features at different resolutions.

[0045] Specifically include: S221. Setting a plurality of convolutional layers and a graph convolutional network (GCNConv) in the encoder to perform convolution processing on the graph structure data at each resolution, and obtaining node features after the single-scale convolution operation on the graph structure data at each resolution.

[0046] First, for the graph data at each resolution (10kb, 50kb, 100kb), let the adjacency matrix of the graph structure be A reso , the input features of each graph structure are X reso , generate hidden layer features through graph convolution operation H reso : H reso = GCNConv(X reso ,A reso ); Hidden layer features H reso As node features after single-scale convolution operation.

[0047] S222. Set up multiple attention mechanism modules, including the first attention mechanism module (S1AttInform), the second attention mechanism module (S2AttInform), and the third attention mechanism module (S3AttInform), which are used to process the node features after single-scale convolution operations at three different resolutions of 10kb, 50kb, and 100kb respectively. H reso , thus obtaining the node potential features of the fusion of nodes and edges at different resolutions.

[0048] (1) First, the node features after the single-scale convolution operation are processed by a multi-layer perceptron (MLP) to obtain the node features after multi-layer perceptron processing.

[0049] MLP is a multi-layer fully connected network that processes the node features after single-scale convolution operations through a series of linear transformations and activation functions. H reso Input the multi-layer perceptron and perform multi-layer perceptual processing, which can be expressed as: x node =MLP ( H reso ); in, x node It is the node feature after multi-layer perception processing.

[0050] (2) Then, the node features after multi-layer perception processing are mapped through the node-edge conversion (node2edge) module to generate edge features of the graph structure.

[0051] In the process of node-edge conversion, the node features after multi-layer perception processing are first mapped to the edges, the relationship between the nodes in the graph structure is calculated, and the edge features are generated by calculating the difference between the receiving node and the sending node.

[0052] Specifically, given a node feature matrix , the receiving relationship matrix rel_rec between nodes and the sending relationship matrix rel_send between nodes, we have: ; ; ; in, It is the node feature matrix formed by the node features after multi-layer perception processing, and its dimension size is , N is the number of nodes, D is the feature dimension; rel_rec is the receiving relationship matrix, which describes which nodes are "receiving" nodes; rel_send is the sending relationship matrix, which describes which nodes are "sending" nodes; receivers is the receiving node feature calculated by the receiving relationship matrix and the node feature matrix; senders is the sending node feature calculated by the sending relationship matrix and the node feature matrix; distance is the difference between the receiving node feature and the sending node feature. Finally, the edge feature edges can be expressed as a new matrix obtained by horizontally splicing the receiving node feature receivers and the sending node feature senders: ; By combining the differences between the receiving node features and the sending node features, graph-structured edge features are generated.

[0053] (3) Next, the edge features edges are processed by a multi-layer perceptron (MLP) to obtain the edge features after multi-layer perceptron processing: ; in, It is the edge feature after multi-layer perception processing.

[0054] (4) Subsequently, the edge features after multi-layer perception processing are mapped through the edge-to-node conversion module to obtain the updated node features.

[0055] The specific operation of edge node conversion is to average the edge features after multi-layer perception processing to obtain new node features, which can help node features aggregate neighbor information. The receiving relationship matrix rel_rec between the nodes is used to map the edges to the nodes as follows: ; in, is the edge feature after multi-layer perception processing, and the dimension size is , N is the number of nodes, is the dimension of edge features; rel_rec is the receiving relationship matrix, rel_rec T is the transpose of the receiving relationship matrix; incoming is the incoming edge feature of each node. All incoming edge features are averaged to obtain the updated node feature for: ; Among them, incoming_size is the number of incoming edges of each node, which is used to average the incoming edge features. This mapping operation generates updated node features, which represent the information obtained by the node through the aggregation of edge features of its neighbors.

[0056] (5) Next, the node features after the single-scale convolution operation are H reso Node features updated by edge node conversion module Splice along the last dimension to get the spliced ​​node features : ; (6) Finally, the concatenated node features are processed by a multi-layer perceptron (MLP). Processing is performed to obtain the node potential features of the nodes and edges after multi-layer perception processing. : ; S223, setting multiple cross-scale fusion modules, performing weighted fusion processing on the node potential features fused by the nodes and edges at different resolutions, and obtaining fused potential space features.

[0057] The cross-scale fusion modules provided in this embodiment include a fusion module S2_to_S1 from the second scale to the first scale, a fusion module S3_to_S1 from the third scale to the first scale, a fusion module S1_to_S2 from the first scale to the second scale, a fusion module S3_to_S2 from the third scale to the second scale, a fusion module S1_to_S3 from the first scale to the third scale, and a fusion module S2_to_S3 from the second scale to the third scale.

[0058] Each fusion module is a cross-scale information fusion module. The fusion module S2_to_S1 fuses the node potential features from the second resolution (e.g., 50kb) with the node potential features from the first resolution (e.g., 10kb). The fusion module S3_to_S1 fuses the node potential features from the third resolution (e.g., 100kb) with the node potential features from the first resolution. The fusion module S1_to_S2 fuses the node potential features from the first resolution with the node potential features from the second resolution. The fusion module S3_to_S2 fuses the node potential features from the third resolution with the node potential features from the second resolution. The fusion module S1_to_S3 fuses the node potential features from the first resolution with the node potential features from the third resolution. The fusion module S2_to_S3 fuses the node potential features from the third resolution with the node potential features from the second resolution.

[0059] The fusion process includes the self-attention calculation of node features and the weighted fusion of features. The specific fusion process includes the following steps: (1) The node features after single-scale convolution operation at the first resolution (10kb), the second resolution (50kb) and the third resolution (100kb) are processed by the first attention mechanism module, the second attention mechanism module and the third attention mechanism module respectively, and the node latent feature attention representation of the node and the edge fused at the first resolution, the node latent feature attention representation of the node and the edge fused at the second resolution and the node latent feature attention representation of the node and the edge fused at the third resolution are obtained.

[0060] For the node features after the single-scale convolution operation at the first resolution , the node potential feature attention representation of the node and edge fusion at the first resolution is obtained through the first attention mechanism module : ; in, represents the receiving relationship matrix between nodes at the first resolution, Represents the sending relationship matrix between nodes at the first resolution.

[0061] For the node features after the single-scale convolution operation at the second resolution , the node potential feature attention representation of the node and edge fusion at the second resolution is obtained through the second attention mechanism module : ; in, represents the receiving relationship matrix between nodes at the second resolution, Represents the sending relationship matrix between nodes at the second resolution.

[0062] For the node features after the single-scale convolution operation at the third resolution , the node potential feature attention representation of the fusion of nodes and edges at the third resolution is obtained through the third attention mechanism module : ; in, represents the receiving relationship matrix between nodes at the third resolution, Represents the sending relationship matrix between nodes at the third resolution.

[0063] (2) By calculating the inner product between the node potential feature attention representations of nodes and edges at different resolutions, the self-attention weight is calculated to obtain the attention weight matrix of each scale fusion, including the attention weight matrix from the first resolution to the second resolution , the attention weight matrix from the first resolution to the third resolution , the attention weight matrix from the second resolution to the first resolution , the attention weight matrix from the second resolution to the third resolution , the attention weight matrix from the third resolution to the first resolution , the attention weight matrix from the third resolution to the second resolution .

[0064] Calculate the inner product between the node potential feature attention representation of the fusion of nodes and edges at the first resolution and the node potential feature attention representation of the fusion of nodes and edges at the second resolution, and use softmax normalization to obtain the attention weight matrix from the first resolution to the second resolution , the attention weight matrix from the second resolution to the first resolution : ; ; Calculate the inner product between the node potential feature attention representation of the fusion of nodes and edges at the first resolution and the node potential feature attention representation of the fusion of nodes and edges at the third resolution, and use softmax normalization to obtain the attention weight matrix from the first resolution to the third resolution , the attention weight matrix from the third resolution to the first resolution : ; ; Calculate the inner product between the node potential feature attention representation of the fusion of nodes and edges at the third resolution and the node potential feature attention representation of the fusion of nodes and edges at the second resolution, and use softmax normalization to obtain the attention weight matrix from the second resolution to the third resolution , the attention weight matrix from the third resolution to the second resolution : ; ; (3) Two-scale feature fusion: Each cross-scale fusion module fuses the attention representation of the node potential features at the corresponding resolution according to the corresponding attention weight matrix to obtain the weighted fusion features of the node potential features at the corresponding resolution.

[0065] The fusion module from the first scale to the second scale uses the attention weight matrix from the first resolution to the second resolution The node potential feature attention representation with the first resolution of 10kb and the node potential feature attention representation with the second resolution of 50kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature of the first resolution. Finally, the weighted fusion feature of the node potential feature of the first resolution is output through the fusion module from the first scale to the second scale, which contains the shared information from the first resolution 10kb and the second resolution 50kb.

[0066] The fusion module from the first scale to the third scale uses the attention weight matrix from the first resolution to the third resolution The node potential feature attention representation of the first resolution 10kb and the node potential feature attention representation of the third resolution 100kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature of the first resolution. Finally, the weighted fusion feature of the node potential feature of the first resolution is output through the fusion module from the first scale to the third scale, which contains the shared information from the first resolution 10kb and the third resolution 100kb.

[0067] The fusion module from the second scale to the first scale uses the attention weight matrix from the second resolution to the first resolution The node potential feature attention representation with the first resolution of 10kb and the node potential feature attention representation with the second resolution of 50kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature of the second resolution. Finally, the weighted fusion feature of the node potential feature of the second resolution is output through the fusion module from the second scale to the first scale, which contains the shared information from the first resolution 10kb and the second resolution 50kb.

[0068] The fusion module from the second scale to the third scale uses the attention weight matrix from the second resolution to the third resolution The node potential feature attention representation at the third resolution of 100kb and the node potential feature attention representation at the second resolution of 50kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature of the second resolution. Finally, the weighted fusion feature of the node potential feature of the third resolution is output through the fusion module from the second scale to the third scale, which contains the shared information from the third resolution 100kb and the second resolution 50kb.

[0069] The third-scale to first-scale fusion module uses the attention weight matrix from the third resolution to the first resolution The node potential feature attention representation of the first resolution 10kb and the node potential feature attention representation of the third resolution 100kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature of the third resolution. Finally, the weighted fusion feature of the node potential feature of the third resolution is output through the fusion module from the third scale to the first scale, which contains the shared information from the first resolution 10kb and the third resolution 100kb.

[0070] The third-scale to second-scale fusion module uses the attention weight matrix from the third resolution to the second resolution The node potential feature attention representation at the third resolution of 100kb and the node potential feature attention representation at the second resolution of 50kb are weighted fused: ; in, It is the weighted fusion feature of the node potential feature at the third resolution. Finally, the weighted fusion feature of the node potential feature at the third resolution is output through the fusion module from the third scale to the second scale, which contains the shared information from the third resolution 100kb and the second resolution 50kb.

[0071] (4) Three-scale feature fusion.

[0072] In the three-scale feature fusion stage, based on the results of the two-scale feature fusion, the VGAE model further weighted-fused the weighted fusion features of the node potential features at different resolutions to obtain the fused latent space features at different resolutions. They represent node features after single-scale convolution operations at three different resolutions (10kb, 50kb, 100kb). The weight w is used to control the fusion ratio and is usually set to 0.3 to balance the contribution of features at different resolutions. The three-scale feature fusion operation is implemented using the following formula: ; ; ; In the formula represents the fused latent space features at the first resolution of 10kb, represents the fused latent space features at the second resolution of 50kb, Represents the fused latent spatial features at the third resolution of 100kb. Through the fusion of three-scale features, the spatial features at different resolutions are weighted to obtain the final fused latent spatial features at different resolutions. .

[0073] S23. Set the decoder of the VGAE model, calculate the similarity between the nodes of the graph structure through the inner product operation of the matrix according to the fused latent space features at different resolutions output by the encoder, so as to reconstruct the predicted values ​​of the Hi-C contact matrix at different resolutions and obtain the predicted contact matrix.

[0074] In the decoding stage, the fused latent space features at different resolutions output by the encoder are , , It is used to reconstruct the contact matrix. After the inner product operation of the matrix, the similarity between the nodes of the graph structure is calculated, and the predicted values ​​of the Hi-C contact matrix at different resolutions are generated: ; ; ; in, represents the predicted value of the Hi-C contact matrix at the first resolution, represents the predicted value of the Hi-C contact matrix at the second resolution, Represents the predicted value of Hi-C contact matrix at the third resolution; fusion of latent space features , represents the fused node features obtained by the encoder; the Sigmoid function is used to ensure that the value of the contact matrix is ​​in the range of [0,1], indicating the connection probability between nodes.

[0075] The inner product operation of the matrix is ​​used to reconstruct the original contact matrix by mapping the input fused latent space features through the low-dimensional space to calculate the similarity between node features (i.e., contact frequency). After processing by the Sigmoid function, a connection probability between nodes is obtained as the final Hi-C contact matrix prediction value, thereby ensuring that the matrix reconstruction result can maintain the spatial structural characteristics of the graph structure data.

[0076] S24. Based on the reconstruction loss and KL divergence loss, construct the loss function of the VGAE model.

[0077] The loss function mainly consists of two parts: one is the reconstruction loss, which is used to measure the difference between the original contact matrix and the predicted contact matrix; ; in is the reconstruction loss, N represents the number of samples, i.e., the total number of elements in the contact matrix; is the element of the original Hi-C contact matrix, indicating the contact frequency between node i and node j; The contact probability value of the Hi-C contact matrix predicted by the model indicates the predicted contact probability between node i and node j. The contact probability value ranges from 0 to 1.

[0078] The other part is the KL divergence loss, which is used to measure the difference between the latent space feature distribution output by the encoder and the standard normal distribution: ; in represents the KL divergence loss, , are the mean and variance of the fused latent space features output by the encoder, respectively.

[0079] The two parts of loss are weightedly added to obtain the final loss function, which is used to optimize the training process of the VGAE model (i.e., the fusion unit model).

[0080] S25. Use the trained fusion unit model to perform feature extraction to obtain the enhanced Hi-C contact matrix.

[0081] After the model training is completed, the original Hi-C data is input into the trained model to obtain a new contact matrix after enhancing each graph structure. Then, the edge information and edge weights in the new contact matrix are extracted, and the new contact matrix is ​​saved as a Hi-C contact matrix format through the pre tool in juicertools. This process provides a basis for the subsequent identification of TADs and Loops domain boundaries, enabling further analysis of the spatial characteristics of the genomic structure.

[0082] S3. Identification of domain boundaries.

[0083] Through the enhanced Hi-C contact matrix, this example further uses arrowhead and hiccups tools in juicertools to accurately identify TADs and Loops domain boundaries in Hi-C data. The identification of these domain boundaries provides important spatial information for studying gene regulation and disease mechanisms.

[0084] Specifically, this embodiment calculates the number of TADs and Loops and the Jaccard index of the boundaries for the identified TADs and Loops to obtain the enhanced contrast effect at the TADs level, and then uses Aggregate Peak Analysis to obtain the corresponding image to represent the enhancement effect of each TAD and Loop in an average sense.

[0085] Based on the same inventive concept, this embodiment also provides a single-cell Hi-C data enhancement system based on multi-resolution and multi-scale, which is used to implement the enhancement method of this embodiment, such as Figure 2 , the enhancement system includes the following modules: A preprocessing module 100 is used to preprocess the multi-resolution Hi-C data, extract the Hi-C contact matrix of each cell at different resolutions, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions; The fusion unit model construction module 200 constructs a variational graph autoencoder VGAE model as the fusion unit model, inputs graph structure data at different resolutions into the fusion unit model, performs feature extraction and matrix reconstruction, realizes cross-scale fusion of graph structure data at different resolutions, obtains fusion latent space features at different resolutions, and reconstructs the contact matrix as the enhanced Hi-C contact matrix.

[0086] Furthermore, the enhanced system of this embodiment may also include an identification module 300 for identifying domain boundaries in Hi-C data through the enhanced Hi-C contact matrix.

[0087] In addition, this embodiment also provides a computer-readable storage medium on which a computer program is stored. When the computer program is run on a computer, the computer executes the enhancement method of this embodiment to achieve enhancement of single-cell Hi-C data.

[0088] Among them, the computer program includes one or more computer instructions. When the computer instructions are loaded and executed on a computer, the process or function described in this embodiment is generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network or other programmable device. The computer instructions can be stored in a computer-readable storage medium, or transferred from one computer-readable storage medium to another computer-readable storage medium. The computer-readable storage medium can be any available medium that can be accessed by a computer, or a data storage device such as a server or a data center that includes one or more available media integrated.

[0089] This embodiment can identify and analyze domains such as TADs and Loops at the single-cell level, which is an innovation compared to traditional batch Hi-C data methods, because traditional batch data usually cannot capture the sparsity and variability of genome structure at the single-cell level. This embodiment can reveal the spatial dynamic characteristics of the genome in a single cell through a multi-resolution graph neural network, especially in different cell types or physiological states, such as changes in TADs and Loops; thereby effectively solving the sparsity problem of single-cell Hi-C data, especially in the identification of TADs and Loops, and providing more detailed and reliable genome structure information through a multi-resolution fusion method, which is helpful for genomic research, especially in providing new technical means for exploring gene regulation and disease mechanisms.

[0090] In summary, the present invention uses a variational graph autoencoder (VGAE) model based on multiple resolution fusion to train the graph structure of single-cell Hi-C data. By introducing the A / B compartments, chromatin loops and topologically related domains TADs information in single-cell Hi-C data, the overall performance of the Hi-C contact matrix is ​​enhanced, so that the Loops and TADs domain boundaries in the contact matrix can be more systematically identified and analyzed. The present invention overcomes the sparsity problem of single-cell Hi-C data through innovative multi-resolution deep learning methods, provides new ideas for three-dimensional genomics research, and effectively enhances the ability to analyze genome structure, especially in the identification of TADs and Loops domains, with significant technical advantages.

[0091] The above embodiment is one of the implementation modes of the present invention, and the implementation modes of the present invention are not limited to the above embodiment. Any other changes, modifications, substitutions, combinations, and simplifications made without departing from the spirit and principle of the present invention shall be equivalent replacement methods and shall be included in the protection scope of the present invention.

Claims

1. A single cell Hi-C data enhancement method, characterized in that: The following steps are involved: S1. Preprocess the multi-resolution Hi-C data, extract the Hi-C contact matrix of each cell at different resolutions, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions; S2. Construct a variational graph autoencoder VGAE model as the fusion unit model, input graph structure data at different resolutions into the fusion unit model, perform feature extraction and matrix reconstruction, realize cross-scale fusion of graph structure data at different resolutions, obtain fusion latent space features at different resolutions, and reconstruct the contact matrix as the enhanced Hi-C contact matrix.

2. The single cell Hi-C data enhancement method according to claim 1, characterized in that: Step S1 includes: S11, preprocessing single-cell Hi-C data at different resolutions to extract Hi-C contact matrices of each cell at different resolutions; the preprocessing includes normalization processing and filling of sparse matrices; S12. Construct the Hi-C contact matrix at each resolution into a graph structure, where the nodes of the graph structure represent the chromosome fragments of the genome, and the edges of the graph structure represent the contact frequencies between different chromosome fragments.

3. The single cell Hi-C data enhancement method according to claim 1, characterized in that: In step S2, the variational graph autoencoder VGAE model adopts an encoder-decoder structure. The encoder is used to extract node features from graph structures of different resolutions, while the decoder reconstructs the contact matrix through the inner product operation of the matrix to learn the genome structure information at different scales.

4. The single cell Hi-C data enhancement method according to claim 1 or 3, characterized in that: Step S2 includes: S21. Design the variational graph autoencoder VGAE model as the fusion unit model. S22, setting the encoder of the VGAE model, inputting multi-resolution graph structure data at the same time so that the VGAE model learns the graph structure node features at each resolution, and extracts the potential structural information through graph convolution operation; then transforming the graph structure node features at different resolutions, enhancing the fusion of the graph structure node features at multiple resolutions, and outputting the fused potential space features at different resolutions; S23, setting a decoder of the VGAE model, calculating the similarity between nodes of the graph structure through an inner product operation of a matrix according to the fused latent space features at different resolutions output by the encoder, so as to reconstruct and generate predicted values ​​of the Hi-C contact matrix at different resolutions, and obtain a predicted contact matrix; S24. Based on the reconstruction loss and the KL divergence loss, a loss function of the VGAE model is constructed to optimize the training process of the fusion unit model; wherein the reconstruction loss is used to measure the difference between the original contact matrix and the predicted contact matrix, and the KL divergence loss is used to measure the difference between the latent space feature distribution output by the encoder and the standard normal distribution; The reconstruction loss is: ； in is the reconstruction loss, N represents the total number of elements in the contact matrix; is the element of the original Hi-C contact matrix, indicating the contact frequency between node i and node j; is the contact probability value of the Hi-C contact matrix predicted by the model, indicating the predicted contact probability between node i and node j; The KL divergence loss is: ； in represents the KL divergence loss, , are the mean and variance of the fused latent space features output by the encoder, respectively; S25. Use the trained fusion unit model to perform feature extraction to obtain the enhanced Hi-C contact matrix.

5. The single cell Hi-C data enhancement method according to claim 4, characterized in that: Step S22 includes: S221, setting a plurality of convolution layers and a graph convolution network in the encoder to perform convolution processing on the graph structure data at each resolution, and obtaining node features after the single-scale convolution operation of the graph structure data at each resolution; S222, setting multiple attention mechanism modules, respectively used to process node features after single-scale convolution operations at different resolutions, so as to obtain node potential features fused with nodes and edges at different resolutions; S223, setting multiple cross-scale fusion modules, performing weighted fusion processing on the node potential features fused by the nodes and edges at different resolutions, and obtaining fused potential space features.

6. The single cell Hi-C data enhancement method according to claim 5, characterized in that: In step S221, for the graph structure data at each resolution, the adjacency matrix of the graph structure is A reso , the input features of each graph structure are X reso , generate hidden layer features through graph convolution operation: < i>H reso = GCNConv(X reso , A reso ); The hidden layer features H reso As node features after single-scale convolution operation; Step S222 includes: (1) The node features after the single-scale convolution operation are processed by a multi-layer perceptron to obtain the node features after multi-layer perceptron processing; (2) The node features after multi-layer perception processing are mapped through the node edge conversion module to generate edge features of the graph structure; In the process of node-edge conversion, the node features after multi-layer perception processing are first mapped to the edges, the relationship between the nodes in the graph structure is calculated, and the edge features are generated by calculating the difference between the receiving node and the sending node; (3) Process the edge features edges through a multi-layer perceptron to obtain edge features after multi-layer perceptron processing; (4) The edge features after multi-layer perception processing are mapped through the edge node conversion module to obtain updated node features; In the process of edge node conversion, given the edge features processed by multi-layer perception And the receiving relationship matrix rel_rec between nodes, perform edge-to-node mapping operation: ； Where rel_rec is the receiving relationship matrix, rel_rec T is the transpose of the receiving relationship matrix; incoming is the incoming edge feature of each node; all incoming edge features are averaged to obtain the updated node features for: ； Among them, incoming_size is the number of incoming edges of each node; (5) Concatenate the node features after the single-scale convolution operation and the node features updated by the edge node conversion module along the last dimension to obtain the concatenated node features; (6) The concatenated node features are processed by a multi-layer perceptron to obtain the node potential features that are a fusion of the nodes and edges after multi-layer perceptron processing.

7. The single cell Hi-C data enhancement method according to claim 5, characterized in that: The attention mechanism module of step S222 includes a first attention mechanism module, a second attention mechanism module, and a third attention mechanism module; The cross-scale fusion module set in step S223 includes a fusion module from the second scale to the first scale, a fusion module from the third scale to the first scale, a fusion module from the first scale to the second scale, a fusion module from the third scale to the second scale, a fusion module from the first scale to the third scale, and a fusion module from the second scale to the third scale; The fusion process of step S223 includes: (1) The first attention mechanism module, the second attention mechanism module, and the third attention mechanism module process the node features after the single-scale convolution operation at the first resolution, the second resolution, and the third resolution respectively, and obtain the node potential feature attention representation of the node and edge fusion at the first resolution , node latent feature attention representation of nodes and edges fused at the second resolution , Node latent feature attention representation of nodes and edges fused at the third resolution ; (2) By calculating the inner product between the node potential feature attention representations of nodes and edges at different resolutions, the self-attention weight is calculated to obtain the attention weight matrix of each scale fusion, including the attention weight matrix from the first resolution to the second resolution , the attention weight matrix from the first resolution to the third resolution , the attention weight matrix from the second resolution to the first resolution , the attention weight matrix from the second resolution to the third resolution , the attention weight matrix from the third resolution to the first resolution , the attention weight matrix from the third resolution to the second resolution ; (3) Two-scale feature fusion: Each cross-scale fusion module fuses the attention representation of the node potential features at the corresponding resolution according to the corresponding attention weight matrix to obtain the weighted fusion features of the node potential features at the corresponding resolution; (4) Three-scale feature fusion: Based on the result of two-scale feature fusion, the weighted fusion features of node potential features at different resolutions are further weighted fused to obtain fused latent space features at different resolutions: ； ； ； In the formula They represent the node features after single-scale convolution operations at three different resolutions, and the weight w is used to control the fusion ratio; represents the fused latent space features of the first resolution, represents the fused latent space features of the second resolution, Represents the fused latent space features at the third resolution.

8. The single cell Hi-C data enhancement method according to claim 1, characterized in that: The enhancement method further comprises the steps of: S3. Identify domain boundaries in Hi-C data using the enhanced Hi-C contact matrix.

9. A single cell Hi-C data enhancement system, characterized in that: For implementing the enhancement method according to any one of claims 1 to 7, the enhancement system comprises the following modules: A preprocessing module is used to preprocess multi-resolution Hi-C data, extract the Hi-C contact matrix at different resolutions for each cell, construct the Hi-C contact matrix at different resolutions into a graph structure, and obtain graph structure data at different resolutions; The fusion unit model construction module constructs a variational graph autoencoder VGAE model as the fusion unit model, inputs graph structure data at different resolutions into the fusion unit model, performs feature extraction and matrix reconstruction, realizes cross-scale fusion of graph structure data at different resolutions, obtains fusion latent space features at different resolutions, and reconstructs the contact matrix as the enhanced Hi-C contact matrix.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and when the computer program is executed on a computer, the computer is enabled to execute the enhancement method according to any one of claims 1 to 8.

Images