Esophageal cancer prognosis scoring method in combination with RSK pathway activity

Abstract

The invention discloses an esophageal cancer prognosis scoring method combined with RSK pathway activity, and belongs to the fields of bioinformatics, molecular prediction models and medical diagnosis assistance. Aiming at clinical requirements proposed in the technical background, the invention mainly aims to provide a risk scoring model constructed based on RSK pathway related genes and an application method thereof, survival prognosis and auxiliary treatment decision can be predicted in esophageal cancer patients, and the model has relatively good generalization ability, stability and mechanism interpretability.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of bioinformatics, molecular prediction model and medical diagnosis assistance, and particularly relates to a risk score model based on RSK (Ribosomal S6 kinase) pathway related genes and transcriptome data and an application method and system thereof in esophageal cancer patient prognosis prediction and individualized treatment decision. BACKGROUND

[0002] Esophageal cancer is one of the malignant tumors with high mortality rate in the world, and the incidence of esophageal squamous cell carcinoma is particularly prominent in some areas of Asia. Although progress has been made in surgery, chemotherapy, radiotherapy and targeted / immunotherapy in recent years, the overall prognosis is still not ideal. Clinically, there is an urgent need for biomarkers or prediction models that can accurately identify high-risk patients to assist individualized treatment strategies.

[0003] In tumor biology, RSK (Ribosomal S6 kinase, RPS6KA family protein kinase) is an important kinase downstream of the MAPK / ERK signaling pathway, which participates in the regulation of cell growth, proliferation, metabolism, apoptosis, migration, etc. by phosphorylating a series of target proteins (including transcription factors, metabolic enzymes, translation regulatory factors, etc.). In recent years, studies have shown that abnormal activation of RSK in various cancers may be related to tumor progression, drug resistance and metastasis. However, there is still a lack of methods that integrate RSK-related activity into the overall transcriptome data and construct quantitative models for patient prognosis prediction.

[0004] Currently, in the direction of cancer transcriptome / prognosis prediction, there are many risk models based on gene signatures (such as differential expression, Cox risk proportion regression model, Least Absolute Shrinkage and Selection Operator, etc.) for patient stratification prediction, but most of these models lack the combination of specific signaling pathways or functional modules, especially the lack of research that takes RSK pathway activity as the core starting point to guide model construction.

[0005] In addition, although there have been some reports on survival prediction models for esophageal cancer, the models are often not perfect in terms of external validation, stability analysis, mechanism explanation and integration with clinical variables. Therefore, there is an urgent need for a prognosis model and method that can combine pathway activity, data-driven and biological priori, with good generalization ability and clinical interpretability.

[0006] Based on the above technical status, the present application proposes a method and system that takes RSK pathway as the core, integrates transcriptome and clinical data, and can be used for esophageal cancer patient prognosis prediction and treatment assistance, in order to achieve a better balance between biological significance, prediction performance and clinical interpretability. SUMMARY

[0007] The present application aims at the clinical needs proposed in the technical background, and the main purpose is to provide a risk score model based on RSK pathway related genes and an application method thereof, which can predict the survival prognosis of esophageal cancer patients and assist in treatment decision making, and has good generalization ability, stability and mechanism explainability.

[0008] To achieve the above-mentioned application purposes, the present application adopts the following technical solutions: an esophageal cancer prognosis scoring method combined with RSK pathway activity, the method comprising:

[0009] Step A: data preparation and preprocessing;

[0010] A1: download the RNA-seq expression matrix and clinical follow-up data of esophageal cancer from TCGA as a data queue, which includes the whole gene sequencing data of the esophageal cancer patient tissue sample, the row name is the corresponding number name of the patient on TCGA, and the column name is the expression amount of tens of thousands of genes of the corresponding tissue sample of the patient, and TCGA represents the Cancer Genome Atlas;

[0011] A2: process the original expression matrix according to the following process: (1) gene ID mapping: if the gene ID is Ensembl ID, map the Ensembl ID to the standard gene symbol; if it is a chip probe ID, map the probe ID to the standard gene symbol. (2) expression standardization: first calculate TPM (transcripts per million, unit: number of transcripts per million) for the RNA-seq original count, and the calculation method is:

