Cancer driver gene identification method based on single-cell pseudo-time sequence

Abstract

The invention discloses a cancer driver gene identification method based on a single cell pseudo time sequence. The method comprises the following steps: 1) extracting a characteristic gene; the method comprises the following steps of 1) establishing a pseudo time sequence, 2) determining the pseudo time sequence, and 3) identifying a significant mutation point and determining the cancer driver gene.According to the method, on the basis of constructed cancer cell development pseudo time sequence information, the change mode of the gene on the pseudo time sequence is revealed, and the cancer driver gene is effectively identified by judging whether the significant mutation point exists on the gene expression quantity pseudo time sequence or not.

Description

Technical Field

[0001] The present invention relates to the fields of statistics and bioinformatics, and in particular to a method for identifying cancer driver genes based on single-cell pseudo-time series. Background Art

[0002] Cancer is a heterogeneous disease characterized by dysregulation of normal cellular growth processes, caused by uncontrolled alterations in multiple genes. Among these altered genes, cancer driver genes are defined as those that directly contribute to the growth and proliferation of cancer cells. [1,2] Therefore, correctly identifying cancer driver genes is crucial for early detection, effective treatment and precision medicine of cancer. A large number of studies have shown that understanding cancer driver genes can help to conduct precision oncology in clinical practice, thereby achieving better treatment responses. [3-5] The recent development of high-throughput technology and international consortium-driven projects has brought researchers a large amount of data, which has brought more information to the discovery of cancer driver genes. The Cancer Genome Atlas (TCGA), which collects multi-omics data of various cancer types, is the most representative example of an international consortium-driven project. [6] In addition, researchers have made efforts to build comprehensive catalogs of cancer driver genes, such as the Network of Cancer Genes (NCG). [7] and COSMIC Cancer Gene Census (CGC) [8] .

[0003] In recent years, due to the advancement of single-cell RNA sequencing (scRNA-seq) technology, researchers can analyze the in vivo transcriptional status of thousands of single cells, which makes more transcriptome studies possible. [9-12] Compared with traditional bulk RNA-seq, which provides average gene expression of a cell population, [13,14] scRNA-seq technology quantifies cell-to-cell variation by analyzing the transcriptome of a single cell. Changes in gene expression in scRNA-seq data reflect the degree of variation of genes from one cell to another, and changes in gene expression play an important role in tissue development. [15-18] Therefore, considering the changes in single-cell gene expression data is important for understanding transcriptional regulation and is more helpful for identifying and analyzing intratumor heterogeneity. [19-21]. Therefore, by analyzing scRNA-seq data and the differences in expression levels between cells, it is more likely to discover precise cancer driver genes. By targeting these driver genes, researchers can develop more effective drugs and treatment strategies and provide personalized treatment and early diagnosis, which is particularly valuable for cancer-related research. In short, scRNA-seq data can study the differences in gene expression between cells and is more likely to discover cancer driver genes that are overlooked in bulk data.

[0004] Detecting cancer driver genes in laboratory experiments is challenging given the complexity of their regulation and the sheer number of genes (over 20,000).

[22] Therefore, it is necessary to develop computational methods that use a variety of genomic data to reveal cancer driver genes and the regulatory mechanisms behind cancer development. In traditional methods for identifying cancer driver genes, researchers mostly focus on mutant genes.

[23] OncodriveFM is a method to identify cancer driver genes based on their mutation frequency.

[24] ActiveDriver is a method that uses the functional impact of genomic mutations to detect cancer driver genes.

[25] Method: OncodriveCLUST detects the enrichment of somatic mutations to discover cancer driver genes

[26] The method can detect cancer genes with a large bias in clustered mutations. However, many mutated genes are not driver genes and may not regulate the progression of cancer, and not all cancer driver genes contain mutations. They may not contain mutations but regulate targets to develop cancer. Genes are also considered cancer drivers. For example, overexpression of KDM5C reduces p54 expression to enhance the proliferation and invasion of gastric cancer cells. KDM5C was identified as a cancer driver gene.

