Cancer prognosis survival analysis method based on collaborative co-expression network

Abstract

The invention relates to a cancer prognosis survival analysis method based on a collaborative co-expression network. The invention relates to the technical field of cancer tumor prognosis analysis, and the method comprises the following steps: collecting cancer multi-omics data based on a public database, and preprocessing the data; based on the processed data, constructing a collaborative co-expression network oriented to the cancer multi-omics data; based on the collaborative co-expression network, a prognosis biomarker is extracted; and according to the extracted prognosis biomarker, carrying out sample clustering for survival analysis, and carrying out prognosis survival analysis. According to the invention, accurate prediction of prognosis categories of cancer patients is realized, and statistical support is provided for biologists to reveal molecular mechanisms causing cancer prognosis survival differences and to assist clinicians to finally realize targeted accurate treatment for prognosis biomarkers.

Description

Technical Field

[0001] The present invention relates to the technical field of cancer tumor prognosis analysis, and is a cancer prognosis survival analysis method based on a collaborative co-expression network. Background Art

[0002] Malignant tumors (commonly known as cancer, hereinafter referred to as such) have the characteristics of high recurrence rate, fast growth rate and low cure rate, which seriously endangers human survival and health. Surgery supplemented by postoperative radiotherapy and chemotherapy is the main method for treating cancer at present. However, the invasion and metastasis characteristics of cancer cells limit the ability of surgery to remove the lesions and cure cancer once, and radiotherapy and chemotherapy bring pain to cancer patients because of their side effects. Prognostic survival analysis can determine the risk of postoperative recurrence, metastasis and even death of patients, provide necessary basis for individualized treatment of cancer patients after surgery, and maximize the survival time of patients. Prognostic survival analysis refers to: constructing an association model between the phenotype or molecular indicators of existing cases and their survival time, and measuring the phenotype or molecular indicators of new patients according to the model to predict their survival time. In fact, cancer recurrence, metastasis, and patient death all have their corresponding time nodes. The present invention refers to the time period from postoperative to the occurrence of the above events as survival time.

[0003] Studies have found that there is no obvious correlation between the phenotype indicators (such as personal information, clinical indicators, etc.) of cancer patients and their survival time. Therefore, constructing a two-dimensional matrix (rows represent different features and columns correspond to different samples) based on multi-omics (such as genomics, transcriptomics, proteomics, etc.) data and extracting features closely related to survival time (i.e. biomarkers, such as DNA, RNA, protein, etc.) from them is a key step in solving cancer prognosis survival analysis. The above process is called feature selection in the field of pattern recognition and machine learning. First, the high-dimensional and small sample characteristics of the constructed two-dimensional matrix lead to the failure of the classic feature selection method in the field of pattern recognition. Limited by the number of cases collected and the cost of collecting data, the number of samples (also known as patients, subjects) is often at least two orders of magnitude different from the number of features (also known as variables), and traditional feature selection methods inevitably face overfitting problems. In addition, taking a certain omics matrix of cancer patients as input and their prognostic survival time as output further increases the difficulty of feature selection. According to the follow-up situation, the patient's survival time corresponds to the time from surgery to the occurrence of an event (such as recurrence, metastasis, death, etc.) or the last follow-up. If the patient is lost to follow-up or is still alive at the end of the experiment, what is recorded is only the patient's follow-up time (right-censored time) rather than the actual survival time.

[0004] In order to extract features that are closely related to the prognosis and survival time of cancer patients from matrices representing different omics, existing research focuses on the following two aspects: ① The idea of ​​screening associated features. The correlation between each feature and the prognosis and survival time is calculated to measure whether the feature is selected, or the noisy features are filtered out using constraints, and only the features that are closely related to the prognosis and survival are retained. ② The idea of ​​combining feature clusters. Closely related features are combined into feature clusters to represent prognostic biomarkers. Both of the above ideas have problems: the former emphasizes the correlation between features and survival time, but ignores the relationship between features; the latter focuses on the correlation between features, ignoring the connection between features and survival time.

[0005] For cancer prognosis survival analysis, there are two urgent problems to be solved: prognostic biomarker extraction and prognostic category discovery. The former aims to select key features closely related to prognostic survival for subsequent survival index construction; the latter aims to discover possible prognostic subtypes in specific cancers or pan-cancers based on different prognostic survival situations, which can be used to predict the patient's survival risk category and provide assistance for individualized clinical treatment. The two are mutually causal: prognostic biomarker extraction can remove redundant features of the two-dimensional matrix, and sample clustering on the two-dimensional matrix of the low-dimensional feature space can help improve the stability and reliability of the discovered prognostic categories; prognostic category discovery divides samples into two or more categories based on the prognostic survival time, and attempts to extract biomarkers with significant differences between the two or more categories as key features of prognostic survival. Research on prognostic category discovery mainly focuses on: univariate sample classification strategy and univariate sample clustering strategy. Research on prognostic biomarker extraction is concentrated on: multivariate feature selection strategy, univariate feature combination strategy, and related feature clustering strategy.

[0006] How to select key biomarkers closely related to cancer prognosis and survival and discover possible prognosis and survival subtypes of specific cancers, these two problems have not been well solved so far, which seriously limits the clinical application of cancer prognosis and survival analysis research. To solve these two problems, existing studies focus on the relationship between features and survival time or features, and focus on the above five strategies respectively. There is an urgent need for an effective method that organically combines the above two focuses and five strategies. First, extract prognosis biomarkers, then cluster samples on the output results of their corresponding feature subspace matrix or regression model, and then predict the survival risk and even survival time of new samples, which conforms to the logical order of cognition. Considering that survival time is continuous, sample clustering strategy is better than sample classification strategy. Compared with univariate or multivariate, it is more important to select prognosis biomarkers that are closely related to prognosis and survival time while paying attention to the correlation between prognosis biomarkers. The relationship between the two has become a "blind spot" in existing research. In addition, no research work has been found that uses cancer multi-omics data as the research object and involves causal relationship inference between prognosis and survival analysis biomarkers. Summary of the invention

[0007] In response to the above problems, the present invention proposes a cancer prognosis survival analysis method based on a collaborative co-expression network. The proposed method is intended to be applied to the prognosis survival analysis of multi-omics data of malignant tumors such as malignant brain glioma, breast cancer, primary liver cancer, gastric adenocarcinoma and esophageal squamous cell carcinoma, select biomarkers closely related to prognosis survival, construct prognosis survival indicators, and use them to discover and compare different prognosis categories of samples, infer the causal relationship of prognosis biomarkers, and achieve accurate prediction of the prognosis category of cancer patients, so as to provide statistical support for biologists to reveal the molecular mechanism that leads to differences in cancer prognosis survival and assist clinicians to ultimately achieve targeted precision treatment for prognosis biomarkers.

[0008] The present invention provides the following technical solutions:

[0009] A cancer prognosis survival analysis method based on a collaborative co-expression network comprises the following steps:

[0010] Step 1: Collect multi-omics data of cancer based on public databases and pre-process the data;

[0011] Step 2: Based on the processed data, construct a collaborative co-expression network for cancer multi-omics data;

[0012] Step 3: Extract prognostic biomarkers based on the synergistic co-expression network;

[0013] Step 4: Based on the extracted prognostic biomarkers, sample clustering for survival analysis is performed, and prognostic survival analysis is performed.

[0014] Preferably, the step 1 is specifically:

[0015] Step 1.1: For the target tumor, collect high-quality RNA-seq data of FPKM expression type from the TCGA public database, convert the expression value to log2 TPM+1, and perform normalization to ensure that the overall expression of the data conforms to the normal distribution;

[0016] Step 1.2: Collect DNA methylation data of target tumor samples in the TCGA database, give priority to high-throughput methylation sequencing data, and use Minfi and ChAMPR packages for background correction, normalization, probe filtering, and batch effect correction;

[0017] Step 1.3: For the target tumor, collect high-quality CNV data of Segment Mean type from the TCGA public database through Genomic Data Commons and convert the data into GISTIC2 format; perform data cleaning methods including removing outliers, correcting errors, and removing redundant information on the CNV data, and then standardize the CNV data to eliminate the impact of experimental batch effects and sample processing differences, and use Z-score standardization to standardize the CNV data;

