A single cell pseudo trajectory identification method based on gene co-expression network

Abstract

The application belongs to the technical field of single cell data analysis, and particularly relates to a single cell pseudo trajectory identification method based on a gene co-expression network. The method comprises the following steps: S1, constructing a single cell transcription data matrix with behavior gene names and cell names; S2, preprocessing the single cell transcription data matrix, constructing a cell expression matrix, performing feature gene screening on the cell expression matrix, and obtaining a single cell gene expression matrix; S3, constructing a gene co-expression network; S4, combining the single cell gene expression matrix and the gene co-expression network, and constructing a gene rank difference expression network; S5, performing integrated clustering on the gene rank difference expression network, obtaining a final clustering result by using similarity network fusion, determining initial clustering, assigning pseudo time to cells in each cluster, and inferring cell pseudo trajectories. The application can identify the differentiation trajectory of cells, and compared with an original gene expression matrix, a more stable identification result can be obtained by using a relative expression order.

Description

Technical Field

[0001] This invention belongs to the field of single-cell data analysis technology, specifically relating to a single-cell pseudo-trajectory recognition method based on gene co-expression networks. 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 algorithms based on scRNA-seq data analysis of dynamic processes such as cell differentiation, cell fate development, and decision-making are mainly divided into three categories.

[0004] The first category is pseudo-trajectory inference algorithms based on tree structures. These algorithms construct minimum spanning trees (MSTs) as the main paths through dimensionality reduction clustering, and then employ various strategies to obtain pseudo-temporal order and phylogenetic structure. For example, TSCAN constructs an MST based on the centroids of cell clusters and then infers the pseudo-temporal order of cells, reducing the complexity of the tree space. Monocle, initially used to reconstruct the differentiation trajectories of human skeletal muscle myoblasts, has a limitation: it cannot predict bifurcation in the phylogenetic tree, thus only analyzing linear differentiation systems. This problem has been addressed in the improved Monocle2. Monocle2 is one of the most mainstream pseudo-trajectory algorithms, iteratively determining the stable positions of cell projections in low-dimensional space using the DDRTree algorithm, generating an MST to represent the inferred phylogenetic structure. Slingshot constructs a cluster-based MST to represent the global phylogenetic tree and uses the synchronous master curve algorithm to fit smooth branch curves to these phylogenetic trees. TiC2D first uses consensus clustering to identify cell clusters and then constructs an MST based on these clusters. The Waterfall algorithm, also based on minimum spanning trees, mitigates the impact of noise by calculating the distance between the minimum spanning trees formed by clusters.

[0005] The second category is pseudo-trace inference algorithms based on graph structures. These algorithms represent data as undirected or directed neighborhood graphs of individual cells, where nodes correspond to cells and edges represent relationships between cells. A suitable path is then selected in the graph to trace biological processes from ancestral cells to different states. Wanderlust constructs an inter-cell graph based on a KNN graph and selects the shortest path in the graph to infer linear trajectories. Wishbone improves upon Wanderlust's method, enabling it to identify bifurcation trajectories. SLICER first constructs an inter-cell k-nearest neighbor graph and then finds the shortest path between the initial and final cells as pseudo-temporal trajectories. GraphDDP starts with user-defined cluster assignments and then uses a force-based graph layout method to infer differentiation trajectories. The PAGA algorithm associates a node with each subpopulation using the Louvain algorithm and connects each node with weighted edges representing statistical measures of connectivity between subpopulations to obtain an undirected graph, thereby inferring cell developmental trajectories. SoptSC uses an inter-cell similarity matrix to infer pseudo-temporal trajectories from an inter-cell graph and uses the minimum spanning tree of the inter-cluster graph to predict cell lineage relationships between clusters.

[0006] The third type of algorithm infers cell differentiation trajectories by analyzing the dynamic changes in cleaved and uncleaved mRNA. Recently, researchers proposed a likelihood-based kinetic model, scVelo, which extends RNA velocity to non-steady-state systems with different dynamics and defines a gene-sharing potential time representing the cell's internal clock. The TIRV algorithm, based on RNA velocity, converts RNA velocity into cosine similarity and finds branch points using branch entropy. CellPath is another method that integrates single-cell gene expression dynamics and RNA velocity information for trajectory inference. However, trajectory inference typically relies on gene expression matrices, which vary across samples, while the relative values ​​between genes remain relatively stable. Therefore, analyzing gene inter-gene relationships not only significantly enhances stability but also allows for the exploration of biological significance in cell differentiation and development based on gene interactions. Summary of the Invention