[0012]

[0013] Wherein, R i is the original sequencing count value of gene i, L i is the transcript length of gene i, and the summation is performed on all detected genes j, and N is the total number of detected genes; then take the logarithmic transformation of TPM:

[0014]

[0015] The expression matrix used for downstream analysis.

[0016] A3: annotate the survival state and survival time corresponding to each tissue sample in the data queue downloaded in A1;

[0017] Step B: construct a RSK related gene set, RSK represents RPS6KA family protein kinase;

[0018] B1: Pathway / target gene set based on prior knowledge: Known genes involved in signal transduction and downstream targets in the RSK pathway, including RPS6KA1 / 3, downstream transcription factors, regulatory proteins, and metabolic targets;

[0019] B2: Data-driven screening of related genes: Given that RPS6KA1 and RPS6KA3 are highly correlated with the RSK pathway, the correlation between RPS6KA1 and RPS6KA3 and all probe genes in the TCGA study cohort was calculated, and genes that are highly positively or negatively correlated with them were screened as a candidate set.

[0020] B3: Combine the pathway / target gene set of prior knowledge extracted from B1 with the data-driven related genes extracted from B2 to form a preliminary RSK-related gene candidate set;

[0021] Step C: Calculate the sample-level RSK activity score;

[0022] C1: The enrichment score (i.e., RSK activity score) of the RSK candidate gene set is calculated at the sample level using single-sample gene set enrichment analysis (ssGSEA) or gene set variation analysis (GSVA). ssGSEA is a method for scoring the enrichment of a single sample on a given gene set, outputting the enrichment score for each sample and each gene set; GSVA is a sample-level gene set activity estimation method based on nonparametric kernel estimation.

[0023] C2: Samples were divided into high-activity and low-activity groups based on the median RSK activity score. Differential expression analysis of each gene between the two groups was performed using the limma package. The log2 fold change (log2FC) and the p-value after multiple test correction (adj.P.Val) were calculated. Simultaneously, the Spearman correlation coefficient r and its p-value between each gene expression and the sample-level RSK activity score were calculated, and the correlation p-value was corrected using Benjamini-Hochberg FDR to obtain the FDR-corrected correlation coefficient rho_fdr. Genes satisfying |log2FC|≥log2(1.5), adj.P.Val≤0.05, |r|≥0.4, and rho_fdr≤0.05 were selected. A composite screening score S was defined. i :

[0024]

[0025] This score considers the intensity of difference, the degree of correlation with pathway activity, and the significance of the difference, facilitating comprehensive ranking. If no gene meets the above strict threshold, the composite screening score S is used. iThe top 50 genes were selected in descending order as backoff candidates to ensure sufficient feature input for subsequent univariate Cox / LASSO modeling.

[0026] Step D: Screening for prognostic-related genes;

[0027] D1: Perform univariate Cox proportional hazards regression analysis on each gene obtained from step C. Model format:

[0028]

[0029] Among them, expr i The expression level of gene i in the sample was then used, and a Wald test was employed with a p-value < 0.05 as the significance screening threshold.

[0030] D2: Subsequently, LASSO-Cox regression (Least Absolute Shrinkage and Selection Operator) was performed on the gene set obtained by single-factor screening, and the penalty parameter λ was selected through cross-validation to obtain the candidate feature subset.

[0031] D3: Based on the dimensionality reduction processing of the above multiple regression analysis, the following risk scoring formula is finally constructed based on the expression of these 11 genes and the regression coefficients:

[0032]

[0033] Where, β i For regression coefficients, expressive measure i This represents the expression level of the corresponding gene in the sample.

[0034] D4: Based on the median risk score, the samples are grouped into high-risk and low-risk groups.

[0035] Furthermore, the following methods were used to evaluate and validate the model performance;

[0036] Step 1: Perform Kaplan-Meier survival analysis within the study cohort, plot the ROC curves at 1 year, 3 years, and 5 years, and calculate the corresponding AUC to verify the model's predictive performance; the ROC curve represents the receiver operating characteristic curve, and the AUC represents the area under the ROC curve and the coordinate axis.

[0037] Step 2: Downstream mechanism and clinical correlation analysis;

