A reeb graph-based single-cell pseudo trajectory topology identification method

Abstract

The application belongs to the technical field of single cell data analysis, and particularly relates to a single cell pseudo trajectory topological identification method based on a Reeb graph. The method comprises the following steps: S1, preprocessing a single cell transcription data set to obtain a single cell gene expression matrix; S2, constructing an optimization model, constructing a single cell connected network, and identifying a cell group; S3, combining a Morse function distribution of the cell and an equivalence class definition to extract a single cell Reeb graph; and S4, combining a maximum value of the Morse function to identify a key node in the Reeb graph, deduce an intermediate stage and a start and end state of cell differentiation, measure a distance influence of the cell in the network and the equivalence class, and identify a cell pseudo trajectory. The application extracts a multi-directional differentiation trajectory of the cell based on a topological structure Reeb graph of the cell connected network. The algorithm framework can convert the cell network into a Reeb graph with less equivalent nodes, and reduces the calculation complexity. Compared with an initial time determination which is completely dependent on prior information, the application provides the properties of the Morse function to determine the key state of cell differentiation.

Description

Technical Field

[0001] This invention belongs to the field of single-cell data analysis technology, specifically relating to a single-cell pseudo-trajectory topology recognition method based on Reeb graphs. Background Technology

[0002] Recent advances in single-cell analysis have enabled the measurement of single-cell transcriptomic profiles at many stages of proliferation and differentiation in organisms. The challenge lies in designing computational techniques and modeling methods capable of automatically reconstructing cell differentiation pathways, sequencing individual cells at the transcriptomic level, and identifying trajectory behaviors along branching lineages.

[0003] Due to the characteristics of single-cell data such as high statistical noise levels, high-dimensional gene expression space, and nonlinearity, pseudo-trajectory inference modeling based on scRNA-seq data analysis of dynamic processes such as cell differentiation, cell fate development, and decision-making is mainly divided into three categories.

[0004] The first type of trajectory inference method infers cell differentiation trajectories based on the minimum spanning tree (MST) algorithm. The Monocle algorithm uses the MST algorithm to find the longest path on a single cell and assigns pseudotime to each cell accordingly. The Monocle2 algorithm uses the MST to learn cell trajectories, updating cell positions by moving cells to the nearest vertex in the MST until the cell trajectory and final position stabilize. Finally, it obtains the cell's pseudotime based on the geodesic distance from the cell to the root vertex along the MST. The TSCAN and Waterfall algorithms construct cell cluster-based MSTs and sort cells using orthogonal projection, reducing the overall computational cost of the MST algorithm. The Slingshot algorithm uses the MST to construct trajectories and proposes a joint master curve algorithm to construct smooth curves. Cells are then projected onto these smooth curves instead of the edges of the MST to calculate cell pseudotime.

[0005] The second category is trajectory inference methods based on graph models. The Diffusion Pseudo-Time (DPT) algorithm uses a weighted k-nearest neighbor (KNN) algorithm to construct cell trajectories and introduces ten diffusion algorithms to calculate the pseudo-time of cells in the "diffusion mapping space." The PAGA algorithm calculates neighborhood relationships by selecting appropriate low-dimensional representations and relevant distance metrics, representing high-dimensional gene expression data as a KNN graph. It uses the Louvain algorithm to partition the vertices of various KNN graphs, combining high-confidence paths generated in the PAGA graph with distance metrics based on random walks to calculate the distance to the root cell and rank the cells. The Monocle3 algorithm uses a graph algorithm to generate trajectories and calculates the geodesic distance from the cell in the trajectory to the selected root node as the pseudo-time. The URD algorithm uses a KNN graph of the transcriptome in the gene expression space to construct trajectories and uses a simulated diffusion algorithm to calculate the pseudo-time of cells.

[0006] The third category is RNA velocity-assisted trajectory inference methods. The VeTra algorithm uses RNA velocity vectors to construct multiple directed graphs and determines the cell's transition state based on the cell's k nearest neighbors. It identifies weakly connected branches by constructing independent cell transition paths and finally clusters these transition paths to obtain the trajectory. The Cytopath algorithm combines the inter-cell transition probability matrix and cell state to construct multiple simulated trajectories and uses RNA velocity to infer the root and terminal states of the cell differentiation trajectory.

[0007] However, in trajectory inference, it is common practice to determine the starting point of a trajectory based on prior information, using the relative relationships between the remaining individual cells or cell clusters and the starting point (such as the shortest path length) as additional information for constructing the trajectory. Using only the starting point information to infer pseudo-times for all cells leads to issues with the rationality and stability of trajectory inference. The Morse function is defined on both the shape and the manifold, and its Hessian matrix at critical points (where the first derivative of the Morse function is 0) degenerates. Reeb plots reveal the topological structure of the manifold and its shape based on the Morse function. Summary of the Invention

[0008] The purpose of this invention is to propose an unsupervised algorithm called SCTRRG (Single-Cell Topology Recognition and Reeb Graph) for reconstructing high-resolution evolutionary paths along topological structures. This method uses multidimensional single-cell data as input to construct network relationships between cells, defines a Morse function to identify key initiation fates and branching cells, and infers trajectories in multi-branched cell fates based on the Reeb graph's topological structure. Two simulated datasets demonstrate that SCTRRG can successfully identify their distinct multi-branched structures with reasonable initiation and branching points. Furthermore, applying SCTRRG to two published single-cell RNA-seq datasets, SCTRRG accurately reconstructs known cell lineages in mouse embryonic and neural stem cell development.

