Reeb graph-based single-cell pseudo-trajectory topology identification method

Abstract

The invention belongs to the technical field of single cell data analysis, and particularly relates to a single cell pseudo-trajectory topology identification method based on a Reeb graph. Comprising 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 performing cell population recognition; s3, extracting a single-cell Reeb map in combination with Morse function distribution and equivalence definition of the cells; and S4, identifying key nodes in the Reeb graph by combining the Morse function extreme value, deducing the middle stage and initial and end states of cell differentiation, measuring the distance influence of the cells in the network and the equivalence class, and identifying the cell pseudo trajectory. According to the method, differentiation tracks of cells in multiple directions are extracted on the basis of a topological structure Reeb graph of a cell connected network. The algorithm framework can convert the cellular network into a Reeb graph with fewer equivalent nodes, so that the calculation complexity is reduced; the nature of the Morse function is provided to determine the critical state of cell differentiation as compared to the initial time determination that completely depends on a priori information.

Description

Technical Field

[0001] The invention belongs to the technical field of single cell data analysis, and in particular relates to a single cell pseudo trajectory topology recognition method based on a Reeb graph. Background Art

[0002] Recent technical advances in single-cell analysis allow for the measurement of single-cell transcriptional profiles during many stages of proliferation and differentiation of an organism. The challenge is to design computational techniques and modeling approaches that can automatically reconstruct cell differentiation pathways, order individual cells based on their transcriptome profile, and identify trajectory behaviors along branching lineages.

[0003] Since single-cell data have characteristics such as high statistical noise level, high-dimensional gene expression space and nonlinearity, pseudo-trajectory inference modeling based on scRNA-seq data to analyze dynamic processes such as cell differentiation, cell fate development and decision-making can be 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 pseudo-time to each cell accordingly; the Monocle2 algorithm uses MST to learn cell trajectories, updates the cell position by moving the cell to the nearest vertex in the MST until the cell trajectory and final position are stable, and finally obtains the pseudo-time of the cell based on the geodesic distance from the cell to the root vertex along the MST. The TSCAN algorithm and the Waterfall algorithm construct MST based on cell clusters, sort cells by orthogonal projection methods, and reduce the total computational complexity of the MST algorithm as a whole. The Slingshot algorithm uses MST to construct trajectories, proposes a joint principal curve algorithm to construct a smooth curve, and then projects the cell onto the smooth curve instead of the edge of the MST to calculate the cell pseudo-time.

[0005] The second category is trajectory inference methods based on graph models. The diffusion pseudo-time (DPT) algorithm uses the weighted k-nearest neighbor (KNN) algorithm to construct the cell's trajectory and introduces ten diffusion algorithms to calculate the pseudo-time of cells in the "diffusion mapping space". The PAGA algorithm calculates the neighborhood relationship by selecting appropriate low-dimensional representations and related distance metrics, represents high-dimensional gene expression data as KNN graphs, uses the Louvain algorithm to divide the vertices of various KNN graphs, combines the high-confidence paths generated in the PAGA graph with the distance metric based on random walks on the graph, and calculates the distance to the root cell to sort 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 pseudo-time; the URD algorithm uses the 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 the RNA velocity-assisted trajectory inference method. The VeTra algorithm uses RNA velocity vectors to construct multiple directed graphs based on the cell's k nearest neighbors to determine the cell's transition state, identifies weakly connected branches by constructing independent cell transition paths, and finally clusters these transition paths to obtain the trajectory together. The Cytopath algorithm combines the inter-cell transition probability matrix and cell state to construct multiple simulation trajectories, and uses RNA velocity to infer the root and terminal states of the cell differentiation trajectory.

[0007] However, when performing trajectory inference, it is usually based on prior information to obtain a certain cell or a certain type of cell as the starting point of the trajectory, and the relative relationship between the remaining single cells or cell clusters and the starting point (such as the shortest path length) is used as additional information to construct the trajectory. Only the starting point information is used to infer all cell pseudo-times, and the rationality and stability of trajectory inference need to be improved. The Morse function is defined on the shape and manifold, and the Hessian matrix of its critical point (the first-order derivative of the Morse function is 0) is degenerate. The Reeb diagram reveals the topological structure of the manifold and the shape based on the Morse function. Summary of the invention

[0008] The purpose of the present invention is to propose an unsupervised algorithm, called SCTRRG (Single Cell Topology Recognition and Reeb Graph), for reconstructing high-resolution developmental paths along topological structures. This method uses multidimensional single-cell data as input, constructs network relationships between cells, defines Morse functions to identify key starting fates and bifurcated cells, and infers trajectories in multi-branched cell fates based on their topological structures through Reeb graphs. Through two simulated data sets, it was demonstrated that SCTRRG can successfully identify their different multi-branched structures with reasonable starting points and branching points. SCTRRG was further applied to two published single-cell RNA-seq data sets, and SCTRRG accurately restored known cell lineages in mouse embryonic and neural stem cell development.

