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

Abstract

The invention 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. Comprising the following steps: S1, constructing a single cell transcription data matrix of behavior gene names and cell names; s2, preprocessing the single-cell transcription data matrix, constructing a cell expression matrix, and performing characteristic gene screening on the cell expression matrix to obtain a single-cell gene expression matrix; s3, constructing a gene co-expression network; s4, constructing a gene rank differential expression network in combination with the single-cell gene expression matrix and the gene co-expression network; and S5, carrying out integrated clustering on the gene rank differential expression network, fusing by using a similarity network to obtain a final clustering result, determining an initial cluster, distributing pseudo time for cells in each cluster, and deducing a cell pseudo trajectory. According to the method, the differentiation track of the cell can be identified, and compared with an original gene expression matrix, a more stable identification result can be obtained by utilizing a relative expression sequence.

Description

Technical Field

[0001] The present invention belongs to the technical field of single-cell data analysis, and in particular relates to a single-cell pseudo-trajectory recognition method based on a gene co-expression network. 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 algorithms 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 category is the pseudo-trajectory inference algorithm based on tree structure. Through dimensionality reduction clustering, the minimum spanning tree (MST) is constructed as the main path, and then various strategies are adopted to obtain the pseudo-time order and lineage structure. For example, TSCAN constructs MST based on the centroid of the cell cluster and then infers the pseudo-time order of the cell. This method reduces the complexity of the tree space. Monocle was originally used to reconstruct the differentiation trajectory of human skeletal muscle myoblasts. One limitation of Monocle is that it cannot predict bifurcations in the lineage, so it can only analyze linear differentiation systems. This problem has been solved in the improved Monocle2. Monocle2 is one of the most mainstream pseudo-trajectory algorithms. It iteratively determines the stable position of cell projections in low-dimensional space based on the DDRTree algorithm, and generates MST to represent the inferred lineage structure. Slingshot constructs a clustering-based MST to represent the global lineage, and uses the synchronized master curve algorithm to fit smooth branch curves to these lineages. TiC2D first uses consensus clustering to identify cell clusters, and then constructs MST based on these clusters. The Waterfall algorithm is also an algorithm based on the minimum spanning tree. It reduces the impact of noise by calculating the minimum spanning tree distance of the cluster composition.

[0005] The second category is the pseudo-trajectory inference algorithms based on graph structures. These algorithms represent data as an undirected or directed neighborhood graph of a single cell, where nodes correspond to cells and edges represent the relationship between cells. Then a suitable path is selected in the graph to track the biological process from the ancestral cell to different states. Wanderlust constructs an intercellular graph based on the KNN graph and selects the shortest path in the graph to infer the linear trajectory. Wishbone improves the Wanderlust method so that Wishbone can identify bifurcated trajectories. SLICER first constructs a k-nearest neighbor graph between cells and then finds the shortest path between the initial cell and the final cell as a pseudo-time trajectory. 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 through the Louvain algorithm and connects each node through a weighted edge representing a statistical measure of connectivity between subpopulations to obtain an undirected graph, thereby inferring the cell development trajectory. SoptSC uses the intercellular similarity matrix to infer pseudo-time from the intercellular graph and uses the minimum spanning tree of the intercluster graph to predict the cell lineage relationship between clusters.

[0006] The third type of algorithm uses the dynamic changes of sheared mRNA and unsheared mRNA to infer the differentiation trajectory of cells. Recently, researchers proposed a likelihood-based kinetic model, scVelo, which can generalize RNA velocity to non-steady states and systems with different dynamics, and define the shared potential time of genes representing the internal clock of the cell. The TIRV algorithm based on RNA velocity converts RNA velocity into cosine similarity and finds branch points through branch entropy. CellPath is another method that integrates single-cell gene expression dynamics and RNA velocity information for trajectory inference. However, when performing trajectory inference, it is usually analyzed based on the gene expression matrix. The gene expression amount varies between different samples, while the relative value between genes is very stable. Therefore, using the correlation between genes for analysis not only greatly enhances stability, but also explores the biological significance of cell differentiation and development based on the interaction between genes. Summary of the invention

