Pseudo time sequence inference method of single cell Hi-C map and related device

Abstract

The invention discloses a pseudo time sequence inference method for a single-cell Hi-C map and a related device. The method comprises the following steps: extracting a single-cell Hi-C contact matrix from an original sequencing library; performing interpolation and feature enhancement on the extracted single cell Hi-C contact matrix to generate an enhanced contact matrix; based on the enhanced contact matrix, extracting multi-order chromatin structural features of the single cells, and embedding each single cell into a low-dimensional space for representation; from the low-dimensional cell representation, multiple discrete pseudo-timings of the cells under multiple k values are calculated, and a continuous development pseudo-timing of each cell is estimated based on these discrete pseudo-timings. The development path of the cells is estimated by adopting the shortest Hamiltonian path, and the final cell development pseudo-timing sequence is estimated based on the discrete pseudo-timing sequence under multiple k values, so that a pseudo-timing sequence inference result is insensitive to the number of cell clusters, and the accuracy is higher.

Description

Technical Field

[0001] This invention belongs to the field of biological computing technology, and specifically relates to a pseudo-time-series inference method and related device for single-cell Hi-C map. Background Technology

[0002] High-throughput chromosome conformation capture (Hi-C) is an emerging genome sequencing technology used to study the three-dimensional structure of chromatin and its gene regulatory functions. Unlike traditional high-throughput sequencing technologies that only obtain genomic sequence information, Hi-C fixes DNA fragments that are in long-range contact through chemical cross-linking. After steps such as enzyme digestion, ligation, labeling, and sequencing, a chromatin interaction matrix, also known as the Hi-C contact matrix or Hi-C map, is finally obtained. The values ​​of the elements in the contact matrix reflect the interaction frequency between chromatin fragments at corresponding locations; fragments with higher contact frequencies are closer in three-dimensional space. Using the Hi-C contact matrix, scientists have discovered many key three-dimensional structures on chromatin, including A / B compartments, topologically associated domains (TADs), chromatin bands, and chromatin loops, from large to small. Early Hi-C sequencing experiments were mainly conducted on population cells, which could not resolve the chromatin structure of individual cells. With the development of sequencing technology, single-cell Hi-C sequencing technology emerged. This technology can reveal the three-dimensional structure of chromatin at the single-cell level, providing a powerful tool for studying the gene regulation role of chromatin three-dimensional structure at the single-cell level.

[0003] While single-cell Hi-C mapping can characterize the spatial interaction frequencies of chromatin fragments in each single cell, similar to other single-cell omics sequencing, the developmental sequence of each single cell in a single-cell Hi-C sequencing library is unknown. This hinders the tracking and study of changes in the three-dimensional structure of chromatin during cell development. Single-cell pseudo-time sequence inference technology aims to infer the dynamic trajectory of cell development during the developmental process based on the inherent differences and connections in single-cell sequencing data, and to rank the chronological order of cell development. Because the inferred developmental time has no units and only represents the sequence, it is called "pseudo" time sequence. Pseudo-time sequence inference based on single-cell Hi-C sequencing data helps reveal changes in the three-dimensional structure of chromatin during cell division, differentiation, and disease development, providing a new perspective for exploring the three-dimensional gene regulatory mechanisms in individual development and disease occurrence.

[0004] Among existing pseudo-temporal inference methods for single-cell Hi-C data, Nagano et al. designed the earliest pseudo-temporal inference method based on biological priors, using several statistical measures to rank cell development time. The HiCRep / MDS method uses HiCRep as a measure of Hi-C contact matrix similarity, employs multidimensional scaling (MDS) dimensionality reduction to embed cells into a two-dimensional space, and calculates cell development time based on projection angles. The Circlet method designs four features on the Hi-C contact matrix, combines these features, and performs dimensionality reduction to map cells onto two semicircles for cell ranking. scHiCPTR employs a strategy combining pre-dimensionality reduction and re-dimensionality reduction to embed cells into a low-dimensional space, performs double graph trimming on the graph of cell developmental proximity relationships, and finally uses diffusion pseudo-time (DPT) to evaluate cell development time.

[0005] Although existing pseudo-temporal inference methods for single-cell Hi-C maps have made significant progress, the following problems still exist: (1) In the cell feature extraction stage, i.e., the cell embedding stage, the higher-order topological structure of chromatin represented by the Hi-C contact matrix is ​​not fully considered, so the embedding vector cannot effectively distinguish cells at different developmental stages; (2) In pseudo-temporal inference, existing methods either do not consider the proximity relationship between cell clusters, or artificially specify the number of cell clusters based on experience, ignoring the fine proximity relationship between cell subpopulations, so the accuracy of pseudo-temporal inference is still not high. Summary of the Invention