[0009] To achieve the aforementioned technical effects and objectives, the first aspect of this invention provides a method for topological recognition of single-cell pseudo-trajectories based on Reeb graphs, comprising:

[0010] S1: Preprocess the single-cell transcription dataset to obtain the single-cell gene expression matrix;

[0011] S2: Based on the single-cell gene expression matrix of S1, an optimized model is constructed, and a single-cell connectivity network is further constructed for cell population identification;

[0012] S3: Based on S2, extract the Reeb map of a single cell by combining the Morse function distribution and equivalence class definition of the cell.

[0013] S4: Building upon S3, this method combines the maximum and minimum values ​​of the Morse function to identify key nodes in the Reeb graph, infer the intermediate stages and initial and final states of cell differentiation, measure the distance influence of cells in the network and equivalence classes, and identify cell pseudotrajectories.

[0014] Preferably, S2 includes:

[0015] S2.1: Constructing an optimization model:

[0016]

[0017] stA T y = y;

[0018]

[0019] y = (1,1,…,1) T

[0020] Where A corresponds to the adjacency matrix of a single-cell connected network, and X is the gene expression matrix of a single cell. * It is the nuclear norm, ||·|| 2,1 It is the sum of the Euclidean norms of all columns; the nonnegative parameter λ is used to balance the nuclear norm and the Euclidean norm; linear constraint A T y = y preserves the invariance of data translation transformations, where y is a vector of length n. This indicates that the i-th cell and the j-th cell are not neighboring structures;

[0021] The Alternating Direction Multiplier Method (ADMM) is used to solve the optimization model of S2.1. The K-Nearest Neighbors (KNN) algorithm is used to determine the nearest neighbor network G of the cell, and the selection of K is based on the BIC principle. It is agreed that non-zero values ​​in Aopt are 1, so that the 0-1 single-cell similarity network Aopt is obtained.

[0022] S2.2: Based on the similarity network Aopt of S2.1, construct a single-cell connected network for cell population identification:

[0023]

[0024] The 0-1 single-cell similarity network Aopt is transformed into a single-cell symmetric network Z, λ using the maximum value method. l Let be the l-th eigenvalue of L arranged in ascending order, and p be the index corresponding to the maximum absolute difference between two adjacent eigenvalues, used to identify the number of cell populations. Combined with spectral clustering algorithms, p cell populations are identified, namely C1, C2, ..., Cp. When Z is a disconnected network, dis represents the Euclidean distance function. Let i represent the i-th cell in cell group u and j-th cell in cell group v, respectively. Calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell groups. Where i * ,j * They are respectively Based on the cell sequence throughout the network, Z is progressively updated into a single-cell connected network; for cases where there are edges between two cell groups, Z... ij =1, then there is an edge between the two cell groups, and the cell group connectivity network CC between the cell groups is further obtained.

[0025] Preferably, S3 includes:

[0026] S3.1: The distribution of Morse function in a single cell is as follows:

[0027]

[0028] Where u(c) i c j The distance f(i) is the geometric distance between the i-th and j-th cells in a disconnected network Z, approximated by Dijkstra's algorithm. i Normalization to To maintain the invariance of the f value relative to the grid scale, where n is the number of cells;

[0029] S3.2: Extracting Reebok images from single cells f :

[0030]

[0031] R f =G(V) r E r )

[0032]

[0033] in, and Let represent the average geodesic distances of the u-th and v-th cells in the i-th cell group, respectively, and d represent the error. An equivalence relation within the same cell population is defined as the difference in the Morse function being within a very small error.

[0034] Based on equivalence classes The definition of Reeb graph R is based on the extraction of Reeb graphs from single-cell connected networks. f R f V r Nodes are determined by representative nodes in equivalence classes, R f Edge set E rThe existence of an edge is determined by the connectivity of cells in the two equivalence classes [p] and [q] in the connected network Z. If there exists a connection between cell s and cell t, then A s，t ≠0, E r There are edges connecting the equivalent nodes [p] and [q].

[0035] Preferably, S4 includes:

[0036] S4.1: Sort the nodes in the Reeb graph of S3:

[0037]

[0038] Among them, V B It is a cell population located at the boundary of the cell population connected network C, with a degree of 1; based on the continuity and equivalence class definition of the Morse function, the largest equivalence class of the Morse function corresponds to R. f The initial equivalent node v on [0] This corresponds to the initial or final state of cell differentiation; the equivalent cell class v with the smallest Morse function. [m] Located in the entire R f The middle position corresponds to the vicinity of the middle position of cell differentiation;

[0039] S4.2: Based on the initial equivalent node v [0] And S3, the entire Reeb diagram R f Calculate the shortest path length from each node in the Reeb graph to the initial equivalent node, denoted as Len. Based on the size of Len, calculate the shortest path length for R. f n r Sort the equivalent nodes:

[0040] T [i] =Len(v [i] ,v [0] );

[0041] S4.3: Measuring cell distance:

[0042]

[0043] Where d(s,j) represents the shortest path length between the s-th and j-th cells in the single-cell connection network Z, D ij The metric is the average distance from cell j to the set of cells in the [i] equivalence class. The metric is the average distance from cell j in cell population c to the cell set of equivalence class [i]. exp is the exponential function, var represents the variance function, and D is the normalized value. ij , It is uniformly denoted as the average distance between cells. The effects of the two distances on cells were characterized from the perspectives of the overall cell network and the local networks of cells within the cell population;

[0044] S4.4: Combining the influence of intercellular distance and the order of cell nodes in the Reeb plot, pseudo-time points for all cells are obtained, thus achieving single-cell pseudo-trajectory inference:

[0045]