[0007] In view of the above problems, the purpose of the present invention is to propose an algorithm based on gene co-expression network, called scRDEN, for finding associated genes and inferring pseudo-trajectories. This method uses single-cell data as input to construct the network relationship between genes. By constructing a gene co-expression network, the gene expression matrix is ​​converted into a gene rank differential expression matrix to infer pseudo-trajectories. Through three data sets, it is proved that scRDEN is superior to other algorithms in trajectory inference and has a good effect on clustering.

[0008] In order to achieve the above technical effects and the above objectives, the present invention provides a single-cell pseudo-trajectory identification method based on a gene co-expression network, comprising:

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

[0010] S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, screen the characteristic genes of the cell expression matrix, and 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: Combine the single-cell gene expression matrix and gene co-expression network to construct a gene rank differential expression network;

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

[0014] Preferably, S3 includes:

[0015] S3.1: Calculate single-cell gene expression matrix A based on Pearson correlation coefficient m,n Gene correlation, get the correlation matrix R m,m :

[0016]

[0017] y k and z k Represents the single cell gene expression matrix A m,n The expression levels of the genes in the yth and zth rows in the kth cell, Represents the single cell gene expression matrix A m,n The average expression of genes in rows y and z in all cells in the matrix, the correlation matrix R m,m Represents the Pearson correlation coefficient between m genes, with a value range of [-1,1];

[0018] S3.2: Correlation matrix R m,m Further correlation tests were performed to retain the relevant genes with significant levels, and the absolute value of the correlation matrix was used as a measure of the correlation strength attribute of the correlation matrix:

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

[0020] Among them, B0 is the correlation matrix after absolute value;

[0021] S3.3: Use the WGCNA method to calculate the threshold b0, retain the gene pairs with correlation greater than the threshold, and then construct the following gene co-expression network B:

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

[0023]

[0024] Among them, G(V B ,E B ) is represented as a gene co-expression network B, V B ,E B are the node set and edge set of gene co-expression network B, E B Contains the element e ij , if the correlation b between the i-th gene and the j-th gene ij If it exceeds the threshold b0, then e ij =1, indicating that there is a direct correlation between two genes in the gene co-expression network B, and the correlation b between the i-th gene and the j-th gene ij Remember B i,j , indicating that the ith gene and the jth gene have been recorded in the node set V of the gene co-expression network B B ; If the correlation b between the i-th gene and the j-th gene ij If it does not exceed the threshold b0, then e ij =0, indicating that there is no direct correlation between the two genes and they are not added to the gene co-expression network B.

[0025] Preferably, S4 includes:

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

[0027]

[0028] (C0) ik A represents the expression value of the ith gene in the mth cell in ascending order Rank()↑. ij Represents the 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 a positive integer;

[0029] S4.2: Gene co-expression network B and rank matrix C0, construct gene rank differential expression network C:

[0030] C=G(Vc ,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 ) is represented as the gene rank differential expression network C, V c ,E c They represent the node set and edge set of the gene rank differential expression network C, respectively. m represents the number of genes, n represents the number of cells, and g represents the number of cells. k ×g l 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 is related to gene l, and the rank of gene k in cell t is Subtract the rank of gene l Get the rank difference value c gt ;

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

[0036]

[0037] c i,j represents 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 jth 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: Use TSNE, PCA, and UMAP to reduce the dimensionality of the adjacency matrix c, and calculate the distance matrix of the adjacency matrix c in the reduced dimensionality space corresponding to each dimensionality reduction method:

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

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

[0042] S5.2: Let the cell similarity network be represented as a graph G(V d ,E d ), vertex V d The corresponding cells {x1,x2,...,x n}, edge E d Represents the similarity of cells, and the edge weight is represented by an n*n similarity matrix W, where W i,j Represents cell x i and x j The similarity between them is calculated by i The average distance between the cells and their neighbors ρ i_ :

[0043]

[0044] N j Represents the cell x i The jth cell with the closest Euclidean distance, k is a hyperparameter set artificially. In the same way, the cell x is calculated j The average distance between the cells and their neighbors ρ j_ ;

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

[0046]

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

[0048]

[0049] Among them, μ is a hyperparameter set empirically;

[0050] S5.5: After the similarity matrix W is normalized, the normalized matrix P is obtained:

[0051]

[0052] S5.6: Constructing Cell x i Its neighboring nodes Ni The local affinity matrix of is denoted as the sparse kernel S:

[0053]

[0054] The normalized matrix P contains all the information about the similarity between each cell and all other cells, while S only contains the information about each cell and 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 the sparse kernel S (v) , the normalized matrix and sparse kernel are iteratively calculated as follows:

[0056] After t iterations:

[0057]

[0058] P t+1 (v) is the state matrix of the vth normalized matrix after t iterations;

[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 combines the features of TSNE, PCA, and UMAP after dimensionality reduction:

[0061]

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

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

[0064] About Eigengap i Sort by size, select the top two eigengap values, and select the i value of the smallest eigengap as the cluster number p for spectral clustering, and get the clustering results CC1, CC2, ..., CC p ; Based on prior knowledge, the starting cluster d1 is established, and the center points d1, d2, ..., d of each cluster in the TSNE dimensionality reduction space are used. p:

[0065]

[0066] Among them, G p represents the number of cells in the pth cluster. When calculating the position of cells in the dimensionality reduction space, the other two dimensionality reduction spaces are not involved in the calculation. Only the TSNE dimensionality reduction space is used to represent the position of cells in the low-dimensional space. i Indicates the position of each cell in the pth cluster in the TSNE corresponding dimensionality reduction space;

[0067] S5.10: Use the Euclidean distance dist to calculate the distance between the starting cluster d1 and other clusters and sort them:

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

[0069]

[0070] Among them, D u Indicates the distance between the center point of the u-th cluster and the starting cluster d1, Indicates that the distance between each cluster center and the starting cluster d1 is ranked in ascending order. The pseudo-time of the f cell in the v-type cell population is the position in ascending order among all center distances, Indicates CC v The pseudo-time values ​​of f cells in the cluster are used to assign the distance between the cluster to which the cell belongs and the starting cluster as the pseudo-time to the cells in the same cluster, and the cell ranking is obtained, thereby realizing the single-cell pseudo-trajectory inference.

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

[0072] 1. The method proposed in the present invention is based on the structure of the gene co-expression network to better extract the characteristics of cells. The unstable gene expression is converted into stable relative gene expression, which provides a basis for analyzing the specificity of cells.

[0073] 2. Use multiple dimensionality reduction methods to reduce the gene rank differential expression matrix, and use similarity network fusion to infer more accurate clustering results and more robust cell trajectories.

[0074] 3. The gene rank differential expression matrix can not only be used for clustering and cell trajectory identification, but also for biological significance analysis based on interrelated genes.

[0075] 4. Compared with other classic methods, this method has better recognition effect and is closer to the real trajectory when inferring differentiation trajectory on the data set. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0077] In the attached picture:

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

[0079] Figure 2 This is a comparison of the POS, Kendall correlation, and Bubble-sort index indicators of scRDEN, Monocle2, and TSCAN pseudo-trajectory inference algorithms in mouse lung epithelial cells.

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

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

[0082] Figure 5 Network topology, clustering coefficient and diversity of scRDEN at different pseudo-times in mouse lung epithelial cells.

[0083] Figure 6 The figure is a flow chart of the method of the present invention. DETAILED DESCRIPTION

[0084] The following is combined with Figure 1-6 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.

[0085] A single-cell pseudo-trajectory identification method based on a gene co-expression network, comprising:

[0086] S1: Construct a single-cell transcriptional data matrix with behavioral gene names as columns and cell names as columns.

[0087] S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, screen the characteristic genes of the cell expression matrix, and obtain the single-cell gene expression matrix; including:

[0088] S2.1: Gene screening:

[0089] In order to better study the expression of related genes in biological development and other processes, mitochondrial genes, spike-in genes and ribosomal genes are screened out according to gene names, usually ranging from dozens to dozens, while retaining genes expressed in more than 10% of cells;

[0090] S2.2: Outlier smoothing: After the initial screening of genes, the abnormally expressed genes are processed. By calculating the average expression μ and variance expression σ of each gene in all cells, check whether the expression range of each gene is within the range of [μ-4σ,μ+4σ] according to the 3σ principle. If it exceeds the boundary, the expression of the gene is smoothed to the gene expression corresponding to the boundary value to obtain the gene expression matrix A0. By smoothing the outliers, the effective gene expression in the single-cell data can be retained.