[0007] To address the aforementioned problems, this invention aims to propose an algorithm based on gene co-expression networks, called scRDEN, for identifying associated genes and inferring pseudo-trajectories. This method uses single-cell data as input to construct a network of relationships between genes. By building a gene co-expression network, the gene expression matrix is ​​transformed into a gene rank difference expression matrix to infer pseudo-trajectories. Using three datasets, it is demonstrated that scRDEN outperforms other algorithms in trajectory inference and also achieves good results in clustering.

[0008] To achieve the aforementioned technical effects and objectives, this invention provides a single-cell pseudo-trajectory recognition method based on gene co-expression networks, comprising:

[0009] S1: Construct a single-cell transcription data matrix with behavioral gene names listed as cell names;

[0010] S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, and screen the cell expression matrix for characteristic genes to obtain the single-cell gene expression matrix;

[0011] S3: Construct a gene co-expression network based on the single-cell gene expression matrix;

[0012] S4: Construct a gene rank differential expression network by combining single-cell gene expression matrix and gene co-expression network;

[0013] S5: Perform integrated clustering on the gene rank differential expression network, use similarity network fusion to obtain the final clustering results, determine the initial clusters, assign pseudo-time to cells in each cluster, and infer cell pseudo-trajectories.

[0014] Preferably, S3 includes:

[0015] S3.1: Calculate the gene expression matrix A of single cells based on Pearson correlation coefficient. m,n Gene correlations were analyzed to obtain the correlation matrix R. m,m :

[0016]

[0017] y k and z k These represent single-cell gene expression matrices A and B, respectively. m,n The expression levels of genes in the y-th and z-th rows in the k-th cell. These represent single-cell gene expression matrices A and B, respectively. m,n The average expression level of genes in all cells in rows y and z of the matrix, and the correlation matrix R m,m represents the Pearson correlation coefficient among m genes, with a numerical range of [-1, 1];

[0018] S3.2: For the correlation matrix R m,m Further correlation testing was conducted, retaining genes with significant correlation levels, and using the absolute value of the correlation matrix as a measure of the correlation strength attribute.

[0019] B 0＝ |(R m,m )|

[0020] Where B0 is the absolute value-based correlation matrix;

[0021] S3.3: The threshold b0 is calculated using the WGCNA method, and gene pairs with correlation greater than the threshold are retained, thereby constructing the following gene co-expression network B:

[0022] B = G(V) B E B )

[0023]

[0024] Among them, G(V) B E B ) represents the gene co-expression network B, V B E B These are the node set and edge set of gene co-expression network B, respectively. B Contains element e ij If the correlation between the i-th gene and the j-th gene is b ij If the threshold b0 is exceeded, then e ij =1 indicates that there is a direct association between two genes in the gene co-expression network B, and the correlation b between the i-th gene and the j-th gene is calculated. ij Write it as B i,j This indicates that the i-th gene and the j-th gene have been recorded in the node set V of the gene co-expression network B. B If the correlation between the i-th gene and the j-th gene is b ij If it does not exceed the threshold b0, then e ij =0 indicates that the two genes are not directly related and will not be added to the gene co-expression network B.

[0025] Preferably, S4 includes:

[0026] S4.1: Constructing A based on cell expression matrix m,n The rank matrix C0:

[0027]

[0028] (C0) ik A represents the rank of the expression value of the i-th gene in the m-th cell, sorted in ascending order by Rank(). ij Gene expression matrix A m,n The expression level of the i-th gene in the j-th cell, (C0) ij ≤m means that the rank of the i-th gene in the j-th cell will not exceed the number of genes m, Z + Represents positive integers;

[0029] S4.2: Construct gene rank differential expression network C from gene co-expression network B and rank matrix C0.

[0030] C = G(V)c E c )

[0031]

[0032] g k ,g l ∈V B ,e kl =1,g k ×g l ∈V c

[0033] 1≤k,l≤n;1≤t≤m

[0034] Among them, G(V) C E C ) represents the gene rank differential expression network C,V c E c Let m and n represent the node set and edge set of the gene rank differential expression network C, respectively, where m represents the number of genes, n represents the number of cells, and g represents the number of nodes and edges. k ×g l This represents a set of gene pairs consisting of gene k and gene l, where all gene pairs constitute the node set in Vc; in the gene co-expression network B, gene k and gene l are related, and the rank of gene k in cell t is... Subtract the rank of gene l Obtain the rank difference value c gt ;