[0006] The purpose of this invention is to provide a pseudo-time-series inference method and related apparatus for single-cell Hi-C atlases, so as to solve the above-mentioned technical problems.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: In a first aspect, the present invention provides a pseudo-time-series inference method for single-cell Hi-C maps, comprising: Single-cell Hi-C contact matrices were extracted from the original sequencing library; The extracted single-cell Hi-C contact matrix is ​​interpolated and feature-enhanced to generate an enhanced contact matrix; Based on the enhanced contact matrix, multi-level chromatin structure features of single cells are extracted, and each single cell is embedded into a low-dimensional space representation. Based on the low-dimensional cell representation, the cell multiple discrete pseudo-temporal sequence under multiple k values ​​is calculated, and the continuous developmental pseudo-temporal sequence of each cell is estimated based on these discrete pseudo-temporal sequences.

[0008] Furthermore, the extraction of the single-cell Hi-C contact matrix from the original sequencing library includes: The process involves acquiring raw high-throughput sequencing reads from single cells in the sequencing library and filtering out low-quality, adapter-contaminated, and PCR-repeated reads to control read quality. Next, a sequence alignment algorithm is used to locate the reads to the reference genome, obtaining the position of each read on the reference genome. Then, the genome is divided into fixed-length bins according to the specified contact matrix resolution, and the number of sequencing reads falling between each bin is counted to form the contact frequencies in the corresponding subscripts of the Hi-C contact matrix.

[0009] Furthermore, the interpolation of the extracted single-cell Hi-C contact matrix includes: Considering the local similarity of interaction frequencies between chromatin segments, the interaction frequencies of target chromosome segment pairs are corrected using the interaction frequencies between coordinately adjacent chromatin segment pairs. Specifically, for a The contact matrix A and the convolution kernel F of size m = 2w + 1 are used, where n is the number of chromatin segments, m is the size of the convolution kernel, and w is a natural number. The values ​​of the elements in the corrected contact matrix B are calculated by scanning point-by-point across the entire matrix using convolution operations. The convolution process is expressed by the formula:

[0010] The interaction frequencies of chromatin fragment pairs that interact at long distances are interpolated using the restart random walk algorithm of the graph. Specifically, first, a probability transition matrix C is constructed, in which the elements... Let C represent the probability of migrating from the i-th chromatin segment to the j-th chromatin segment. C is calculated according to the following formula:

[0011] Based on the probability transition matrix C, the random walk process is calculated iteratively:

[0012] in, Let represent the matrix after the t-th random walk, with initial values ​​of . , It is an identity matrix, where p represents the probability of returning to the starting node in a random walk, and the condition for termination of the iteration is... The number of iterations is less than the specified threshold, or the number of iterations reaches the specified maximum number.

[0013] Furthermore, the feature enhancement of the extracted single-cell Hi-C contact matrix includes: Sort all elements in the contact matrix from largest to smallest, set elements less than the 20th percentile to 0, and set elements greater than or equal to the 20th percentile to 1, to obtain a binarized contact matrix. An undirected, unweighted graph is constructed using a binarized contact matrix as the adjacency matrix, and the local and global importance of edges is considered for scoring. The local importance of edges is scored using a dice coefficient (DC). The dice coefficient between the i-th node and the j-th node is calculated as follows:

[0014] in, Let represent the set of neighboring nodes of the i-th node. Indicates the number of elements in the set; Global importance is scored using edge betweenness centrality (EB). The method for calculating edge betweenness centrality is as follows:

[0015] in, Let C represent the number of shortest paths passing through nodes i and j, C represent the total number of shortest paths in the graph, and V be the set of nodes in the graph. The importance score of an edge is a weighted sum of the normalized dice coefficient DC and the edge betweenness centrality EB, expressed as:

[0016] in, Let be the maximum value of the dice coefficient of the edge in the graph. This represents the maximum value of the edge betweenness centrality in the graph; The importance of edges is used to measure the frequency and intensity of interactions between chromatin segments corresponding to nodes. As elements in the contact matrix after feature enhancement.

[0017] Furthermore, the step of extracting multi-level chromatin structure features of single cells based on the enhanced contact matrix and embedding each single cell into a low-dimensional representation includes: The HiCS single-cell hierarchical TADs domain identification package in Python is used to extract multi-level TADs hierarchical features of chromatin domains. Specifically, based on the hierarchical TADs domain identification results of the HiCS method, a three-dimensional feature is encoded for each bin, with each dimension of the feature representing the hierarchical TADs domain information of each bin: for non-TADs boundary points, the feature is encoded as [0,0,0]; for first-order boundary points, the feature is encoded as [1,0,0]; for second-order boundary points, the feature is encoded as [1,1,0]; and for third-order boundary points, the feature is encoded as [1,1,1]. The specific method of embedding each cell into a representation is as follows: Dimensionality reduction is performed on the TADs structural domain features of each cell hierarchy; let the library contain m cells, and each cell contain n bins; for first-order boundary points, the first elements of the feature encoding vectors of all bins of each cell are sequentially arranged to form an n-dimensional vector, and the vectors of the m cells in the library are combined to form an n-dimensional vector. The initial feature matrix of the first-order boundary points is obtained by using principal component analysis (PCA) to linearly reduce the dimensionality of the n-dimensional column vectors in the initial feature matrix to 1000 dimensions, resulting in a... The first-order TADs structural domain boundary point feature dimensionality reduction matrix is ​​obtained by sequentially reducing the dimensionality of the second-order and third-order boundary point features, respectively. The dimensionality reduction matrices of the boundary point features of first-order, second-order, and third-order TADs structural domains are stacked row by row to form a... The intermediate temporary matrix is ​​used to reduce the 3000-dimensional column vectors in the intermediate temporary matrix to 5 dimensions using the nonlinear dimensionality reduction method UMAP, resulting in a matrix with dimension 1. The cell embedding matrix, where each column vector of the cell embedding matrix represents the embedding domain representation of each cell.