[0046] If we consider the influence of the overall cell network on cell v, then i = 1, 2, ..., n r v = 1, 2, ..., n; if we consider the local networks of cells in the cell population, then

[0047] A second aspect of the present invention provides a single-cell pseudo-trajectory topology recognition device based on Reeb graphs. The device is constructed based on the single-cell pseudo-trajectory topology recognition method of Reeb graphs to realize single-cell pseudo-trajectory topology recognition.

[0048] A third aspect of the present invention provides a single-cell pseudo-trajectory topology recognition device based on Reeb graphs. The control device includes a processor and a memory. The memory stores a computer program, which is loaded and executed by the processor to implement the single-cell pseudo-trajectory topology recognition method based on Reeb graphs.

[0049] A fourth aspect of the present invention provides a computer-readable storage medium storing a computer program, which is loaded and executed by a processor to implement a single-cell pseudo-trajectory topology recognition method based on Reeb graphs.

[0050] Compared with the prior art, the beneficial effects of the present invention are:

[0051] 1. The method proposed in this invention is based on the Reeb graph of the cellular connection network topology to extract the differentiation trajectories of cells in different directions. On the one hand, the algorithm framework can transform the cellular network into a Reeb graph with fewer equivalent nodes, thereby reducing computational complexity; on the other hand, compared with the initial time determination that relies entirely on prior information, it provides the properties of the Morse function to determine the initial time.

[0052] 2. Based on existing networks or trees to measure inter-cell relationships, we define equivalence classes for cells and propose REEB diagrams to simulate positional relationships during cell differentiation.

[0053] 3. Using the Morse function distribution of cells to quantitatively characterize the initial or terminal state and intermediate stage of cell differentiation, it was found that cells with the highest Morse function value are often at the initial or terminal moment of differentiation;

[0054] 4. SCTRRG can accurately identify the lineage structure of cell differentiation corresponding to the four datasets and infer its pseudo-time and branch trajectories. Compared with other classic methods, it has better accuracy in identifying pseudo-trajectories. Attached Figure Description

[0055] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used together with the embodiments of the invention to explain the invention and do not constitute a limitation thereof.

[0056] In the attached diagram:

[0057] Figure 1 This is a framework diagram of the SCTRRG method of this invention; (A) The input of SCTRRG is a single-cell matrix of n cells and m genes. (B) The network is clustered using a nonlinear optimization model, and 6 cell clusters (C0, C1, C2, C3, C4, and C5) are obtained through spectral clustering. (C) Distribution of six Morse functions (MF). (D) Cell lineage. The initial cluster (C0) is identified based on the connected network structure defined by MF. (E) The Reeb graph reflects the basic topological structure of the connected cell network based on MF. In the two-dimensional spectral embedding space, the pseudo-time of the method (F) (normalized to 0-1) has an optional initial state (maximum number of boundary component nodes) and a median state (minimum number of branch component nodes) in the evolution trajectory.

[0058] Figure 2 This is a schematic diagram illustrating the transition from single-cell similar networks to sparsely connected networks.

[0059] Figure 3 Calculate the equivalence class for Morse functions.

[0060] Figure 4 The results are as follows: (a) Cell distribution in the low-dimensional space composed of second and third dimensionality reduction features is labeled with initial classification labels (1, 2, 3) based on the diffusion mapping dimensionality reduction method. (b) Six cell subgroups and cell lineage structures are identified using the SCTRRG algorithm. (c) Cell distribution in the spectral embedding space composed of two-dimensional and three-dimensional dimensionality reduction features is colored using the Morse function values ​​of the cells based on the spectral embedding dimensionality reduction method. (d) Comparison of Morse function distributions of the six cell groups. The label colors are consistent with those in Figure b. (e) Reeb plot of the connected cell network in the S1 data. The largest labeled point represents the initial equivalent node of unit group C1, and the smallest labeled point represents the intermediate equivalent node of unit group C2. The label colors are consistent with those of C1, C2, C3, C4, and C5 in Figure B, with orange representing C6. (f) Cell distribution in the spectral embedding space, stained and labeled by pseudo-time points of the cells.

[0061] Figure 5To identify cell lineages and pseudotimes in the multi-branch simulation data S2 using the SCTRRG algorithm: (a) Based on the diffusion mapping dimensionality reduction method, the cell distribution in the low-dimensional space composed of the second and third dimensionality reduction features is labeled with initial classification labels (1,2,3,4,5). (b) Nine cell subpopulations and cell lineage structures are identified using the SCTRRG algorithm. (c) Based on the spectral embedding dimensionality reduction method, the cell distribution in the spectral embedding space composed of the second and third dimensionality reduction features is colored using the Morse function values ​​of the cells. (d) Comparison of the Morse function distributions of the nine groups of cells. The label colors are consistent with Figure b. (e) Reeb plot of the multi-branch simulation data (S2) connecting cell networks. The largest labeled point represents the initial equivalent node belonging to C1, and the smallest labeled point represents the intermediate equivalent node belonging to C2. The label colors are consistent with Figure b. (f) Cell distribution in the spectral embedding space, colored and labeled by the pseudotimes of the cells.