[0091] S2.3: log transformation: Before selecting characteristic genes, in order to eliminate the dimensional differences in gene expression, the gene expression matrix A0 obtained by initial screening and abnormal processing is log2 transformed:

[0092] A1=log2(A0+1)

[0093] A1 represents the normalized expression data. In order to keep the expression value consistent with 0, the expression values ​​of all genes were added with 1 and then log-processed.

[0094] S2.4: Feature gene screening: In order to effectively avoid the interference of external factors and capture the stable performance of the cell population, the intrinsic entropy model was selected to select feature genes. This model is based on information theory and believes that the intrinsic entropy and extrinsic entropy of a gene actually represent the inherent fluctuation and "noisy" fluctuation of the gene, respectively. Intrinsic entropy can capture the degree of regulatory fluctuation of each gene in the cell process, and can identify information genes from single-cell transcription data in order to accurately cluster and classify cells. We selected high intrinsic entropy genes with an average expression greater than I after log processing (I is a threshold selected artificially) to obtain a high-quality expression matrix A. m,n (m is the number of genes, n is the number of cells) This completes the processing work for the next step of network construction.

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

[0096] S3.1: Calculate single-cell gene expression matrix A based on Pearson correlation coefficient m,n Gene correlation, get the correlation matrix R m,m :

[0097]

[0098] y k and z kRepresents the single cell gene expression matrix A m,n The expression levels of the genes in the yth and zth rows in the kth cell, Represents the single cell gene expression matrix A m,n The average expression of genes in rows y and z in all cells in the matrix, the correlation matrix R m,m Represents the Pearson correlation coefficient between m genes, with a value range of [-1,1];

[0099] S3.2: Correlation matrix R m,m Further correlation tests were performed to retain the relevant genes with significant levels, and the absolute value of the correlation matrix was used as a measure of the correlation strength attribute of the correlation matrix:

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

[0101] Among them, B0 is the correlation matrix after absolute value;

[0102] S3.3: Use the WGCNA method to calculate the threshold b0, retain the gene pairs with correlation greater than the threshold, and then construct the following gene co-expression network B:

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

[0104]

[0105] Among them, G(V B ,E B ) is represented as a gene co-expression network B, V B ,E B are the node set and edge set of gene co-expression network B, E B Contains the element e ij , if the correlation b between the i-th gene and the j-th gene ij If it exceeds the threshold b0, then e ij =1, indicating that there is a direct correlation between two genes in the gene co-expression network B, and the correlation b between the i-th gene and the j-th gene ij Remember B i,j , indicating that the ith gene and the jth gene have been recorded in the node set V of the gene co-expression network B B ; If the correlation b between the i-th gene and the j-th gene ij If it does not exceed the threshold b0, then e ij =0, indicating that there is no direct correlation between the two genes and they are not added to the gene co-expression network B.

[0106] S4: Combine the single-cell gene expression matrix and gene co-expression network to construct a gene rank differential expression network; including:

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

[0108]

[0109] (C0) ik A represents the expression value of the ith gene in the mth cell in ascending order Rank()↑. ij Represents the 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 a positive integer;

[0110] S4.2: Gene co-expression network B and rank matrix C0, construct gene rank differential expression network C:

[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 ) is represented as the gene rank differential expression network C, V c ,E c They represent the node set and edge set of the gene rank differential expression network C, respectively. m represents the number of genes, n represents the number of cells, and g represents the number of cells. k ×g l 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 is related to gene l, and the rank of gene k in cell t is Subtract the rank of gene l Get the rank difference value c gt ;

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

[0117]

[0118] c i,j represents 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 jth 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 integrated clustering on the gene rank differential expression network, use similarity network fusion to obtain the final clustering result, determine the starting cluster, assign pseudo time to cells in each cluster, and infer cell pseudo trajectory. Including:

[0120] S5.1: Use TSNE, PCA, and UMAP to reduce the dimensionality of the adjacency matrix c, and calculate the distance matrix of the adjacency matrix c in the reduced dimensionality space corresponding to each dimensionality reduction method:

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

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