[0038] F1: Gene differential expression analysis was performed on the high / low risk groups separately, and genes with significant differences were selected for enrichment analysis to elucidate the potential biological differential pathways between the two groups.

[0039] F2: Using the drug sensitivity prediction database GDSC, predict the difference in response to SL0101 between high- and low-risk groups, and plot a risk group vs. drug sensitivity comparison graph; GDSC represents anticancer drug sensitivity genomics, and SL0101 is a drug that acts on the RSK protein pathway.

[0040] F3: Construct a nomograph that integrates risk scores with several clinical variables to predict an individual’s survival probability at 1 / 3 / 5 years to assess predictive consistency.

[0041] The method of this invention can achieve the following technical effects or advantages:

[0042] 1. Quantify RSK pathway activity into sample-level scores to characterize the variation levels of this signaling pathway in different patients;

[0043] 2. Based on this, by combining differential expression, correlation analysis, univariate Cox and LASSO methods, genes that are highly associated with prognosis are screened out to form a robust risk score model (RSK-related risk score).

[0044] 3. Provide a systematic model evaluation process (KM analysis, ROC curve, etc.);

[0045] 4. Through enrichment analysis, drug sensitivity prediction, and clinical correlation analysis, reveal the biological mechanisms behind risk grouping and their relationship with clinical characteristics;

[0046] 5. Construct a nomogram to link risk scores with routine clinical variables, providing personalized survival probability prediction and decision support. Attached Figure Description

[0047] Figure 1 This is a scatter plot of bar graphs showing the high-risk and low-risk groups enriched on the RSK pathway / candidate gene set using ssGSEA.

[0048] Figure 2 This is a heatmap of expression on all samples in the study cohort across the 11 genes ultimately included in the risk model.

[0049] Figure 3 Survival analysis was performed after the samples were divided into high-risk and low-risk groups based on the constructed RSK-related Risk Score.

[0050] Figure 4 This is the ROC curve.

[0051] Figure 5 This is a graph comparing the predicted response of SL0101 drug between the high-risk and low-risk groups.

[0052] Figure 6 It is called Nomogram. Implementation

[0053] The method of the present invention will be described in detail below with reference to specific embodiments, but it should not be construed as limiting the scope of protection of the present invention.

[0054] Example 1: Data Preparation and Preprocessing

[0055] Gene expression levels and clinical follow-up information from TCGA-ESCA were downloaded from TCGA (https: / / portal.gdc.cancer.gov / ) to form a study cohort (n = 196).

[0056] Furthermore, the original TCGA FPKM was converted to TPM form, and log2(TPM + 1) was calculated to obtain the normalized expression matrix; for all samples, low-expression genes with TPM < 1 in more than 50% of the samples were removed.

[0057] Furthermore, ComBat was used to perform batch correction on multiple datasets to make the expression levels of each dataset more consistent, and each sample was annotated with its OS status and OS time.

[0058] Example 2: Constructing a candidate set of RSK-related genes

[0059] Searching databases such as KEGG, Reactome, and MSigDB, we identified 16 known target genes of the RSK signaling pathway and its downstream components, which are as follows:

[0060] RPS6KA1, RPS6KA2, RPS6KA3, RPS6KA6, MAPK3, YBX1, GAB2, PDCD4, BAD, ESR1, IRS2, SORBS2, CIC, DUSP6, RPS6, MTOR.

[0061] Furthermore, in the study cohort, the Pearson correlation coefficients between RPS6KA1 and all genes were calculated, and genes with |r| ≥ 0.5 and p < 0.05 were retained as candidates; similarly, similarity calculations were performed on RPS6KA3.

[0062] Furthermore, the two sets of candidate genes are combined with the gene set of the pathway knowledge base to obtain the final RSK-related candidate gene set.

[0063] Example 3: Calculation of sample-level RSK activity score and correlation analysis

[0064] The R package GSVA was used to calculate the enrichment score of the candidate gene set in each sample, named the RSK activity score, and the score of each sample in the study cohort was visualized.

[0065] The samples were divided into high-activity and low-activity groups based on the median RSK score; differential expression analysis was performed on the high vs. low-activity groups using the limma package to screen for up- and down-regulated genes.