[0035] S4.3: The adjacency matrix c of the gene rank differential expression network C represents the rows of associated gene pairs and the columns of cells. The adjacency matrix c is normalized as follows:

[0036]

[0037] c i,j Let represent the rank difference value of the i-th row in the j-th cell of the adjacency matrix c. min(c j ) represents the minimum rank difference of the j-th column in the adjacency matrix c, max(c j ) represents the maximum rank difference of the j-th column in the adjacency matrix.

[0038] Preferably, S5 includes:

[0039] S5.1: Apply three dimensionality reduction methods—TSNE, PCA, and UMAP—to the adjacency matrix c, and calculate the distance matrix of c in the dimensionality reduction space corresponding to each method:

[0040] ρ (v) (x i ,x j ) = dist(x i ,xj ), v = 1, 2, 3.

[0041] Where, ρ (v) (x i ,x j ) represents cell x i and cell x j The Euclidean distance in the three dimensionality reduction spaces, where v represents the three different dimensionality reduction methods;

[0042] S5.2: Let the cell similarity network be represented as graph G(V) d E d Vertex V d Corresponding cells {x1,x2,...,x n}, edge E d The similarity of cells is represented by an n*n similarity matrix W, where Wi represents the edge weights. i,j Represents cell x i and x j The similarity between cells is calculated. i The average distance ρ between it and its neighboring cells i_ :

[0043]

[0044] N j Indicates the relationship with cell x i The j-th cell with the closest Euclidean distance is found, where k is a manually set hyperparameter. Similarly, cell x is calculated. j The average distance ρ between it and its neighboring cells j_ ;

[0045] S5.3: Calculate cell x j The average distance ε between it and its neighboring cells i,j :

[0046]

[0047] S5.4: Based on S5.1-S5.3, transform the distance matrix into a similarity matrix W:

[0048]

[0049] Where μ is a hyperparameter set empirically;

[0050] S5.5: The normalized matrix P is obtained after standardizing the similarity matrix W:

[0051]

[0052] S5.6: Constructing Cells x i Its neighboring node Ni The local affinity matrix, denoted as the sparse kernel S:

[0053]

[0054] The normalized matrix P contains all the information about the similarity of each cell to all other cells, while S only contains the information about each cell to the k cells with the highest similarity.

[0055] S5.7: Calculate the normalized matrix P of the similarity matrix after dimensionality reduction using TSNE, PCA, and UMAP. (v) and sparse kernel S (v) The normalized matrix and sparse kernel are then iteratively calculated as follows:

[0056] After t iterations:

[0057]

[0058] P t+1 (v) It is the state matrix after the v-th normalized matrix has been iterated t times;

[0059] When the maximum number of iterations is reached or the state matrix no longer changes, the matrix is ​​no longer updated.

[0060] S5.8: Calculate the overall similarity matrix E, which integrates the features after dimensionality reduction from TSNE, PCA, and UMAP.

[0061]

[0062] P (i) This represents the final state matrix of the three dimensionality reduction methods, TSNE, PCA, and UMAP, after the previous iteration.

[0063] S5.9: Calculate the Laplacian matrix L = DE, where D represents the degree matrix and E represents the overall similarity matrix; obtain the eigenvalues ​​λ and eigenvectors x through Lx = λx, and define the eigenvalue difference eigengap. i =λ i+1 -λ i , where λ i It is the i-th eigenvalue;

[0064] eigengap i Sort by size, select the top two eigengap values, and then select the smaller eigengap value as the cluster number p for spectral clustering, resulting in clustering results CC1, CC2, ..., CC p Based on prior knowledge, an initial cluster d1 is established, and the centroids d1, d2, ..., d1 of each cluster in the TSNE dimensionality reduction space are used. p:

[0065]

[0066] Among them, G p This represents the number of cells in the p-th cluster. When calculating the cell's position in the reduced-dimensional space, the other two reduced-dimensional spaces are not used in the calculation; only the TSNE reduced-dimensional space is used to represent the cell's position in the low-dimensional space. X i This represents the position of each cell in the p-th cluster in the TSNE-reduced space;

[0067] S5.10: Calculate the distances between the initial cluster d1 and other clusters using Euclidean distance dist, and then sort them:

[0068] D u =dist(d u ,d1),

[0069]