[0062] Figure 6 SCTRRG is used to infer cell lines and pseudotime during the differentiation of mouse neural stem cells. (a) Based on the diffusion mapping dimensionality reduction method, the distribution of mouse neural stem cells in the low-dimensional space composed of one-dimensional and two-dimensional dimensionality reduction features is labeled with the first known label (1,2,3,4). (b) Cell lineages formed by four cell subpopulations (qNSC, NB, aNSC and Oli) identified by the SCTRRG algorithm. The distribution of 272 neural stem cells in the spectral embedding space composed of one-dimensional and two-dimensional dimensionality reduction features is stained using the Morse function value of the cells. The label color is consistent with Figure b. (c) Based on the spectral embedding method, the distribution of 272 neural stem cells in the spectral embedding space composed of one-dimensional and two-dimensional dimensionality reduction features is stained using the Morse function value of the cells. (d) Comparison of Morse function distribution of the four groups of cells. The label color is consistent with Figure (b). (e) Reeb plot of cellular connectivity network extracted from MNS data

[125] . The largest labeled point represents the initial equivalent node of qNSC, and the smallest labeled point represents the intermediate equivalent node of NB. (f) Cell distribution in the spectral embedding space, stained and labeled by pseudo-times of cells.

[0063] Figure 7To infer cell lines and pseudo-times during mouse embryo implantation cell differentiation using SCTRRG. (a) Based on diffusion mapping, the distribution of mouse embryonic cells in the low-dimensional space constitutes the second and fourth dimensionality reduction features in the initial stage of implantation-labeled cell differentiation (1c, 2c, 4c, 8c, 16c, 32c and 64c). (b) Cell lineage subgroups formed by 9 cell types were analyzed by the SCTRRG algorithm (1c+2c, 4c, 8c, 16c, 32c-TE, 64c-TE, ICM, EPI and PE). (c) Distribution of 442 mouse embryonic cells stained with Morse function values ​​in the low-dimensional space of the spectral embedding. (d) Comparison of Morse function distributions of the 9 cell groups. Label colors are consistent with Figure b. (e) Reeb graphs were extracted from the connection network of the MLP data

[16] . The largest labeled point represents the initial equivalent node belonging to the 1C (1-cell) stage, and the smallest labeled point represents the intermediate equivalent node belonging to the 16C (16-cell) stage. Label colors are consistent with Figure b. (f) Cell distribution in the spectral embedding space, stained and labeled by pseudo-times of cells.

[0064] Figure 8 Comparison of pseudo-temporal inference algorithms;

[0065] Figure 9 This is a flowchart of the method of the present invention. Detailed Implementation

[0066] The following is in conjunction with the appendix Figure 1-9 The preferred embodiments of the present invention will be described herein. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0067] Example 1:

[0068] A single-cell pseudo-trajectory topology recognition method based on Reeb graphs, comprising:

[0069] S1: Preprocess the single-cell transcription dataset to obtain the single-cell gene expression matrix;

[0070] S2: Based on the single-cell gene expression matrix of S1, an optimized model is constructed, and a single-cell connectivity network is further built for cell population identification. This includes:

[0071] S2.1: Constructing an optimization model:

[0072]

[0073] stA T y = y;

[0074]

[0075] y = (1,1,…,1)T

[0076] Where A corresponds to the adjacency matrix of a single-cell connected network, and X is the gene expression matrix of a single cell. * It is the nuclear norm, ||·|| 2,1 It is the sum of the Euclidean norms of all columns; the nonnegative parameter λ is used to balance the nuclear norm and the Euclidean norm; linear constraint A T y = y preserves the invariance of data translation transformations, where y is a vector of length n. This indicates that the i-th cell and the j-th cell are not neighboring structures;

[0077] The Alternating Direction Multiplier Method (ADMM) is used to solve the optimization model of S2.1. The K-Nearest Neighbors (KNN) algorithm is used to determine the nearest neighbor network G of the cell, and the selection of K is based on the BIC principle. It is agreed that non-zero values ​​in Aopt are 1, so that the 0-1 single-cell similarity network Aopt is obtained.

[0078] S2.2: Based on the similarity network Aopt of S2.1, construct a single-cell connected network for cell population identification:

[0079]

[0080] The 0-1 single-cell similarity network Aopt is transformed into a single-cell symmetric network Z, λ using the maximum value method. l Let be the l-th eigenvalue of L arranged in ascending order, and p be the index corresponding to the maximum absolute difference between two adjacent eigenvalues, used to identify the number of cell populations. Combined with spectral clustering algorithms, p cell populations are identified, namely C1, C2, ..., Cp. When Z is a disconnected network, dis represents the Euclidean distance function. Let i represent the i-th cell in cell group u and j-th cell in cell group v, respectively. Calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell groups. Where i * ,j * They are respectively Based on the cell sequence throughout the network, Z is progressively updated into a single-cell connected network; for cases where there are edges between two cell groups, Z... ij =1, which means there is a connection between the two cell groups, and we can further obtain the cell group connectivity network CC between the cell groups.

[0081] S3: Building upon S2, and combining the Morse function distribution and equivalence class definition of cells, extract the Reeb map for single cells. This includes:

[0082] S3.1: The distribution of Morse function in a single cell is as follows:

[0083]

[0084] Where u(c) i c j The distance f(i) is the geometric distance between the i-th and j-th cells in a disconnected network Z, approximated by Dijkstra's algorithm. i Normalization to To maintain the invariance of the f value with respect to the grid scale, where n is the number of cells; the literature shows that the f function can be approximated as a Morse function.

[0085] S3.2: Extracting Reebok images from single cells f :

[0086]

[0087] R f =G(V) r E r )

[0088]

[0089] in, and Let represent the average geodesic distances of the u-th and v-th cells in the i-th cell group, respectively, and d represent the error. An equivalence relation within the same cell population is defined as the difference in the Morse function being within a very small error.