[27] .

[0005] The functions of molecular components in human cells are interdependent. Diseases are rarely caused by abnormalities in a single gene, but are reflected in the perturbations of complex biological networks. When genes interact in complex biological networks, the regulatory effects of mutant genes on gene products may be greater than the effects they produce themselves.

[28] With this in mind, some cancer driver gene identification methods have begun to integrate biological network information.

[29] The method integrates different biological data, including genomic data, influence diagrams of biological pathway information, and gene expression data, and constructs a bipartite graph of genes to detect the effects of mutant genes on genes with peripheral expression.

[30] Methods analyze the controllability of gene regulatory networks to guide key nodes in the network to discover cancer driver genes. Although these methods can elucidate the molecular mechanisms of cancer development at the network level, they require large data sets to produce reliable results.

[0006] It has been documented that cancer is driven by irregularities (including mutations) at different molecular levels, such as genetic, epigenetic, and transcriptomic levels, which can alter normal biological processes in a dynamic manner over time, giving rise to different types of cancer.

[31] Andres et al. recognized the dynamic characteristics of cancer and proposed the dynamic cancer driver gene (DCD)

[32] The DCD method assumes a causal relationship between driver genes and cancer progression and uses causal inference techniques to identify cancer driver genes. However, the intricate interactions between driver genes and cancer progression cannot be oversimplified as causal relationships. This complexity may hinder the accurate identification of cancer driver genes, resulting in a relatively low identification rate of the DCD method. Summary of the invention

[0007] The purpose of the present invention is to propose a method for identifying cancer driver genes based on single-cell pseudo-time series, in response to the fact that the existing technology ignores the dynamic information of cancer. This method is based on the constructed pseudo-time series information of cancer cell development, reveals the change pattern of genes in the pseudo-time series, and effectively identifies cancer driver genes by determining whether there are significant mutation points in the pseudo-time series of gene expression.

[0008] The technical solution for achieving the purpose of the present invention is:

[0009] The method for identifying cancer driver genes based on single-cell pseudo-time series includes the following steps:

[0010] 1) Extract characteristic genes:

[0011] Screening out important genes from scRNA-seq data, that is, extracting characteristic genes: Given a connected PPI network P = (V, E), each node v in the network P i ∈V represents the gene g i The protein generated by the expression, each undirected edge (v i ,v j )∈E represents the gene g i and g j The interaction relationship between proteins, where i≠j, so assuming g iIt represents genes and nodes in network P at the same time. Given a gene set G, a subnetwork of network P is constructed by taking the intersection of genes in network P and gene set G. For convenience of representation, the new subnetwork is still represented by P = (V, E). The matrix A |v|×|v| Represents the adjacency matrix of the subnetwork P, where element a ij See formula (1) for definition. Element a ij where i,j=1,2,...,|V|:

[0012]

[0013] First, the importance of each node in the sub-network P is ranked based on the eigenvector centrality. The higher the eigenvector centrality of the node, the greater its influence in the sub-network P. i ,i=1,2,...,|V|,v i The eigenvector centrality C EI [v i ] is defined as v i The weighted cumulative sum of the feature vector centrality of neighbor nodes is shown in formula (2):

[0014]

[0015] Where λ is a non-zero constant, formula (2) can be converted into matrix formula (3):

[0016] λC EI =C EI A, (3)

[0017] C EI is the eigenvector of the adjacency matrix A with eigenvalue λ, and C EI The top N genes in terms of eigenvector centrality are taken as the feature gene set;

[0018] 2) Determine the pseudo time series:

[0019] Given the characteristic gene expression matrix X m×N , X m×N The middle row represents a set of cell samples C = {c i |i=1,2,…,m}, the column represents a set of genes G={g j |j=1,2,…N}, using the PhenoPath framework