[0070] Among them, D u This represents the distance between the cluster center of cluster u and the initial cluster d1. This indicates that the distances between each cluster center and the starting cluster d1 are ranked in ascending order. The pseudo-time of cell f in cell population v is its position in ascending order of all center distances. CC v The pseudo-time values ​​of f cells in a cluster are used to assign pseudo-time values ​​to cells in the same cluster, with the distance between the cell's cluster and the starting cluster as the pseudo-time value. This results in the cell ranking and enables single-cell pseudo-trajectory inference.

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

[0072] 1. The method proposed in this invention is based on the structure of gene co-expression networks, which better extracts cellular characteristics. It transforms unstable gene expression levels into stable relative gene expression, providing a basis for analyzing cell specificity.

[0073] 2. Dimensionality reduction is performed on the gene rank differential expression matrix using multiple dimensionality reduction methods, and similarity network fusion is used to infer more accurate clustering results and more robust cell trajectories.

[0074] 3. Gene rank differential expression matrix can not only be used for clustering and cell trajectory recognition, but also for analyzing the biological significance of interrelated genes.

[0075] 4. Compared with other classic methods, this method achieves better recognition results and is closer to the real trajectory when inferring differentiation trajectories on the dataset. Attached Figure Description

[0076] 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.

[0077] In the attached diagram:

[0078] Figure 1 This describes the scRDEN algorithm flow.

[0079] Figure 2 Comparison of POS, Kendall correlation, and Bubble-sort index in mouse lung epithelial cells between scRDEN, Monocle2, and TSCAN pseudo-trajectory inference algorithms.

[0080] Figure 3 Comparison of scRDEN with other clustering methods (SC3, Kmeans, Spectral) in mouse lung epithelial cells for ARI.

[0081] Figure 4 Comparison of scRDEN with other clustering methods (SC3, Kmeans, Spectral) in mouse lung epithelial cells using NMI.

[0082] Figure 5 The network topology, clustering coefficients, and diversity of scRDEN at different pseudotimes in mouse lung epithelial cells are shown.

[0083] Figure 6 This is a flowchart of the method of the present invention. Detailed Implementation

[0084] The following is in conjunction with the appendix Figure 1-6 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.

[0085] A single-cell pseudo-trajectory recognition method based on gene co-expression networks includes:

[0086] S1: Construct a single-cell transcription data matrix with behavioral gene names listed as cell names.

[0087] S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, and perform feature gene screening on the cell expression matrix to obtain the single-cell gene expression matrix; including:

[0088] S2.1: Initial gene screening:

[0089] To better study the expression of related genes in biological development and other processes, mitochondrial genes, spike-in genes, and ribosome genes are screened out based on gene names, usually by several dozen, while retaining genes expressed in more than 10% of cells.

[0090] S2.2: Outlier Smoothing: This process handles outlier genes after initial screening. The mean expression level μ and variance σ of each gene across all cells are calculated. The expression range of each gene is checked against the 3σ rule to ensure it falls within the range [μ-4σ, μ+4σ]. If the expression exceeds the boundary, its level is smoothed to the boundary value, resulting in the gene expression matrix A0. Outlier smoothing preserves valid gene expression in single-cell data.

[0091] S2.3: Log Transformation: Before selecting feature genes, in order to eliminate the dimensional differences in gene expression, the gene expression matrix A0 obtained from the initial screening and anomaly treatment is transformed using log2:

[0092] A1 = log2(A0 + 1)

[0093] A1 represents the standardized expression data. To maintain consistency with expression values ​​of 0, the expression levels of all genes were incremented by one before logarithmic processing.

[0094] S2.4: Feature Gene Screening: To effectively avoid interference from external factors and capture the stability of cell populations, an intrinsic entropy model was chosen to select feature genes. This model, based on information theory, posits that a gene's intrinsic entropy and extrinsic entropy represent its inherent and "noisy" fluctuations, respectively. Intrinsic entropy can capture the degree of regulatory fluctuation of each gene during cellular processes and can identify informative genes from single-cell transcription data for accurate cell clustering and classification. We selected high intrinsic entropy genes with an average expression level greater than I (where I is a manually selected threshold) after log processing, obtaining a high-quality expression matrix A. m,n (m is the number of genes, n is the number of cells) This completes the processing for the next step of network construction.

[0095] S3: Construct a gene co-expression network based on the single-cell gene expression matrix. This includes:

[0096] S3.1: Calculate the gene expression matrix A of single cells based on Pearson correlation coefficient. m,n Gene correlations were analyzed to obtain the correlation matrix R. m,m :

[0097]

[0098] y k and z kThese represent single-cell gene expression matrices A and B, respectively. m,n The expression levels of genes in the y-th and z-th rows in the k-th cell. These represent single-cell gene expression matrices A and B, respectively. m,n The average expression level of genes in all cells in rows y and z of the matrix, and the correlation matrix R m,m represents the Pearson correlation coefficient among m genes, with a numerical range of [-1, 1];

[0099] S3.2: For the correlation matrix R m,m Further correlation testing was conducted, retaining genes with significant correlation levels, and using the absolute value of the correlation matrix as a measure of the correlation strength attribute.

[0100] B0 = |(R m,m )|

[0101] Where B0 is the absolute value-based correlation matrix;

[0102] S3.3: The threshold b0 is calculated using the WGCNA method, and gene pairs with correlation greater than the threshold are retained, thereby constructing the following gene co-expression network B:

[0103] B = G(V) B E B )

[0104]

[0105] Among them, G(V) B E B ) represents the gene co-expression network B, V B E B These are the node set and edge set of gene co-expression network B, respectively. B Contains element e ij If the correlation between the i-th gene and the j-th gene is b ij If the threshold b0 is exceeded, then e ij =1 indicates that there is a direct association between two genes in the gene co-expression network B, and the correlation b between the i-th gene and the j-th gene is calculated. ij Write it as B i,j This indicates that the i-th gene and the j-th gene have been recorded in the node set V of the gene co-expression network B. B If the correlation between the i-th gene and the j-th gene is b ij If it does not exceed the threshold b0, then e ij =0 indicates that the two genes are not directly related and will not be added to the gene co-expression network B.

[0106] S4: Construct a gene rank-differential expression network by combining single-cell gene expression matrices and gene co-expression networks; including:

[0107] S4.1: Constructing A based on cell expression matrix m,n The rank matrix C0:

[0108]

[0109] (C0) ik A represents the rank of the expression value of the i-th gene in the m-th cell, sorted in ascending order by Rank(). ij Gene expression matrix A m,n The expression level of the i-th gene in the j-th cell, (C0) ij ≤m means that the rank of the i-th gene in the j-th cell will not exceed the number of genes m, Z + Represents positive integers;

[0110] S4.2: Construct gene rank differential expression network C from gene co-expression network B and rank matrix C0.

[0111] C = G(V) c E c )

[0112]

[0113] g k ,g l ∈V B ,e kl =1,g k ×g l ∈V c

[0114] 1≤k,l≤n;1≤t≤m

[0115] Among them, G(V) C E C ) represents the gene rank differential expression network C,V c E c Let m and n represent the node set and edge set of the gene rank differential expression network C, respectively, where m represents the number of genes, n represents the number of cells, and g represents the number of nodes and edges. k ×g l This represents a set of gene pairs consisting of gene k and gene l, where all gene pairs constitute the node set in Vc; in the gene co-expression network B, gene k and gene l are related, and the rank of gene k in cell t is... Subtract the rank of gene l Obtain the rank difference value c gt ;

[0116] S4.3: The adjacency matrix c of the gene rank differential expression network C represents the rows of associated gene pairs and the columns of cells. The adjacency matrix c is normalized as follows:

[0117]

[0118] c i,j Let represent the rank difference value of the i-th row in the j-th cell of the adjacency matrix c. min(c j ) represents the minimum rank difference of the j-th column in the adjacency matrix c, max(c j ) represents the maximum rank difference of the j-th column in the adjacency matrix.

[0119] S5: Perform ensemble clustering on the gene rank differential expression network, use similarity network fusion to obtain the final clustering results, determine the initial clusters, assign pseudo-times to cells in each cluster, and infer cell pseudo-trajectories. This includes:

[0120] S5.1: Apply three dimensionality reduction methods—TSNE, PCA, and UMAP—to the adjacency matrix c, and calculate the distance matrix of c in the dimensionality reduction space corresponding to each method:

[0121] ρ (v) (x i ,x j ) = dist(x i ,x j ), v = 1, 2, 3.

[0122] Where, ρ (v) (x i ,x j ) represents cell x i and cell x j The Euclidean distance in the three dimensionality reduction spaces, where v represents the three different dimensionality reduction methods;

