Colorectal cancer drug relocation method based on multi-omics integration

Abstract

The invention discloses a colorectal cancer drug relocation method based on multi-omics integration. The system comprises a multi-omics data acquisition and preprocessing module, a tumor microenvironment analysis module, a specific disease network construction module, a multi-dimensional drug relocation module and a result evaluation module. And the tumor microenvironment analysis module comprises cell heterogeneity identification, cell map construction, cell annotation and tumor cell subset annotation. The specific disease network construction module comprises tumor feature expression program extraction, expression program screening, meta-program construction, clinical related meta-program recognition and specific disease protein interaction network construction. And the multi-dimensional drug relocation module comprises a module for identifying diseases by using a random walk algorithm, carrying out drug screening based on disturbance data, carrying out drug screening based on network proximity and carrying out comprehensive drug relocation. From the perspective of single cell data, element programs related to colorectal cancer survival are excavated, corresponding modules are designed, and the efficiency and precision of colorectal cancer targeted drug screening are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of computational biology, and particularly relates to a colorectal cancer drug repositioning method based on multi-omics integration. BACKGROUND

[0002] Colorectal cancer (CRC) is a malignant tumor of the digestive system with high morbidity and mortality worldwide. Despite various treatment methods such as surgery, radiotherapy, chemotherapy, and targeted therapy, the overall prognosis of colorectal cancer patients is still not ideal due to strong tumor cell heterogeneity, unclear biomarker identification, and rapid development of drug resistance.

[0003] With the development of single-cell sequencing technology (scRNA-seq), researchers can analyze the functional state, spatial distribution, and interaction characteristics of various cells in the tumor microenvironment at single-cell resolution, significantly improving the accuracy of understanding disease development. At the same time, complex network analysis methods can be used to construct disease-related gene regulatory networks or protein interaction networks, providing new theoretical foundations and technical means for identifying key modules of diseases and discovering drug targets.

[0004] The theory of "disease-related genes interacting in an interaction group" believes that disease phenotypes are caused by complex gene interaction networks. In this process, genes related to disease pathophysiology usually cluster together in a non-random manner. This non-random clustering sub-network is called "disease module". They participate in the formation of disease phenotypes as a whole. Previous studies have found that in such networks, disease-related genes not only tend to be adjacent to each other, but also have high connectivity between them. Further research suggests that even potential drug targets may be clustered in disease modules.

[0005] Traditional drug development has a long cycle and high cost, while drug repositioning technology (i.e., old drugs for new uses) as a low-cost and efficient drug development strategy has received widespread attention in cancer treatment research in recent years. The core of this technology is to use existing drugs for new disease indications, thereby saving time for early screening and toxicity verification.

[0006] Currently, although some studies have attempted to combine single-cell data with protein networks or drug perturbed cell line data to reveal disease mechanisms. However, there is still a lack of a perfect method system for systematically integrating single-cell data with these two to achieve high-precision drug repositioning prediction. SUMMARY

[0007] The purpose of the present application is to provide a colorectal cancer drug repositioning method based on multi-omics integration.

[0008] The present application comprises a multi-omics data acquisition and preprocessing module, a tumor microenvironment analysis module, a specific disease network construction module, a multi-dimensional drug repositioning module and a result evaluation module.

[0009] The multi-omics data acquisition and preprocessing module collects multi-omics data and performs standardized preprocessing on these data, providing a basis for the construction of subsequent modules; the multi-omics data includes single-cell sequencing data of colorectal cancer patients disclosed by the GEO database, batch transcriptome sequencing data and clinical data of colorectal cancer patients downloaded from the TCGA website, drug data and drug target gene data downloaded from the DrugBank database, three-stage data of drug perturbed colorectal cancer cell lines obtained from the LINCS database, and human protein-protein interaction network data.

[0010] The tumor microenvironment analysis module includes cell heterogeneity recognition, construction of cell atlas, cell annotation and tumor cell subpopulation annotation. The cell heterogeneity recognition is based on single-cell transcriptome sequencing data, and uses unsupervised analysis method to recognize cell subpopulations with heterogeneity in the tissue, depicts the composition and state distribution of different cell types in the tumor microenvironment, including cell dimension reduction and visualization, cell clustering. The construction of cell atlas is based on UMAP nonlinear dimension reduction algorithm, and the results are mapped to the graph for two-dimensional visualization, forming an intuitive cell atlas, each point representing a single cell, and its position in the graph reflecting its expression similarity with other cells. The cell annotation is based on the clustering results, and each cell cluster is functionally annotated. The tumor cell subpopulation annotation is further in-depth mining of the functional heterogeneity on the basis of the tumor cells identified in the preliminary annotation.