[0018] Step 1.3: Collect ChIP-Seq FASTQ format data of target tumor samples in the TCGA database, use FastQC to perform quality control on the raw sequencing data to identify and deal with low-quality reads, contamination or other potential problems, and use Trimmomatic or Cutadapt tools to remove adapter sequences and low-quality reads; use alignment tools to align the cleaned reads to the reference genome, use SAMtools to filter out reads with low alignment quality, multiple-point alignments, and non-specific alignments, and use Picard Tools or SAMtools to remove PCR duplications to ensure data accuracy.

[0019] Preferably, the step 2 is specifically:

[0020] Cox regression based on the joint enumeration of two features is used to achieve the fitting of the prognosis survival time, and Pearson correlation based on the joint enumeration of two features is used to measure the correlation between features. The above two metrics are unified to mark the weights of any two node edges in the fully connected network; the cluster center automatic search clustering algorithm based on density descending is used to delete the edges with lower weights in the network, so as to realize network construction and feature cluster mining. The cluster center automatic search algorithm based on density descending is used to cluster the weights of all edges in the network in density descending order, and delete the edges with lower weight density from the network until a disconnected feature cluster appears;

[0021] Step 2.1: Based on the Cox regression, the Cox proportional hazard regression model is used to fit the survival time by enumerating the pairwise features, so as to reveal the significant relationship between any pair of features and prognosis survival. By performing the maximum likelihood estimation on the logarithmic partial likelihood function of the hazard function, the regression coefficient β corresponding to any pair of features is obtained. i ,β j ) T , when it is assumed that the features are uncorrelated, any component of the regression coefficient obeys the Gaussian distribution, and the Wald test is used to determine whether the regression coefficient component deviates significantly from zero;

[0022] The sample size can be further expanded by reordering the samples and their corresponding survival results. The reordered p-value corresponding to the Wald statistic can be generated on each feature. All pairwise features are enumerated to obtain paired p-values, which are used to select significant feature pairs closely related to survival from all pairwise feature combinations. In order to deal with overfitting, a penalty constraint is added to the logarithmic partial likelihood function corresponding to the risk function. Considering that the hazard ratio may change with survival time, a weighted Cox regression model is used instead.

[0023] Step 2.2: Pearson correlation-based joint enumeration of two features. To measure the correlation between two features, the Pearson correlation coefficient r is used to calculate the correlation between any two features, and the t statistic of n samples is constructed based on this. The statistic is expressed as: By using the sample reordering method, the order of sample values ​​of one feature in any pair of features is disrupted, and the above t statistic is calculated repeatedly. The corresponding reordered p value is expressed as:

[0024]

[0025] Among them, t(i,j) is the statistic corresponding to the correlation coefficient r of features i and j, t b (i,j) represents the statistic generated by randomly reordering the values ​​of any feature, and B is the number of random reorderings;

[0026] Step 2.3: Determine the association relationship of network nodes. The collaborative co-expression network uses a single feature as a network node and the relationship between two features as an edge. The weight of the edge between nodes is calculated by comprehensively considering the significant relationship between the feature and the survival time and the correlation between the features. The following metric is constructed to measure the weight of the edge between nodes:

[0027]

[0028] Where p2(i) and p2(j) represent the binary Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively, and p1(i) and p1(j) represent the univariate Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively. represents rounding up; the first term of m(i,j) is the "synergistic" part that measures the contribution of features i and j to survival regression, and the second term is derived from the reordered p-value obtained by the joint enumeration of two features based on Pearson correlation, corresponding to the "co-expression" part; m(i,j) is used to measure the association between nodes i and j: the larger its value, the more likely it is that nodes i and j are derived from the same feature cluster, and vice versa. When the value of m(i,j) is less than zero, its value is set to the maximum value of the weight in the network, and the weighted undirected graph constructed is called a synergistic-co-expression network;

[0029] Step 2.4: Based on the feature cluster mining of the collaborative-co-expression network, for the weight of the collaborative-co-expression network edge, which is a measure of the node association relationship, an automatic cluster center search clustering algorithm based on density descending is introduced. The measurement threshold of the edge in the network is limited by the optimal automatic estimation of the number of categories and cluster boundaries to determine whether there is an associated edge between any two nodes in the network, thereby degenerating the network from a weighted undirected graph to an undirected graph.

[0030] Preferably, the step 3 is specifically:

[0031] Step 3.1: Bottom-up “cross-cluster enumeration” biomarker selection. According to the “cross-cluster hypothesis”, enumerate the univariate within each feature cluster, combine the univariate from c different feature clusters into a multivariate, and calculate its regression coefficient on the Cox survival hazard ratio function;

[0032] Step 3.2: Top-down "cross-cluster random" biomarker selection, introducing a sample resampling voting strategy, randomly selecting a single variable within each feature cluster, combining the single variables from c different feature clusters into multivariates, and calculating their regression coefficients for the Cox survival hazard ratio function; assuming that the features are uncorrelated, analogous to the pairwise feature joint enumeration method based on Cox regression, each group of "cross-cluster" multivariates can obtain corresponding c re-ranked p-values, and each round of multivariate voting is performed on the multivariates with significant re-ranked p-values; after l rounds of voting, the features with the highest votes are selected from each feature cluster to form a significant prognostic feature set;

[0033] Step 3.3: Unsupervised feature selection of joint multivariates; single variables need to be extracted from different feature clusters to form joint multivariates for regression analysis. For bottom-up "cross-cluster enumeration" biomarker selection, c re-ranked p-values ​​corresponding to each enumerated multivariate are recorded; for top-down "cross-cluster random" biomarker selection, multiple groups of c-dimensional p-value vectors corresponding to multiple randomly extracted multivariates can be recorded by setting a threshold. In the c-dimensional feature space, all enumerated or re-ranked p-value vectors that meet the threshold limit are automatically searched and clustered based on the descending density of the cluster center, and the cluster closest to the origin is selected to form a c-dimensional feature set. The Mahalanobis distance from the scattered points to the origin is calculated. When the bottom-up "cross-cluster enumeration" biomarker selection is adopted, the multidimensional variable with the closest distance to the origin is selected as the significant prognostic feature set; if the top-down "cross-cluster random" biomarker selection is adopted, it is still necessary to determine whether the multidimensional variable with the closest distance to the origin corresponds to the feature combination with the highest cumulative votes;

[0034] Step 3.4: Construction of prognostic survival analysis indicators. Taking the Cox proportional hazard regression model as an example, since it is difficult to estimate the baseline risk function or baseline survival function, the linear part of the risk score is used to construct the prognostic survival indicator. The risk score of the sample is a linear combination of the corresponding sample values ​​of the prognostic characteristics, and the combination coefficient is the regression coefficient.

[0035] Preferably, the step 4 is specifically:

[0036] Step 4.1: Sample clustering for prognostic survival analysis: To discover the prognostic survival category, the key biomarkers from different omics or the prognostic survival analysis indicators built by different omics data are combined as the input of the convolutional neural network, and its output is regarded as a new feature subspace. The samples are automatically searched and clustered based on the cluster center in descending density, so as to discover the prognostic category to which the samples belong;

[0037] Step 4.2: Draw Kaplan-Meier survival analysis curves for samples of different categories to observe the changes in survival time within the category and the differences in survival time between categories. Use the log-rank test to measure whether there are significant quantitative differences in survival time between different categories. Draw a risk score survival analysis graph to observe the changes in survival time between samples of different categories.

[0038] Step 4.3: Analysis of significant differences in prognosis between classes. Class comparison is achieved using statistics based on sample distribution differences to evaluate whether there are significant differences in risk scores for different categories. On this basis, attempts are made to construct different types of classifiers using training samples, and the accuracy and effectiveness of the classification results are verified using the prognosis survival results of the test samples.

[0039] Preferably, the cross-group assumption is as follows:

[0040] The significant feature set closely related to the prognosis survival time is not entirely composed of a single significant variable related to the prognosis survival time, and there is no obvious correlation between the variables in the significant feature set.

[0041] A cancer prognosis survival analysis system based on a collaborative co-expression network, the system comprising:

[0042] A data acquisition module, which collects multi-omics data of cancer based on a public database and pre-processes the data;

[0043] A network construction module, wherein the network construction module constructs a collaborative co-expression network for cancer multi-omics data based on the processed data;

[0044] An extraction module, wherein the extraction module extracts prognostic biomarkers based on a collaborative co-expression network;