[33] Obtaining single-cell pseudo-time series: First, select a gene g that plays an important role in the process of cell carcinogenesis, and use the m-dimensional vector composed of the expression of the gene g in m cell samples as the time covariate, denoted by P g , based on the covariate P g , using the PhenoPath framework to analyze each cell sample ci Assign a pseudo-time score, denoted as Ph(c i ), i = 1, 2, ..., m. Rearrange the cell samples according to the score Ph(·), that is, if the cell sample coordinates k < r (k ≠ r), then Ph(c k ) ≤ Ph(c r ). Thus, a single-cell gene expression pseudo-time series matrix is obtained, denoted as X T×N , where T = m;

[0020] 3) Identify significant mutation points:

[0021] Given the single-cell gene expression pseudo-time series matrix X T×N , in X T×N , the columns represent a set of gene sets G = {g j | j = 1, 2, … N}, and the rows represent a set of cell samples C = {c i | i = 1, 2, …, T} sorted by pseudo-time. The matrix element x ij , i = 1, 2, ..., T, j = 1, 2, ..., N represents the expression level of gene g j in the c i cell sample. The column vector X[,j] = {x 1j , …, x Tj} represents the gene expression values of gene g j sorted by pseudo-time. X[,j] = {x 1j , …, x Tj} is called the pseudo-time series of the expression level of g j ;

[0022] Assume that a set of pseudo-time series of gene g expression levels x1, x2, …, x T is given. Randomly divide x1, x2, …, x T into and two parts. Let the distribution functions of the random variables in these two parts be F1(x) and F2(x) respectively. If F1(x) ≠ F2(x), it means that a significant mutation has occurred in the gene expression level at the time point. To identify this mutation point, design the following statistic:

[0023]

[0024] where e x is the natural exponential function, θ is a scaling factor used to adjust the function gradient and control the change rate of the natural exponential function when (x k - x i ) is close to 0. When the pseudo-time series follows a continuous distribution, the statistic U t,T can be obtained through the following recursive relationship:

[0025] U t,T =U t-1,T +V t,T , (5)

[0026] When t=2,3,…,T, we have:

[0027]

[0028] AndU 1,T =V 1,T , then the point τ that satisfies the following relationship is identified as the mutation point:

[0029] K τ =|U τ,T |=max|U t,T |, (7)

[0030] For the mutation point τ, the bootstrap resampling method is used to determine whether τ is a significant mutation point, as shown in formula (8):

[0031]

[0032] in is the bootstrap resampling sample According to the K of the mutation point τ Statistics, p boot The calculation is based on the original data K τ Statistics and statistics obtained by resampling The comparison shows that if p boot <0.05, then the mutation point τ is a significant mutation point;

[0033] By identifying the pseudo-time matrix X of single-cell gene expression T×N Each gene g j ,j=1,2,...,N, whether there is a significant mutation point in the pseudo time series, that is, if p boot <0.05, then g j Identified as a cancer driver gene;

[0034] The following is an algorithm for identifying significant mutation points, as shown in the following algorithm:

[0035]

[0036]

[0037] This technical solution uses scRNA-seq data, which can more accurately capture differential expression between cells compared to traditional methods based on batch data, introduces pseudo-time analysis, and better reflects the dynamic evolution characteristics of cancer. In addition, this technical solution proposes an innovative time series change point recognition algorithm to identify change points that have changed significantly in single-cell pseudo-time series, and uses bootstrap resampling to test the significance of the identified change points to determine cancer driver genes. The method proposed in this technical solution has been effectively verified in real single-cell cancer datasets and can accurately identify driver genes related to the cancer development process. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 is a flow chart of an embodiment;

[0039] Figure 2 This is a schematic diagram of the gene expression change of the gene ERBB2 on the pseudo time axis in the embodiment;

[0040] Figure 3 is a comparison diagram of the results obtained by the embodiment and the DCD method;

[0041] Figure 4 This is a schematic diagram of the gene expression level change of the gene ESR1 in the embodiment on the pseudo time axis;

[0042] Figure 5 Schematic diagram of the change in gene expression level of gene BTG1 on the pseudo time axis in the example. DETAILED DESCRIPTION