[0011] The specific disease network construction module comprises tumor characteristic expression program extraction, expression program screening, meta-program construction, clinical relevant meta-program identification, and specific disease protein interaction network construction. The tumor characteristic expression program extraction is to obtain tumor cells in cell annotation results, divide the single cell expression matrix according to sample sources, filter out samples with small cell numbers, select genes with high expression amounts, and calculate gene variation score weights to obtain a weight matrix, which is used for weighted non-negative matrix factorization method to factorize the expression matrix, extract potential tumor characteristic expression programs, and obtain two matrices of cell x expression program and expression program x gene. The expression program screening is to screen a representative expression characteristic set in different tumor samples to construct a robust expression program. The meta-program construction is to cluster the robust expression program to construct a meta-program. First, the gene overlap ratio between the robust expression programs is calculated, and a similarity matrix is constructed according to the gene overlap ratio. Then, the similarity is converted into a distance value to form a distance matrix. The ward method is used for hierarchical clustering of the robust expression program to obtain the meta-program. The clinical relevant meta-program identification uses the ssGSEA algorithm to score the bulk RNA-seq data to evaluate the enrichment degree of the meta-program in the TCGA-COAD and TCGA-READ samples. According to the median of the meta-program score, the samples are grouped into enrichment and non-enrichment groups. The Kaplan-Meier curve is used to evaluate the correlation between the enrichment degree of each meta-program and the survival time of the samples, and the meta-program with P-Value≤0.05 is selected as the clinical relevant meta-program. The specific disease protein interaction network construction uses the AUCell scoring function to evaluate the expression of the clinical relevant meta-program in each tumor cell subpopulation in the single cell data of the tumor cell subpopulation. The tumor cell subpopulation with significant expression of the clinical relevant meta-program is extracted to construct a cell-specific protein interaction network, and the protein-protein interaction with a Pearson correlation coefficient P-Value≤0.05 is retained.

[0012] The multi-dimensional drug repositioning module comprises disease module identification using random walk algorithm, drug screening based on perturbation data, drug screening based on network proximity, and comprehensive drug repositioning. The disease module identification using random walk algorithm is to take the number of occurrences of each gene in the meta-program as the node weight to construct a cell-specific protein interaction network with weight, and perform random walk algorithm on the network to identify the disease. The drug screening based on perturbation data is to calculate the perturbation score of all drug perturbation experiments on the disease module, and the perturbation score <0, significance level <0.05 is the predicted potential target drug. The drug screening based on network proximity is to calculate the network proximity of the disease module and each drug, and the network proximity less than a set value, a significance level A drug with a value less than 0.05 is a predicted potential targeted drug. The comprehensive drug repositioning is to select the intersection of drug screening based on perturbation data and drug screening based on network proximity as the final predicted potential targeted drug.

[0013] The result evaluation module includes result evaluation of drug prediction results using AUC index of ROC curve, F1-score index and biological experiment verification. Among them, the AUC index of ROC curve result evaluation: the abscissa is the false positive rate , the ordinate is the true positive rate , draw the ROC curve and calculate the area AUC under the ROC curve, find the model parameters under the optimal prediction result, AUC is used to evaluate the overall classification ability of the model, the larger the AUC value, the closer to 1 indicates that the model performance is better. F1-score index result evaluation: the accuracy, precision, recall and F1-score index are used to quantitatively evaluate whether the predicted drug has clinical application potential. , , , ; wherein, represents the number of non-clinical drugs predicted as non-clinical drugs, represents the number of clinical drugs predicted as non-clinical drugs, represents the number of non-clinical drugs predicted as clinical drugs, represents the number of clinical drugs predicted as clinical drugs; and the F1-score index represents the harmonic mean of the model accuracy and recall rate.

[0014] The present application fuses single-cell sequencing data, protein interaction network and drug perturbation data, and proposes a systematic and repeatable drug repositioning method, which improves the efficiency and accuracy of colorectal cancer targeted drug screening, and realizes systematic drug repositioning from molecular mechanism analysis to clinical application. The method can be applied to other tumor or complex disease data-driven drug discovery research with expression heterogeneity. BRIEF DESCRIPTION OF DRAWINGS

[0015] Figure 1 It is a whole schematic diagram of the method of the present application. DETAILED DESCRIPTION

[0016] As shown in Figure 1 , the multi-omics integrated colorectal cancer precision drug repositioning method includes a multi-omics data acquisition and preprocessing module, a tumor microenvironment analysis module, a specific disease network construction module, a multi-dimensional drug repositioning module and a result evaluation module.

[0017] I. Multi-omics Data Acquisition and Preprocessing Module: This module collects data from multiple omics disciplines, including single-cell transcriptomics, batch transcriptomics, pharmacomics, and proteomics, and performs standardized preprocessing on this data, providing a foundation for the construction of subsequent modules and algorithm design. This module includes the following sub-tasks: (1-1) Multi-omics data collection: The module's data includes: single-cell sequencing data (scRNA-seq) of colorectal cancer patients published in the GEO database; bulk transcriptome sequencing data (bulk RNA-seq) and clinical data (TCGA-COAD / READ) of colorectal cancer patients downloaded from the TCGA website; drug data and drug target gene data downloaded from the DrugBank database; data on three phases (Phase I, II, and III) of drug-perturbed colorectal cancer cell lines (HT29, HT116, etc.) obtained from the LINCS database; and human protein-protein interaction network data.