[0090] Based on equivalence classes The definition of Reeb graph R is based on the extraction of Reeb graphs from single-cell connected networks. f R f V r Nodes are determined by representative nodes in equivalence classes, R f Edge set E r The existence of an edge is determined by the connectivity of cells in two equivalence classes ([p], [q]) within the connected network Z. If there exists a connection between cell s and cell t, then A s，t ≠0, E r There are edges connecting the equivalent nodes [p] and [q].

[0091] S4: Building upon S3, this section combines Morse function extrema to identify key nodes in the Reeb graph, infer intermediate stages and initial / final states of cell differentiation, measure the distance influence of cells in the network and equivalence classes, and identify cell pseudotrajectories. This includes:

[0092] S4.1: Sort the nodes in the Reeb graph of S3:

[0093]

[0094] Among them, V B It is a cell population located at the boundary of the cell population connected network C, with a degree of 1; based on the continuity and equivalence class definition of the Morse function, the largest equivalence class of the Morse function corresponds to R. f The initial equivalent node v on [0] This corresponds to the initial or final state of cell differentiation; the equivalent cell class v with the smallest Morse function. [m] Located in the entire R f The middle position corresponds to the vicinity of the middle position of cell differentiation, as a possible branching state;

[0095] S4.2: Based on the initial equivalent node v [0] And S3, the entire Reeb diagram R f Calculate the shortest path length from each node in the Reeb graph to the initial equivalent node, denoted as Len. Based on the size of Len, calculate the shortest path length for R. f n r Sort the equivalent nodes:

[0096] T [i] =Len(v [i] ,v [0] )

[0097] S4.3: Measuring cell distance:

[0098]

[0099] Where d(s,j) represents the shortest path length between the s-th and j-th cells in the single-cell connection network Z, D ij The metric is the average distance from cell j to the set of cells in the [i] equivalence class. The metric is the average distance from cell j in cell population c to the cell set of equivalence class [i]. exp is the exponential function, var represents the variance function, and D is the normalized value. ij , It is uniformly denoted as the average distance between cells. The effects of the two distances on cells were characterized from the perspectives of the overall cell network and the local networks of cells within the cell population;

[0100] S4.4: Combining the influence of intercellular distance and the order of cell nodes in the Reeb plot, pseudo-time points for all cells are obtained, thus achieving single-cell pseudo-trajectory inference:

[0101]

[0102] If we consider the influence of the overall cell network on cell v, then i = 1, 2, ..., n r v = 1, 2, ..., n; if we consider the local networks of cells in the cell population, then v∈C u .

[0103] Example 2:

[0104] A single-cell pseudo-trajectory topology recognition device based on Reeb graphs is provided. The device is constructed based on the above-mentioned single-cell pseudo-trajectory topology recognition method based on Reeb graphs, and realizes single-cell pseudo-trajectory topology recognition.

[0105] Example 3:

[0106] A single-cell pseudo-trajectory topology recognition device based on Reeb graphs, the control device including a processor and a memory, wherein the memory stores instructions executable by the processor, the instructions being executed by the processor to enable the processor to perform the control method described in Embodiment 1.

[0107] Example 4:

[0108] A computer-readable storage medium storing a computer program, the computer program being loaded by a processor and executed by the Reeb graph-based single-cell pseudo-trajectory topology recognition method described in Example 1.

[0109] Simulation experiment:

[0110] Obtain single-cell datasets. These include two simulated datasets from PROSSTT, a real scRNA-seq dataset, and the q-PCR dataset in Table 1. PROSSTT simulates single-cell RNA-seq transcriptomes of cells differentiating along a complex phylogenetic tree structure. The first simulation had two branches with 80 cells, and the second simulation had three branches with 150 cells.

[0111] The scRNA-seq dataset contains 272 mouse neural stem (MNS) cells with two branches. The single-cell q-PCR dataset includes 442 late-stage MLP mouse cells with three known branching structures. A branch represents a pathway from the initial cell to one of the endpoints determining the cell's fate.

[0112] Table 1 Single-cell dataset

[0113]

[0114] The effectiveness of the SCTRRG algorithm in identifying single-cell differentiation trajectories was investigated. The SCTRRG algorithm was applied to two simulated single-cell datasets and two real single-cell datasets. Specifically, the effectiveness of the SCTRRG algorithm was validated using two simulated single-cell datasets, and its effectiveness was validated using two real datasets. The single-cell experimental datasets were used to verify that the SCTRRG algorithm, combined with prior gene expression information, could identify different cell lines and branches.

[0115] 1. Simulation of two branch campuses

[0116] First, we used the single-cell pseudotrack simulation software PROSSTT to generate simulation data S1 of a bibranch structure, which mainly contained the expression of 500 genes in 80 cells.

[0117] exist Figure 4 In (a), the S1 data mainly includes three cell types (labeled as 1, 2, and 3) and two differentiation directions. Based on the spectral clustering method of cell connectivity networks, the SCTRRG algorithm consensus distinguishes six cell types (C1, C2, C3, C4, C5, and C6) and accurately identifies the cell lineages of two branches (e.g., ...). Figure 4 (b) is shown. Figure 4 (c) shows the positional distribution and Morse function distribution of the 80 units in the spectral embedding space.

[0118] It can be observed that the S1 data in the spectral embedding dimensionality reduction space is closer to several lines than in the diffusion mapping dimensionality reduction space, without significant scattering fluctuations. Furthermore, based on the distribution of Morse function values ​​for each monomer in the spectral embedding low-dimensional space, it can be seen that almost all monomers with Morse function values ​​close to 1 (the yellow portion) are distributed at the boundary positions of the spectral embedding space.