[0123] S5.2: Let the cell similarity network be represented as graph G(V) d E d Vertex V d Corresponding cells {x1,x2,...,x n}, edge E d The similarity of cells is represented by an n*n similarity matrix W, where Wi represents the edge weights. i,j Represents cell x i and x j The similarity between cells is calculated. i The average distance ρ between it and its neighboring cells i_ :

[0124]

[0125] N j Indicates the relationship with cell x iThe j-th cell with the closest Euclidean distance is found, where k is a manually set hyperparameter. Similarly, cell x is calculated. j The average distance ρ between it and its neighboring cells j_ ;

[0126] S5.3: Calculate cell x j The average distance ε between it and its neighboring cells i,j :

[0127]

[0128] S5.4: Based on S5.1-S5.3, transform the distance matrix into a similarity matrix W:

[0129]

[0130] Where μ is a hyperparameter set empirically;

[0131] S5.5: The normalized matrix P is obtained after standardizing the similarity matrix W:

[0132]

[0133] S5.6: Constructing Cells x i Its neighboring node N i The local affinity matrix, denoted as the sparse kernel S:

[0134]

[0135] The normalized matrix P contains all the information about the similarity of each cell to all other cells, while S only contains the information about each cell to the k cells with the highest similarity.

[0136] S5.7: Calculate the normalized matrix P of the similarity matrix after dimensionality reduction using TSNE, PCA, and UMAP. (v) and sparse kernel S (v) The normalized matrix and sparse kernel are then iteratively calculated as follows:

[0137] After t iterations:

[0138]

[0139] P t+1 (v) It is the state matrix after the v-th normalized matrix has been iterated t times;

[0140] When the maximum number of iterations is reached or the state matrix no longer changes, the matrix is ​​no longer updated.

[0141] S5.8: Calculate the overall similarity matrix E, which integrates the features after dimensionality reduction from TSNE, PCA, and UMAP.

[0142]

[0143] P (i) This represents the final state matrix of the three dimensionality reduction methods, TSNE, PCA, and UMAP, after the previous iteration.

[0144] S5.9: Calculate the Laplacian matrix L = DE, where D represents the degree matrix and E represents the overall similarity matrix; obtain the eigenvalues ​​λ and eigenvectors x through Lx = λx, and define the eigenvalue difference eigengap. i =λ i+1 -λ i , where λ i It is the i-th eigenvalue;

[0145] eigengap i Sort by size, select the top two eigengap values, and then select the smaller eigengap value as the cluster number p for spectral clustering, resulting in clustering results CC1, CC2, ..., CC p Based on prior knowledge, an initial cluster d1 is established, and the centroids d1, d2, ..., d1 of each cluster in the TSNE dimensionality reduction space are used. p :

[0146]

[0147] Among them, G p This represents the number of cells in the p-th cluster. When calculating the cell's position in the reduced-dimensional space, the other two reduced-dimensional spaces are not used in the calculation; only the TSNE reduced-dimensional space is used to represent the cell's position in the low-dimensional space. X i This represents the position of each cell in the p-th cluster in the TSNE-reduced space;

[0148] S5.10: Calculate the distances between the initial cluster d1 and other clusters using Euclidean distance dist, and then sort them:

[0149] D u =dist(d u ,d1),

[0150]

[0151] Among them, D u This represents the distance between the cluster center of cluster u and the initial cluster d1. This indicates that the distances between each cluster center and the starting cluster d1 are ranked in ascending order. The pseudo-time of cell f in cell population v is its position in ascending order of all center distances. CC v The pseudo-time values ​​of f cells in a cluster are used to assign pseudo-time values ​​to cells in the same cluster, with the distance between the cell's cluster and the starting cluster as the pseudo-time value. This results in the cell ranking and enables single-cell pseudo-trajectory inference.

[0152] Simulation experiment:

[0153] Data from mouse lung epithelial cells were analyzed. The data included 101 cells and 23,363 genes. After preprocessing the gene expression matrix of the mouse lung epithelial cells, expression matrices for 101 cells and 3,768 genes were obtained, which were used to construct a gene co-expression network. Then, a gene rank differential expression network was constructed using the processed gene expression matrix combined with the gene co-expression network. The gene rank differential expression network employed three dimensionality reduction methods: PCA, TSNE, and UMAP. A similarity matrix was constructed in the dimensionality-reduced space for network fusion. Based on the fused network, three cell clusters were identified. Initial clusters were determined based on prior information, and pseudo-time was allocated to cells within each cluster based on the distance between each cluster and the initial cluster.