[0043] The content of the present invention is further described below in conjunction with the drawings and embodiments, but the present invention is not limited thereto.

[0044] Example:

[0045] 1. Data preparation for this example:

[0046] Download the real single-cell breast cancer transcriptomics data from the public database (National Center for Biotechnology Information Gene Expression Omnibus, referred to as NCBI GEO), the accession number is GSE75688

[34] . The expression value TPM of the downloaded single-cell gene expression matrix is ​​converted using log(TPM+1). In the method steps, the converted log(TPM+1) value is called the gene expression value. The single-cell breast cancer gene expression data contains 317 cancer cell samples and 10029 genes. The data is preprocessed to delete genes that are not expressed in more than 20% of the cells. The preprocessed single-cell gene expression data contains 317 cancer cell samples and 9551 genes. Download Vinayagam et al.

[35] The human protein-protein interaction network provided, referred to as the PPI network, contains 5436 nodes and 34814 edges.

[0047] 2. Setting parameters for this example:

[0048] Number of characteristic genes N = 1381;

[0049] Time covariate gene g = 'ERBB2';

[0050] Scaling factor parameter

[0051] parameter

[0052] Reference Figure 1 , a cancer driver gene identification method based on single-cell pseudo-time series, comprising the following steps:

[0053] 1) Extract characteristic genes:

[0054] Screen out important genes from scRNA-seq data, that is, extract characteristic genes. Given a connected PPI network P = (V, E), each node v in the network P i ∈V represents the gene g i The protein generated by the expression, each undirected edge (v i ,v j )∈E represents the gene g i and g j The interaction relationship between proteins, i≠j, so assuming g i It represents both genes and nodes in network P. Given a gene set G, a subnetwork of network P is constructed by taking the intersection of genes in network P and gene set G. For convenience, the new subnetwork is still represented by P = (V, E). Subnetwork P contains 3453 nodes and 16996 edges. Matrix A |V|×|V| Represents the adjacency matrix of the subnetwork P, where element a ij See formula (1) for definition. Element a ij where i,j=1,2,...,|V|:

[0055]

[0056] First, based on eigenvector centrality, the importance of each node in sub-network P is ranked. The higher the value of the eigenvector centrality of a node, the greater its influence in sub-network P. Given node v i , i = 1, 2, ..., |V|, v i 's eigenvector centrality C EI [v i is defined as the weighted sum of the eigenvector centralities of v i 's neighbor nodes, as shown in formula (2):

[0057]

[0058] where λ is a non-zero constant, and formula (2) can be transformed into matrix formula (3):

[0059] λC EI = C EI A, (3)

[0060] C EI is the eigenvector corresponding to the eigenvalue λ of the adjacency matrix A. The top 1381 genes ranked by eigenvector centrality in C EI are used as the feature gene set;

[0061] 2) Determine the pseudo-time series:

[0062] Given the obtained feature gene expression matrix X m×N , in this example, m = 317 and N = 1381. The rows in X m×N represent a set of cell sample collections C = {c i |i = 1, 2, …, m}, and the columns represent a set of gene collections G = {g j |j = 1, 2, … N}. Use the PhenoPath framework

[33] to obtain the single-cell pseudo-time series: First, select a gene ERBB2 that plays an important role in the process of cell carcinogenesis. The m-dimensional vector composed of the expression levels of ERBB2 in m cell samples is used as the time covariate, denoted as P ERBB2 . Based on the covariate P ERBB2 , the PhenoPath framework scores each cell sample c i , denoted as Ph(c i ), i = 1, 2, ..., m. Rearrange the cell samples according to the scores Ph(·), that is, if the cell sample coordinates k < r (k ≠ r), then Ph(c k ) ≤ Ph(c r ). Thus, the single-cell gene expression pseudo-temporal matrix is obtained, denoted as X T×N , where T = m. Figure 2 is the gene expression change of ERBB2 on the pseudo time axis;

[0063] 3) Identify significant mutation points:

[0064] Given a single-cell gene expression pseudo-time matrix X T×N In this example, T = 317, N = 1381. T×N The middle column represents a set of genes G = {g j |j=1,2,…N}, and the rows represent the set of cell samples C={c i |i=1,2,…,T}, matrix element x ij , i=1,2,...,T, j=1,2,...,N represents gene g j In c i The expression level on the cell sample, column vector X[,j] = {x 1j ,…,x Tj} represents gene g j Gene expression values ​​sorted by pseudo time, X[,j] = {x 1j ,…,x Tj} is called g j Pseudo time series of expression;

[0065] Assume that a set of pseudo time series x1, x2, …, x2 of gene expression is given. T , x1,x2,…,x T Arbitrary and Let the distribution functions of the random variables in these two parts be F1(x) and F2(x) respectively. If F1(x)≠F2(x), it means that the gene expression level is A significant mutation occurred at the time point. To identify the mutation point, the following statistics are designed:

[0066]

[0067] Where e is the natural exponential function, is a scaling factor used to adjust the function gradient and control the natural exponential function in (x k -x i ) is close to 0. In this example, Set to 0.5. When the pseudo time series follows a continuous distribution, the statistic U t,T It can be obtained through the following recursive relationship:

[0068] U t,T =U t-1,T +V t,T , (5)

[0069] When t=2,3,…,T, we have:

[0070]

[0071] AndU 1,T =V 1,T , then the point τ that satisfies the following relationship is identified as the mutation point:

[0072] K τ =|U τ,T |=max|U t,T |, (7)

[0073] For the mutation point τ, the bootstrap resampling method is used to determine whether τ is a significant mutation point, as shown in formula (8):

[0074]

[0075] in is the bootstrap resampling sample According to the K of the mutation point τ Statistics, In this example Set to 1000. boot The calculation is based on the original data K τ Statistics and statistics obtained by resampling The comparison shows that if p boot <0.05, then the mutation point τ is a significant mutation point;

[0076] By identifying the pseudo-time matrix X of single-cell gene expression T×N Each gene g j ,j=1,2,...,N, whether there is a significant mutation point in the pseudo time series to determine g j Is it a cancer driver gene? boot <0.05, then g j Identified as a cancer driver gene;

[0077] The following is an algorithm for identifying significant mutation points, as shown in the following algorithm:

[0078]

[0079]

[0080] The experimental results obtained in this example are compared with the COSMIC Cancer Gene Census (CGC) benchmark database for the accuracy of cancer driver gene identification, and with the DCD

[32] Compared with the results obtained by the method, the experimental results of this example have better recognition accuracy. Figure 3 As shown in Figure 2, the CGC cancer driver gene identified in this example, gene ESR1 is a cancer driver gene that can be identified by both the DCD method and this example method. The gene expression level change of gene ESR1 on the pseudo time axis is shown in Figure 2. Figure 4 As shown; the gene BTG1 is not identified by the DCD method, but it is a cancer driver gene identified by this method. The gene expression changes of the gene BTG1 on the pseudo time axis are shown as follows Figure 5 As shown, this example can not only discover CGC cancer driver genes that can be identified by other methods, but also identify CGC cancer driver genes that cannot be identified by other methods, and Figure 4 and Figure 5 The middle dotted line shows that this example can find cell sample points where gene expression values ​​mutate over pseudo time, and can effectively identify cancer driver genes in combination with cancer dynamic characteristics.

[0081] References:

[0082] [1] Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz Jr LA, Kinzler KW. Cancergenome landscapes. science. 2013;339(6127):1546–1558.

[0083] [2] Paik S, Shak S, Tang G, Kim C, Baker J, Cronin M, et al. A multigene assay to predict recurrence of tamoxifen-treated, node-negative breastcancer. New England Journal of Medicine. 2004; 351(27): 2817–2826.

[0084] [3]Gray RG,Quirke P,Handley K,Lopatin M,Magill L,Baehner FL,etal.Validation study of aquantitative multigene reverse transcriptase–polymerase chain reaction assay forassessment of recurrence risk in patientswith stage II colon cancer.Journal of clinicaloncology.2011;29(35):4611–4619.