[0119] From the perspective of Morse distribution in cells, cells located in the center or branching stages of the cell line (C2, C3, C4) have higher Morse values ​​than those in the root (C1) and some leaves (C5, C6). Cell branch C2 contains the minimum Morse function value for its corresponding single cell, while leaf cell C5 contains the maximum Morse function value for its corresponding single cell. Figure 4 (d)).

[0120] 2. Multi-branch simulation

[0121] Then, S2 was generated using the single-cell pseudotrack simulation software PROSSTT. S2 is a simulation data set from PROSSTT containing three branching structures, mainly including the expression of 500 genes in 150 cells. Figure 5 As shown in (a), the S2 data mainly includes 5 cell types (denoted as: 1, 2, 3, 4, 5) and 3 differentiation directions. Based on the spectral clustering method, the SCTRRG algorithm identified 9 cell types (represented as C1, C2, C3, C4, C5, C6, C7, C8, and C9), and accurately identified cell lineages containing 3 branches (e.g., Figure 5 (b) is shown). In Figure 5 (c) and Figure 5In (d), the spatial distribution of the reproductive embedding locations of 150 cells can be compared with the Morse function distribution of single cells and 9 cells. It was found that the Morse function values ​​of the reproductive embedding boundary locations of almost all cells were close to 1 (yellow portion), mainly including the root (C1) and leaf (C6, C8, C9) cell lines. Morse function values ​​close to 0 (red) were located in the middle or branching stages of the cell lines (C3, C4, C5, and C7). Among these, the branching cell populations C3 and C5 contained single cells corresponding to the minimum Morse function value, and the leaf cell population C8 contained single cells corresponding to the maximum Morse function value. Figure 5 (d) Furthermore, the SCTRRG algorithm extracts the topological Reeb graph of the cellular network in S2 and marks the distribution of the initial and intermediate nodes of the Reeb graph as the maximum (green) and minimum (blue) nodes of the Morse function. Figure 5 (e)).

[0122] In the absence of prior biological information, the initial node for cell differentiation can be selected in VB. The node with the highest cell equivalence class value in the C9Morse function can be chosen as the initial node, using the cell point at the end of the branch as the initial point, combined with the distance of a single cell connected to the network path; single-cell pseudomoments are determined, and dimensionality reduction through downward embedding is rendered in space. Figure 5 (f)). With Figure 5 (d) and Figure 5 (f) Comparison: Along the trajectory of cell differentiation, from the initial point to the intermediate stage (corresponding to cell C2 in the spatial location of cell distribution), the average Morse function value of the corresponding cell population changes from 0.60 to 0.10 (C9→C5). From the intermediate stage to the three final stages, the median Morse function value changes from 0.10 to 0.81 (C5→C1). The variation within the ranges of 0.10–0.80 (C5→C6) and 0.1–0.84 (C5→C8) indicates that the average change in Morse function within the cell population increases with the early stage of differentiation and decreases with the increase of the overall differentiation process.

[0123] 3. Bifurcation recognition of mouse neural stem cell differentiation

[0124] Furthermore, the SCTRRG algorithm was used to identify and infer cell lineages and pseudotimes during mouse neural stem cell differentiation. The MNS dataset mainly contains a single-cell dataset of 43,309 genes expressed in 272 neural stem cells. After preprocessing with highly variable genes (Fano factor > 0.5, mean value between 0.01 and 3.5), 8,846 genes were selected. Figure 6As shown in (a), the MNS dataset mainly contains four labeled cell populations (denoted as 1, 2, 3, and 4, respectively). In the low-dimensional space of the diffusion map, the partial distributions of labeled 3-cell and labeled 4-cell populations overlap in the low-dimensional space, which increases the difficulty of inferring the direction of cell differentiation. Based on the differential expression and prior information of eight marker genes (Cd24a, Tubb3, Sox9, Id2, Sox10, Mbp, Egr1, Fos) in mouse neural stem cells, the initial cell population should be qNSC (high expression of Sox9). Based on the different expressions of the eight marker genes, the SCTRRG algorithm identified four cell populations, labeled as cell populations qNSC, NB, Oli, and aNSC, and accurately identified the cell lineages of two branches. Figure 6 (b). In Figure 6 (c) and Figure 6 In (d), the spatial distribution, single cell, and cell population Morse function distribution of the germinal embedding sites of 272 neural stem cells can be compared. It can be observed that cells with Morse function values ​​close to 1 (yellow portion) are mostly distributed in the spatial location on the left (Oli cell group). From an overall distribution perspective, the root cell lineage (cell qNSC) and leaf cell lineages (cell populations Oli and aNSC) are located below the Morse function value in the clade (NB) cell group. The leaf cell population Oli contains [data missing]. The single cell corresponding to the maximum Morse function value ([data missing]). Figure 6 (d) Furthermore, the SCTRRG algorithm extracted the Reeb graph of the mouse neural stem cell connectivity network. The initial and intermediate nodes of the Reeb graph can be labeled as the points of maximum (pink) and minimum (brown) values ​​of the Morse function, respectively. Figure 6 (e)). Comparison Figure 6 (a) and Figure 6 (e) It can be observed that some low-dimensional features of neural stem cells overlap in the third and fourth label categories, making it impossible to distinguish the trends of cell differentiation. In the Reeb plot extracted by the SCTRRG algorithm, the connections between cell groups and the overall trend of cell differentiation are clearly visible. Finally, the pseudo-motion of individual cells is determined by combining the path distance of the mouse neural stem cell connectivity network and presented in the spectral embedding space, as shown below. Figure 6 As shown in (f). Comparison Figure 6 (d) and Figure 6(f) Along the cell differentiation trajectory, the mean Morse function of mouse neural stem cells changed from the initial qNSC phase to the two phases (NB and aNSC cell groups) from 0.46 to 0.37 (qNSC→NB) and 0.46 to 0.62 (qNSC→aNSC), respectively. From the NB phase to the final stage, the median Morse function changed from 0.37 to 0.99 (NB→Oil), respectively. This indicates that along the branching movement of mouse neural stem cells, the mean Morse function underwent a process of first decreasing and then increasing.