[0045] A clustering module, wherein the clustering module performs sample clustering for survival analysis based on extracted prognostic biomarkers;

[0046] An analysis module is used to perform prognosis survival analysis based on the results of sample clustering.

[0047] A computer-readable storage medium stores a computer program, which is executed by a processor to implement a cancer prognosis survival analysis method based on a collaborative co-expression network.

[0048] A computer device comprises a memory and a processor, wherein the memory stores a computer program, and the processor implements a cancer prognosis survival analysis method based on a collaborative co-expression network when executing the computer program.

[0049] The present invention has the following beneficial effects:

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

[0051] The present invention is applied to the prognosis and survival analysis of cancer. On the basis of realizing the selection of prognostic biomarkers, the establishment of survival indicators, the discovery and comparison of prognostic categories, and the prediction of categories, it provides statistical support for revealing the molecular mechanisms that lead to differences in cancer prognosis and survival and ultimately realizing targeted precision treatment for key biomarkers. The innovation of the present invention is mainly reflected in the following two aspects:

[0052] (1) The concept of synergy-coexpression is proposed, and the correlation metrics based on the survival regression model and the coexpression network are combined. The cluster center automatic search clustering technology based on density descending is introduced to construct the synergy-coexpression network.

[0053] (2) Based on the assumption that significant biomarkers are cross-cluster, a bottom-up biomarker enumeration method and a top-down biomarker random selection method are designed to achieve joint multivariate unsupervised biomarker selection;

[0054] The present invention uses mature theories such as Cox regression, Pearson correlation, unsupervised clustering, supervised classification, differential expression analysis, Bayesian network, etc. to provide strong theoretical support for the implementation of the present invention. The technical route is reasonably laid out, the logical relationship between each part is close, and it is advanced layer by layer. First, a metric based on Cox regression and Pearson correlation is constructed to measure the weight of the edge between nodes, and a clustering algorithm is used to degenerate the synergistic-coexpression network from a weighted undirected graph to an undirected graph; then, for the bottom-up or top-down feature selection results, a clustering algorithm is introduced to realize unsupervised feature selection of joint multivariates, and a prognostic survival analysis indicator is constructed; then, the survival analysis indicators of the constructed prognostic features are clustered to complete the prognostic category discovery and class comparison; finally, a Bayesian network is used to realize the causal relationship inference of the prognostic features, and the survival risk category to which the new sample belongs is predicted. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0056] Figure 1 This is a framework diagram for the research content of cancer multi-omics prognostic survival analysis based on synergy-co-expression network;

[0057] Figure 2 It is a technical route diagram of the method of the present invention;

[0058] Figure 3 Schematic diagram for constructing synergy-coexpression networks for cancer multi-omics data;

[0059] Figure 4 Schematic diagram for prognostic biomarker extraction;

[0060] Figure 5 Schematic diagram for prognostic category discovery in cancer multi-omics data;

[0061] Figure 6 The result of hierarchical clustering on all pairwise significant miRNA expression levels;

[0062] Figure 7 The results of hierarchical clustering on all pairwise significant miRNA correlation coefficients. DETAILED DESCRIPTION

[0063] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0064] In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc., indicating the orientation or positional relationship, are based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention. In addition, the terms "first", "second", and "third" are used for descriptive purposes only, and cannot be understood as indicating or implying relative importance.

[0065] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "installed", "connected", and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0066] In addition, the technical features involved in the different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0067] The present invention is described in detail below in conjunction with specific embodiments. Specific embodiment one:

[0069] according to Figure 1-Figure 7 As shown, the specific optimization technical solution adopted by the present invention to solve the above technical problems is: the present invention relates to a cancer prognosis survival analysis method based on a collaborative co-expression network.

[0070] A cancer prognosis survival analysis method based on a collaborative co-expression network comprises the following steps:

[0071] Step 1: Collect multi-omics data of cancer based on public databases and pre-process the data;

[0072] Step 2: Based on the processed data, construct a collaborative co-expression network for cancer multi-omics data;

[0073] Step 3: Extract prognostic biomarkers based on the synergistic co-expression network;

[0074] Step 4: Based on the extracted prognostic biomarkers, sample clustering for survival analysis is performed to conduct prognostic survival analysis. Specific embodiment 2:

[0076] The difference between the second embodiment of the present invention and the first embodiment is that:

[0077] The step 1 is specifically as follows:

[0078] Step 1.1: For the target tumor, collect high-quality RNA-seq data of FPKM expression type from the TCGA public database, convert the expression value to log2 TPM+1, and perform normalization to ensure that the overall expression of the data conforms to the normal distribution;

[0079] Step 1.2: Collect DNA methylation data of target tumor samples in the TCGA database, give priority to high-throughput methylation sequencing data, and use Minfi and ChAMPR packages for background correction, normalization, probe filtering, and batch effect correction;

[0080] Step 1.3: For the target tumor, collect high-quality CNV data of Segment Mean type from the TCGA public database through Genomic Data Commons and convert the data into GISTIC2 format; perform data cleaning methods including removing outliers, correcting errors, and removing redundant information on the CNV data, and then standardize the CNV data to eliminate the impact of experimental batch effects and sample processing differences, and use Z-score standardization to standardize the CNV data;

[0081] Step 1.3: Collect ChIP-Seq FASTQ format data of target tumor samples in the TCGA database, use FastQC to perform quality control on the raw sequencing data to identify and deal with low-quality reads, contamination or other potential problems, and use Trimmomatic or Cutadapt tools to remove adapter sequences and low-quality reads; use alignment tools to align the cleaned reads to the reference genome, use SAMtools to filter out reads with low alignment quality, multiple-point alignments, and non-specific alignments, and use Picard Tools or SAMtools to remove PCR duplications to ensure data accuracy. Specific embodiment three:

[0083] The difference between the third embodiment of the present invention and the second embodiment is that:

[0084] The step 2 is specifically as follows:

[0085] Cox regression based on the joint enumeration of two features is used to achieve the fitting of the prognosis survival time, and Pearson correlation based on the joint enumeration of two features is used to measure the correlation between features. The above two metrics are unified to mark the weights of any two node edges in the fully connected network; the cluster center automatic search clustering algorithm based on density descending is used to delete the edges with lower weights in the network, so as to realize network construction and feature cluster mining. The cluster center automatic search algorithm based on density descending is used to cluster the weights of all edges in the network in density descending order, and delete the edges with lower weight density from the network until a disconnected feature cluster appears;

[0086] Step 2.1: Based on the Cox regression, the Cox proportional hazard regression model is used to fit the survival time by enumerating the pairwise features, so as to reveal the significant relationship between any pair of features and prognosis survival. By performing the maximum likelihood estimation on the logarithmic partial likelihood function of the hazard function, the regression coefficient β corresponding to any pair of features is obtained. i ,β j ) T , when it is assumed that the features are uncorrelated, any component of the regression coefficient obeys the Gaussian distribution, and the Wald test is used to determine whether the regression coefficient component deviates significantly from zero;

[0087] The sample size can be further expanded by reordering the samples and their corresponding survival results. The reordered p-value corresponding to the Wald statistic can be generated on each feature. All pairwise features are enumerated to obtain paired p-values, which are used to select significant feature pairs closely related to survival from all pairwise feature combinations. In order to deal with overfitting, a penalty constraint is added to the logarithmic partial likelihood function corresponding to the risk function. Considering that the hazard ratio may change with survival time, a weighted Cox regression model is used instead.

[0088] Step 2.2: Pearson correlation-based joint enumeration of two features. To measure the correlation between two features, the Pearson correlation coefficient r is used to calculate the correlation between any two features, and the t statistic of n samples is constructed based on this. The statistic is expressed as: By using the sample reordering method, the order of sample values ​​of one feature in any pair of features is disrupted, and the above t statistic is calculated repeatedly. The corresponding reordered p value is expressed as:

[0089]

[0090] Among them, t(i,j) is the statistic corresponding to the correlation coefficient r of features i and j, t b (i,j) represents the statistic generated by randomly reordering the values ​​of any feature, and B is the number of random reorderings;

[0091] Step 2.3: Determine the association relationship of network nodes. The collaborative co-expression network uses a single feature as a network node and the relationship between two features as an edge. The weight of the edge between nodes is calculated by comprehensively considering the significant relationship between the feature and the survival time and the correlation between the features. The following metric is constructed to measure the weight of the edge between nodes:

[0092]

[0093] Where p2(i) and p2(j) represent the binary Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively, and p1(i) and p1(j) represent the univariate Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively. represents rounding up; the first term of m(i,j) is the "synergistic" part that measures the contribution of features i and j to survival regression, and the second term is derived from the reordered p-value obtained by the joint enumeration of two features based on Pearson correlation, corresponding to the "co-expression" part; m(i,j) is used to measure the association between nodes i and j: the larger its value, the more likely it is that nodes i and j are derived from the same feature cluster, and vice versa. When the value of m(i,j) is less than zero, its value is set to the maximum value of the weight in the network, and the weighted undirected graph constructed is called a synergistic-co-expression network;

[0094] Step 2.4: Based on the feature cluster mining of the collaborative-co-expression network, for the weight of the collaborative-co-expression network edge, which is a measure of the node association relationship, an automatic cluster center search clustering algorithm based on density descending is introduced. The measurement threshold of the edge in the network is limited by the optimal automatic estimation of the number of categories and cluster boundaries to determine whether there is an associated edge between any two nodes in the network, thereby degenerating the network from a weighted undirected graph to an undirected graph. Specific embodiment four:

[0096] The difference between the fourth embodiment of the present invention and the third embodiment is that:

[0097] The step 3 is specifically as follows:

[0098] Step 3.1: Bottom-up “cross-cluster enumeration” biomarker selection. According to the “cross-cluster hypothesis”, enumerate the univariate within each feature cluster, combine the univariate from c different feature clusters into a multivariate, and calculate its regression coefficient on the Cox survival hazard ratio function;

[0099] Step 3.2: Top-down "cross-cluster random" biomarker selection, introducing a sample resampling voting strategy, randomly selecting a single variable within each feature cluster, combining the single variables from c different feature clusters into multivariates, and calculating their regression coefficients for the Cox survival hazard ratio function; assuming that the features are uncorrelated, analogous to the pairwise feature joint enumeration method based on Cox regression, each group of "cross-cluster" multivariates can obtain corresponding c re-ranked p-values, and each round of multivariate voting is performed on the multivariates with significant re-ranked p-values; after l rounds of voting, the features with the highest votes are selected from each feature cluster to form a significant prognostic feature set;

[0100] Step 3.3: Unsupervised feature selection of joint multivariates; single variables need to be extracted from different feature clusters to form joint multivariates for regression analysis. For bottom-up "cross-cluster enumeration" biomarker selection, c re-ranked p-values ​​corresponding to each enumerated multivariate are recorded; for top-down "cross-cluster random" biomarker selection, multiple groups of c-dimensional p-value vectors corresponding to multiple randomly extracted multivariates can be recorded by setting a threshold. In the c-dimensional feature space, all enumerated or re-ranked p-value vectors that meet the threshold limit are automatically searched and clustered based on the descending density of the cluster center, and the cluster closest to the origin is selected to form a c-dimensional feature set. The Mahalanobis distance from the scattered points to the origin is calculated. When the bottom-up "cross-cluster enumeration" biomarker selection is adopted, the multidimensional variable with the closest distance to the origin is selected as the significant prognostic feature set; if the top-down "cross-cluster random" biomarker selection is adopted, it is still necessary to determine whether the multidimensional variable with the closest distance to the origin corresponds to the feature combination with the highest cumulative votes;

[0101] Step 3.4: Construction of prognostic survival analysis indicators. Taking the Cox proportional hazard regression model as an example, since it is difficult to estimate the baseline risk function or baseline survival function, the linear part of the risk score is used to construct the prognostic survival indicator. The risk score of the sample is a linear combination of the corresponding sample values ​​of the prognostic characteristics, and the combination coefficient is the regression coefficient. Specific embodiment five:

[0103] The difference between the fifth embodiment of the present invention and the fourth embodiment is that:

[0104] The step 4 is specifically as follows:

[0105] Step 4.1: Sample clustering for prognostic survival analysis: To discover the prognostic survival category, the key biomarkers from different omics or the prognostic survival analysis indicators built by different omics data are combined as the input of the convolutional neural network, and its output is regarded as a new feature subspace. The samples are automatically searched and clustered based on the cluster center in descending density, so as to discover the prognostic category to which the samples belong;

[0106] Step 4.2: Draw Kaplan-Meier survival analysis curves for samples of different categories to observe the changes in survival time within the category and the differences in survival time between categories. Use the log-rank test to measure whether there are significant quantitative differences in survival time between different categories. Draw a risk score survival analysis graph to observe the changes in survival time between samples of different categories.

[0107] Step 4.3: Analysis of significant differences in prognosis between classes. Class comparison is achieved using statistics based on sample distribution differences to evaluate whether there are significant differences in risk scores for different categories. On this basis, attempts are made to construct different types of classifiers using training samples, and the accuracy and effectiveness of the classification results are verified using the prognosis survival results of the test samples. Specific embodiment six:

[0109] The difference between the sixth embodiment of the present invention and the fifth embodiment is that:

[0110] Synergy represents the correlation between the feature pair composed of two biomarkers and prognostic survival time; co-expression represents the correlation between two biomarkers. Specific embodiment seven:

[0112] The difference between the seventh embodiment of the present invention and the sixth embodiment is that:

[0113] The cross-group assumptions are as follows:

[0114] The significant feature set closely related to the prognosis survival time is not entirely composed of a single significant variable related to the prognosis survival time, and there is no obvious correlation between the variables in the significant feature set. Specific embodiment eight:

[0116] The difference between the eighth embodiment of the present invention and the seventh embodiment is that:

[0117] The present invention provides a cancer prognosis survival analysis system based on a collaborative co-expression network, the system comprising:

[0118] A data acquisition module, which collects multi-omics data of cancer based on a public database and pre-processes the data;

[0119] A network construction module, wherein the network construction module constructs a collaborative co-expression network for cancer multi-omics data based on the processed data;

[0120] An extraction module, wherein the extraction module extracts prognostic biomarkers based on a collaborative co-expression network;

[0121] A clustering module is used to cluster samples for survival analysis based on extracted prognostic biomarkers and perform prognostic survival analysis. Specific embodiment nine:

[0123] The difference between the ninth embodiment of the present invention and the eighth embodiment is that:

[0124] The present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the program is executed by a processor to implement a cancer prognosis survival analysis method based on a collaborative co-expression network. Specific embodiment ten:

[0126] The difference between the tenth embodiment of the present invention and the ninth embodiment is that:

[0127] The present invention provides a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and the processor implements a cancer prognosis survival analysis method based on a collaborative co-expression network when executing the computer program. Specific embodiment eleven:

[0129] The difference between the eleventh embodiment of the present invention and the tenth embodiment is that:

[0130] The present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0131] The previous research of this invention found that the significant feature set closely related to the prognosis survival time is not composed entirely of a single significant variable related to the prognosis survival time; that is, there is no obvious correlation between the variables in the significant feature set. (In fact, this hypothesis comes from a series of research works such as the covariate selection method related to tumor prognosis survival proposed previously). This phenomenon is called . This shows that the features in the significant feature set come from different feature clusters in the co-expression network model, rather than the same feature cluster where a single significant feature is located. This is completely contrary to the practice of selecting the cluster where a single significant variable is located from the co-expression network model to achieve feature selection. Take the two features that constitute a significant feature pair (i.e., a bigram) as an example: one and only one of the features is closely related to the prognosis survival time; when the other feature is added, the fitting performance of the feature bigram for the prognosis survival time is further enhanced. In fact, the existing co-expression network model is a tree structure, and focusing on the correlation between a pair of features and the prognosis survival time is precisely the key to graphing (i.e., networking) the feature set.

[0132] Based on the above considerations, the cross-cluster hypothesis is regarded as a collaborative expression of the binary group of prognostic significant features, and a prognostic survival analysis method based on the collaborative-coexpression network is proposed to achieve the integration of the five types of strategies in existing research (such as Figure 1 ). Among them, The correlation between the feature pair consisting of two biomarkers and the prognostic survival time; It represents the correlation between two biomarkers.

[0133] For RNA-seq data, Chip-seq data, DNA methylation data, copy number variation data and other multi-omics data, we will build synergy-coexpression networks and select prognostic biomarkers from them; combine prognostic biomarkers from different omics to generate prognostic survival analysis indicators, and build Bayesian networks to achieve causal relationship inference of different prognostic biomarkers. Specifically, we plan to conduct research in the following three aspects.