[0085] [4]Weinstein JN,Collisson EA,Mills GB,Shaw KR,Ozenberger BA,EllrottK,et al.The cancergenome atlas pan-cancer analysis project.Naturegenetics.2013;45(10):1113–1120.

[0086] [5]Dressler L,Bortolomeazzi M,Keddar MR,Misetic H,Sartini G,Acha-Sagredo A,et al.

[0087] Comparative assessment of genes driving cancer and somatic evolutionin non-cancertissues:an update of the Network of Cancer Genes(NCG)resource.Genome biology. 2022;23(1):35.

[0089] [6]Sondka Z,Bamford S,Cole CG,Ward SA,Dunham I,Forbes SA.The COSMICCancerGene Census:describing genetic dysfunction across all humancancers.Nature ReviewsCancer.2018;18(11):696–705.

[0090] [7]Nawy T.Single-cell sequencing.Nature methods.2014;11(1):18–18.

[0091] [8]Jaitin DA,Kenigsberg E,Keren-Shaul H,Elefant N,Paul F,Zaretsky I,et al.Massivelyparallel single-cell RNA-seq for marker-free decomposition oftissues into cell types.

[0092] Science.2014;343(6172):776–779.

[0093] [9]Gierahn TM,Wadsworth MH,Hughes TK,Bryson BD,Butler A,Satija R,etal.Seq-Well:

[0094] portable,low-cost RNA sequencing of single cells at highthroughput.Nature methods.

[0095] 2017;14(4):395–398.

[0096]

[10] Zheng GX,Terry JM,Belgrader P,Ryvkin P,Bent ZW,Wilson R,etal.Massively paralleldigital transcriptional profiling of single cells.Naturecommunications.2017;8(1):14049.

[11] Oesper L,Mahmoody A,Raphael BJ.THetA:inferring intra-tumor heterogeneity fromhigh-throughput DNA sequencingdata.Genome biology.2013;14:1–21.

[0097]

[12] Roth A,Khattra J,Yap D,Wan A,Laks E,Biele J,et al.PyClone:statistical inference ofclonal population structure in cancer.Naturemethods.2014;11(4):396–398.

[0098]

[13] Navin N,Kendall J,Troge J,Andrews P,Rodgers L,McIndoo J,etal.Tumour evolutioninferred by single-cell sequencing.Nature.2011;472(7341):90–94.

[0099]

[14] Pollen AA,Nowakowski TJ,Chen J,Retallack H,Sandoval-Espinosa C,Nicholas CR,et al.Molecular identity of human outer radial glia duringcortical development.Cell.2015;163(1):55–67.

[0100]

[15] Qi R,Ma A,Ma Q,Zou Q.Clustering and classification methods forsingle-cell RNA-sequencing data.Briefings in bioinformatics.2020;21(4):1196–1208.

[0101]

[16] Lake BB,Ai R,Kaeser GE,Salathia NS,Yung YC,Liu R,et al.Neuronalsubtypes anddiversity revealed by single-nucleus RNA sequencing of the humanbrain.Science.2016;352(6293):1586–1590.

[0102]

[17] Patel AP,Tirosh I,Trombetta JJ,Shalek AK,Gillespie SM,Wakimoto H,et al.Single-cellRNA-seq highlights intratumoral heterogeneity in primaryglioblastoma.Science.2014;344(6190):1396–1401.

[0103]

[18] Guo X,Zhang Y,Zheng L,Zheng C,Song J,Zhang Q,et al.Globalcharacterization of Tcells in non-small-cell lung cancer by single-cellsequencing.Nature medicine.2018;24(7):978–985.

[0104]

[19] Peng J,Sun BF,Chen CY,Zhou JY,Chen YS,Chen H,et al.Single-cellRNA-seqhighlights intra-tumoral heterogeneity and malignant progression inpancreatic ductaladenocarcinoma.Cell research.2019;29(9):725–738.