[0123] S5.2: Let the cell similarity network be represented as a graph G(V d ,E d ), vertex V d The corresponding cells {x1,x2,...,x n}, edge E d Represents the similarity of cells, and the edge weight is represented by an n*n similarity matrix W, where W i,j Represents cell x i and x j The similarity between them is calculated by i The average distance between the cells and their neighbors ρ i_ :

[0124]

[0125] N j Represents the cell x iThe jth cell with the closest Euclidean distance, k is a hyperparameter set artificially. In the same way, the cell x is calculated j The average distance between the cells and their neighbors ρ j_ ;

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

[0127]

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

[0129]

[0130] Among them, μ is a hyperparameter set empirically;

[0131] S5.5: After the similarity matrix W is normalized, the normalized matrix P is obtained:

[0132]

[0133] S5.6: Constructing Cell x i Its neighboring nodes N i The local affinity matrix of is denoted as the sparse kernel S:

[0134]

[0135] The normalized matrix P contains all the information about the similarity between each cell and all other cells, while S only contains the information about each cell and 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 the sparse kernel S (v) , the normalized matrix and sparse kernel are iteratively calculated as follows:

[0137] After t iterations:

[0138]

[0139] P t+1 (v) is the state matrix of the vth normalized matrix after t iterations;

[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 combines the features of TSNE, PCA, and UMAP after dimensionality reduction:

[0142]

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

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

[0145] About Eigengap i Sort by size, select the top two eigengap values, and select the i value of the smallest eigengap as the cluster number p for spectral clustering, and get the clustering results CC1, CC2, ..., CC p ; Based on prior knowledge, the starting cluster d1 is established, and the center points d1, d2, ..., d of each cluster in the TSNE dimensionality reduction space are used. p :

[0146]

[0147] Among them, G p represents the number of cells in the pth cluster. When calculating the position of cells in the dimensionality reduction space, the other two dimensionality reduction spaces are not involved in the calculation. Only the TSNE dimensionality reduction space is used to represent the position of cells in the low-dimensional space. i Indicates the position of each cell in the pth cluster in the TSNE corresponding dimensionality reduction space;

[0148] S5.10: Use the Euclidean distance dist to calculate the distance between the starting cluster d1 and other clusters and sort them:

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

[0150]

[0151] Among them, D u Indicates the distance between the center point of the u-th cluster and the starting cluster d1, Indicates that the distance between each cluster center and the starting cluster d1 is ranked in ascending order. The pseudo-time of the f cell in the v-type cell population is the position in ascending order among all center distances, Indicates CC v The pseudo-time values ​​of f cells in the cluster are used to assign the distance between the cluster to which the cell belongs and the starting cluster as the pseudo-time to the cells in the same cluster, and the cell ranking is obtained, thereby realizing the single-cell pseudo-trajectory inference.

[0152] Simulation experiment:

[0153] Mouse lung epithelial cell data were analyzed. Mouse lung epithelial cell data included 101 cells and 23,363 genes. After preprocessing the gene expression matrix of mouse lung epithelial cells, 101 cells and 3,768 gene expression matrices were obtained to create a gene co-expression network. Then, the 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 used three dimensionality reduction methods: PCA, TSNE, and UMAP, and the similarity matrix was constructed in the dimensionality reduction space for network fusion. Three cell clusters were determined based on the fused network, and the starting cluster was determined based on prior information. The distance between each cluster and the starting cluster was used as pseudo time to assign to the cells in the cluster.

[0154] In trajectory inference analysis, we compared scRDEN with Monocle2 and TSCAN algorithms in terms of POS, Kendall correlation, and Bubble-sort index. Figure 2 It can be found that the scRDEN method is superior to the other two trajectory inference algorithms in all three indicators, proving that the scRDEN method can infer more accurate cell trajectories.

[0155] In cluster analysis, we compared scRDEN with SC3 clustering, Kmeans clustering, and Spectralclustering in terms of ARI and NMI. Figure 3 , Figure 4 It can be found that the ARI and NMI indicators of the scRDEN method in cluster analysis are only lower than those of the SC3 clustering method, but better than Kmeans clustering and spectral clustering. We selected 12 genes to display the gene rank difference network under different pseudo-times. The darker the edge color between genes, the greater the rank difference between genes.

[0156] The network structure under different pseudo-times is quantified by clustering coefficient and diversity index. Figure 5 It can be found 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 late stages of development.

[0157] Figure 1 scRDEN algorithm process and clustering and trajectory inference analysis. Figure 2 This is a comparison of the POS, Kendall correlation, and Bubble-sort index indicators of scRDEN, Monocle2, and TSCAN trajectory inference algorithms in mouse lung epithelial cell data. Figure 3 ARI comparison of scRDEN and other clustering methods (SC3, Kmeans, Spectral) on mouse lung epithelial cell data. Figure 4 Comparison of NMI between scRDEN and other clustering methods (SC3, Kmeans, Spectral) on mouse lung epithelial cell data. Figure 5 Represents the network topology, clustering coefficient and diversity of mouse lung epithelial cells at different pseudo-times.

[0158] 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 identification method based on gene co-expression network, characterized in that: include: S1: Construct a single-cell transcriptional data matrix with behavioral gene names listed as cell names; S2: Preprocess the single-cell transcription data matrix, construct the cell expression matrix, screen the characteristic genes of the cell expression matrix, and obtain the single-cell gene expression matrix; S3: Construct a gene co-expression network based on the single-cell gene expression matrix; S4: Combine the single-cell gene expression matrix and gene co-expression network to construct a gene rank differential expression network; S5: Perform integrated clustering on the gene rank differential expression network, use similarity network fusion to obtain the final clustering result, determine the starting cluster, assign pseudo time to the cells in each cluster, and infer cell pseudo trajectory.

2. The single-cell pseudo-trajectory identification method based on gene co-expression network according to claim 1, characterized in that: S3 includes: S3.1: Calculate single-cell gene expression matrix A based on Pearson correlation coefficient m,n Gene correlation, get the correlation matrix R m,m : y k and z k Represents the single cell gene expression matrix A m,n The expression levels of the genes in the yth and zth rows in the kth cell, Represents the single cell gene expression matrix A m,n The average expression of genes in rows y and z in all cells in the matrix, the correlation matrix R m,m Represents the Pearson correlation coefficient between m genes, with a value range of [-1,1]; S3.2: Correlation matrix R m,m Further correlation tests were performed to retain the relevant genes with significant levels, and the absolute value of the correlation matrix was used as a measure of the correlation strength attribute of the correlation matrix: B 0＝ |(R m,m )| Among them, B0 is the correlation matrix after absolute value; S3.3: Use the WGCNA method to calculate the threshold b0, retain the gene pairs with correlation greater than the threshold, and then construct the following gene co-expression network B: B=G(V B ,E B ) Among them, G(V B ,E B ) is represented as a gene co-expression network B, V B ,E B are the node set and edge set of gene co-expression network B, E B Contains the element e ij , if the correlation b between the i-th gene and the j-th gene ij If it exceeds the threshold b0, then e ij =1, indicating that there is a direct correlation between two genes in the gene co-expression network B, and the correlation b between the i-th gene and the j-th gene ij Remember B i,j , indicating that the ith gene and the jth gene have been recorded in the node set V of the gene co-expression network B B ; If the correlation b between the i-th gene and the j-th gene ij If it does not exceed the threshold b0, then e ij =0, indicating that there is no direct correlation between the two genes and they are not added to the gene co-expression network B.

3. The single-cell pseudo-trajectory identification method based on gene co-expression network according to claim 2, characterized in that: S4 include: S4.1: Constructing a cell-based expression matrix A m,n The rank matrix C0 is: (C0) ik A represents the expression value of the ith gene in the mth cell in ascending order Rank()↑. ij Represents the 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 a positive integer; S4.2: Gene co-expression network B and rank matrix C0, construct gene rank differential expression network C: C=G(V c ,E c ) g k ,g l ∈V B ,e kl =1.g k ×g l ∈V c 1≤k,l≤n;1≤t≤m Among them, G(V C ,E C ) is represented as the gene rank differential expression network C, V c ,E c They represent the node set and edge set of the gene rank differential expression network C, m represents the number of genes, n represents the number of cells, and g k ×g l 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 is related to gene l, and the rank of gene k in cell t is Subtract the rank of gene l Get the rank difference value c gt ; S4.3: The rows of the adjacency matrix c of the gene rank differential expression network C represent the associated gene pairs, and the columns represent the cells. The adjacency matrix c is normalized: c i,j represents 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 jth column in the adjacency matrix c, max(c j ) represents the maximum rank difference of the j-th column in the adjacency matrix.

4. The single-cell pseudo-trajectory identification method based on gene co-expression network according to claim 3, characterized in that: S5 includes: S5.1: Use TSNE, PCA, and UMAP to reduce the dimensionality of the adjacency matrix c, and calculate the distance matrix of the adjacency matrix c in the reduced dimensionality space corresponding to each dimensionality reduction method: ρ (v) (x i ,x j )=dist(x i ,x j ),v=1,2,3. Among them, ρ (v) (x i ,x j ) represents cell x i and cell x j The Euclidean distance in three dimensionality reduction spaces, v represents three different dimensionality reduction methods; S5.2: Let the cell similarity network be represented as a graph G(V d ,E d ), vertex V d The corresponding cells {x1,x2,...,x n }, edge E d Represents the similarity of cells, and the edge weight is represented by an n*n similarity matrix W, where W i,j Represents cell x i and x j The similarity between them is calculated by i The average distance between the cells and their neighbors ρ i_ : N j Represents the cell x i The jth cell with the closest Euclidean distance, k is a hyperparameter set artificially. In the same way, the cell x is calculated j The average distance between the cells and their neighbors ρ j_ ; S5.3: Calculate cell x j The average distance between its neighboring cells ε i,j : S5.4: Based on S5.1-S5.3, convert the distance matrix into a similarity matrix W: Among them, μ is a hyperparameter set empirically; S5.5: After the similarity matrix W is normalized, the normalized matrix P is obtained: S5.6: Constructing Cell x i Its neighboring nodes N i The local affinity matrix of is denoted as the sparse kernel S: The normalized matrix P contains all the information about the similarity between each cell and all other cells, while S only contains the information about each cell and the k cells with the highest similarity; S5.7: Calculate the normalized matrix P of the similarity matrix after dimensionality reduction using TSNE, PCA, and UMAP (v) and the sparse kernel S (v) , the normalized matrix and sparse kernel are iteratively calculated as follows: After t iterations: P t+1 (v) is the state matrix of the vth normalized matrix after t iterations; 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 combines the features of TSNE, PCA, and UMAP after dimensionality reduction: P (i) 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 L = DE, where D represents the degree matrix; E represents the overall similarity matrix; obtain the eigenvalue λ and eigenvector x through Lx = λx, and define the eigenvalue difference eigengap i =λ i+1 -λ i , where λ i is the i-th eigenvalue; About Eigengap i Sort by size, select the top two eigengap values, and select the i value of the smallest eigengap as the cluster number p for spectral clustering, and get the clustering results CC1, CC2, ..., CC p ; Based on prior knowledge, the starting cluster d1 is established, and the center points d1, d2, ..., d of each cluster in the TSNE dimensionality reduction space are used. p : Among them, G p represents the number of cells in the pth cluster. When calculating the position of cells in the dimensionality reduction space, the other two dimensionality reduction spaces are not involved in the calculation. Only the TSNE dimensionality reduction space is used to represent the position of cells in the low-dimensional space. i Indicates the position of each cell in the pth cluster in the TSNE corresponding dimensionality reduction space; S5.10: Use the Euclidean distance dist to calculate the distance between the starting cluster d1 and other clusters and sort them: D u =dist(d u ,d1), Among them, D u Indicates the distance between the center point of the u-th cluster and the starting cluster d1, Indicates that the distance between each cluster center and the starting cluster d1 is ranked in ascending order. The pseudo-time of the f cell in the v-type cell population is the position in ascending order among all center distances, Indicates CC v The pseudo-time values ​​of f cells in the cluster are used to assign the distance between the cluster to which the cell belongs and the starting cluster as the pseudo-time to the cells in the same cluster, and the cell ranking is obtained, thereby realizing the single-cell pseudo-trajectory inference.

Images