[0125] In summary, numerical results demonstrate that the SCTRRG algorithm can infer the two branch trajectories during mouse neural stem cell differentiation and determine the pseudomoments for each cell. The changing trend of the average degree of the cell population Morse function is related to the cell differentiation process.

[0126] 4. Multibranch recognition before mouse embryo implantation

[0127] Finally, the SCTRRG algorithm can be used to identify and infer multi-branch trajectories and pseudo-times during preimplantation cell differentiation in mouse embryos. The MLP dataset mainly consists of a single-cell dataset containing 48 genes expressed in 442 cells. Since the number of cells is greater than the number of genes, no gene selection was performed. The MLP dataset mainly contains cell populations at seven stages (labeled as 1C, 2C, 4C, 8C, 16C, 32C, and 64C).

[0128] exist Figure 7 In the low-dimensional space of the diffusion map shown in (a), cells exhibit multiple differentiation directions. Cells at the intermediate 8C, 16C, and 32C stages are relatively disordered, and the connections between cells are indistinguishable. Figure 7 (b) and Figure 7 In (c), based on the maximum value of the Morse function corresponding to the initial cell differentiation, the initial cells of the cell group can be at the 1c+2c stage. The SCTRRG algorithm can obtain nine cell types: 1c+2c, 4c, 8c, 16c, 32c-TE, 64c-TE, ICM, EPI, and PE stages, accurately identifying the three branching lineages in the cell differentiation process. The initial cell population during mouse embryonic development is the initial stage of cell differentiation. Figure 7 (c) and Figure 7 In (d), 442 units of the Morse function distribution can be compared. The Morse function value is close to 1 (yellow part). Most units decrease by three units in the embedding space: 1c+2c, EPI and PE. From the overall distribution, the cells located in the cell line branch (16c) with a Morse average of close to 0.1 are much lower than the cell line root (64c-TE, EPI, PE) and leaf cells (1c+2c).

[0129] Furthermore, the SCTRRG algorithm extracted the topological Reeb graph of the connectivity network during the preimplantation phase of mouse embryos. The nodes that correspond to the initial and intermediate nodes of the Reeb graph, representing the nodes with the largest Morse function (red) and the smallest Morse function (green), can be labeled as follows: Figure 7 (e)). With Figure 7 (a) In contrast, the Reeb plot clearly shows the connections between the 16C and 32C-TE cells and the ICM.

[0130] Finally, by combining the path distances of cell communication networks during the preimplantation phase of mouse embryos, the pseudomoments of individual cells can be determined, and... Figure 7 (f) is presented in the spectral embedding space. (Comparison) Figure 7 (d) and Figure 7 (f) Following the trajectory of preimplantation cell differentiation in mice, from the initial stage 1c+2c to the intermediate stage (16c), the mean Morse function of the cells changes from 0.80 to 0.11. From the intermediate stage to the three endpoints (PE, EPI, 64c-TE), the median Morse function changes from 0.11 to 0.64 (16C→PE), from 0.11 to 0.64 (16C→EPI), and from 0.11 to 0.44 (16C→64C-TE), respectively. This indicates that the mouse embryo is in the process of cell differentiation before implantation. The average Morse function of the cell population as a whole first decreases and then increases along each branch.

[0131] In summary, the SCTRRG algorithm can infer three complex branching trajectories during the preimplantation phase of mouse embryos and determine the pseudomoments for each cell. The variation trend of the Morse function is related to the cellular evolutionary process of different branching trajectories.

[0132] 5. Comparison of Pseudo-trajectory algorithms

[0133] To verify the accuracy of the algorithm, the classic single-cell pseudo-trajectory recognition algorithms DPT and Monocle2 were compared with the SCTRRG algorithm. Simultaneously, the Kendall rank correlation coefficient between the cell's prior order (simulated data) and cell differentiation stage (real data) and the pseudo-time point was calculated as an evaluation index for the accuracy of pseudo-trajectory judgment. Figure 8 As shown, at pseudo-time points in the three datasets, the SCTRRG algorithm outperforms DPT and Monocle2 inference results. In simulated data, the results of the SCTRRGp and SCTRRGcp algorithms are not significantly different, while in real-world data, the SCTRRGp algorithm outperforms the other algorithms.

[0134] Therefore, when using the SCTRRG algorithm to identify single-cell pseudosites, it is recommended to infer the pseudomoments of cells by considering path distances within the cell network. For single-cell data lacking information on cell differentiation stages, Kendall's rank correlation coefficient cannot be calculated, such as in the mouse neural stem cell dataset. In summary, the proposed method can effectively infer two- and multi-branch lineage processes in single cells across various simulated and real-world datasets. Furthermore, its inference accuracy on these three datasets outperforms two other classic algorithms, further demonstrating the rationality and effectiveness of the proposed method.