[0134] (1) Research on the construction method of synergistic-co-expression network for cancer multi-omics data

[0135] The two-dimensional matrices corresponding to the multi-omics data are used as input to construct the corresponding synergy-coexpression network. Figure 1As shown in the red box in , it is divided into two parts: "synergy" and "co-expression". On the one hand, focus on the connection between biomarkers and survival time: adopt a bottom-up feature combination strategy to design the association measurement between pairwise biomarkers and survival time as the "synergy" part of the constructed network; different from the univariate feature combination strategy, the "synergy" part of the constructed network only focuses on the joint enumeration of pairwise features. On the other hand, explore the correlation between biomarkers: adopt a related feature clustering strategy to traverse the correlation between pairwise biomarkers as the "co-expression" part of the constructed network; different from the related feature clustering strategy, the "co-expression" part of the constructed network focuses on the correlation between pairwise features, rather than the correlation clustering of features. The synergy-co-expression network regards features as network nodes and the association between pairwise features as edges. The key to constructing the above network is to determine whether there is an associated edge between any two nodes. The present invention attempts to unify the "synergy" and "co-expression" information of pairwise features into the weight of the edge connecting any two nodes, and uses a cluster center automatic search clustering algorithm based on density descending to screen out edges with smaller weights, so as to obtain a synergy-co-expression network for cancer multi-omics data.

[0136] (2) Research on the selection method of prognostic biomarkers by combining multiple variables

[0137] Multivariate prognostic biomarker selection Figure 1 As shown in the green box in , it is divided into two parts: "cross-cluster enumeration" and "cross-cluster random". The above two prognostic biomarker selection methods are both based on the cross-cluster hypothesis, that is, the significant features related to survival time come from different feature clusters (or feature clusters, for the convenience of network representation, they are all referred to as feature clusters below). For the constructed synergistic-coexpression network, the present invention will try two feature selection methods, namely "cross-cluster enumeration" from bottom to top and "cross-cluster random" from top to bottom, based on the idea of ​​"cross-cluster", to obtain a set of key biomarkers closely related to survival time. It should be further explained that when constructing a synergistic-coexpression network, only two features need to be concerned about; the "cross-cluster enumeration" and "cross-cluster random" here are both carried out between different feature clusters in the network. If the constructed synergistic-coexpression network has c feature clusters, the dimension of each round of enumeration or random selection of features is c.

[0138] (3) Research on prognostic category discovery methods for cancer multi-omics data

[0139] Prognostic survival categories were found to be Figure 1As shown in the blue box in , for the prognostic biomarker sets selected from different omics data, the present invention intends to carry out the following research in sequence: construct a prognostic survival analysis indicator by combining prognostic biomarkers from different omics; use a bottom-up sample clustering strategy on the constructed prognostic survival analysis indicator to discover the prognostic survival category of the sample; in order to verify the reliability of the selected biomarkers and sample categories, use a bottom-up sample classification strategy to determine the multivariate statistical hypothesis test results and classification results of the sample prognostic category; and realize the prediction of the category to which the new sample belongs.

[0140] Research objectives

[0141] In cancer multi-omics analysis, by taking into account the correlation between features and between features and prognostic survival time, selecting key biomarkers closely related to prognostic survival, and discovering possible prognostic survival subtypes of specific cancers, have become two urgent problems to be solved. The present invention will focus on cancer multi-omics data with survival time and survival event occurrence status (such as death or not), and construct a synergistic-coexpression network based on the cross-cluster hypothesis; it is divided into "cross-cluster enumeration" and "cross-cluster random" feature selection methods to achieve joint multivariate prognostic biomarker selection; for the submatrix corresponding to the selected prognostic biomarker and its prognostic survival analysis index, an unsupervised clustering method is used to find different risk categories of sample prognostic survival; the present invention realizes accurate prediction of the prognostic category to which cancer patients belong. In this way, more accurate statistical models and computational tools are provided for biologists to study prognostic survival analysis using cancer multi-omics data.

[0142] Key issues to be addressed

[0143] In order to achieve the above research objectives, the present invention intends to focus on solving the following three problems.

[0144] (1) Determination of the association relationship between nodes in the synergistic-coexpression network

[0145] The synergy-coexpression network is a small-world network that combines pairwise features. The "cross-cluster hypothesis" states that a pair of features that are significantly correlated with survival time originate from different feature clusters, while related pairwise features come from the same feature cluster. How to unify the measures of "synergy" and "coexpression" on the associated edges of the synergy-coexpression network, calculate their weights and automatically remove the edges with lower weights to obtain feature clusters that are not connected to each other is a key problem that needs to be solved urgently.

[0146] (2) Unsupervised feature selection of multiple variables

[0147] After establishing the synergistic-coexpression network, it is necessary to select a set of prognostic biomarkers from the network for a series of subsequent applications (such as designing prognostic survival analysis indicators, etc.). Based on the "cross-cluster hypothesis", two feature selection methods, bottom-up "cross-cluster enumeration" and top-down "cross-cluster randomization", can be used. Regardless of the method, key features need to be selected from different feature clusters to form a joint multivariate for Cox regression analysis. How to achieve the optimal selection of prognostic features in an unsupervised manner is a key issue that needs to be solved urgently.

[0148] The present invention insists on using the expression method of "unsupervised feature selection", and its viewpoint is as follows: in the process of constructing the synergistic-coexpression network, survival information is indeed used when calculating the weights of the edges between two nodes; but in the step of realizing the joint multivariate feature selection, only the proposed unsupervised clustering algorithm is used to perform unsupervised clustering on the feature scores or p-values ​​enumerated across clusters or randomly selected across clusters. The entire clustering process does not involve the prognostic survival time or survival status at all. For a detailed description, please see the "Proposed Research Plan" section below.

[0149] (3) Prognostic category discovery for cancer multi-omics data

[0150] In the past decade, there have been various debates on whether to conduct univariate sample classification or univariate sample clustering. The focus of the conflict is: how to divide samples into different survival categories. For key biomarkers selected from different omics data that are closely related to prognostic survival time, how to construct a cross-omics prognostic survival model and connect the sample clustering strategy with the sample classification strategy is a key issue that needs to be solved urgently.

[0151] Study plan

[0152] In view of the above research content and research objectives, this paper takes the cancer multi-omics public data (TCGA) with prognosis and survival information as the research object, and will cover the collection and analysis of transcriptomics (RNA-seq), genomics (DNA methylation, CNV), and proteomics (Chip-seq) data. For the multi-omics data of a certain cancer, the specific data preprocessing scheme is as follows:

[0153] ① For the target tumor, we plan to collect high-quality RNA-seq data of FPKM expression type from the TCGA public database, and convert the expression value into log2(TPM+1) for normalization, aiming to ensure that the overall expression of the data conforms to the normal distribution and ensure the accuracy of the experiment.

[0154] ②Collect DNA methylation data of target tumor samples in the TCGA database, and give priority to high-throughput methylation sequencing (such as Illumina 450K or EPIC BeadChip) data. Use R packages such as Minfi and ChAMP for background correction, normalization, probe filtering, and batch effect correction.

[0155] ③ For target tumors, high-quality CNV data of SegmentMean type from the TCGA public database were collected through Genomic Data Commons (GDC). For further analysis, the data was converted to GISTIC2 format; to ensure data quality, the CNV data was first cleaned by methods including removing outliers, correcting errors, and removing redundant information, and then the CNV data was standardized to eliminate the effects of experimental batch effects and sample processing differences. Z-score standardization is proposed to standardize CNV data.

[0156] ④ FASTQ format data of chromatin immunoprecipitation sequencing (ChIP-Seq) of target tumor samples were collected in the TCGA database. FastQC was used to perform quality control on the raw sequencing data to identify and deal with low-quality reads, contamination or other potential problems. Tools such as Trimmomatic or Cutadapt were used to remove adapter sequences and low-quality reads. Alignment tools (such as BWA, Bowtie2 or STAR) were used to align the cleaned reads to the reference genome. SAMtools was used to filter out reads with low alignment quality, multi-point alignment, and non-specific alignment. Further experiments were conducted using tools such as Picard Tools or SAMtools to remove PCR duplications to ensure data accuracy.