[0018] (1-2) Preprocessing of single-cell sequencing data: (1-2-1) Quality Control: In single-cell sequencing data, the number of expressed genes, total expression level (UMI counts), number of gene types, and the proportion of mitochondrial and ribosomal genes in each cell are first counted. Based on these indicators, low-quality cells with abnormal gene expression levels or excessively high mitochondrial / ribosomal ratios are removed to ensure the accuracy and reliability of downstream analysis.

[0019] (1-2-2) Standardization and Normalization: To eliminate systematic bias caused by sequencing depths between different cells, the LogNormalize method was used to standardize the original expression matrix. Specifically, the normalization process was as follows: To achieve relative uniformity in expression levels, `count` represents the raw UMI count of a gene in a cell, and `total` represents the total expression level of all genes detected in that cell. Then, the `ScaleData` function is used to normalize the gene expression values, correcting for differences in order of magnitude between different genes.

[0020] (1-2-3) High-variability gene screening and feature normalization: High-variability genes with large expression differences (variance, standard deviation) are extracted from the whole gene expression matrix and selected as high-information feature genes, representing potential differences in cell state or functional changes. Subsequently, these feature genes are normalized to a uniform scale to provide representative input features for subsequent analysis steps such as dimensionality compression and cell clustering.

[0021] (1-3) Preprocessing of TCGA clinical data: The bulk RNA-seq gene expression matrix data of TCGA-COAD and TCGA-READ were converted into TPM (Transcripts Per Million) format, and each value in the matrix was added by 1 and then taken the natural logarithm for standardization processing. The clinical data were obtained for age, gender, survival status, survival time, etc. The gene expression matrix data and clinical data were matched according to the sample name, and non-tumor samples and samples with unknown survival status were filtered out.

[0022] (1-4) Drug and protein network data preprocessing: (1-4-1) Drug target data preprocessing: The pairing information of drugs and their target genes was extracted from the DrugBank database, the drug name format was unified, and the DrugBank ID was used as the unique identifier. The target gene symbol was standardized to HGNC official naming and verified by Entrez ID. After eliminating the entries with missing information or mapping failure, a high-quality correspondence table of drugs and their target genes was constructed to provide standardized target set input for subsequent network analysis.

[0023] (1-4-2) Drug perturbation data preprocessing: Drug treatment experiment data of colorectal cancer related cell lines were obtained from the LINCS database. The screening fields included drug name, cell line number, dose information, treatment time, treatment type (such as DMSO or Drug), and corresponding gene expression matrix. First, the DMSO control group and drug treatment group of the same cell line and treatment time were matched in pairs to construct a standardized drug perturbation experiment instance. Then, experimental records with missing information or less than three times of repetition were eliminated to improve data reliability. The gene expression matrix of each sample was converted to TPM and standardized, and then the expression fold change of the drug treatment group compared with the corresponding control group at each gene was calculated. According to the expression change, all genes were sorted in descending order to construct a sorted list of perturbation response: where g represents the gene, and the sequence order reflects the regulation strength of the drug on the gene expression.

[0024] (1-4-3) Protein network data preprocessing: High-confidence human protein-protein interaction (PPI) data were downloaded from authoritative protein interaction databases (such as STRING, BioGRID). The interaction pairs with experimental verification support were retained, and the protein naming format was unified (such as Gene Symbol) to construct a directed network graph structure, where the node represents the protein and the edge represents the functional interaction between proteins. This network will serve as the basic framework for network distance calculation and module identification of disease modules and drug targets.

[0025] II. Tumor microenvironment analysis module, including cell heterogeneity recognition, cell atlas construction and cell annotation, and tumor cell subpopulation annotation. In addition to tumor cells, the tumor microenvironment also includes immune cells, stromal cells, etc. To perform dimensionality reduction and clustering analysis on single-cell gene expression matrix data, obtain a two-dimensional spatial position map of the cells, and according to the cell marker genes, different cells are annotated into cell lineages, a cell atlas is constructed, and the data of tumor cells is separately taken out for further process analysis, and according to the annotation of cell subpopulations according to marker genes, this process is called tumor microenvironment analysis.

[0026] (2-1) Cell heterogeneity recognition: This module aims to identify cell subpopulations with heterogeneity within the tissue based on single-cell transcriptome sequencing data using unsupervised analysis methods, and to characterize the composition and state distribution of different cell types in the tumor microenvironment.

[0027] (2-1-1) Cell dimensionality reduction and visualization: Due to the high dimensionality of single-cell transcriptome data, it is easily affected by dimensionality and noise interference. Therefore, this method uses unsupervised dimensionality reduction to compress the features of high-dimensional expression matrix and extract the most representative expression patterns. PCA (Principal Component Analysis) is used for linear dimensionality reduction of high-variable genes, and the principal components are extracted as the main expression characteristics of cells; and on this basis, further through nonlinear mapping UMAP (Uniform Manifold Approximation and Projection), the local structure between cells is preserved, which helps to reveal the potential cell clustering structure and improve the effect of subsequent clustering and visualization.