[0018] Furthermore, the calculation of the cell multiple discrete pseudo-temporal sequence under multiple k-values ​​based on the low-dimensional cell representation includes: Iterate through the range of k values ​​specified by the user [k min k max For each integer k value within the range, calculate a set of discrete pseudo-time series. For the selected k value, the discrete pseudo-temporal sequence is calculated as follows: k-Means clustering is performed on the low-dimensional embedding domain representation of the obtained single cell, and the discrete pseudo-temporal sequence of each cluster is calculated based on the shortest Hamiltonian distance path. Specifically: For a given number of cell clusters k, the k-means clustering method is used to cluster the cells into k clusters, and an undirected weighted graph G=(V) is constructed with the cells of the k clusters as nodes. G E G ), as a proximity diagram between cell clusters, where V G and E G Let N be the set of nodes and the set of edges in the graph, respectively. The number of edges is N = k(k-1) / 2. Define the weight D of the edge between node i and node j. ij Let be the Euclidean distance between the centers of the i-th cluster and the j-th cluster; Sort the edges of graph G in ascending order according to their weights to obtain the edge sequence arranged in ascending order of weights, denoted as Ea = {e1, e2, ..., e...}. N}, where e1 is the edge with the smallest weight, e N The edge with the largest weight; Initialize a new unweighted graph S=(V S E S ), where V S = V G E S ={e1}, update the unweighted graph S according to the following rules: starting from e2, sequentially from E... a Take one e from the middle i , will e i When adding graph S, first determine if the degree of graph S is no greater than 2, then keep e. i Otherwise remove e i Secondly, determine if E is satisfied. S The number of elements in the middle reaches k-1, or e i =e N If the unweighted graph S is not updated, then the graph S obtained is the shortest Hamiltonian path through the centers of k cell clusters. Based on the order of the shortest Hamiltonian paths represented in Figure S, the k cell clusters are renamed sequentially as 1, 2, 3, ..., k, serving as the discrete pseudo-time for the entire cell cluster. Cells within each cluster have the same discrete pseudo-time. Each cell can obtain (k... max -k min +1) discrete pseudo-time labels.

[0019] Furthermore, the estimation of the continuous developmental pseudo-timeline for each cell based on these discrete pseudo-timelines includes: For all m cells (k max -k min +1) discrete pseudo-temporal labels are used for feature standardization, and the first two principal components are calculated using PCA; Based on the variable loadings in PCA, all pseudo-time series with the first principal component greater than the second principal component are screened out to obtain the largest set of interrelated pseudo-time series, while excluding pseudo-time series that are not related to them. The selected discrete pseudo-time series are averaged to generate an initial coarse continuous pseudo-time series; The initial coarse continuous pseudo-timeline was smoothed using locally weighted regression scatter smoothing, resulting in the final developmental pseudo-timeline for each cell.

[0020] Secondly, the present invention provides a pseudo-time-series inference system for single-cell Hi-C maps, comprising: The data extraction module is used to extract single-cell Hi-C contact matrices from the raw sequencing library; The enhancement module is used to interpolate and enhance the features of the extracted single-cell Hi-C contact matrix to generate an enhanced contact matrix. The low-dimensional representation module is used to extract multi-level chromatin structure features of single cells based on the enhanced contact matrix and embed each single cell into a low-dimensional space representation. The computation output module is used to compute multiple discrete pseudo-temporal sequences of cells under multiple k values ​​based on low-dimensional cell representations, and to estimate the continuous developmental pseudo-temporal sequence of each cell based on these discrete pseudo-temporal sequences.

[0021] Thirdly, the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the pseudo-time-series inference method for single-cell Hi-C atlases.

[0022] Fourthly, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the pseudo-time-series inference method for a single-cell Hi-C map.

[0023] Compared with the prior art, the present invention has the following technical effects: This invention extracts multi-level chromatin structural features (such as hierarchical information of TADs) to more finely distinguish cell developmental states, making it particularly suitable for analyzing cellular heterogeneity in complex biological processes such as embryonic development or disease progression. It employs multiple k-value calculations to compute discrete pseudo-time series and reduces dependence on the number of cell clusters (k-value) through averaging and smoothing, thereby reducing subjective bias and improving the stability and reproducibility of inference results. From data preprocessing to feature extraction and pseudo-time series computation, each step is interconnected, ensuring the method's applicability at single-cell resolution and providing a reliable computational tool for developmental biology and disease research.