[0009] In order to achieve the above technical effects and the above objectives, the present invention provides a single cell pseudo trajectory topology recognition method based on a Reeb graph in a first aspect, comprising:

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

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

[0012] S3: Based on S2, the single-cell Reeb graph is extracted by combining the cell’s Morse function distribution and equivalence class definition;

[0013] S4: Based on S3, combined with the Morse function maximum, the key nodes in the Reeb graph are identified, the intermediate stages and initial and final states of cell differentiation are inferred, the distance influence of cells in the network and equivalence class is measured, and cell pseudo-trajectories are identified.

[0014] Preferably, S2 comprises:

[0015] S2.1: Build optimization model:

[0016]

[0017] sA T y=y;

[0018]

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

[0020] Among them, A corresponds to the adjacency matrix of the single-cell connectivity network, X is the single-cell gene expression matrix, ||·|| * is the nuclear norm, ||·|| 2,1 is the sum of the Euclidean norms of all columns, and the non-negative parameter λ is used to balance the nuclear norm and the Euclidean norm; the linear constraint A T y=y maintains the invariance of data translation transformation, y is a vector of length n, Indicates that the i-th cell and the j-th cell are not neighbor structures;

[0021] The alternating direction multiplier method ADMM was used to solve the optimization model of S2.1, and the K-nearest neighbor algorithm KNN was used to determine the cell's neighbor network G. The selection of K was based on the BIC principle. The non-zero value in Aopt was agreed to be 1, and a 0-1 single-cell similarity network Aopt was obtained.

[0022] S2.2: Based on the similarity network Aopt of S2.1, a single-cell connectivity network is constructed to identify cell populations:

[0023]

[0024] The 0-1 single-cell similarity network Aopt is converted into a single-cell symmetric network Z,λ by the maximum method l is the lth eigenvalue of L in order of size, p is the subscript corresponding to the maximum absolute value difference between two adjacent eigenvalues, which is used to identify the number of cell groups. Combined with the spectral clustering algorithm, p cell groups are identified, namely C1, C2, ... Cp; when Z is a non-connected network, dis represents the Euclidean distance function, Respectively represent the i-th cell in the u-type cell group and the j-th cell in the v-type cell group, and calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell group where i * ,j * They are In the corresponding cell order of the entire network, Z is gradually updated to a single-cell connected network; for the case where there is an edge between two cell groups, Z ij =1, then there is an edge between the two cell groups, and further a cell group connectivity network CC between the cell groups is obtained.

[0025] Preferably, S3 includes:

[0026] S3.1: The single-cell Morse function distribution is:

[0027]

[0028] Where u(c i , c j ) is the geometric distance between the i-th cell and the j-th cell in the non-connected network Z, which is approximated by the Dijkstra algorithm. i Normalized to To maintain the invariance of the f value with respect to the grid scale, n is the number of cells;

[0029] S3.2: Extracting single-cell Reeb maps f :

[0030]

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

[0032]

[0033] in, and They represent the average geodesic distances of the uth and vth cells in the i-th cell population, d represents the error, The equivalence relationship within the same cell population is defined as the difference of the Morse function within a small error;

[0034] Based on equivalence class According to the definition of single cell connectivity network, Reeb graph R is extracted f , R f V r The nodes are determined by the representative nodes in the equivalence class, R f Edge Set E rWhether the edge exists is determined by the connectivity of cells in the two equivalence classes [p] and [q] in the connected network Z. If there is a connection between s cell and t cell, A s，t ≠0, E r There are edges between 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 is a cell cluster located at the boundary of the cell cluster connected network CC, 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] , corresponding to the initial state or terminal state of the cell differentiation process; the equivalent cell class v with the smallest Morse function [m] Located throughout R f The middle position of the cell corresponds to the middle position of cell differentiation.

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

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

[0041] S4.3: Measuring cell distances:

[0042]

[0043] Where d(s,j) represents the shortest path length between the sth and jth cells in the single-cell connection network Z, and D ij It measures the average distance from the j cell to the [i] equivalent class cell set. It measures the average distance from cell j in cell population c to the equivalent cell set [i]. exp is the exponential function, var is the variance function, and the normalized D ij , The average distance between cells The effects of the two distances on cells are characterized from the perspective of the overall cell network and the local network of cells in the cell population;

[0044] S4.4: Combining the influence of inter-cell distance and the order of cell nodes in the Reeb graph, the pseudo-times of all cells are obtained, that is, the single-cell pseudo-trajectory inference is realized:

[0045]

[0046] If the influence of the whole cell network on v cell is considered, then i=1,2,...,n r ,v=1,2,...,n; if we consider the local network of cells in the cell group, then

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