[0066] The correlation between the expression of each candidate gene and the RSK score was calculated, and genes that were highly correlated with the score were selected as further modeling candidates.

[0067] Example 4: Screening of prognosis-related genes and construction of risk models

[0068] Univariate Cox regression was performed on the above candidate genes, and 13 genes with p < 0.05 were selected as preliminary candidates in the study cohort.

[0069] Furthermore, these 13 genes were regressed using the LASSO-Cox model, and the optimal λ value was selected through 10-fold cross-validation, resulting in 11 genes with non-zero coefficients as stable genes. Their regression coefficient β i The risk score model (RSK-related Risk Score) is constructed using the final regression coefficients for each gene.

[0070] The score is calculated as follows: RSK-related Risk Score = (0.215864776692421 × L) MID1IP1 ) + (-0.323696299 × L ACOT11 ) + (0.108927691570885 × L SLC17A1 ) + (-0.012259331 ×L GNL3L ) + (0.0742212902027045 × L DLL4 ) + (-0.206205801 × L JMJD8 ) +(0.0185269553552937 × L HEPACAM2 ) + (0.306009510784265 × L GXYLT1 )+(0.0489925075169666 × L CLCA1 )+ (0.140785855752085 × LMID1 )+ (0.162550339236486× L GK (The data is presented in a table format, where each characteristic gene name corresponds to its gene expression level in the sample.)

[0071] Furthermore, the samples are divided into high-risk and low-risk groups based on the median of the RiskScore or other split points.

[0072] Example 5: Model Evaluation and Validation

[0073] Kaplan-Meier survival curves were compared between the high- and low-risk groups in the study cohort, and the Log-rank p-value was calculated. ROC curves for 1 year, 3 years, and 5 years were plotted, and AUC values ​​were calculated to test the consistency between the predicted probability and the actual observation.

[0074] Example 6: Downstream Mechanisms and Clinical Association Analysis

[0075] After differential expression analysis between high- and low-risk groups in Example 3, significantly upregulated / downregulated genes (|log2FC| ≥ 1, adjusted p < 0.05) were identified. GO / KEGG / GSEA enrichment analysis was performed using tools such as clusterProfiler and fgsea to identify significantly enriched biological processes and pathways.

[0076] Using drug-sensitivity data from the GDSC database, combined with the R-package oncoPredict method, we predicted the half-maximal inhibitory concentration (IC50) or response value of each sample to several commonly used anticancer drugs, and compared the box plots of the predicted response differences between high- and low-risk groups on these drugs.

[0077] A nomogram model was built using risk scores as variables to provide an individual's 1-year / 3-year / 5-year survival probability, and a nomogram was plotted.

Claims

1. A prognostic scoring method for esophageal cancer that incorporates RSK pathway activity, the method comprising: Step A: Data preparation and preprocessing; A1: Download the RNA-seq expression matrix and clinical follow-up data of esophageal cancer from TCGA as a data cohort. This cohort includes all gene sequencing data of tissue samples from esophageal cancer patients. The row name of this cohort is the corresponding number name of the patient on TCGA, and the column name is the expression level of tens of thousands of genes in the corresponding tissue sample of the patient. TCGA stands for Cancer Genome Atlas. A2: Annotate the original representation matrix; A3: Annotates the corresponding survival status and survival time for each tissue sample in the data queue downloaded from A1; Step B: Construct the RSK-related gene set, where RSK stands for RPS6KA family protein kinases; B1: Pathway / target gene set based on prior knowledge: Known genes involved in signal transduction and downstream targets in the RSK pathway, including RPS6KA1 / 3, downstream transcription factors, regulatory proteins, and metabolic targets; B2: Data-driven screening of related genes: Given that RPS6KA1 and RPS6KA3 are highly correlated with the RSK pathway, the correlation between RPS6KA1 and RPS6KA3 and all probe genes in the TCGA study cohort was calculated, and genes that are highly positively or negatively correlated with them were screened as a candidate set. B3: Combine the pathway / target gene set of prior knowledge extracted from B1 with the data-driven related genes extracted from B2 to form a preliminary RSK-related gene candidate set; Step C: Calculate the sample-level RSK activity score; C1: The enrichment score of the RSK candidate gene set, i.e., the RSK activity score, is calculated at the sample level using single-sample gene set enrichment analysis or gene set variation analysis; among them, ssGSEA is a method for enriching a single sample on a given gene set, outputting the enrichment score for each sample and each gene set; GSVA is a sample-level gene set activity estimation method based on nonparametric kernel estimation. C2: The samples were divided into high-activity and low-activity groups based on the median RSK activity score. Differential expression analysis of each gene between the two groups was performed using the limma package. The log2 fold change of the gene and the p-value after multiple test correction were calculated. At the same time, the Spearman correlation coefficient r and its p-value between the expression of each gene and the sample-level RSK activity score were calculated. The correlation p-value was corrected by Benjamini-Hochberg FDR to obtain the FDR-corrected correlation coefficient rho_fdr. Genes that met the following conditions were selected: |log2FC|≥log2(1.5), adj.P.Val≤0.05, |r|≥0.4, and rho_fdr≤0.

05. If no gene meets the above strict threshold, then the composite screening score S is used. i The top 50 genes were selected in descending order as backoff candidates to ensure sufficient feature input for subsequent univariate Cox / LASSO modeling; Step D: Screening for prognostic-related genes; D1: Perform univariate Cox proportional hazards regression analysis on each gene obtained from step C. Model format: ； Among them, expr i The expression level of gene i in the sample was used, and a Wald test was then performed with a p-value < 0.05 as the significance screening threshold. D2: Subsequently, LASSO-Cox regression (Least Absolute Shrinkage and Selection Operator) was performed on the gene set obtained by single-factor screening, and the penalty parameter λ was selected through cross-validation to obtain the candidate feature subset; D3: Based on the dimensionality reduction processing of the above multiple regression analysis, the following risk scoring formula is finally constructed based on the expression of these 11 genes and the regression coefficients: ； Where, β i For regression coefficients, expressive measure i This represents the expression level of the corresponding gene in the sample. D4: Based on the median risk score, the samples are grouped into high-risk and low-risk groups.

2. The esophageal cancer prognostic scoring method based on RSK pathway activity as described in claim 1, characterized in that, The specific method for step A2 is as follows: (1) Gene ID mapping: If the gene ID is Ensembl ID, then Ensembl ID is mapped to a standard gene symbol; if it is a microarray probe ID, then probe ID is also mapped to a standard gene symbol. (2) Expression level standardization: For the raw RNA-seq count, the TPM (transcripts per million) is calculated first. The calculation method is as follows: ； Among them, R i L represents the raw sequencing count of gene i. i Let TPM be the transcript length of gene i, summed over all detected genes j, where N is the total number of detected genes; then, a logarithmic transformation is performed on the TPM: ； The expression matrix used for downstream analysis.

3. The esophageal cancer prognostic scoring method based on RSK pathway activity as described in claim 1, characterized in that, In step C2, the composite screening score S is defined. i : 。 4. The esophageal cancer prognostic scoring method based on RSK pathway activity as described in claim 1, characterized in that, The following methods were used for model performance evaluation and validation; Step 1: Perform Kaplan-Meier survival analysis within the study cohort, plot the ROC curves at 1 year, 3 years, and 5 years, and calculate the corresponding AUC to verify the model's predictive performance; the ROC curve represents the receiver operating characteristic curve, or ROC curve, and the AUC represents the area under the ROC curve and the coordinate axis. Step 2: Downstream mechanism and clinical correlation analysis; F1: Gene differential expression analysis was performed on the high / low risk groups separately, and genes with significant differences were selected for enrichment analysis to elucidate the potential biological differential pathways between the two groups. F2: Using the drug sensitivity prediction database GDSC, predict the difference in response to SL0101 between high- and low-risk groups, and plot a risk group vs. drug sensitivity comparison graph; GDSC represents anticancer drug sensitivity genomics, and SL0101 is a drug that acts on the RSK protein pathway. F3: Construct a nomograph that integrates risk scores with several clinical variables to predict an individual’s survival probability at 1 / 3 / 5 years to assess predictive consistency.

Images