[0024] Using the hierarchical features of TAD domains to represent cells can more effectively distinguish cells at different developmental stages. In the pseudo-temporal inference step, the proximity relationship between cell clusters is considered, the shortest Hamiltonian path is used to estimate the cell development path, and the final cell development pseudo-temporal is estimated based on discrete pseudo-temporal estimation under multiple k values. This makes the pseudo-temporal inference results insensitive to the number of cell clusters and more accurate. Attached Figure Description

[0025] Figure 1 This is a flowchart of a pseudo-time-series inference method for single-cell Hi-C maps provided in Embodiment 1 of the present invention.

[0026] Figure 2 This is an architecture diagram of the pseudo-time-series inference system for single-cell Hi-C mapping provided in Embodiment 2 of the present invention.

[0027] Figure 3This is a scatter plot of cell embeddings on the cell cycle dataset provided in Embodiment 3 of the present invention.

[0028] Figure 4 This is a scatter plot of pseudo-time series inference on the cell cycle dataset provided in Embodiment 3 of the present invention. Detailed Implementation

[0029] The present invention will be further described below with reference to the accompanying drawings: Example 1, please refer to Figure 1 This invention provides a pseudo-time-series inference method for single-cell Hi-C maps, comprising: Single-cell Hi-C contact matrices were extracted from the original sequencing library; The extracted single-cell Hi-C contact matrix is ​​interpolated and feature-enhanced to generate an enhanced contact matrix; Based on the enhanced contact matrix, multi-level chromatin structure features of single cells are extracted, and each single cell is embedded into a low-dimensional space representation. Based on the low-dimensional cell representation, the cell multiple discrete pseudo-temporal sequence under multiple k values ​​is calculated, and the continuous developmental pseudo-temporal sequence of each cell is estimated based on these discrete pseudo-temporal sequences.

[0030] This invention extracts multi-level chromatin structural features, enabling more precise differentiation of cell developmental states. It is particularly suitable for analyzing cellular heterogeneity in complex biological processes such as embryonic development or disease progression. Multiple k-values ​​are used to calculate discrete pseudo-time series, and averaging and smoothing processes reduce dependence on cell cluster size, thereby reducing subjective bias and improving the stability and reproducibility of the inference results.

[0031] Example 2: This invention provides a pseudo-time-series inference method for single-cell Hi-C maps, comprising: Step 1: Extract the single-cell Hi-C contact matrix from the original sequencing library.

[0032] In this implementation, firstly, all reads from the Hi-C paired-end sequencing library are acquired, and HiC-Pro software is used to filter out low-quality, adapter-contaminated, and PCR-repeated reads to control read quality. Secondly, the reads are located to the reference genome using BWA to obtain the position of each read on the reference genome. Then, the genome is divided into fixed-length bins according to the specified contact matrix resolution, and the number of sequencing reads falling between each pair of bins is counted to form the contact frequency in the corresponding subscript of the Hi-C contact matrix.

[0033] Step 2: Interpolate and enhance the features of the single-cell Hi-C contact matrix.

[0034] In this implementation, the specific method for interpolating the single-cell Hi-C contact matrix is ​​as follows: (1) Considering the local similarity of the interaction frequencies between chromatin segments, the interaction frequencies of the target chromosome segment pairs are corrected by using the interaction frequencies between chromatin segments that are adjacent in coordinates.

[0035] Specifically, this process is achieved using convolution operations, for a... The contact matrix A and the convolution kernel F of size m = 2w + 1 are given, where n is the number of chromatin segments, m is the size of the convolution kernel, and w is a natural number. In this embodiment, the size of the convolution kernel F is... The elements in the convolution kernel are all 1s. The convolution operation calculates the values ​​of the elements in the corrected contact matrix B by scanning point-by-point across the entire matrix. This convolution process can be expressed by the formula:

[0036] (2) Considering the global similarity of the interaction frequency patterns between nodes with large linear distances in the genome, the map restart random walk algorithm is used to interpolate the interaction frequencies of chromatin fragment pairs that interact at long distances.

[0037] Specifically, first, a probability transition matrix C is constructed, in which the elements... Let C represent the probability of migrating from the i-th chromatin segment to the j-th chromatin segment. C can be calculated using the following formula:

[0038] Based on the probability transition matrix C, the process of restarting the random walk can be calculated iteratively:

[0039] in, Let represent the matrix after the t-th random walk, with initial values ​​of . , It is an identity matrix, where p represents the probability of returning to the starting node in a random walk, and the condition for termination of the iteration is... The number of iterations is less than the specified threshold, or the number of iterations reaches the specified maximum number.

[0040] In this implementation, the restart probability p is set to 0.2, the iteration termination threshold is set to 10⁻⁶, and the maximum number of iterations is set to 20.

[0041] In this implementation, the specific method for feature enhancement of the single-cell Hi-C contact matrix is ​​as follows: (1) Binarization of the contact matrix. Sort all elements in the contact matrix from largest to smallest, set elements less than the 20th percentile to 0, and set elements greater than or equal to the 20th percentile to 1, to obtain the binarized contact matrix.