[0048] The third aspect of the present invention provides a single-cell pseudo-trajectory topology recognition device based on a Reeb graph, wherein the control device comprises a processor and a memory, wherein a computer program is stored in the memory, and the computer program is loaded and executed by the processor to implement the single-cell pseudo-trajectory topology recognition method based on a Reeb graph.

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

[0050] Compared with the prior art, the present invention has the following beneficial effects:

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

[0052] 2. Based on the existing network or tree to measure the relationship between cells, the equivalence class of cells is defined, and the REEB diagram is proposed to simulate the positional relationship during cell differentiation.

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

[0054] 4. SCTRRG can accurately identify the lineage structure of cell differentiation corresponding to the four data sets, infer its pseudo-time and each branch trajectory. Compared with other classic methods, the pseudo-trajectory accuracy recognition effect is better. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0056] In the attached picture:

[0057] Figure 1 The framework diagram of the SCTRRG method of the present 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) Six types of Morse function (MF) distribution. (D) Cell lineage. MF is defined based on the connected network structure to identify the initial cluster (C0). (E) The Reeb diagram reflects the basic topological structure of the connected cell network based on MF. In the two-dimensional spectral embedding space, the method (F) unit pseudo-time (normalized to 0-1) has an optional initial state (the number of boundary component nodes is the largest) and a median state (the number of branch component nodes is the smallest) in the development trajectory.

[0058] Figure 2 Schematic diagram of the transition from a single-cell similarity network to a sparsely connected network.

[0059] Figure 3 Compute equivalence classes for the Morse function.

[0060] Figure 4 As shown in the figure: (a) Based on the diffusion mapping dimensionality reduction method, the distribution of cells in the low-dimensional space composed of the second and third dimensionality reduction features is marked by the initial classification labels (1, 2, 3). (b) Six cell subgroups and cell lineage structures are identified by the SCTRRG algorithm. (c) Based on the spectral embedding dimensionality reduction method, the distribution of cells in the spectral embedding space composed of two-dimensional dimensionality reduction features and three-dimensional dimensionality reduction features is colored by the Morse function value of the cell. (d) Comparison of the Morse function distribution of 6 groups of cells. The label color is consistent with Figure b. (e) Reeb diagram of the connected cell network in the S1 data. The maximum marked point represents the initial equivalent node of the cell group C1, and the minimum marked point represents the intermediate equivalent node of the cell group C2. The marking color is consistent with C1, C2, C3, C4 and C5 in Figure B, and orange represents C6. (f) The distribution of cells in the spectral embedding space is stained and marked by the pseudo-moments of the cells.

[0061] Figure 5Identification of cell lineages and pseudo-times of multi-branch simulated data S2 for 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 the initial classification labels (1, 2, 3, 4, 5). (b) Nine cell subpopulations and cell lineage structures are identified by 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 by the Morse function value of the cell. (d) Comparison of Morse function distributions of 9 groups of cells. The label color is consistent with Figure b. (e) Reeb graph of multi-branch simulated data (S2) connecting the cell network. The maximum labeled point represents the initial equivalent node belonging to C1, and the minimum labeled point represents the intermediate equivalent node belonging to C2. The label color is consistent with Figure b. (f) Cell distribution in the spectral embedding space, colored by the pseudo-time of the cell.

[0062] Figure 6 SCTRRG is used to infer cell lineages and pseudo-times during the differentiation of mouse neural stem cells. (a) Based on the diffusion map 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 marked with the first known label (1, 2, 3, 4). (b) The cell lineages formed by the 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 by the Morse function value of the cell. 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 with the Morse function value of the cell. (d) Comparison of the Morse function distribution of the four groups of cells. The label color is consistent with Figure (b). (e) Reeb graph of the cellular connection network extracted from the MNS data

[125] . The maximum labeled point represents the initial equivalent node of the qNSC, and the minimum labeled point represents the intermediate equivalent node of the NB. (f) Distribution of cells in the spectral embedding space, labeled by the pseudo-time of the cells.

[0063] Figure 7To infer cell lineages and pseudo-times during differentiation of mouse embryonic implantation cells using SCTRRG. (a) Based on the diffusion map-based dimensionality reduction method, the distribution of mouse embryonic cells in the low-dimensional space consists of the second and fourth dimensionality reduction features at the initial stages of implantation cell differentiation (1c, 2c, 4c, 8c, 16c, 32c and 64c). (b) The cell lineage subsets formed by the 9 cell types were analyzed by the SCTRRG algorithm (1c+2c, 4c, 8c, 16c, 32c-TE, 64c-TE, ICM, EPI and PE). (c) The distribution of 442 mouse embryonic cells in the spectral embedding low-dimensional space colored by the Morse function value. (d) Comparison of the Morse function distribution of the 9 groups of cells. The label colors are consistent with those in Figure b. (e) Reeb graph extracted from the connection network of the MLP data