[0105]

[20] Pham VVH,Lin L,Bracken C,et al.Computational methods for cancerdriver discovery:A survey.Theranostics 2021;11:5553–68.

[0106]

[21] Dimitrakopoulos CM,Beerenwinkel N.Computational approaches forthe identificationof cancer genes and pathways.Wiley Interdiscip Rev SystBiol Med.2017;9(1):e1364

[0107]

[22] Linehan WM,Srinivasan R,Schmidt LS.The genetic basis of kidneycancer:a metabolicdisease.Nat Rev Urol.2010;7:277

[0108]

[23] Lawrence MS,Stojanov P,Polak P,Kryukov GV,Cibulskis K,SivachenkoA.et al.Mutational heterogeneity in cancer and the search for new cancer-associated genes.Nature.2013;499(7457):214-218

[0109]

[24] Gonzalez-Perez A,Lopez-Bigas N.Functional impact bias revealscancer drivers.Nucleic Acids Res.2012;40(21):e169-e169

[0110]

[25] Vandin F.Computational methods for characterizing cancermutational heterogeneity.Front Genet.2017;8:83

[0111]

[26] Tamborero D,Gonzalez-Perez A,Lopez-Bigas N.OncodriveCLUST:exploiting thepositional clustering of somatic mutations to identify cancergenes.Bioinformatics.2013;29(18):2238-44

[0112]

[27] Xu L,Wu W,Cheng G,Qian M,Hu K,Yin G.et al.Enhancement ofproliferation and invasion of gastric cancer cell by KDM5C via decrease inp53 expression.Technol CancerRes Treat.2017;16(2):141-149

[0113]

[28] Barabási A-L,Gulbahce N,Loscalzo J.Network medicine:a network-based approach tohuman disease.Nat Rev Genet 2011;12:56–68.

[0114]

[29] Bashashati A,Haffari G,Ding J,Ha G,Lui K,Rosner J.etal.DriverNet:uncovering theimpact of somatic driver mutations ontranscriptional networks in cancer.Genome Biol.

[0115] 2012;13(12):R124-R124

[0116]

[30] Pham VVH,Liu L,Bracken CP,Goodall GJ,Long Q,Li J.et al.CBNA:Acontrol theorybased method for identifying coding and non-coding cancerdrivers.PLoS Comput Biol.

[0117] 2019;15(12):e1007538

[0118]

[31] Fan J,Slowikowski K,Zhang F.Single-cell transcriptomics incancer:computationalchallenges and opportunities.Exp Mol Med 2020;52(9):1452–65.

[0119]

[32] Cifuentes-Bernal AM,Pham VVH,Li X et al.Dynamic cancer drivers:acausal approachfor cancer driver discovery based on bio-pathologicaltrajectories.Brief Funct Genomics

[0120] 2022;21:455–65.

[0121]

[33] Campbell KR,Yau C.Uncovering pseudotemporal trajectories withcovariates fromsingle cell and bulk expression data.Nat Commun 2018;9(1):1–12.

[0122]

[34] Chung W,Eum HH,Lee H-O,et al.Single-cell rna-seq enablescomprehensive tumourand immune cell profiling in primary breast cancer.NatCommun 2017;8(1):15081.

[0123]

[35] Vinayagam A,Stelzl U,Foulle R,et al.A directed proteininteraction network forinvestigating intracellular signal transduction.SciSignal 2011;4(189):rs8–8。

Claims