[0042] (2) Construct a graph and score the importance of edges. Construct an undirected, unweighted graph using the binarized contact matrix as the adjacency matrix, and score the local and global importance of edges.

[0043] The local importance of edges is scored using a dice coefficient (DC). The dice coefficient between the i-th node and the j-th node is calculated as follows:

[0044] in, Let represent the set of neighboring nodes of the i-th node. This indicates the number of elements in the set.

[0045] Global importance is scored using edge betweenness centrality (EB). The method for calculating edge betweenness centrality is as follows:

[0046] in, Let C represent the number of shortest paths passing through nodes i and j, C represent the total number of shortest paths in the graph, and V be the set of nodes in the graph.

[0047] The importance score of an edge is a weighted sum of the normalized dice coefficient DC and the betweenness centrality EB, which can be expressed as:

[0048] in, Let be the maximum value of the dice coefficient of the edge in the graph. It represents the maximum betweenness centrality of the edges in the graph.

[0049] Finally, the importance of edges is used to measure the frequency and intensity of interactions between chromatin segments corresponding to nodes. As elements in the contact matrix after feature enhancement.

[0050] Step 3: Extract the multi-level chromatin structure features of single cells on the contact matrix after interpolation, and embed each single cell into a low-dimensional space representation based on the multi-level chromatin structure features.

[0051] In this implementation, the specific method for extracting multi-level chromatin structure features of single cells is as follows: The multi-level TADs hierarchical features of chromatin domains are extracted using Python's HiCS single-cell hierarchical TADs domain recognition package. Specifically, based on the hierarchical TADs domain recognition results of the HiCS method, a three-dimensional feature is encoded for each bin. Each dimension of the feature represents the hierarchical TADs domain information of each bin: for non-TADs boundary points, the feature is encoded as [0,0,0]; for first-order boundary points, the feature is encoded as [1,0,0]; for second-order boundary points, the feature is encoded as [1,1,0]; and for third-order boundary points, the feature is encoded as [1,1,1].

[0052] In this implementation, the specific method for embedding each cell into a representation is as follows: (1) Dimensionality reduction is performed on the TADs structural domain features of each cell hierarchy. Let there be m cells in the library, and each cell have n bins. For first-order boundary points, the first element of the feature encoding vector of all bins in each cell is sequentially arranged to form an n-dimensional vector. The vectors of the m cells in the library are then combined to form an n-dimensional vector. The initial feature matrix of the first-order boundary points is linearly reduced to 1000 dimensions by principal component analysis (PCA), resulting in a... The dimensionality reduction matrix of the first-order TADs structural domain boundary point features is obtained. Similarly, the dimensionality reduction of the second-order and third-order boundary point features is performed sequentially to obtain the second-order and third-order TADs structural domain boundary point feature dimensionality reduction matrices, respectively.

[0053] (2) Stack the first-order, second-order, and third-order TADs structural domain boundary point feature dimension reduction matrices row by row to form a The intermediate temporary matrix is ​​further reduced to 5 dimensions using the nonlinear dimensionality reduction method UMAP, resulting in a matrix with dimension 1. The cell embedding matrix, where each column vector of the cell embedding matrix represents the embedding domain representation of each cell.

[0054] Step 4: Calculate the discrete pseudo-time series under multiple k values.

[0055] In this implementation, the method for calculating discrete pseudo-time series under multiple k-values ​​is as follows: traverse the user-specified range of k-values ​​[k min k max Given all integer k values ​​in the range, calculate a set of discrete pseudo-time series for each k value.

[0056] For the selected k value, the method for calculating the discrete pseudo-temporal sequence is as follows: perform k-Means clustering on the low-dimensional embedding domain representation of the single cell obtained in step 3, and calculate the discrete pseudo-temporal sequence of each cluster based on the shortest Hamiltonian distance path.

[0057] Specifically: (1) For a given number of cell clusters k, the k-means clustering method is used to cluster the cells into k clusters, and an undirected weighted graph G=(V) is constructed with the cells of the k clusters as nodes. G E G ), as a proximity diagram between cell clusters, where V G and E G Let N be the set of nodes and the set of edges in the graph, respectively. The number of edges is N = k(k-1) / 2. Define the weight D of the edge between node i and node j. ij Let be the Euclidean distance between the centers of the i-th cluster and the j-th cluster.

[0058] (2) Sort the edges of graph G in ascending order according to their weights to obtain the edge sequence arranged in ascending order of weights, denoted as E. a ={e1, e2, ..., e N}, where e1 is the edge with the smallest weight, e N The edge with the largest weight.

[0059] (3) Initialize a new unweighted graph S=(V S E S ), where V S = V G E S ={e1}, update the unweighted graph S according to the following rules: starting from e2, sequentially from E... a Take one e from the middle i , will e i When adding graph S, first determine if the degree of graph S is no greater than 2, then keep e. i Otherwise remove e i Secondly, determine if E is satisfied. S The number of elements in the middle reaches k-1, or e i =e N If the unweighted graph S is not updated, then the graph S is stopped. The final graph S is the shortest Hamiltonian path through the centers of the k cell clusters.