[0154] In trajectory inference analysis, we compared scRDEN with Monocle2 and TSCAN algorithms on three metrics: POS, Kendall correlation, and Bubble-sort index. Figure 2 It can be observed that the scRDEN method outperforms the other two trajectory inference algorithms in all three metrics, proving that the scRDEN method can infer more accurate cell trajectories.

[0155] In cluster analysis, we compared scRDEN with SC3 clustering, K-means clustering, and Spectralclustering on the ARI and NMI metrics. Figure 3 , Figure 4 It can be observed that the scRDEN method has lower ARI and NMI in cluster analysis than the SC3 clustering method, but is superior to K-means clustering and Spectral clustering. We selected 12 genes to illustrate the gene rank difference network at different pseudo-times, where darker edge colors between genes represent greater rank differences.

[0156] The network structure under different pseudo-times was quantified using clustering coefficients and diversity indices. Figure 5 It can be observed that the clustering coefficient and diversity show a trend of first increasing and then decreasing, which is consistent with the fact that cell activity is intense in the early and middle stages of development and tends to be stable in the later stages of development.

[0157] Figure 1 This document outlines the scRDEN algorithm process and clustering and trajectory inference analysis. Figure 2 This study compares the POS, Kendall correlation, and Bubble-sort index of three trajectory inference algorithms—scRDEN, Monocle2, and TSCAN—in mouse lung epithelial cell data. Figure 3 A comparison of ARI in mouse lung epithelial cell data between scRDEN and other clustering methods (SC3, Kmeans, Spectral). Figure 4 Comparison of NMI between scRDEN and other clustering methods (SC3, Kmeans, Spectral) in mouse lung epithelial cell data. Figure 5 The network topology, clustering coefficients, and diversity of mouse lung epithelial cells at different pseudo-time points are represented.

[0158] 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 identifying pseudotrajectories in single cells based on gene co-expression networks, characterized in that: include: S1: Construct a single-cell transcription data matrix with behavioral gene names listed as cell names; S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, and screen the cell expression matrix for characteristic genes to obtain the single-cell gene expression matrix; S3: Construct a gene co-expression network based on the single-cell gene expression matrix; S3 includes: S3.1: Calculation of single-cell gene expression matrix based on Pearson correlation coefficient Gene correlations were analyzed to obtain a correlation matrix. : ； and These represent single-cell gene expression matrices. The expression levels of genes in the y-th and z-th rows in the k-th cell. These represent single-cell gene expression matrices. The average expression level of genes in all cells in rows y and z of the matrix, and the correlation matrix. The Pearson correlation coefficient represents the correlation coefficient among m genes, with a numerical range of 1. ; S3.2: For the correlation matrix Further correlation testing was conducted, retaining genes with significant correlation levels, and using the absolute value of the correlation matrix as a measure of the correlation strength attribute. ； in, B 0 It is the correlation matrix after absolute value transformation; S3.3: Calculate the threshold using the WGCNA method b 0. Gene pairs with correlation greater than a threshold are retained, and the following gene co-expression network B is constructed: ； in, Represented as gene co-expression network B, These are the node set and edge set of the gene co-expression network B, respectively. Contains elements If the correlation between the i-th gene and the j-th gene is... b ij Exceeding the threshold b 0, then This indicates that there is a direct association between two genes in gene co-expression network B, and the correlation between the i-th gene and the j-th gene is expressed as follows. Record This indicates that the i-th gene and the j-th gene have been recorded in the node set of the gene co-expression network B. If the correlation between the i-th gene and the j-th gene is... Not exceeding the threshold b 0, then This indicates that the two genes are not directly related and will not be included in gene co-expression network B; S4: Construct a gene rank differential expression network by combining single-cell gene expression matrix and gene co-expression network; S5: Perform integrated clustering on the gene rank differential expression network, use similarity network fusion to obtain the final clustering results, determine the initial clusters, assign pseudo-time to cells in each cluster, and infer cell pseudo-trajectories.