1. A method for identifying cancer driver genes based on single-cell pseudo-time series, characterized in that: The steps include: 1) Extract characteristic genes: Screening out important genes from scRNA-seq data, that is, extracting characteristic genes: Given a connected PPI network P = (V, E), each node v in the network P i ∈V represents the gene g i The protein generated by the expression, each undirected edge (v i , v j )∈E represents the gene g i and g j The interaction relationship between proteins, where i≠j, so assuming g i It represents genes and nodes in network P at the same time. Given a gene set G, a subnetwork of network P is constructed by taking the intersection of genes in network P and gene set G. For convenience of representation, the new subnetwork is still represented by P = (V, E). The matrix A |V|×|V| Represents the adjacency matrix of the subnetwork P, where element a ij See formula (1) for definition, element a ij where i, j = 1, 2, ..., |V|: First, the importance of each node in the sub-network P is ranked based on the eigenvector centrality. The higher the eigenvector centrality of the node, the greater its influence in the sub-network P. Given a node v i , i = 1, 2, ..., |V|, v i The eigenvector centrality C EI [v i ] is defined as v i The weighted cumulative sum of the feature vector centrality of neighbor nodes is shown in formula (2): Where λ is a non-zero constant, formula (2) can be converted into matrix formula (3): λC EI =C EI A, (3) C EI is the eigenvector of the adjacency matrix A with eigenvalue λ, and C EI The top N genes in terms of eigenvector centrality are taken as the feature gene set; 2) Determine the pseudo time series: Given the obtained characteristic gene expression matrix X m×N , X m×N where the rows represent a set of cell sample collections C = {c i | i = 1, 2,..., m}, and the columns represent a set of characteristic gene sets G = {g j | j = 1, 2,..., N}, use the PhenoPath framework to obtain single-cell pseudotime series: First, select a gene g that plays an important role in the process of cell carcinogenesis, and use the m-dimensional vector composed of the expression levels of this gene g on m cell samples as the time covariate, denoted as P g . Based on the covariate P g , use the PhenoPath framework to assign a pseudotime score to each cell sample c i , denoted as Ph(c i ), i = 1, 2,..., m. Rearrange the cell samples according to the score Ph(·), that is, if the cell sample coordinates k < r, k ≠ r, then Ph(c k ) ≤ Ph(c r ). Thus, a single-cell gene expression pseudotemporal matrix is obtained, denoted as X T×N , where T = m; 3) Identify significant mutation points: Given a single-cell gene expression pseudo-time matrix X T×N , X T×N The middle column represents a set of genes G = {g j |j=1, 2, ...N}, and the rows represent the set of cell samples sorted by pseudo time C={c i |i=1,2,...,T},matrix element x ij , i = 1, 2, ..., T, j = 1, 2, ..., N represents gene g j In c i The expression level on the cell sample, column vector X[,j] = {x 1j , ..., x Tj } represents gene g j Gene expression values ​​sorted by pseudo time, X[,j] = {x 1j , ..., x Tj } is called gene g j Pseudo time series of expression; Assume that a set of pseudo time series x1, x2, ..., x2 of gene expression is given. T , x1, x2, …, x T Arbitrary and Let the distribution functions of the random variables in these two parts be F1(x) and F2(x) respectively. If F1(x)≠F2(x), it means that the gene expression level is A significant mutation occurred at the time point. To identify the mutation point, the following statistics are designed: where e x is the natural exponential function, is a scaling factor used to adjust the function gradient and control the natural exponential function in (x k -x i ) is close to 0. When the pseudo time series follows a continuous distribution, the statistic U t，T It can be obtained through the following recursive relationship: U t，T =U t-1，T +V t，T , (5) When t=2,3,...,T,we have: And U 1，T =V 1，T , then the point τ that satisfies the following relationship is identified as the mutation point: K τ =|U t，T |=max|U t，T |, (7) For the mutation point τ, the bootstrap resampling method is used to determine whether τ is a significant mutation point, as shown in formula (8): in is the bootstrap resampling sample According to the K of the mutation point τ Statistics, p boot The calculation is based on the original data K τ Statistics and statistics obtained by resampling The comparison shows that if p boot <0.05, then the mutation point τ is a significant mutation point; By identifying the pseudo-time matrix X of single-cell gene expression T×N Each gene g j , j = 1, 2, ..., N, whether there is a significant mutation point in the pseudo time series, that is, if p boot <0.05, then g j Identified as a cancer driver gene; The following is an algorithm for identifying significant mutation points, as shown in the following algorithm:

Images