[0060] (4) According to the order of the shortest Hamiltonian path shown in Figure S, the k cell clusters are renamed as 1, 2, 3, ..., k, which serve as the discrete pseudo-time of the entire cell cluster. Each cell in the cluster has the same discrete pseudo-time.

[0061] After step 4, each cell can obtain (k max -k min+1) discrete pseudo-time labels.

[0062] Step 5: Estimate the developmental pseudo-timing of each cell based on the multiple discrete pseudo-timing.

[0063] In this implementation, the specific method for estimating the pseudo-timeline of development for each cell is as follows: (1) For all m cells (k max -k min +1) discrete pseudo-time series labels are used for feature standardization, and the first two principal components are calculated using PCA.

[0064] (2) Based on the variable loadings in PCA, filter out all pseudo-time series where the first principal component is greater than the second principal component to obtain the largest set of interrelated pseudo-time series, and exclude pseudo-time series that are not related to it.

[0065] (3) Average the selected discrete pseudo-time series to generate the initial coarse continuous pseudo-time series.

[0066] (4) The initial coarse continuous pseudo-time series was smoothed using Locally Weighted Regression Scatterplot Smoothing (LOESS) as the final developmental pseudo-time series for each cell.

[0067] Example 3: This example uses the Cell-cycle dataset released by Nagano et al. in 2017 to conduct experiments. The dataset contains four states of the cell cycle. After filtering out single cells with a contact number of less than 5k, the number of cells in each state of the dataset is as follows: 280 cells in G1 stage, 303 cells in ES stage, 262 cells in MS stage, and 326 cells in G2 stage.

[0068] First, the results of cell embedding are visualized, such as... Figure 3 The figure shows a scatter plot of the first two dimensions of the cell embedding vector. It can be seen that the cell embedding scatter plot of the present invention presents a clear ring-shaped topological structure, and cells at different stages can be clearly distinguished.

[0069] Secondly, the pseudo-timing information inferred by this invention is displayed in pseudo-color on a scatter plot in the embedded space, such as... Figure 4 As shown in the figure, the pseudo-chronology of cell development increases along the circular trajectory shown in the figure. The cell goes through four stages, G1, ES, MS and G2, to complete a complete cell cycle, which is consistent with the biological process.

[0070] Finally, the pseudo-temporal inference results of the present invention on cell cycle datasets were quantitatively compared with those of five existing methods, namely HiCRep / MDS, Nagano et al., Circlet / MCM, Circlet / CDD, and Circlet / MCM+CDD, using the Kendall rank correlation coefficient index. As shown in Table 1, it can be seen that the pseudo-temporal inference method of the single-cell Hi-C map of the present invention achieved the optimal Kendall rank correlation coefficient, which demonstrates the effectiveness and advancement of the present invention.

[0071] Table 1. Comparison of Kendall rank correlation coefficients for pseudo-time series inference using different methods on the cell cycle dataset.

[0072] In another embodiment of the present invention, a pseudo-temporal inference system for single-cell Hi-C atlases is provided, which can be used to implement the above-mentioned pseudo-temporal inference method for single-cell Hi-C atlases. Specifically, the system includes: The data extraction module is used to extract single-cell Hi-C contact matrices from the raw sequencing library; The enhancement module is used to interpolate and enhance the features of the extracted single-cell Hi-C contact matrix to generate an enhanced contact matrix. The low-dimensional representation module is used to extract multi-level chromatin structure features of single cells based on the enhanced contact matrix and embed each single cell into a low-dimensional space representation. The computation output module is used to compute multiple discrete pseudo-temporal sequences of cells under multiple k values ​​based on low-dimensional cell representations, and to estimate the continuous developmental pseudo-temporal sequence of each cell based on these discrete pseudo-temporal sequences.

[0073] The module division in this embodiment of the invention is illustrative and represents only one logical functional division. In actual implementation, other division methods may be used. Furthermore, the functional modules in the various embodiments of the invention can be integrated into a single processor, exist as separate physical entities, or be integrated into a single module. The integrated modules described above can be implemented in hardware or as software functional modules.

[0074] In another embodiment of the present invention, a computer device is provided, comprising a processor and a memory. The memory stores a computer program, which includes program instructions. The processor executes the program instructions stored in the computer storage medium. The processor may be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing and control core of the terminal, suitable for implementing one or more instructions, specifically suitable for loading and executing one or more instructions from the computer storage medium to achieve a corresponding method flow or corresponding function. The processor described in this embodiment of the present invention can be used for the operation of a pseudo-timing inference method for single-cell Hi-C mapping.