2. The single-cell pseudo-trajectory recognition method based on gene co-expression networks according to claim 1, characterized in that: S4 includes: S4.1: Constructing a cell expression matrix rank matrix C 0 : ； This indicates that the expression values ​​of the i-th gene in the m-th cell are sorted in ascending order. The ranking Gene expression matrix The expression level of the i-th gene in the j-th cell. This means that the rank of the i-th gene in the j-th cell will not exceed the number of genes m. Represents positive integers; S4.2: Gene Co-expression Network B and Rank Matrix C 0 Construct a gene rank differential expression network C: ； in, Represented as gene rank differential expression network C, Let m and n represent the node set and edge set of the gene rank differential expression network C, respectively, where m represents the number of genes and n represents the number of cells. This represents a set of gene pairs consisting of gene k and gene l, where all gene pairs constitute the node set in Vc; in the gene co-expression network B, gene k and gene l are related, and the rank of gene k in cell t is... Subtract the rank of gene l Obtain the rank difference value ; S4.3: The adjacency matrix c of the gene rank differential expression network C represents the rows of associated gene pairs and the columns of cells. The adjacency matrix c is normalized as follows: ； This represents the rank difference value of the i-th row in the j-th cell of the adjacency matrix c. This represents the minimum rank difference in the j-th column of the adjacency matrix c. This represents the maximum rank difference in the j-th column of the adjacency matrix.

3. The single-cell pseudo-trajectory recognition method based on gene co-expression networks according to claim 2, characterized in that: S5 includes: S5.1: Apply three dimensionality reduction methods—TSNE, PCA, and UMAP—to the adjacency matrix c, and calculate the distance matrix of c in the dimensionality reduction space corresponding to each method: ； in, Represents cells x i and cells x j The Euclidean distance in the three dimensionality reduction spaces, where v represents the three different dimensionality reduction methods; S5.2: Let the cell similarity network be represented as a graph. ,vertex Corresponding cells ,side The similarity of cells is represented by a single edge weight. Similarity matrix W It means that among them Represents cells x i and x j Calculate the similarity between cells x i Average distance to its neighboring cells : ； Indicates with cells x i The closest European-style distance j Each cell, k These are artificially set hyperparameters; the same method is used to calculate cells. x j The average distance between it and its neighboring cells ; S5.3: Calculating Cells x j The average distance between it and its neighboring cells : ； S5.4: Based on S5.1-S5.3, transform the distance matrix into a similarity matrix W: ； in, It is a hyperparameter set based on experience; S5.5: The normalized matrix P is obtained after standardizing the similarity matrix W: ； S5.6: Constructing Cells x i Its neighboring nodes N i The local affinity matrix, denoted as the sparse kernel S: ； The normalized matrix P contains all the information about the similarity of each cell to all other cells, while S only contains the information about the cell with the highest similarity. k Information about each cell; S5.7: Calculate the normalized matrix of the similarity matrix after dimensionality reduction using TSNE, PCA, and UMAP. and sparse kernel The normalized matrix and sparse kernel are then iteratively calculated as follows: After t iterations: ； It is the state matrix after the v-th normalized matrix has been iterated t times; When the maximum number of iterations is reached or the state matrix no longer changes, the matrix is ​​no longer updated. S5.8: Calculate the overall similarity matrix E, which integrates the features after dimensionality reduction from TSNE, PCA, and UMAP. ； This represents the final state matrix of the three dimensionality reduction methods, TSNE, PCA, and UMAP, after the previous iteration. S5.9: Calculate the Laplace matrix , D Degree matrix; E Represents the overall similarity matrix; through Obtain eigenvalues and eigenvectors Define the eigenvalue difference ,in It is the i-th eigenvalue; right Sort by size and select the top two. Values, select the smaller one. of i Value as cluster number p Perform spectral clustering to obtain the clustering results. Based on prior knowledge, the initial cluster d1 is established, and the centroids of each cluster in the TSNE dimensionality reduction space are used. : ； in, Indicates the first p When calculating the cell position in a cluster, the other two dimensionality reduction spaces are not used in the calculation; only the TSNE dimensionality reduction space is used to represent the cell position in the low-dimensional space. Indicates the first p The position of each cell in each cluster in the TSNE-reduced space; S5.10: Calculate the distances between the initial cluster d1 and other clusters using Euclidean distance dist, and then sort them: ； in, This represents the distance between the cluster center of cluster u and the initial cluster d1. This indicates that the distances between each cluster center and the starting cluster d1 are ranked in ascending order. The pseudo-time of cell f in cell population v is its position in ascending order of all center distances. express The pseudo-time values ​​of f cells in a cluster are used to assign pseudo-time values ​​to cells in the same cluster, with the distance between the cell's cluster and the starting cluster as the pseudo-time value. This results in the cell ranking and enables single-cell pseudo-trajectory inference.

Images