[16] . The maximum labeled point represents the initial equivalent node belonging to the 1C (1-cell) stage, and the minimum labeled point represents the intermediate equivalent node belonging to the 16C (16-cell) stage. The label colors are consistent with those in Figure b. (f) Distribution of cells in the spectral embedding space, labeled by the pseudo-time of the cells.

[0064] Figure 8 For pseudo-temporal reasoning algorithm comparison;

[0065] Fig. 9 The present invention is a flow chart of the method. DETAILED DESCRIPTION

[0066] The following is combined with Figure 1-9 The preferred embodiments of the present invention are described. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0067] Embodiment 1:

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

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

[0070] S2: Based on the single-cell gene expression matrix of S1, an optimization model is constructed to further construct a single-cell connectivity network for cell population identification. Including:

[0071] S2.1: Build optimization model:

[0072]

[0073] sA T y=y;

[0074]

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

[0076] Among them, A corresponds to the adjacency matrix of the single-cell connectivity network, X is the single-cell gene expression matrix, ||·|| * is the nuclear norm, ||·|| 2,1 is the sum of the Euclidean norms of all columns, and the non-negative parameter λ is used to balance the nuclear norm and the Euclidean norm; the linear constraint A T y=y maintains the invariance of data translation transformation, y is a vector of length n, Indicates that the i-th cell and the j-th cell are not neighbor structures;

[0077] The alternating direction multiplier method ADMM was used to solve the optimization model of S2.1, and the K-nearest neighbor algorithm KNN was used to determine the cell's neighbor network G. The selection of K was based on the BIC principle. The non-zero value in Aopt was agreed to be 1, and a 0-1 single-cell similarity network Aopt was obtained.

[0078] S2.2: Based on the similarity network Aopt of S2.1, a single-cell connectivity network is constructed to identify cell populations:

[0079]

[0080] The 0-1 single-cell similarity network Aopt is converted into a single-cell symmetric network Z,λ by the maximum method l is the lth eigenvalue of L in order of size, p is the subscript corresponding to the maximum absolute value difference between two adjacent eigenvalues, which is used to identify the number of cell groups. Combined with the spectral clustering algorithm, p cell groups are identified, namely C1, C2, ... Cp; when Z is a non-connected network, dis represents the Euclidean distance function, Respectively represent the i-th cell in the u-type cell group and the j-th cell in the v-type cell group, and calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell group where i * ,j * They are In the corresponding cell order of the entire network, Z is gradually updated to a single-cell connected network; for the case where there is an edge between two cell groups, Z ij =1, it can be considered that there is an edge between the two cell groups, and the cell group connectivity network CC between the cell groups is further obtained.

[0081] S3: Based on S2, combined with the cell's Morse function distribution and equivalence class definition, extract the single-cell Reeb graph. Including:

[0082] S3.1: The single-cell Morse function distribution is:

[0083]

[0084] Where u(c i , c j ) is the geometric distance between the i-th cell and the j-th cell in the non-connected network Z, which is approximated by the Dijkstra algorithm. i Normalized to To keep the f value invariant with respect to the grid scale, n is the number of cells; the literature shows that the f function can be approximately regarded as a Morse function.

[0085] S3.2: Extracting single-cell Reeb maps f :

[0086]

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

[0088]

[0089] in, and They represent the average geodesic distances of the uth and vth cells in the i-th cell population, d represents the error, The equivalence relationship within the same cell population is defined as the difference of the Morse function within a small error;

[0090] Based on equivalence class According to the definition of single cell connectivity network, Reeb graph R is extracted f , R f V r The nodes are determined by the representative nodes in the equivalence class, R f The edge set E r Whether the edge exists is determined by the connectivity of cells in the two equivalence classes ([p], [q]) in the connected network Z. If there is a connection between s cell and t cell, A s，t ≠0, E r There are edges between the equivalent nodes [p] and [q].

[0091] S4: Based on S3, combined with the Morse function maximum value, identify the 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 class, and identify cell pseudo-trajectories. Including:

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

[0093]

[0094] Among them, V B is a cell cluster located at the boundary of the cell cluster connected network CC, 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] , corresponding to the initial state or terminal state of the cell differentiation process; the equivalent cell class v with the smallest Morse function [m] Located throughout R f The middle position of , corresponding to near the middle position of cell differentiation, serves as a possible branching state;

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

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

[0097] S4.3: Measuring cell distances:

[0098]

[0099] Where d(s,j) represents the shortest path length between the sth and jth cells in the single-cell connection network Z, and D ij It measures the average distance from the j cell to the [i] equivalent class cell set. It measures the average distance from cell j in cell population c to the equivalent cell set [i]. exp is the exponential function, var is the variance function, and the normalized D ij , The average distance between cells The effects of the two distances on cells are characterized from the perspective of the overall cell network and the local network of cells in the cell population;