[0028] (2-1-2) Cell clustering: In the reduced space, use unsupervised clustering method to identify cell subpopulations. The specific steps include: Adjacency graph construction: based on the PCA low-dimensional representation, the similarity between any pair of cells is calculated, and a KNN graph is constructed to represent the relationship network between cells.

[0029] Unsupervised clustering analysis: Louvain community discovery algorithm is used to mine the modular structure on the KNN graph, and the cells are automatically divided without human label intervention, realizing data-driven subpopulation identification.

[0030] This process can efficiently identify the clustering characteristics of main cell types such as immune cells, stromal cells and tumor cells.

[0031] (2-2) Cell atlas construction and annotation: To fully characterize the composition and distribution characteristics of various cells in the tumor microenvironment, based on dimensionality reduction and clustering analysis, a cell atlas is constructed and cell lineages are annotated.

[0032] (2-2-1) Cell atlas construction: On the basis of the UMAP nonlinear dimensionality reduction algorithm, the results are mapped to a two-dimensional visual map. While maintaining the similarity of the local structure, the global topological information is maximally preserved, the results are mapped to the graph to form an intuitive cell atlas, each point represents a single cell, and its position in the graph reflects its expression similarity with other cells. Provide a structural basis for subsequent cell type annotation and microenvironment analysis.

[0033] (2-2-2) Cell annotation: Based on the clustering results, each cell cluster is functionally annotated, and the specific process is as follows: Marker gene identification: statistics of significantly up-regulated genes in each cell cluster; Reference comparison: match these genes with the cell type marker gene set defined in the literature (such as PTPRC representing immune cells, COL1A1 representing fibroblasts, EPCAM representing epithelial cells, etc.); Automatic annotation and classification output: according to the expression pattern and reference database (such as PanglaoDB, CellMarker) to annotate the cell major category; the cell output is classified into main types such as tumor cells, immune cells (T cells, B cells, macrophages) and stromal cells (fibroblasts, endothelial cells) and so on.

[0034] This process combines several methods to annotate cells, laying the foundation for subsequent tumor cell refinement analysis, immune atlas reconstruction, and cell-cell interaction research.

[0035] (2-3) Tumor cell subpopulation annotation On the basis of tumor cells identified in the preliminary annotation, further explore their functional heterogeneity. The specific method is as follows: Functional gene set scoring: reference to the functionally related gene set defined in previous literature (such as proliferation, EMT, angiogenesis, etc.), using the AUCell algorithm (SCENIC package) gene set scoring method to score each tumor cell; Subgroup functional annotation: according to the scores of different functional modules, identify tumor cell subgroups with specific biological states, and annotate them with functional labels such as "proliferative", "stemness", "immune escape" and so on; Visualization display: map different functional subgroups to the UMAP atlas to show their differences in spatial distribution and expression characteristics.

[0036] This process helps to reveal the diversity of functional states of tumor cells in the microenvironment, providing support for disease progression interpretation and targeted therapy.

[0037] III. Specific disease network construction module, after identifying different cell types in single cells, this module focuses on tumor cell subgroups, and further explores key expression programs closely related to disease progression. This module includes tumor characteristic expression program extraction, expression program screening, construction of sub-programs, identification of clinically relevant sub-programs, and construction of specific disease protein interaction networks.

[0038] (3-1). Tumor characteristic expression program extraction: Obtain tumor cells in cell annotation results and divide single cell expression matrix according to sample source, filter out samples with less than 50 cells, select the top 7000 genes with the highest expression in each sample, and calculate CNV (gene variation score) weight to obtain a weight matrix, which is used for weighted non-negative matrix factorization (wNMF) method to factorize the expression matrix, extract potential tumor characteristic expression programs, and obtain cell x expression program and expression program x gene two matrices.

[0039] Mathematically, the problem can be represented by the following formula: ; Where C is the weight matrix, V is the single cell matrix of cell x gene, W is the cell x expression program matrix, and H is the expression program x gene matrix. Generally, the top 50 genes with the highest score in the gene sequence corresponding to each expression program are selected as the tumor characteristic expression program (Gene Signature) of the expression program.

[0040] (3-2). Expression program screening: To screen out a set of representative expression characteristics in different tumor samples, construct a robust expression program (Robust Program) that meets the following standards: a. Intra-sample stability: If the gene overlap ratio between two expression programs in the same tumor sample reaches or exceeds 70%, it is considered that the expression program shows consistency within the sample and has intra-sample stability.

[0041] b. Inter-sample consistency: If the gene overlap rate of two expression programs from different tumor samples is not less than 20%, it means that they have shared characteristics across samples and are considered to have inter-sample consistency.

[0042] c. Intra-sample redundancy identification: If two robust expression programs in the same tumor sample have a gene overlap ratio of more than 20%, it is determined that there is redundancy and needs to be removed or merged in subsequent processing.