[0157] Based on the above data, the present invention intends to carry out prognostic survival analysis research in four parts. The overall technical roadmap is as follows: Figure 2 As shown in the figure. The solid arrows indicate the order of the research plan. For the multi-omics data of a certain cancer, the processing plan is as follows: after preprocessing the data of different omics, the synergistic-coexpression network construction and the extraction of prognostic significant biomarkers are realized respectively; the prognostic significant biomarkers from different omics are combined as the input of the convolutional neural network, and its output is used for prognostic category discovery; on this basis, the causal relationship inference of the prognostic biomarker set from different omics is performed.

[0158] (1) Construction of collaborative-co-expression networks for cancer multi-omics data

[0159] This part serves as the basis for the entire project. The specific research work includes: how to combine pairwise features to achieve the fitting of prognosis survival; how to calculate the correlation between pairwise features; on this basis, consider how to unify the fitting results and correlation metrics into the weight metrics of the edges between network nodes and delete the key issue of edges with lower weights.

[0160] like Figure 3 As shown in the figure, for a certain group data of cancer, Cox regression based on the joint enumeration of two features is proposed to achieve the fitting of the prognosis survival time. Pearson correlation based on the joint enumeration of two features is used to measure the correlation between features. The above two metrics are unified to mark the weights of any two node edges in the fully connected network; the proposed cluster center automatic search clustering algorithm based on density descending is used to delete the edges with lower weights in the network, so as to achieve network construction and feature cluster mining. Here, the proposed cluster center automatic search algorithm based on density descending is used to cluster the weights of all edges in the network in density descending order, and delete the edges with lower weight density (corresponding to several Gaussian distributions with lower density where the weight is located) from the network until a disconnected feature cluster appears. For the specific content of the proposed clustering algorithm, please refer to the "Research Foundation" section below and its corresponding references.

[0161] Pairwise feature joint enumeration based on Cox regression. Research on survival analysis involves a variety of models. The Cox proportional hazard regression model has been integrated into various application statistical software because it is easy to understand and implement, and has been widely used. The present invention adopts the Cox proportional hazard regression model to fit the survival time by enumerating pairwise features, thereby revealing the significant relationship between any pair of features and prognosis survival. By performing maximum likelihood estimation on the logarithmic partial likelihood function of the risk function, the regression coefficient β corresponding to any pair of features can be obtained. i ,β j ) T . If the features are assumed to be uncorrelated, any component of the regression coefficient follows a Gaussian distribution. The Wald test can be used to determine whether the regression coefficient component deviates significantly from zero. The sample size can be further expanded by reordering the samples and their corresponding survival results, and a reordered p-value corresponding to the Wald statistic can be generated for each feature. Enumerating all pairwise features to obtain paired p-values ​​can be used to select significant feature pairs that are closely related to survival from all pairwise feature combinations. To deal with overfitting, it is possible to consider adding a penalty constraint to the logarithmic partial likelihood function corresponding to the risk function; considering that the hazard ratio may change with survival time, a weighted Cox regression model can be used instead.

[0162] Pairwise feature joint enumeration based on Pearson correlation. To measure the correlation between two features, the Pearson correlation coefficient r is used to calculate the correlation between any two features, and the t statistic of n samples is constructed based on this. The statistic can be expressed as: We still use the sample reordering method to disrupt the order of sample values ​​of any two features and repeat the calculation of the above t statistic. The corresponding reordered p value can be expressed as:

[0163]

[0164] Among them, t(i,j) is the statistic corresponding to the correlation coefficient r of features i and j, t b (i, j) represents the statistic generated by randomly reordering the values ​​of any feature, and B is the number of random reorderings.

[0165] Determination of the association relationship of network nodes. The collaborative-coexpression network uses a single feature as a network node and the relationship between two features as an edge. The weight of the edge between nodes is calculated by comprehensively analyzing the significant relationship between the feature and the survival time and the correlation between the features. On the one hand, the Cox regression model based on the joint enumeration of two features identifies the association relationship between paired network nodes and the survival time. The smaller the reranking p-value corresponding to any pair of feature nodes, the greater the contribution of these two nodes to the regression of the survival risk proportional function. On the other hand, the Pearson correlation coefficient based on the joint enumeration of two features identifies the correlation between two nodes. The smaller the reranking p-value corresponding to any pair of feature nodes, the more correlated the two nodes are. According to the "cross-cluster hypothesis", two features that are significantly correlated with survival time come from different feature clusters; and the more correlated the two features are, the more likely they are to come from the same feature cluster. Therefore, the present invention intends to construct the following metric for measuring the weights of the edges between nodes:

[0166]

[0167] Where p2(i) and p2(j) represent the binary Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively, and p1(i) and p1(j) represent the univariate Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively. Represents rounding up. The first term of m(i,j) is the "synergistic" part that measures the contribution of features i and j to survival regression, and the second term is derived from the reordered p-value obtained by the joint enumeration of two features based on Pearson correlation, which corresponds to the "co-expression" part. m(i,j) is used to measure the association between nodes i and j: the larger its value, the more likely it is that nodes i and j are derived from the same feature cluster, and vice versa. In addition, when the value of m(i,j) is less than zero, its value is set to the maximum value of the weight in the network. The weighted undirected graph constructed in this way is called .

[0168] Feature cluster mining based on synergy-coexpression network. Aiming at the node association relationship measurement of the edge weight of the synergy-coexpression network, the proposed cluster center automatic search clustering algorithm based on density descending is introduced. The measurement threshold of the edge in the network is limited by the optimal automatic estimation of the number of categories and cluster boundaries to determine whether there is an associated edge between any two nodes in the network, so as to degenerate the network from a weighted undirected graph to an undirected graph. Among them, the unconnected nodes are regarded as originating from different feature clusters, and the synergy-coexpression network is degenerated into several non-overlapping feature clusters.

[0169] (2) Extraction of prognostic biomarkers

[0170] This part is based on the "cross-cluster hypothesis" and aims to extract a set of significant biomarkers (i.e., prognostic feature set) that are closely related to prognosis and survival. Specific research work includes: how to enumerate all cross-cluster feature combinations from the bottom up to select a prognostic biomarker set; how to randomly extract cross-cluster feature combinations from the top down to select a prognostic biomarker set; for any of the above feature selection schemes, consider the key issue of how to automatically select a prognostic feature set from a large number of multidimensional features, that is, the problem of unsupervised feature selection of joint multiple variables; how to construct clinical indicators for survival analysis based on the selected prognostic feature set.

[0171] like Figure 4 As shown in the figure, the obtained synergistic-coexpression network can clearly mark different feature clusters. According to the "cross-cluster hypothesis", bottom-up "cross-cluster" feature enumeration and top-down random "cross-cluster" feature selection are respectively adopted on a certain omics matrix. Ignoring the correlation of significant features, an unsupervised feature selection scheme for measuring multivariate can be designed based on the automatic search clustering algorithm of cluster center in descending density. On this basis, the selected prognostic features are used to calculate the risk score based on the Cox proportional hazard model to realize the construction of prognostic survival analysis indicators.

[0172] Bottom-up "cross-cluster enumeration" biomarker selection. According to the "cross-cluster hypothesis", single variables are enumerated within each feature cluster, and the single variables from c different feature clusters are combined into multivariate variables to calculate their regression coefficients for the Cox survival hazard ratio function. This method is essentially a multivariate joint enumeration based on Cox regression. Assuming that the features are uncorrelated, analogous to the pairwise feature joint enumeration method based on Cox regression, each group of "cross-cluster" multivariate can obtain corresponding c re-ranked p-values, which can be used to select a set of significant prognostic feature sets that are most closely related to survival from all c-dimensional multivariate variables. Similarly, it is possible to consider adding penalty constraints to the log-likelihood function corresponding to the risk function or using a weighted Cox regression model to improve the existing model.

[0173] Top-down "cross-cluster random" biomarker selection. If the number of feature clusters is too large, consider using a top-down "cross-cluster random" feature selection method. Introduce a sample resampling voting strategy, randomly select a single variable within each feature cluster, combine the single variables from c different feature clusters into multivariates, and calculate their regression coefficients for the Cox survival hazard ratio function; assuming that the features are uncorrelated, analogous to the pairwise feature joint enumeration method based on Cox regression, each group of "cross-cluster" multivariates can obtain corresponding c re-ranked p values, and each round counts the multivariates with significant re-ranked p values; after l rounds of voting, select the features with the highest votes from each feature cluster to form a significant prognostic feature set.