[0100] S4.4: Combining the influence of inter-cell distance and the order of cell nodes in the Reeb graph, the pseudo-times of all cells are obtained, that is, the single-cell pseudo-trajectory inference is realized:

[0101]

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

[0103] Embodiment 2:

[0104] A single-cell pseudo-trajectory topology recognition device based on a Reeb graph is constructed based on the single-cell pseudo-trajectory topology recognition method of the Reeb graph, and realizes single-cell pseudo-trajectory topology recognition.

[0105] Embodiment 3:

[0106] A single-cell pseudo-trajectory topology recognition device based on a Reeb graph, the control device comprising a processor and a memory, wherein the memory stores instructions executable by the processor, and the instructions are executed by the processor so that the processor can execute the control method described in Example 1.

[0107] Embodiment 4:

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

[0109] Simulation experiment:

[0110] Obtain single-cell datasets. Including two simulated datasets from PROSSTT, a real scRNA-seq dataset, and the q-PCR datasets in Table 1. PROSSTT simulates single-cell RNA-seq transcriptomes of cells that differentiate along a complex lineage tree structure. The first simulation has two branches and 80 cells, and the second simulation has three branches and 150 cells.

[0111] The scRNA-seq dataset contains 272 cells with two branches of mouse neural stem (MNS) cell differentiation. The single-cell q-PCR dataset includes 442 late-stage mouse cells with three known branch structures of the MLP. A branch represents a path from the initial cell to one of the end points of the entire cell fate decision.

[0112] Table 1 Single cell datasets

[0113]

[0114] The effectiveness of the SCTRRG algorithm in identifying single-cell differentiation trajectories was studied. The SCTRRG algorithm was applied to two simulated single-cell datasets and two real single-cell datasets. The effectiveness of the SCTRRG algorithm was verified using two single-cell simulated datasets, and the effectiveness of the algorithm was verified using two real datasets. The single-cell experimental dataset was used to verify that the SCTRRG algorithm combined with previous gene expression information to identify different cell lineages and branches.

[0115] 1. Simulation of two campuses

[0116] First, we used the single-cell pseudo-track simulation software PROSSTT to generate the simulation data S1 of the double-branch structure, which mainly contains the expression of 500 genes in 80 cells.

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

[0118] It can be found that the S1 data is closer to several lines in the spectrum embedding reduced dimensionality space than in the diffusion mapping reduced dimensionality space, and there is no large scatter fluctuation. At the same time, according to the distribution of the Morse function values ​​of each monomer in the spectrum embedding low-dimensional space, it can be seen that almost all monomers with Morse function values ​​close to 1 (yellow part) are distributed at the boundary position of the spectrum embedding space.

[0119] From the Morse distribution of cells, the Morse values ​​of cells located in the center or branch stage of the cell line (C2, C3, C4) are higher than those of the root (C1) and some (C5, C6) leaves of the cell line. The cell branch C2 contains the minimum Morse function value corresponding to a single cell, and the cell leaf C5 contains the maximum Morse function corresponding to a single cell ( Figure 4 (d)).

[0120] 2. Multi-branch simulation

[0121] Then, the single-cell pseudo-track simulation software PROSSTT was used to generate S2, which is the simulation data containing three branch structures in PROSSTT, mainly including the expression of 500 genes in 150 cells. Figure 5 As shown in (a), the S2 data mainly contains 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 types of cells (denoted as C1, C2, C3, C4, C5, C6, C7, C8 and C9), and accurately identified the cell lineage containing 3 branches (such as Figure 5 (b)). Figure 5 (c) and Figure 5In (d), the spatial distribution of the reproductive embedding positions of 150 cells can be compared with the Morse function distribution of a single cell and 9 cells. It is found that the Morse function values ​​of the reproductive embedding boundary positions of almost all cells in the spatial distribution are close to 1 (yellow part), mainly including the cell lineage root (C1) and leaves (C6, C8, C9); Morse function values ​​close to 0 (red) are located in the middle or branching stage of the cell lineage (C3, C4, C5 and C7), among which the branching cell populations C3 and C5 contain single cells corresponding to the minimum Morse function value, and the leaf cell population C8 contains the single cell corresponding to the maximum Morse function value ( Figure 5 (d)). Further, the SCTRRG algorithm extracts the topological Reeb graph of the cellular connection network in S2 and marks the distribution of the maximum node (green) and minimum node (blue) of the Morse function corresponding to the initial and intermediate nodes of the Reeb graph ( Figure 5 (e)).