[0043] (3-3) Constructing meta-programs: This part clusters the robust expression programs to construct meta-programs. First, the gene overlap ratio between robust expression programs is calculated, and a similarity matrix is constructed accordingly. Then, the similarity is converted into distance values to form a distance matrix. The robust expression programs are hierarchically clustered using the ward method to obtain meta-programs. Meta-programs consisting of 10 or more robust expression programs are retained. Then, the frequency of all genes in the meta-programs in their robust expression programs is counted, and genes with a frequency higher than 25% are retained as the gene characteristics of the meta-programs. Finally, meta-programs with more than 25 genes are retained for subsequent analysis.

[0044] (3-4) Identifying clinically relevant meta-programs: The ssGSEA algorithm is used to score the bulk RNA-seq data to evaluate the enrichment of meta-programs in TCGA-COAD and TCGA-READ samples. According to the median of the meta-program scores, the samples are grouped into enrichment and non-enrichment groups. The correlation between the enrichment degree of each meta-program and the survival time of the samples is evaluated using the Kaplan-Meier curve, and meta-programs with P-Value≤0.05 (Log-rank test) are selected as clinically relevant meta-programs.

[0045] (3-5) Constructing specific disease protein interaction networks: The AUCell scoring function is used to evaluate the expression of clinically relevant meta-programs in each tumor cell subpopulation in single-cell data. The tumor cell subpopulation with significant expression of clinically relevant meta-programs is extracted to construct cell-specific protein interaction networks, and protein-protein interactions with a Pearson correlation coefficient P-Value≤0.05 are retained.

[0046] Four, multi-dimensional drug repositioning module, including using random walk algorithm to identify diseases, drug screening based on perturbation data, drug screening based on network proximity, and comprehensive drug repositioning.

[0047] (4-1) Random walk algorithm to identify disease module: The number of occurrences of each gene (one-to-one corresponding to protein) in the meta-program is taken as the node weight to construct a weighted cell-specific protein interaction network. The random walk algorithm is performed on this network to identify disease modules, a. Randomly select a (protein) node from the network as the initial node of the network module.

[0048] b. In each iteration, select a candidate node i from the first-order neighbors of the current selected node. When the candidate node i meets both connectivity significance and module score increase, add the candidate node i to the network module.

[0049] ; denotes the connectivity significance of node i, is the number of nodes connected to node i in the network module, is the degree of node i, m is the number of nodes in the network module, and N is the number of nodes in the protein-protein interaction network. When i ≤ 0.05, it indicates that node i satisfies the connectivity significance condition.

[0050] denotes the network module score before adding node i, denotes the network module score after adding node i, s(i) is the weight of node i, M represents the network module after adding the node, and μ is the average node weight of the protein-protein interaction network. When , it indicates that node i satisfies the module score increase condition.

[0051] Repeat steps a and b until no node can satisfy the connectivity significance and module score increase conditions.

[0052] To generate the final disease module, a large number of network modules are usually generated first, and then network modules with less than 10 nodes are removed. The network modules are sorted in descending order according to the module score defined above, and the network modules with scores in the top 5% are selected. The nodes with a frequency greater than 10% are retained to construct the final disease module.

[0053] (4-2). Drug screening based on perturbation data: Calculate the perturbation score of all drug perturbation experiments on the disease module. The drugs with high and significant ( <0.05) perturbation scores are the predicted potential targeted drugs, where denotes the significance level after multiple test correction, and less than 0.05 is considered significant. For a given disease module gene set S and the gene differential fold change ranking list before and after each drug perturbation Perform GSEA enrichment analysis.

[0054] ; wherein, is the cumulative weight score of hit items (weighted hits), is the cumulative penalty score of miss items, is the differential fold score of gene in the gene ranking list, and is the weight parameter. denotes the enrichment score, denotes the significance of the permutation test, and S denotes the given disease module protein set. ​​​​​

[0055] <0.05、

[0056] (4-3). Drug screening based on network proximity: Calculate the network proximity between disease module and each drug (target), drugs with small and significant network proximity (p<0.05) are predicted potential targeted drugs. Given a disease module protein set S and drug target set T, the network proximity is calculated as .

[0057] ; where is the number of elements in set S, is the number of elements in set T, is the shortest distance between node s and node t in the network, calculated using the permutation test .

[0058] (4-4). Comprehensive drug repositioning: i.e. drugs with close network proximity to disease module and down-regulate the gene expression of disease module, such drugs are candidate drugs. Select the intersection of the candidate drugs of (4-2) and (4-3) as the final predicted potential targeted drugs.

[0059] Five, the result evaluation module, including the drug prediction results using ROC curve AUC index, F1-score index and biological experiment verification result evaluation.

[0060] Network proximity less than the set value, drug perturbation score is significant (p<0.05) and perturbation score <0.05) or perturbation score ≥0 drugs have clinical value, set as positive samples; network proximity greater than or equal to the set value, drug perturbation score is not significant (p≥0.05) or perturbation score ≥0 drugs are set as negative samples.