[0174] Unsupervised feature selection of joint multivariates. Regardless of the above methods, it is necessary to extract single variables from different feature clusters to form joint multivariates for regression analysis. For the bottom-up "cross-cluster enumeration" biomarker selection, c reordered p-values ​​corresponding to each enumerated multivariate are recorded; for the top-down "cross-cluster random" biomarker selection, multiple groups of c-dimensional p-value vectors corresponding to multiple randomly extracted multivariates can be recorded by setting a threshold. In the c-dimensional feature space, all enumerated or reordered p-value vectors that meet the threshold limit are automatically clustered based on density descending cluster center search, and the cluster closest to the origin is selected to form a c-dimensional feature set, and the Mahalanobis distance from the scattered points to the origin is calculated. If the bottom-up "cross-cluster enumeration" biomarker selection is used, the multidimensional variable with the closest distance to the origin is selected as the significant prognostic feature set; if the top-down "cross-cluster random" biomarker selection is used, it is still necessary to determine whether the multidimensional variable with the closest distance to the origin corresponds to the feature combination with the highest cumulative votes.

[0175] Construction of prognostic survival analysis indicators. Taking the Cox proportional hazard regression model as an example, since the baseline risk function or baseline survival function is difficult to estimate, the risk score (risk score) - the linear part of this risk function is often used to construct the prognostic survival indicator. The risk score of the sample is a linear combination of the sample values ​​corresponding to the prognostic feature, and the combination coefficient is the regression coefficient. If other models are used, it is possible to consider directly using the submatrix element value corresponding to the prognostic feature set as the indicator of sample prognostic survival analysis.

[0176] (3) Prognostic category discovery for cancer multi-omics data

[0177] This part aims to discover different categories of prognostic samples. Specific research work includes: how to cluster samples based on the established prognostic survival analysis indicators; how to verify the reliability of prognostic categories based on sample distribution differences. It is planned to use a clustering algorithm based on density descending cluster center automatic search clustering to cluster samples of prognostic survival analysis indicators, and to use hypothesis testing technology to analyze the degree of expression differences between categories. The schematic diagram of prognostic category discovery is shown in the figure below. Figure 5 shown.

[0178] Sample clustering for prognostic survival analysis. In order to discover the category of prognostic survival, the key biomarkers from different omics or the prognostic survival analysis indicators built by different omics data are combined as the input of the convolutional neural network, and its output is regarded as a new feature subspace. The samples are automatically searched and clustered based on the cluster center in descending density, so as to realize the discovery of the prognostic category to which the samples belong. Kaplan-Meier survival analysis curves are drawn between samples of different categories to observe the changes in survival time within the class and the differences in survival time between classes, and the log-rank test is used to measure whether there are significant differences in survival time between different categories. Furthermore, the risk score survival analysis graph is drawn to observe the changes in survival time between samples of different categories.

[0179] Analysis of significant differences in prognosis between categories. To verify the reliability of the prognostic categories found, a series of statistical methods based on sample distribution differences (such as Welch two-sample t test, one-way analysis of variance, Hotelling two-sample t test, and 2 Tests, etc.) are used to realize class comparison and evaluate whether there are significant differences in risk scores of different categories. On this basis, attempts are made to construct different types of classifiers using training samples, and the accuracy and effectiveness of the classification results are verified using the prognosis survival results of test samples.

[0180] The present invention discovered that the two microRNAs significantly correlated with survival time come from different characteristic clusters. This can be used as the most direct experimental verification of the "cross-cluster hypothesis". In addition, the applicant has independently developed 5 practical original differential expression and prognostic survival analysis tools (including: JSDSA, JCD-DEA, ECFS-DEA, AFS-DEA and IOFS-SA). The specific content of the research work that has been carried out and is being carried out is detailed in the research basis section.

[0181] Research foundation;

[0182] Since December 2012, the applicant has been engaged in the research of differential expression analysis and prognosis survival analysis of cancer multi-omics data. In recent years, the applicant has done a lot of preliminary research on this topic, mainly including: (1) Research on covariate selection methods related to cancer prognosis survival; (2) Research on clustering methods based on automatic search of cluster centers based on density descending; (3) Research on re-clustering methods for prognosis survival category discovery; (4) Development of related differential expression analysis and prognosis survival analysis tools. The above research has laid the foundation for the development of this project.

[0183] Study on covariate methods related to cancer prognosis and survival

[0184] This research work discovered a pair of microRNAs that lead to significant differences in the prognosis of malignant gliomas through bottom-up multivariate joint enumeration, and further experimental results verified the "cross-cluster hypothesis". The applicant, together with the relevant medical team, conducted statistical experiments and pathway verification on the microRNA expression spectrum of the TCGA public data set, designed a roadmap based on covariate selection, and discovered a pair of key microRNAs that lead to differences in the prognosis and survival of malignant gliomas. While verifying the idea of ​​bottom-up multivariate joint enumeration, this research work discovered that the two microRNAs that are significantly correlated with survival time may come from different gene clusters. Figure 6 It was revealed that the selected two significant microRNAs (miR-222 and miR-10b) came from different gene clusters. Figure 7 The solid black boxes in the figure represent the gene clusters, and the dashed black boxes represent the pairwise significant genes. This indicates that there is a significant positive correlation between genes within the gene clusters, and pairwise significant genes cross clusters. The above results can be used as the most direct experimental verification of the "cross-cluster hypothesis".

[0185] The applicant has also independently developed five practical original differential expression and prognostic survival analysis tools, including: JSDSA, JCD-DEA, ECFS-DEA, AFS-DEA and IOFS-SA. JSDSA is a bottom-up prognostic survival analysis tool based on Cox regression. The development of this tool provides firm technical support for the realization of the joint enumeration of two features based on Cox regression and the bottom-up "cross-group enumeration" feature selection. JCD-DEA is a bottom-up differential expression analysis tool that combines sample distribution and classification differences. ECFS-DEA is a top-down differential expression analysis tool based on integrated classification. AFS-DEA is its network version. The development of the above tools can provide strong technical support for the analysis of significant prognostic differences between classes. IOFS-SA is a special tool for the construction of prognostic survival analysis indicators and the discovery of prognostic categories. It realizes interactive re-clustering of samples for the constructed prognostic survival analysis indicators.

[0186] The above is only a preferred implementation of a cancer prognosis survival analysis method based on a collaborative co-expression network. The protection scope of a cancer prognosis survival analysis method based on a collaborative co-expression network is not limited to the above embodiments. All technical solutions under this idea belong to the protection scope of the present invention. It should be pointed out that for those skilled in the art, several improvements and changes without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A cancer prognosis survival analysis method based on a collaborative co-expression network, characterized by: The following steps are involved: Step 1: Collect multi-omics data of cancer based on public databases and pre-process the data; Step 2: Based on the processed data, construct a collaborative co-expression network for cancer multi-omics data; Step 3: Extract prognostic biomarkers based on the synergistic co-expression network; Step 4: Based on the extracted prognostic biomarkers, sample clustering for survival analysis is performed, and prognostic survival analysis is performed.

2. The method according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1: For the target tumor, collect high-quality RNA-seq data of FPKM expression type from the TCGA public database, convert the expression value to log2 TPM+1, and perform normalization to ensure that the overall expression of the data conforms to the normal distribution; Step 1.2: Collect DNA methylation data of target tumor samples in the TCGA database, give priority to high-throughput methylation sequencing data, and use Minfi and ChAMPR packages for background correction, normalization, probe filtering, and batch effect correction; Step 1.3: For the target tumor, collect high-quality CNV data of SegmentMean type from the TCGA public database through Genomic Data Commons and convert the data into GISTIC2 format; perform data cleaning methods including removing outliers, correcting errors, and removing redundant information on the CNV data, and then standardize the CNV data to eliminate the impact of experimental batch effects and sample processing differences, and use Z-score standardization to standardize the CNV data; Step 1.3: Collect ChIP-Seq FASTQ format data of target tumor samples in the TCGA database, use FastQC to perform quality control on the raw sequencing data to identify and deal with low-quality reads, contamination or other potential problems, and use Trimmomatic or Cutadapt tools to remove adapter sequences and low-quality reads; use alignment tools to align the cleaned reads to the reference genome, use SAMtools to filter out reads with low alignment quality, multiple-point alignments, and non-specific alignments, and use Picard Tools or SAMtools to remove PCR duplications to ensure data accuracy.