[0122] In the absence of prior biological information, the initial node of cell differentiation can be selected in VB. The cell with the highest equivalence class in the C9Morse function value can be selected as the initial node, and the cell point at the end of the branch is used as the initial point, combined with the distance of the single cell connection on the network path; determine the single cell pseudo moment, and render it in space to embed the dimension reduction ( Figure 5 (f)). Figure 5 (d) and Figure 5 (f) Comparison: along the trajectory of cell differentiation, from the initial point to the middle stage (in the spatial position of cell distribution, corresponding to cell C2), the average value of the Morse function of the corresponding cell group changes from 0.60 to 0.10 (C9→C5), and from the middle stage to the three end stages, the median of the cell Morse function changes from 0.10 to 0.81 (C5→C1). The Morse function changes in the range of 0.10-0.80 (C5→C6) and 0.1-0.84 (C5→C8), indicating that the average change of the Morse function in the cell group increases with the early stage of the differentiation process and decreases with the increase of the overall differentiation process.

[0123] 3. Bifurcation Identification of Mouse Neural Stem Cell Differentiation

[0124] In addition, the SCTRRG algorithm was used to identify and infer cell lineages and pseudotimes during the differentiation of mouse neural stem cells. The MNS dataset mainly contains a single-cell dataset of 43,309 genes expressed in 272 neural stem cells. After preprocessing of highly variable genes (Fano factor > 0.5, average value between [0.01 and 3.5]), 8846 genes were screened. Figure 6As shown in (a), the MNS dataset mainly contains 4 labeled cell populations (respectively represented as: 1, 2, 3, and 4). In the low-dimensional space of the diffusion map, the distribution of labeled 3 cells and labeled 4 cells partially overlaps in the low-dimensional space, which increases the difficulty of inferring the direction of cell differentiation. Based on the differential expression of 8 marker genes (Cd24a, Tubb3, Sox9, Id2, Sox10, Mbp, Egr1, and Fos) in mouse neural stem cells and prior information, the initial cell population should be qNSC (Sox9 high expression). Based on the different expressions of the 8 marker genes, the SCTRRG algorithm identified 4 cell populations, which were labeled as cell populations qNSC, NB, Oli, and aNSC, and accurately identified the cell lineages of the two branches. Figure 6 (b) In Figure 6 (c) and Figure 6 In (d), the spatial distribution of the germline embedding positions of 272 neural stem cells, the Morse function distribution of individual cells and cell groups can be compared. It can be found that cells with Morse function values ​​close to 1 (yellow part) are mostly distributed in the spatial position of the germline embedding on the left (Oli cell group). From the overall distribution, the root cell line (cell qNSC) and leaf cell line (cell group Oli and aNSC) located below the Morse value of the cell line branch phase (NB) cells, and the leaf cell group Oli contains. The single cell corresponding to the maximum Morse function value ( Figure 6 (d)). Furthermore, the SCTRRG algorithm extracted the topological Reeb graph of the mouse neural stem cell connection network. The initial node and intermediate node of the Reeb graph can be marked, corresponding to the maximum point (pink) and minimum point (brown) of the Morse function, respectively ( Figure 6 (e)). Compare Figure 6 (a) and Figure 6 (e) It can be found that in the third and fourth category labels, some low-dimensional features of neural stem cells overlap, making it impossible to distinguish the trend of cell differentiation. In the Reeb diagram extracted by the SCTRRG algorithm, the connection between cell groups and the overall trend of cell differentiation can be clearly seen. 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 in Figure 6 (f) shows. Figure 6 (d) and Figure 6(f) Along the cell differentiation trajectory, the mean values ​​of the Morse function of the corresponding cell groups of mouse neural stem cells from the initial moment of cell qNSC to the two phases (NB and aNSC cell groups) changed from 0.46 to 0.37 (qNSC→NB) and 0.46 to 0.62 (qNSC→aNSC), and the median values ​​of the cell Morse function changed from 0.37 to 0.99 (NB→Oil) from the NB stage to the end stage. This indicates that the mean values ​​of the corresponding Morse function of mouse neural stem cells move along the branches, and the mean values ​​of the corresponding Morse functions experience a process of first decreasing and then increasing.

[0125] In summary, the numerical results show that the SCTRRG algorithm can infer two branching trajectories in the differentiation process of mouse neural stem cells and determine the pseudo-moment of each cell. The trend of the average degree of the Morse function of the cell population is related to the cell differentiation process.

[0126] 4. Identification of multiple branches in mouse embryos before implantation

[0127] Finally, the SCTRRG algorithm can be used to identify and infer multi-branch trajectories and pseudo-time during cell differentiation before mouse embryo implantation. The MLP dataset mainly contains single-cell datasets of 48 genes expressed in 442 cells. Since the number of cells is greater than the number of genes, the genes here are not filtered. The MLP dataset mainly contains cell populations in seven stages (labeled 1C, 2C, 4C, 8C, 16C, 32C, and 64C).