[0061] ​​​​​​​​​The AUC index of the ROC curve result evaluation: the accuracy of the model is measured by calculating the area AUC under the ROC curve, the ROC curve is a tool for evaluating the performance of a binary classification model, which intuitively reflects the discriminant ability of the model by plotting the relationship between the false positive rate ( ) and the true positive rate ( ) of the model at different thresholds. The horizontal axis false positive rate represents the proportion of all actual negative samples that are incorrectly predicted as positive. The vertical axis true positive rate represents the proportion of all actual positive samples that are correctly predicted as positive: , . By plotting curve and calculating , and finding the model parameters under the optimal prediction result, where The area under the curve is defined as to evaluate the overall classification ability of the model, and The larger the value, the closer the value to 1 indicates the better performance of the model.

[0062] F1-score index result evaluation: the accuracy accuracy, precision, recall and F1-score indicators are used to quantitatively evaluate whether the predicted drug has clinical application potential: , , , ; wherein, The number of non-clinical drugs predicted as non-clinical drugs, The number of clinical drugs predicted as non-clinical drugs, The number of non-clinical drugs predicted as clinical drugs, The number of clinical drugs predicted as clinical drugs. The F1-score index represents the harmonic mean of the accuracy and recall of the model.

[0063] Biological experiment verification: In addition to evaluating the prediction results by the above two methods, the predicted drugs are verified by literature related to colorectal cancer drug prediction, and the drugs not reported in the literature are verified by biological experiments.

Claims

1. A multi-omics-integrated drug repositioning method for colorectal cancer, characterized by: It includes modules for multi-omics data acquisition and preprocessing, tumor microenvironment analysis, specific disease network construction, multi-dimensional drug repositioning, and result evaluation. The multi-omics data acquisition and preprocessing module acquires multi-omics data and performs standardized preprocessing on these data to provide a foundation for the construction of subsequent modules. The tumor microenvironment analysis module includes cell heterogeneity identification, cell atlas construction, cell annotation, and tumor cell subset annotation. The specific disease network construction module includes tumor feature expression program extraction, expression program screening, metaprogram construction, identification of clinically relevant metaprograms, and construction of specific disease protein interaction network; The multi-dimensional drug relocation module includes a disease identification module using a random walk algorithm, drug screening based on perturbation data, drug screening based on network proximity, and comprehensive drug relocation. The result evaluation module includes evaluating the drug prediction results using the AUC index of the ROC curve, the F1-score index, and biological experimental verification.

2. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 1, characterized in that, In the multi-omics data acquisition and preprocessing module, the multi-omics data includes: single-cell sequencing data of colorectal cancer patients published in the GEO database, batch transcriptome sequencing data and clinical data of colorectal cancer patients downloaded from the TCGA website, drug data and drug target gene data downloaded from the DrugBank database, data on three stages of drug-induced perturbation of colorectal cancer cell lines obtained from the LINCS database, and human protein-protein interaction network data. The specific preprocessing steps for single-cell sequencing data are as follows: (A) Quality control: First, count the number of expressed genes, total expression level, number of gene types, and the proportion of mitochondrial genes and ribosome genes in each cell; based on the above indicators, remove low-quality cells with abnormal gene expression levels or mitochondrial / ribosome ratios greater than a set threshold. (B) Standardization and Normalization: The original representation matrix is ​​standardized using the LogNormalize method, in the form of... 'count' represents the raw UMI count of a gene in a cell, and 'total' represents the total expression level of all genes detected in the cell. Then, the ScaleData function is used to normalize the gene expression values ​​and correct for differences in the order of magnitude between different genes. (C) High-mutation gene screening and feature normalization: High-mutation genes with expression differences greater than a set value are extracted from the whole gene expression matrix and used as high-information feature genes, representing potential differences in cell state or functional changes. The expression differences include variance and standard deviation. Then, these feature genes are normalized to a uniform scale and used as representative input features. The specific preprocessing of TCGA bulk transcriptome sequencing data and clinical data is as follows: the bulk RNA-seq gene expression matrix data of TCGA-COAD and TCGA-READ are converted into TPM format, and each value in the matrix is ​​incremented by 1 and the natural logarithm is taken for standardization. Information on age, sex, survival status and survival time is obtained from the clinical data. The gene expression matrix data and clinical data are matched according to the sample name, and non-tumor samples and samples with unknown survival status are filtered out. The preprocessing of drug data and drug target gene data specifically involves: extracting the pairing information between drugs and their target genes from the DrugBank database, standardizing the drug name format and using DrugBank ID as the unique identifier, standardizing the target gene symbols to the official HGNC name, and verifying them through Entrez ID; after removing entries with missing information or failed mapping, constructing a correspondence table between drugs and their target genes to provide a standardized set of target genes for subsequent network analysis. Preprocessing of drug perturbation data: Drug treatment experimental data of colorectal cancer-related cell lines were obtained from the LINCS database. Screening fields included: drug name, cell line number, dosage information, treatment time, treatment type, and corresponding gene expression matrix. First, DMSO control groups and drug treatment groups with the same cell line and treatment time were paired to construct standardized drug perturbation experimental instances. Then, experimental records with missing experimental information or fewer than three repetitions were removed to improve data reliability. TPM transformation and standardization were performed on the gene expression matrix of each sample, and the fold change in expression of each gene in the drug treatment group compared to the corresponding control group was calculated. All genes were sorted in descending order based on expression changes to construct a perturbation response ranking list, where the sequence order reflects the intensity of drug regulation of gene expression. Protein network data preprocessing: Human protein-protein interaction data were downloaded from the protein interaction database; interaction pairs with experimental validation were retained, and the naming format of proteins was standardized. An undirected network graph structure was constructed, in which nodes represent proteins and edges represent functional interactions between proteins; this network serves as the basic framework for calculating network distances between disease modules and drug targets and for module identification.

3. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 1, characterized in that: In the tumor microenvironment analysis module, the cellular heterogeneity identification is based on single-cell transcriptome sequencing data. Unsupervised analysis methods are used to identify heterogeneous cell subpopulations within the tissue, characterizing the composition and state distribution of different cell types in the tumor microenvironment; including: (1) Cell dimensionality reduction and visualization: Unsupervised dimensionality reduction is used to compress the features of the high-dimensional expression matrix and extract representative expression patterns; Principal component analysis is used to linearly reduce the dimensionality of highly variable genes and extract principal components as the main expression features of cells; Local structures between cells are preserved through nonlinear mapping UMAP. (2) Cell clustering; In the dimensionality-reduced space, unsupervised clustering methods are used to identify cell subpopulations; First, an adjacency graph is constructed, and based on the low-dimensional representation of principal component analysis, the similarity between any pair of cells is calculated, and a KNN graph is constructed to represent the relationship network between cells; Then, unsupervised clustering analysis is performed: the Louvain community detection algorithm is used to mine modular structures on the KNN graph, and cells are automatically divided to achieve data-driven subpopulation identification; The construction of the cell atlas is specifically based on the UMAP nonlinear dimensionality reduction algorithm, and the results are visualized and mapped in two dimensions. The results are mapped onto the graph to form an intuitive cell atlas, where each point represents a single cell, and its position in the graph reflects its expression similarity with other cells. The cell annotation mentioned above is specifically based on clustering results, performing functional annotation on each cell cluster. The specific process is as follows: Marker gene identification: Identify significantly upregulated genes in each cell cluster; Reference comparison: Match significantly upregulated genes with the set of cell type marker genes defined in the literature; Annotation and classification output: Cell categories are annotated based on expression patterns and reference databases; cell outputs are classified into major types, including tumor cells, immune cells, and stromal cells; The tumor cell subset annotation described above is based on the tumor cells identified in the initial annotation, and further explores their functional heterogeneity; the specific method is as follows: Functional gene set scoring: Referring to the functionally relevant gene sets defined in previous literature, each tumor cell was scored using the AUCell algorithm gene set scoring method; Subpopulation functional annotation: Based on the scores of different functional modules, tumor cell subpopulations with specific biological states are identified and annotated with functional labels; Visualization: Different functional subgroups are mapped onto the UMAP map to show their differences in spatial distribution and expression characteristics.

4. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 1, characterized in that: In the specific disease network construction module, the extraction of tumor feature expression program specifically involves: obtaining tumor cells from the cell annotation results, dividing the single-cell expression matrix according to the sample source, filtering out samples with fewer than A cells, selecting the top B genes with the highest expression levels in each sample and calculating the gene variation score weights to obtain a weight matrix, using a weighted nonnegative matrix factorization method to factorize the expression matrix, extracting potential tumor feature expression programs, and obtaining two matrices: cell × expression program and expression program × gene, where A = 30–80 and B = 5000–10000. The expression program screening involves selecting a representative set of expression features from different tumor samples and constructing a robust expression program. The construction of the metaprogram involves clustering robust expression programs to construct the metaprogram. First, the gene overlap ratio between robust expression programs is calculated, and a similarity matrix is ​​constructed accordingly. Then, the similarity is converted into distance values ​​to form a distance matrix. The Ward method is used to perform hierarchical clustering of the robust expression programs to obtain the metaprogram. The method for identifying clinically relevant metaprograms involves using the ssGSEA algorithm to score bulk RNA-seq data and assess the enrichment of metaprograms in TCGA-COAD and TCGA-READ samples. Samples are then grouped into enriched and non-enriched groups based on the median metaprogram score. Kaplan-Meier curves are used to evaluate the correlation between the enrichment of each metaprogram and sample survival time, and metaprograms with a P-value ≤ 0.05 are selected as clinically relevant metaprograms. The aforementioned construction of a disease-specific protein-protein interaction network involves using the AUCell scoring function to evaluate the expression of clinically relevant meta-programs in single-cell data of tumor cell subpopulations; extracting tumor cell subpopulations with significantly expressed clinically relevant meta-programs to construct a cell-specific protein-protein interaction network, while retaining protein-protein interactions with a Pearson correlation coefficient (P-Value) ≤ 0.