3. The method according to claim 2, characterized in that: The step 2 is specifically as follows: Cox regression based on the joint enumeration of two features is used to achieve the fitting of the prognosis survival time, and Pearson correlation based on the joint enumeration of two features is used to measure the correlation between features. The above two metrics are unified to mark the weights of any two node edges in the fully connected network; the cluster center automatic search clustering algorithm based on density descending is used to delete the edges with lower weights in the network, so as to realize network construction and feature cluster mining. The cluster center automatic search algorithm based on density descending is used to cluster the weights of all edges in the network in density descending order, and delete the edges with lower weight density from the network until a disconnected feature cluster appears; Step 2.1: Based on the Cox regression, the Cox proportional hazard regression model is used to fit the survival time by enumerating the pairwise features, so as to reveal the significant relationship between any pair of features and prognosis survival. By performing the maximum likelihood estimation on the logarithmic partial likelihood function of the hazard function, the regression coefficient β corresponding to any pair of features is obtained. i ,β j ) T , when it is assumed that the features are uncorrelated, any component of the regression coefficient obeys the Gaussian distribution, and the Wald test is used to determine whether the regression coefficient component deviates significantly from zero; The sample size can be further expanded by reordering the samples and their corresponding survival results. The reordered p-value corresponding to the Wald statistic can be generated on each feature. All pairwise features are enumerated to obtain paired p-values, which are used to select significant feature pairs closely related to survival from all pairwise feature combinations. In order to deal with overfitting, a penalty constraint is added to the logarithmic partial likelihood function corresponding to the risk function. Considering that the hazard ratio may change with survival time, a weighted Cox regression model is used instead. Step 2.2: Pearson correlation-based joint enumeration of two features. To measure the correlation between two features, the Pearson correlation coefficient r is used to calculate the correlation between any two features, and the t statistic of n samples is constructed based on this. The statistic is expressed as: By using the sample reordering method, the order of sample values ​​of one feature in any pair of features is disrupted, and the above t statistic is calculated repeatedly. The corresponding reordered p value is expressed as: Among them, t(i,j) is the statistic corresponding to the correlation coefficient r of features i and j, t b (i,j) represents the statistic generated by randomly reordering the values ​​of any feature, and B is the number of random reorderings; Step 2.3: Determine the association relationship of network nodes. The collaborative co-expression network uses a single feature as a network node and the relationship between two features as an edge. The weight of the edge between nodes is calculated by comprehensively considering the significant relationship between the feature and the survival time and the correlation between the features. The following metric is constructed to measure the weight of the edge between nodes: Where p2(i) and p2(j) represent the binary Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively, and p1(i) and p1(j) represent the univariate Cox regression reranking p-values ​​obtained by combining the two features i and j, respectively. represents rounding up; the first term of m(i,j) is the "synergistic" part that measures the contribution of features i and j to survival regression, and the second term is derived from the reordered p-value obtained by the joint enumeration of two features based on Pearson correlation, corresponding to the "co-expression" part; m(i,j) is used to measure the association between nodes i and j: the larger its value, the more likely it is that nodes i and j are derived from the same feature cluster, and vice versa. When the value of m(i,j) is less than zero, its value is set to the maximum value of the weight in the network, and the weighted undirected graph constructed is called a synergistic-co-expression network;. Step 2.4: Based on the feature cluster mining of the collaborative-co-expression network, for the weight of the collaborative-co-expression network edge, which is a measure of the node association relationship, an automatic cluster center search clustering algorithm based on density descending is introduced. The measurement threshold of the edge in the network is limited by the optimal automatic estimation of the number of categories and cluster boundaries to determine whether there is an associated edge between any two nodes in the network, thereby degenerating the network from a weighted undirected graph to an undirected graph.

4. The method according to claim 3, characterized in that: The step 3 is specifically as follows: Step 3.1: Bottom-up "cross-cluster enumeration" biomarker selection. According to the "cross-cluster hypothesis", single variables are enumerated within each feature cluster, and the single variables from c different feature clusters are combined into multivariate variables to calculate their regression coefficients for the Cox survival hazard ratio function; Step 3.2: Top-down "cross-cluster random" biomarker selection, introducing a sample resampling voting strategy, randomly selecting a single variable within each feature cluster, combining the single variables from c different feature clusters into multivariates, and calculating their regression coefficients for the Cox survival hazard ratio function; assuming that the features are uncorrelated, analogous to the pairwise feature joint enumeration method based on Cox regression, each group of "cross-cluster" multivariates can obtain corresponding c re-ranked p-values, and each round counts the multivariates with significant re-ranked p-values; after l rounds of voting, the features with the highest votes are selected from each feature cluster to form a significant prognostic feature set; Step 3.3: Unsupervised feature selection of joint multivariates; single variables need to be extracted from different feature clusters to form joint multivariates for regression analysis. For bottom-up "cross-cluster enumeration" biomarker selection, c re-ranked p-values ​​corresponding to each enumerated multivariate are recorded; for top-down "cross-cluster random" biomarker selection, multiple groups of c-dimensional p-value vectors corresponding to multiple randomly extracted multivariates can be recorded by setting a threshold. In the c-dimensional feature space, all enumerated or re-ranked p-value vectors that meet the threshold limit are automatically searched and clustered based on the descending density of the cluster center, and the cluster closest to the origin is selected to form a c-dimensional feature set. The Mahalanobis distance from the scattered points to the origin is calculated. When the bottom-up "cross-cluster enumeration" biomarker selection is adopted, the multidimensional variable with the closest distance to the origin is selected as the significant prognostic feature set; if the top-down "cross-cluster random" biomarker selection is adopted, it is still necessary to determine whether the multidimensional variable with the closest distance to the origin corresponds to the feature combination with the highest cumulative votes; Step 3.4: Construction of prognostic survival analysis indicators. Taking the Cox proportional hazard regression model as an example, since it is difficult to estimate the baseline risk function or baseline survival function, the linear part of the risk score is used to construct the prognostic survival indicator. The risk score of the sample is a linear combination of the corresponding sample values ​​of the prognostic characteristics, and the combination coefficient is the regression coefficient.

5. The method according to claim 4, characterized in that: The step 4 is specifically as follows: Step 4.1: Sample clustering for prognostic survival analysis: To discover the prognostic survival category, the key biomarkers from different omics or the prognostic survival analysis indicators built by different omics data are combined as the input of the convolutional neural network, and its output is regarded as a new feature subspace. The samples are automatically searched and clustered based on the cluster center in descending density, so as to discover the prognostic category to which the samples belong; Step 4.2: Draw Kaplan-Meier survival analysis curves for samples of different categories to observe the changes in survival time within the category and the differences in survival time between categories. Use the log-rank test to measure whether there are significant quantitative differences in survival time between different categories. Draw a risk score survival analysis graph to observe the changes in survival time between samples of different categories. Step 4.3: Analysis of significant differences in prognosis between classes. Class comparison is achieved using statistics based on sample distribution differences to evaluate whether there are significant differences in risk scores for different categories. On this basis, attempts are made to construct different types of classifiers using training samples, and the accuracy and effectiveness of the classification results are verified using the prognosis survival results of the test samples.

6. The method according to claim 5, characterized in that: Synergy represents the correlation between the feature pair composed of two biomarkers and prognostic survival time; co-expression represents the correlation between two biomarkers.

7. The method according to claim 4, characterized in that: The cross-group assumptions are as follows: The significant feature set closely related to the prognosis survival time is not entirely composed of a single significant variable related to the prognosis survival time, and there is no obvious correlation between the variables in the significant feature set.

8. A cancer prognosis survival analysis system based on a collaborative co-expression network, characterized by: The system comprises: A data acquisition module, which collects multi-omics data of cancer based on a public database and pre-processes the data; A network construction module, wherein the network construction module constructs a collaborative co-expression network for cancer multi-omics data based on the processed data; An extraction module, wherein the extraction module extracts prognostic biomarkers based on a collaborative co-expression network; A clustering module is used to cluster samples for survival analysis based on extracted prognostic biomarkers, and to perform prognostic survival analysis.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: The program is executed by a processor to implement the method according to claims 1-7.

10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the method of claims 1-7 is implemented.

Images