[0135] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A method for topological recognition of single-cell pseudo-trajectories based on Reeb graphs, characterized in that: include: S1: Preprocess the single-cell transcription dataset to obtain the single-cell gene expression matrix; S2: Based on the single-cell gene expression matrix of S1, an optimized model is constructed, and a single-cell connectivity network is further constructed for cell population identification; S2.1: Constructing an optimization model: ； Where A corresponds to the adjacency matrix of a single-cell connected network. X This is a single-cell gene expression matrix. It is the nuclear norm. It is the sum of the Euclidean norms of all columns, a non-negative parameter. Used to balance the nuclear norm and Euclidean norm; linear constraints A T y = y preserves the invariance of data translation transformation. y It is a vector of length n. Indicates the first i The cell and the first j The individual cells are not neighboring structures; The Alternating Direction Multiplier Method (ADMM) is used to solve the optimization model of S2.

1. The K-Nearest Neighbors (KNN) algorithm is used to determine the nearest neighbor network G of the cell, and the selection of K is based on the BIC principle. It is agreed that non-zero values ​​in Aopt are 1, so that the 0-1 single-cell similarity network Aopt is obtained. S2.2: Based on the similarity network Aopt of S2.1, construct a single-cell connected network for cell population identification: ； The 0-1 single-cell similarity network Aopt is transformed into a single-cell symmetric network Z using the maximum value method. It is the first in L arranged in order of size. l 1 eigenvalue, p The index corresponding to the maximum absolute difference between two adjacent feature values ​​is used to identify the number of cell populations. Combined with spectral clustering algorithms, p cell populations are identified, namely C1, C2, ..., Cp. When Z is a disconnected network, dis represents the Euclidean distance function. Let i represent the i-th cell in cell group u and j-th cell in cell group v, respectively. Calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell groups. ,in They are respectively Based on the cell order throughout the network, Z is progressively updated into a single-cell connected network; this applies to cases where there are edges between two cell groups. Then there are edges connecting the two cell groups, which further leads to the cell group connectivity network between cell groups. CC ; S3: Based on S2, extract the Reeb map of a single cell by combining the Morse function distribution and equivalence class definition of the cell. S4: Building upon S3, this method combines the maximum and minimum values ​​of the Morse function to identify key nodes in the Reeb graph, infer the intermediate stages and initial and final states of cell differentiation, measure the distance influence of cells in the network and equivalence classes, and identify cell pseudotrajectories.

2. The method for single-cell pseudo-trajectory topology recognition based on Reeb graphs according to claim 1, characterized in that: S3 includes: S3.1: The distribution of Morse function in a single cell is as follows: ； in, The geometric distance between the i-th and j-th cells in a disconnected network Z is approximated using Dijkstra's algorithm. The function... Normalization to In order to maintain f The value is invariant relative to the grid scale, where n is the number of cells; S3.2: Extracting Reebok images from single cells R f : ； in, and Let u and v represent the average geodesic distances of the u-th and v-th cells in the i-th cell group, respectively. d Indicates error. An equivalence relation within the same cell population is defined as the difference in the Morse function being within a very small error. Based on equivalence classes The definition of Reeb graph extraction based on single-cell connected networks. R f , R f V r Nodes are determined by representative nodes in the equivalence class. R f Connecting edges E r The existence of an edge is determined by the connectivity of cells in the two equivalence classes [p] and [q] in the connected network Z. If there exists a connection between cell s and cell t, then... , E r There are edges connecting the equivalent nodes [p] and [q].

3. The method for single-cell pseudo-trajectory topology recognition based on Reeb graphs according to claim 2, characterized in that: S4 includes: S4.1: Sort the nodes in the Reeb graph of S3: ； in, V B It is located in the cell population connectivity network CC The cell population at the boundary has a degree of 1; based on the continuity and equivalence class definition of the Morse function, the largest equivalence class of the Morse function corresponds to... R f Initial equivalent nodes on v [0] This corresponds to the initial or final state of cell differentiation; the equivalent cell class with the smallest Morse function. v [m] Located in the whole R f The middle position corresponds to the vicinity of the middle position of cell differentiation; S4.2: Based on initial equivalent nodes v [0] And the entire Reeb map of S3 R f Calculate the shortest path length from each node in the Reeb graph to the initial equivalent node, denoted as Len. Based on the size of Len, [the following steps are performed]. R f n r Sort the 10 equivalent nodes: ； S4.3: Measuring cell distance: ； in, d ( s,j () represents the shortest path length between the s-th and j-th cells in a single-cell connection network Z. D ij The metric is the average distance from cell j to the set of cells in the [i] equivalence class. The metric measures the average distance from cell j in cell population c to the cell set of equivalence class [i]. exp is the exponential function, var represents the variance function, and normalization is used. It is uniformly denoted as the average distance between cells. The effects of the two distances on cells are characterized from the perspectives of the overall cell network and the local network of cells within the cell population. S4.4: Combining the influence of intercellular distance and the order of cell nodes in the Reeb plot, pseudo-time points for all cells are obtained, thus achieving single-cell pseudo-trajectory inference: ； If we consider the effect of the overall cellular network on V cells, then If we consider the local networks of cells within a cell population, then .

4. A single-cell pseudo-trajectory topology recognition device based on Reeb graphs, characterized in that: This method is constructed based on the Reeb graph-based single-cell pseudo-trajectory topology recognition method according to any one of claims 1-3.

5. A single-cell pseudo-trajectory topology recognition device based on Reeb graphs, characterized in that: The identification device includes a processor and a memory, wherein the memory stores a computer program, which is loaded and executed by the processor to implement the single-cell pseudo-trajectory topology identification method based on Reeb graphs as described in any one of claims 1-3.

6. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which is loaded and executed by a processor to implement the single-cell pseudo-trajectory topology recognition method based on Reeb graphs as described in any one of claims 1-3.

Images