05.

5. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 4, characterized in that: The representative set of expression features meets the following criteria: a. Intrasample stability: If the gene overlap between two expression programs in the same tumor sample reaches or exceeds 70%, the expression program is considered to show consistency within the sample and has the characteristic of intrasample stability. b. Inter-sample consistency: If two expression programs from different tumor samples have a gene overlap rate of not less than 20%, it indicates that they share characteristics across samples and are considered to have inter-sample consistency. c. Intrasample redundancy identification: If two expression programs that have been identified as robust in the same tumor sample have a gene overlap ratio of more than 20%, they are considered to have redundancy and need to be deduplicated or merged in subsequent processing.

6. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 4, characterized in that: In the construction of the metaprogram, a metaprogram consisting of 10 or more robust expression programs is retained. The frequency of occurrence of all genes in the metaprogram in its robust expression program is counted. Genes with a frequency higher than 25% are retained as gene features of the metaprogram. Metaprograms with more than 25 genes are retained for subsequent analysis.

7. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 1, characterized in that: In the multi-dimensional drug relocation module, the random walk algorithm for disease identification module uses the occurrence frequency of each gene in the metaprogram as the node weight to construct a weighted cell-specific protein interaction network, and executes the random walk algorithm on the network to identify diseases. The aforementioned drug screening based on perturbation data calculates the perturbation score of all drug perturbation experiments on the disease module. <0, significance level Drugs with a value <0.05 are considered potential targeted drugs; The aforementioned drug screening based on network proximity calculates the network proximity between the disease module and each drug. Less than the set value, significance level Drugs with a value <0.05 are considered potential targeted drugs; The aforementioned integrated drug relocation selects the intersection of drug screening based on perturbation data and drug screening based on network proximity as the final predicted potential targeted drug.

8. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 7, characterized in that: The specific details of the random walk algorithm for disease identification module are as follows: ① Randomly select a node from the network as the initial node of the network module; ② In each iteration, candidate node i is selected from the first-order neighbors of the currently selected node. When candidate node i satisfies both connectivity significance and module score increase, candidate node i is added to the network. Connectivity saliency of node i ; This represents the number of nodes in the network module that are connected to node i. It is a node i The degree, m is the number of nodes in the network module, and N is the number of nodes in the protein interaction network; when When ≤0.05, it means that node i satisfies the connectivity significance condition; Network module score before adding node i ; Network module score after adding node i ; s(i) is the weight of node i, M represents the network module after adding the node, and μ is the average node weight of the protein interaction network; when When, it means that node i satisfies the condition for increasing module score; Repeat steps ① and ② until no node satisfies the connectivity significance and module score increase conditions.

9. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 7, characterized in that: The aforementioned drug screening based on perturbation data specifically involves: sorting a given set of genes for a disease module, S, and a list of gene differences before and after each drug perturbation. Perform GSEA enrichment analysis; , , ;in, The cumulative weighted score for the hit items. The cumulative penalty score for each missed item. Genes in the gene sorting list Difference multiple score, These are weight parameters. Represents enrichment fractions. The significance of the permutation test is indicated by S, where S represents the set of proteins in a given disease module. right Standardization is performed to obtain the drug perturbation score. ,right Multiple hypothesis testing correction was performed to obtain Repeat the above steps until the drug perturbation score and significance result are calculated for each drug experiment, and retain them. <0.05、 For drugs with a perturbation score <0, the perturbation score is used according to the drug perturbation fraction. The candidates were sorted in ascending order. The aforementioned drug screening based on network proximity specifically refers to: Given a disease module protein set S and a drug target set T, calculate the network proximity. ; ;in, This represents the number of elements in set S. This represents the number of elements in set T. The shortest distance between nodes s and t in the network is represented by the permutation test. .

10. The method for drug repositioning in colorectal cancer based on multi-omics integration as described in claim 1, characterized in that: In the results evaluation module, the AUC index of the ROC curve is evaluated: the horizontal axis represents the false positive rate. The vertical axis represents the true positive rate. Plot the ROC curve and calculate the area under the ROC curve (AUC) to find the model parameters with the best prediction results. AUC is used to evaluate the overall classification ability of the model. The larger the AUC value and the closer it is to 1, the better the model performance. F1-score evaluation: Accuracy, precision, recall, and F1-score are used to quantitatively assess the clinical application potential of the predicted drugs. , , , ;in, This indicates the number of non-clinical drugs predicted as non-clinical drugs. This indicates the number of clinical drugs predicted as non-clinical drugs. This indicates the number of non-clinical drugs predicted to be clinical drugs. This indicates the number of clinical drugs predicted as clinical drugs; the F1-score represents the harmonic mean of the model's precision and recall.

Images