[0128] exist Figure 7 In the low-dimensional space of the diffusion map shown in (a), it can be seen that cells have multiple differentiation directions. The cells in the intermediate 8C, 16C, and 32C stages are relatively chaotic, and the connections between cells cannot be distinguished. Figure 7 (b) and Figure 7 In (c), according to the maximum value of the Morse function corresponding to the initial cell differentiation, the initial cells of the cell group can be in the 1c+2c stage. The SCTRRG algorithm can obtain 9 types of cells, namely 1c+2c, 4c, 8c, 16c, 32c-TE, 64c-TE, ICM, EPI and PE stages, accurately identifying the three branch lineages in the cell differentiation process. The initial cell group in the mouse embryonic development process is the initial stage of cell differentiation. Figure 7 (c) and Figure 7 In (d), the Morse function distribution of 442 units can be compared, and the Morse function value is close to 1 (yellow part). Most of the units fall three units lower 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 value close to 0.1 are much lower than those in the cell line root (64c-TE, EPI, PE) and leaf cells (1c+2c).

[0129] In addition, the SCTRRG algorithm extracted the topological Reeb graph of the connection network during the preimplantation process of mouse embryos. The initial nodes and intermediate nodes of the Reeb graph can be marked as the maximum node (red) and minimum node (green) of the Morse function corresponding to them ( Figure 7 (e)). Figure 7 (a) Compared with the Reeb diagram, the connection between 16C and 32C-TE and ICM of the cell population can be clearly seen.

[0130] Finally, combined with the path distances of the cellular communication network during the mouse preimplantation process, the pseudo-moments of individual cells can be determined and Figure 7 (f) The spectrum is presented in the embedding space. Figure 7 (d) and Figure 7 (f) Along the trajectory of cell differentiation before mouse implantation, the mean value of the corresponding cell Morse function changes from 0.80 to 0.11 from the initial moment of 1c+2c to the middle stage (16c). From the middle stage to the three ends (PE, EPI, 64c-TE), the median of the cell 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). This indicates that the mouse embryo is in the process of cell differentiation before implantation. The average Morse function of the cell group 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 process of mouse embryos and determine the pseudo-moment of each cell. The changing trend of the Morse function is related to the cell evolution process of different branching trajectories.

[0132] 5. Pseudo-trajectory algorithm comparison

[0133] In order to test the accuracy of the algorithm, two classic single-cell pseudo-trajectory recognition algorithms, DPT and Monocle2, were selected to compare with the SCTRRG algorithm. At the same time, the Kendall rank correlation coefficient between the cell prior order (simulated data) and the cell differentiation stage (real data) and the pseudo-time was calculated as an evaluation index for the accuracy of pseudo-trajectory judgment. Figure 8 As shown in the figure, at the pseudo-time of the three data sets, the inference results of the SCTRRG algorithm are better than those of DPT and Monocle2. In the simulated data, the results of the SCTRRGp algorithm and the SCTRRGcp algorithm are not much different, while in the actual data, the results of the SCTRRGp algorithm are better than those of other algorithms.

[0134] Therefore, when using the SCTRRG algorithm to identify single-cell pseudo-sites, it is recommended to infer the pseudo-moments of cells by considering the path distance in the cell network. For single-cell data without information on the cell differentiation stage, the Kendall rank correlation coefficient cannot be calculated, such as the mouse neural stem cell dataset. In summary, the method proposed in this paper can well infer the two-branch and multi-branch lineage processes of single cells in different simulated data and real data, and the inference accuracy on these three datasets can be better than that of two other very classic algorithms, which can further prove the rationality and effectiveness of this method.

[0135] The above shows and describes the basic principles, main features and advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention to be protected. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A single-cell pseudo-trajectory topology recognition method based on Reeb graph, characterized in that: include: S1: Preprocess the single-cell transcriptional dataset to obtain the single-cell gene expression matrix; S2: Based on the single-cell gene expression matrix of S1, an optimization model is constructed to further construct a single-cell connectivity network for cell population identification; S3: Based on S2, the single-cell Reeb graph is extracted by combining the cell’s Morse function distribution and equivalence class definition; S4: Based on S3, combined with the Morse function maximum, the key nodes in the Reeb graph are identified, the intermediate stages and initial and final states of cell differentiation are inferred, the distance influence of cells in the network and equivalence class is measured, and cell pseudo-trajectories are identified.