[0075] In another embodiment of the present invention, a storage medium is provided, specifically a computer-readable storage medium (Memory), which is a memory device in a computer device used to store programs and data. It is understood that the computer-readable storage medium here can include both the built-in storage medium in the computer device and extended storage media supported by the computer device. The computer-readable storage medium provides storage space that stores the terminal's operating system. Furthermore, the storage space also stores one or more instructions suitable for loading and execution by a processor. These instructions can be one or more computer programs (including program code). It should be noted that the computer-readable storage medium here can be high-speed RAM or non-volatile memory, such as at least one disk storage device. The processor can load and execute one or more instructions stored in the computer-readable storage medium to implement the corresponding steps of the pseudo-timing inference method for single-cell Hi-C mapping in the above embodiments.

[0076] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0077] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0078] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0079] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0080] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A pseudo-time-series inference method for single-cell Hi-C maps, characterized in that, include: Single-cell Hi-C contact matrices were extracted from the original sequencing library; The extracted single-cell Hi-C contact matrix is ​​interpolated and feature-enhanced to generate an enhanced contact matrix; Based on the enhanced contact matrix, multi-level chromatin structure features of single cells are extracted, and each single cell is embedded into a low-dimensional space representation. Based on the low-dimensional cell representation, the cell multiple discrete pseudo-temporal sequence under multiple k values ​​is calculated, and the continuous developmental pseudo-temporal sequence of each cell is estimated based on these discrete pseudo-temporal sequences.

2. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The extraction of the single-cell Hi-C contact matrix from the original sequencing library includes: The process involves acquiring raw high-throughput sequencing reads from single cells in the sequencing library and filtering out low-quality, adapter-contaminated, and PCR-repeated reads to control read quality. Next, a sequence alignment algorithm is used to locate the reads to the reference genome, obtaining the position of each read on the reference genome. Then, the genome is divided into fixed-length bins according to the specified contact matrix resolution, and the number of sequencing reads falling between each bin is counted to form the contact frequencies in the corresponding subscripts of the Hi-C contact matrix.

3. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The interpolation of the extracted single-cell Hi-C contact matrix includes: Considering the local similarity of interaction frequencies between chromatin segments, the interaction frequencies of target chromosome segment pairs are corrected using the interaction frequencies between coordinately adjacent chromatin segment pairs. Specifically, for a The contact matrix A and the convolution kernel F of size m = 2w + 1 are used, where n is the number of chromatin segments, m is the size of the convolution kernel, and w is a natural number. The values ​​of the elements in the corrected contact matrix B are calculated by scanning point-by-point across the entire matrix using convolution operations. The convolution process is expressed by the formula: The interaction frequencies of chromatin fragment pairs that interact at long distances are interpolated using the graph restart random walk algorithm. Specifically, first, a probability transition matrix C is constructed, in which the elements... Let C represent the probability of migrating from the i-th chromatin segment to the j-th chromatin segment. C is calculated according to the following formula: Based on the probability transition matrix C, the random walk process is calculated iteratively: in, Let represent the matrix after the t-th random walk, with initial values ​​of . , It is an identity matrix, where p represents the probability of returning to the starting node in a random walk, and the condition for termination of the iteration is... The number of iterations is less than the specified threshold, or the number of iterations reaches the specified maximum number.

4. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The feature enhancement of the extracted single-cell Hi-C contact matrix includes: Sort all elements in the contact matrix from largest to smallest, set elements less than the 20th percentile to 0, and set elements greater than or equal to the 20th percentile to 1, to obtain a binarized contact matrix. An undirected, unweighted graph is constructed using a binarized contact matrix as the adjacency matrix, and the local and global importance of edges is considered for scoring. The local importance of edges is scored using a dice coefficient (DC). The dice coefficient between the i-th node and the j-th node is calculated as follows: in, Let represent the set of neighboring nodes of the i-th node. Indicates the number of elements in the set; Global importance is scored using edge betweenness centrality (EB). The method for calculating edge betweenness centrality is as follows: in, Let C represent the number of shortest paths passing through nodes i and j, C represent the total number of shortest paths in the graph, and V be the set of nodes in the graph. The importance score of an edge is a weighted sum of the normalized dice coefficient DC and the edge betweenness centrality EB, expressed as: in, Let be the maximum value of the dice coefficient of the edge in the graph. This represents the maximum value of the edge betweenness centrality in the graph; The importance of edges is used to measure the frequency and intensity of interactions between chromatin segments corresponding to nodes. As elements in the contact matrix after feature enhancement.

5. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The process of extracting multi-level chromatin structure features of single cells based on the enhanced contact matrix and embedding each single cell into a low-dimensional representation includes: The HiCS single-cell hierarchical TADs domain identification package in Python is used to extract multi-level TADs hierarchical features of chromatin domains. Specifically, based on the hierarchical TADs domain identification results of the HiCS method, a three-dimensional feature is encoded for each bin, with each dimension of the feature representing the hierarchical TADs domain information of each bin: for non-TADs boundary points, the feature is encoded as [0,0,0]; for first-order boundary points, the feature is encoded as [1,0,0]; for second-order boundary points, the feature is encoded as [1,1,0]; and for third-order boundary points, the feature is encoded as [1,1,1]. The specific method of embedding each cell into a representation is as follows: Dimensionality reduction is performed on the TADs structural domain features of each cell hierarchy; let the library contain m cells, and each cell contain n bins; for first-order boundary points, the first elements of the feature encoding vectors of all bins of each cell are sequentially arranged to form an n-dimensional vector, and the vectors of the m cells in the library are combined to form an n-dimensional vector. The initial feature matrix of the first-order boundary points is obtained by using principal component analysis (PCA) to linearly reduce the dimensionality of the n-dimensional column vectors in the initial feature matrix to 1000 dimensions, resulting in a... The first-order TADs structural domain boundary point feature dimensionality reduction matrix is ​​obtained by sequentially reducing the dimensionality of the second-order and third-order boundary point features, respectively. The dimensionality reduction matrices of the boundary point features of first-order, second-order, and third-order TADs structural domains are stacked row by row to form a... The intermediate temporary matrix is ​​used to reduce the 3000-dimensional column vectors in the intermediate temporary matrix to 5 dimensions using the nonlinear dimensionality reduction method UMAP, resulting in a matrix with dimension 1. The cell embedding matrix, where each column vector of the cell embedding matrix represents the embedding domain representation of each cell.

6. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The calculation of the cell multiple discrete pseudo-temporal sequence under multiple k-values ​​based on the low-dimensional cell representation includes: Iterate through the range of k values ​​specified by the user [k min k max For each integer k value within the range, calculate a set of discrete pseudo-time series. For the selected k value, the discrete pseudo-temporal sequence is calculated as follows: k-Means clustering is performed on the low-dimensional embedding domain representation of the obtained single cell, and the discrete pseudo-temporal sequence of each cluster is calculated based on the shortest Hamiltonian distance path. Specifically: For a given number of cell clusters k, the k-means clustering method is used to cluster the cells into k clusters, and an undirected weighted graph G=(V) is constructed with the cells of the k clusters as nodes. G E G ), as a proximity diagram between cell clusters, where V G and E G Let N be the set of nodes and the set of edges in the graph, respectively. The number of edges is N = k(k-1) / 2. Define the weight D of the edge between node i and node j. ij Let be the Euclidean distance between the centers of the i-th cluster and the j-th cluster; Sort the edges of graph G in ascending order according to their weights to obtain the edge sequence arranged in ascending order of weights, denoted as Ea = {e1, e2, ..., e...}. N }, where e1 is the edge with the smallest weight, e N The edge with the largest weight; Initialize a new unweighted graph S=(V S E S ), where V S = V G E S ={e1}, update the unweighted graph S according to the following rules: starting from e2, sequentially from E... a Take one e from the middle i , will e i When adding graph S, first determine if the degree of graph S is no greater than 2, then keep e. i Otherwise remove e i Secondly, determine if E is satisfied. S The number of elements in the middle reaches k-1, or e i =e N If the unweighted graph S is not updated, then the graph S obtained is the shortest Hamiltonian path through the centers of k cell clusters. Based on the order of the shortest Hamiltonian paths represented by Figure S, the k cell clusters are renamed sequentially as 1, 2, 3, ..., k, serving as the discrete pseudo-time for the entire cell cluster. Cells within each cluster have the same discrete pseudo-time; each cell can obtain (k... max -k min +1) discrete pseudo-time labels.

7. The pseudo-time-series inference method for single-cell Hi-C maps according to claim 1, characterized in that, The estimation of the continuous developmental pseudo-timing of each cell based on these discrete pseudo-timings includes: For all m cells (k max -k min +1) discrete pseudo-temporal labels are used for feature standardization, and the first two principal components are calculated using PCA; Based on the variable loadings in PCA, all pseudo-time series with the first principal component greater than the second principal component are screened out to obtain the largest set of interrelated pseudo-time series, while excluding pseudo-time series that are not related to them. The selected discrete pseudo-time series are averaged to generate an initial coarse continuous pseudo-time series; The initial coarse continuous pseudo-timeline was smoothed using locally weighted regression scatter smoothing, resulting in the final developmental pseudo-timeline for each cell.

8. A pseudo-time-series inference system for single-cell Hi-C maps, characterized in that, include: The data extraction module is used to extract single-cell Hi-C contact matrices from the raw sequencing library; The enhancement module is used to interpolate and enhance the features of the extracted single-cell Hi-C contact matrix to generate an enhanced contact matrix. The low-dimensional representation module is used to extract multi-level chromatin structure features of single cells based on the enhanced contact matrix and embed each single cell into a low-dimensional space representation. The computation output module is used to compute multiple discrete pseudo-temporal sequences of cells under multiple k values ​​based on low-dimensional cell representations, and to estimate the continuous developmental pseudo-temporal sequence of each cell based on these discrete pseudo-temporal sequences.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the pseudo-time-series inference method for a single-cell Hi-C map as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the pseudo-time-series inference method for a single-cell Hi-C map as described in any one of claims 1 to 7.

Images