2. The single cell pseudo trajectory topology recognition method based on Reeb graph according to claim 1, characterized in that: S2 includes: S2.1: Build optimization model: stA T y=y; y＝(1,1…,1) T Among them, A corresponds to the adjacency matrix of the single-cell connectivity network, X is the single-cell gene expression matrix, ||·|| * is the nuclear norm, ||·|| 2,1 is the sum of the Euclidean norms of all columns, and the non-negative parameter λ is used to balance the nuclear norm and the Euclidean norm; the linear constraint A T y=y maintains the invariance of data translation transformation, y is a vector of length n, Indicates that the i-th cell and the j-th cell are not neighbor structures; The alternating direction multiplier method ADMM was used to solve the optimization model of S2.1, and the K-nearest neighbor algorithm KNN was used to determine the cell's neighbor network G. The selection of K was based on the BIC principle. The non-zero value in Aopt was agreed to be 1, and a 0-1 single-cell similarity network Aopt was obtained. S2.2: Based on the similarity network Aopt of S2.1, a single-cell connectivity network is constructed to identify cell populations: The 0-1 single-cell similarity network Aopt is converted into a single-cell symmetric network Z,λ by the maximum method l is the lth eigenvalue of L in order of size, p is the subscript corresponding to the maximum absolute value difference between two adjacent eigenvalues, which is used to identify the number of cell groups. Combined with the spectral clustering algorithm, p cell groups are identified, namely C1, C2, ... Cp; when Z is a non-connected network, dis represents the Euclidean distance function, Respectively represent the i-th cell in the u-type cell group and the j-th cell in the v-type cell group, and calculate the two cells corresponding to the minimum Euclidean distance between cells in the cell group where i * ,j * They are In the corresponding cell order of the entire network, Z is gradually updated to a single-cell connected network; for the case where there is an edge between two cell groups, Z ij =1, then there is an edge between the two cell groups, and further a cell group connectivity network CC between the cell groups is obtained.

3. The single cell pseudo trajectory topology recognition method based on Reeb graph according to claim 2, characterized in that: S3 includes: S3.1: The single-cell Morse function distribution is: Where u(c i , c j ) is the geometric distance between the i-th cell and the j-th cell in the non-connected network Z, which is approximated by the Dijkstra algorithm. i Normalized to To maintain the invariance of the f value with respect to the grid scale, n is the number of cells; S3.2: Extracting single-cell Reeb maps f : R f =G(V r ,E r ) in, and They represent the average geodesic distances of the uth and vth cells in the i-th cell population, d represents the error, The equivalence relationship within the same cell population is defined as the difference of the Morse function within a small error; Based on equivalence class According to the definition of single cell connectivity network, Reeb graph R is extracted f , R f V r The nodes are determined by the representative nodes in the equivalence class, R f Edge Set E r Whether the edge exists is determined by the connectivity of cells in the two equivalence classes [p] and [q] in the connected network Z. If there is a connection between s cell and t cell, A s，t ≠0, E r There are edges between the equivalent nodes [p] and [q].

4. The single cell pseudo trajectory topology recognition method based on Reeb graph according to claim 3, characterized in that: S4 includes: S4.1: Sort the nodes in the Reeb graph of S3: Among them, V B is a cell cluster located at the boundary of the cell cluster connected network CC, 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] , corresponding to the initial state or terminal state of the cell differentiation process; the equivalent cell class v with the smallest Morse function [m] Located throughout R f The middle position of the cell corresponds to the middle position of cell differentiation. S4.2: Based on the initial equivalent node v [0] and S3 entire Reeb graph R f , calculate the shortest path length from each node in the Reeb graph to the initial equivalent node, denoted as Len. According to the size of Len, R f n r Sort equivalent nodes: T [i] =Len(v [i] ,v [0] ); S4.3: Measuring cell distances: Where d(s,j) represents the shortest path length between the sth and jth cells in the single-cell connection network Z, and D ij It measures the average distance from the j cell to the [i] equivalent class cell set. It measures the average distance from cell j in cell population c to the equivalent cell set [i]. exp is the exponential function, var is the variance function, and the normalized D ij , The average distance between cells The effects of the two distances on cells are characterized from the perspective of the overall cell network and the local network of cells in the cell population; S4.4: Combining the influence of inter-cell distance and the order of cell nodes in the Reeb graph, the pseudo-times of all cells are obtained, that is, the single-cell pseudo-trajectory inference is realized: If the influence of the whole cell network on v cell is considered, then i=1,2,...,n r ,v=1,2,...,n; if we consider the local network of cells in the cell group, then 5. A single-cell pseudo-trajectory topology recognition device based on Reeb graph, characterized in that: The method is constructed based on the single-cell pseudo-trajectory topology recognition method based on the Reeb graph as described in any one of claims 1 to 4.

6. A single-cell pseudo-trajectory topology recognition device based on Reeb graph, characterized in that: The control device includes a processor and a memory, wherein a computer program is stored in the memory, and the computer program is loaded and executed by the processor to implement the single-cell pseudo-trajectory topology recognition method based on the Reeb graph as described in any one of claims 1 to 4.

7. 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 a Reeb graph as described in any one of claims 1 to 4.

Images