Hepatocellular carcinoma prognosis model construction method based on immunogen cell death related gene and application

Abstract

The invention relates to the technical field of hepatocellular carcinoma, in particular to a hepatocellular carcinoma prognosis model construction method based on immunogen cell death related genes and application, and an accurate prognosis prediction model is constructed by integrating the immunogen cell death related genes (IRGs) and molecular characteristics of hepatocellular carcinoma. A training set and a verification set provided by TCGA and ICGC databases are utilized, so that the model can perform effective sample analysis under a large-scale data background. Through consistency clustering analysis, the optimal clustering number is determined, the samples are orderly divided into different molecular subtypes, and it is ensured that the samples in each subtype have similar molecular characteristics. According to the invention, the capability of capturing liver cancer heterogeneity on the molecular level of the model is increased, so that the accuracy of prognosis prediction is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of hepatocellular carcinoma, in particular to a hepatocellular carcinoma prognosis model construction method based on immunogenic cell death related genes and application. BACKGROUND

[0002] Hepatocellular carcinoma, also known as liver cancer, is the most common type of primary liver cancer, with a global incidence ranking fifth among the most common cancers and the fourth leading cause of cancer-related deaths. Men are at higher risk of liver cancer than women. There are mainly two types of primary liver cancer - liver hepatocellular carcinoma (LIHC) and intrahepatic cholangiocarcinoma, and less common cancers such as angiosarcoma and hepatoblastoma. LIHC originates from hepatocytes, the main parenchymal cells of the liver. LIHC accounts for more than 80% of primary liver cancer cases worldwide. The pathogenesis of LIHC is complex and involves multiple factors, such as cell cycle dysregulation, DNA methylation changes, chromosomal instability, immune regulation, and epithelial-mesenchymal transition. Although different serum or tissue prognostic genes have been studied, their clinical utility has not been widely accepted. In addition, 80% of patients are already in the advanced stage at the first visit, losing the opportunity for radical surgery. Despite radical resection of LIHC in some patients, the high recurrence and metastasis rates after surgery remain a major challenge for patient survival. In addition, despite the emergence of new tumor treatment methods in recent years, such as immunotherapy using checkpoint inhibitors, new prognostic genes have not yet entered routine practice. Therefore, we urgently need an effective prediction model to accurately assess patient prognosis and screen for prognostic markers that reflect tumor progression at the molecular level, which may help to more accurately achieve individualized survival prediction and treatment.

[0003] Immunogenic cell death (ICD) is a form of regulated cell death that is induced by certain chemotherapeutic drugs, oncolytic viruses, physical and chemical therapies, photodynamic therapy, and radiotherapy. It can activate the immune system in immunocompetent hosts to fight cancer. ICD can release damage-associated molecular patterns (DAMPs) from dying tumor cells, thereby activating tumor-specific immune responses and eliciting long-term efficacy of anticancer drugs through a combination of direct killing of cancer cells and antitumor immunity. Studies have shown that inhibition of STAT3 can induce ICD in LIHC cells. Mecheliolide (a natural guaiacol lactone) can induce ICD in LIHC cells in vitro and in vivo through reactive oxygen species-mediated endoplasmic reticulum stress. FOLFOX (a combination therapy) can significantly inhibit LIHC tumor growth through ICD-related antitumor immunity.

[0004] Therefore, the present application aims to identify prognostic genes in LIHC patients based on immunogenic cell death related genes (IRGs), and perform molecular subtype typing, construct and verify the prognostic model of IRGs, so as to provide certain support for the prognosis and treatment of LIHC. SUMMARY

[0005] The present application aims to provide a liver cancer prognosis model construction method and application based on immunogenic cell death related genes, identify prognostic genes in LIHC patients based on immunogenic cell death related genes (IRGs), and perform molecular subtype typing, construct and verify the prognostic model of IRGs, so as to provide support for the prognosis and treatment of LIHC.

[0006] To achieve the above technical purposes and achieve the above technical effects, the present application is implemented by the following technical solutions:

[0007] A liver cancer prognosis model construction method based on immunogenic cell death related genes (IRGs), comprising the following steps:

[0008] S1: Obtain the gene expression data of liver cancer tissue samples and normal tissues in the TCGA-LIHC database as a training set, and obtain the gene expression data of liver cancer tissue samples in the LIRI-JP project of the ICGC database as a validation set;

[0009] S2: Perform consensus clustering analysis on the training set samples based on immunogenic cell death related genes using consistency clustering analysis, determine the optimal cluster number according to the cumulative distribution function, and classify the samples into subtype 1 and subtype 2;

[0010] S3: Use the DEseq2 package to screen the differentially expressed genes DEGs1 between HCC tissues and normal tissues, threshold: |log2FC|>1.5, adjust.P<0.05; screen the differentially expressed genes DEGs2 between subtype 1 and subtype 2, threshold: |log2FC|>1.5 and adj.P<0.05; take the intersection of DEGs1 and DEGs2, and obtain IRGs-related differentially expressed genes (IRGs-DEGs);

[0011] S4: Perform univariate Cox regression analysis (R package survival, threshold P<0.01) on the IRGs-DEGs obtained in step S3; screen out significant genes for LASSO regression analysis (R package glmnet, λmin=0.00972); and determine the prognostic genes;

[0012] S5: Establishing the risk score formula:

[0013] Risk score = 0.03210648*EXP(KRT20) + 0.02462363*EXP(MMP12) + 0.02566009*EXP(AGR2) + 0.02882947*EXP(CXCL5) - 0.02414141*EXP(CYP3A4) + 0.06227324*EXP(MAGEA8);

[0014] S6: Using the median of the risk score of the training set (TCGA-LIHC) as the threshold to divide the patients into high-risk and low-risk groups, and grouping the patients in the validation set according to the median risk score, and verifying the survival difference between the high-risk and low-risk groups through the Kaplan-Meier survival curve.

[0015] Further, the prognostic genes include: KRT20, MMP12, AGR2, CXCL5, CYP3A4, MAGEA8.

[0016] Further, the immunogenic cell death-related genes IRGs in step S2 include ENTPD1, NT5E, CALR, HMGB1, HSP90AA1, ATG5, BAX, CASP8, PDIA3, EIF2AK3, PIK3CA, CXCR3, IFNA1, IFNB1, IL10, IL6, TNF, CASP1, IL1R1, IL1B, NLRP3, P2RX7, LY96, MYD88, TLR4, CD4, CD8B, CD8A, FOXP3, IFNG, IFNGR1, IL17A, IL17RA, PRF1, a total of 34 genes.

[0017] In another aspect, the present application provides a biomarker combination for prognosis evaluation of hepatocellular carcinoma, which includes the following 6 genes: KRT20, MMP12, AGR2, CXCL5, CYP3A4, MAGEA8.

[0018] Further, the biomarker combination calculates the prognosis risk by detecting the gene expression amount in the liver tissue sample and substituting into the risk score formula:

[0019] Risk score = 0.03210648*EXP(KRT20) + 0.02462363*EXP(MMP12) + 0.02566009*EXP(AGR2) + 0.02882947*EXP(CXCL5) - 0.02414141*EXP(CYP3A4) + 0.06227324*EXP(MAGEA8).

[0020] RiskScore > 0.1427 indicates that the patient belongs to the high-risk group (significantly reduced overall survival);

[0021] RiskScore ≤ 0.1427 indicates that the patient belongs to the low-risk group (increased response rate to immunotherapy).

[0022] In another aspect, the present application provides a kit for detecting the above-mentioned biomarker combination, comprising: primer pairs or probes specific for amplifying KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8 genes; and reagent components for qRT-PCR detection thereof.

[0023] The beneficial effects of the present application are:

[0024] The present application integrates immunogenic cell death-related genes (IRGs) and molecular characteristics of hepatocellular carcinoma to construct a precise prognostic prediction model. The training set and validation set provided by the TCGA and ICGC databases enable the model to effectively analyze samples in a large-scale data context. Through consistency clustering analysis, the optimal cluster number is determined, and the samples are sequentially divided into different molecular subtypes, ensuring that the samples within each subtype have similar molecular characteristics. The present application increases the ability of the model to capture the heterogeneity of liver cancer at the molecular level, thereby improving the accuracy of prognosis prediction. The model further optimizes gene selection through LASSO regression, enabling it to retain the most valuable gene combination in a complex biological context, greatly improving the accuracy of individual patient prognosis prediction.

[0025] The present application divides patients into high-risk and low-risk groups by constructing a risk score formula. By using the mRNA expression levels of key prognostic genes such as KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8, individual risk can be quantified, and personalized clinical decisions can be guided accordingly. Not only can the overall survival of patients be predicted, but the predictive power of the model can also be evaluated through ROC curve analysis, showing a significant correlation between risk score and individual prognosis. This process combines gene expression differences, clinical pathological feature analysis and survival data to provide a comprehensive tool to optimize treatment plans, reduce ineffective medical expenses and side effect risks caused by inappropriate treatment.

[0026] The present application deeply analyzes the differences in immune cell infiltration between high and low risk groups, and finds that there is significantly enhanced immune cell activity in the high risk group, especially in activated CD4 T cells, dendritic cells, etc. However, this phenomenon may suggest the existence of tumor immune escape mechanism, thereby affecting the treatment effect. By quantifying the abundance of various immune cell infiltrations through the ssGSEA algorithm, combined with Spearman correlation analysis, the complex interaction between the immune features in the prognostic model and the specific immune microenvironment of liver cancer is revealed. Not only does it provide new insights into the immune escape mechanism of liver cancer, but also provides data support for achieving more effective immunotherapy, which helps to design more targeted immunotherapy strategies.

[0027] The present application enables KRT20, MMP12 and other genes to enter clinical detection as reliable biomarkers by screening specific prognostic related genes, helping to estimate the treatment response and prognosis of patients. The mRNA expression levels of these genes have been verified to have differences between normal and disease states, providing evidence for their effectiveness as biomarkers. Combined with drug sensitivity analysis, this scheme also provides a new tool for evaluating potential chemotherapy efficacy and immunotherapy response, demonstrating the application potential of these genes in treatment strategy optimization. This not only accelerates the clinical transformation of biomarkers, but also provides basic support for the development of targeted therapeutic drugs for liver cancer in the future.

[0028] Of course, implementing any product of the present application does not necessarily require all the advantages described above to be achieved at the same time. BRIEF DESCRIPTION OF DRAWINGS

[0029] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed for the embodiment description will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor.

[0030] Figure 1 Heat map of expression values of IRGs in normal and LIHC samples;

[0031] Figure 2 Schematic diagram of ConsensusClusterPlus analysis; Left: consensus matrix heat map, the consensus matrix is a sample square matrix, the numerical value represents the possibility of two belonging to the same cluster, the value range is from 0 to 1, and the color is from white to dark blue, Right: Consensus cumulative distribution function diagram, for each K corresponding consensus matrix, a 100 bin column chart is used to calculate the cumulative distribution;

[0032] Figure 3 Figure 7 is a schematic diagram of KM survival curves plotted based on different sample subtypes; Top: horizontal axis represents survival time of each patient, vertical axis represents survival probability, blue is subtype 2, red is subtype 1; Bottom: number of samples of different subtypes at different times;

[0033] Figure 4 Figure 8 is a volcano plot and a heat map of differentially expressed genes between normal and disease; A: volcano plot; horizontal axis represents differential fold log2(fold change), vertical axis represents -log10(FDR). Each point in the figure represents a gene, the upper orange triangle represents up-regulated differential genes, the lower blue triangle represents down-regulated differential genes, gray represents non-significant differential genes, and Noregulated represents unknown differential expression genes; B: heat map, each small square represents the expression value of each gene, and the color represents the expression amount of the gene, the larger the expression amount, the redder the color (red for high expression, blue for low expression), and different colors represent different groups;

[0034] Figure 5 Figure 9 is a volcano plot and a heat map of differentially expressed genes between different subtypes; A: volcano plot; horizontal axis represents differential fold log2(fold change), vertical axis represents -log10(FDR). Each point in the figure represents a gene, the upper orange triangle represents up-regulated differential genes, the lower blue triangle represents down-regulated differential genes, gray represents non-significant differential genes, and Noregulated represents unknown differential expression genes; B: heat map, each small square represents the expression value of each gene, and the color represents the expression amount of the gene, the larger the expression amount, the redder the color (red for high expression, blue for low expression), and different colors represent different groups;

[0035] Figure 6 Figure 10 is a Venn diagram of differentially expressed genes between normal and disease and differentially expressed genes between different subtypes; Green DEGs1 are differentially expressed genes between normal and disease, and red DEGs2 are differentially expressed genes between different subtypes;

[0036] Figure 7Figure 1. Forest plot for univariate Cox regression analysis. The first column is the gene, the second column is the statistical P value of the corresponding gene, the third column is the numerical change of hazard ratio (HR) and 95% confidence interval, where the hazard ratio is used to estimate the fold change of the risk of outcome events due to the presence of a certain factor, and the fourth column is the graphical representation of HR and its confidence interval. The vertical line in the middle of the figure is the null line, i.e. HR = 1, indicating that there is no statistical association between the studied factor and the outcome. Each horizontal line is the 95% confidence interval of the study, and the small square (or other graphics) in the center of the horizontal line is the point estimate of the HR value. If the horizontal line of the 95% confidence interval of a study intersects with the null vertical line, it is considered that there is no statistical association between the studied factor and the outcome. If the horizontal line falls to the left of the null vertical line, it is considered that the studied factor is beneficial to the occurrence of the outcome. If the horizontal line falls to the right of the null line, it is considered that the studied factor is not conducive to the occurrence of the outcome. Red represents HR greater than 1, and green represents HR less than 1.

[0037] Figure 8 Figure 2. LASSO regression schematic diagram. Left: The graph of the penalty term parameter, the horizontal coordinate is the log(lambda) value, and the vertical coordinate is the degree of freedom, representing the cross-validation error. In actual analysis, the position where the cross-validation error is the smallest is expected, and in the figure, the left dashed line position is the position where the cross-validation error is the smallest. According to this position (lambda.min), the top horizontal coordinate log(Lambda) is determined, and the number of characteristic genes is displayed on the upper side. Find the optimal log(Lambda) value, and find the corresponding genes and their coefficients in the right graph. Right: The horizontal coordinate is log(lambda), and the vertical coordinate is the coefficient of the gene. Different variables change their coefficients after being penalized by λ;

[0038] Figure 9 Figure 3. KM survival curve schematic diagram based on high and low risk groups according to risk score. Upper: The horizontal coordinate is the survival time of each patient, and the vertical coordinate is the survival probability. Red is the high risk group, and green is the low risk group. Lower: The horizontal coordinate is the survival time of each patient, and the vertical coordinate is the high risk group and the low risk group.

[0039] Figure 10 Figure 4. Risk factor association diagram of the training set. A: Risk score curve, horizontal coordinate for each patient, vertical coordinate for risk score, red for high risk group, blue for low risk group. B: Shows the survival status of patients in the training set, horizontal coordinate for each patient, vertical coordinate for survival time, red for survival, blue for death. The dashed line represents the median risk score, and the patients are divided into low risk group and high risk group. C: Expression distribution of prognostic genes. The redder the color, the higher the expression value, and vice versa. The horizontal coordinate is the patient, and the vertical coordinate is each gene.

[0040] Figure 11Figure 1 is a schematic diagram of a ROC curve for a training set; the horizontal X-axis is the false positive rate, the closer the X-axis is to zero, the higher the accuracy; the vertical Y-axis is called the true positive rate, the larger the Y-axis, the true positive and false positive can be constructed by evaluating different thresholds, a curve is formed, which extends from the lower left to the upper right, and bends to the upper left. A classifier with no discrimination between positive and negative classes will form a diagonal line, at false positive rate 0 and true positive rate 0 and false positive rate 1 and true positive rate. Different colors represent different times; AUC is the area under the curve, the larger, the higher the accuracy;

[0041] Figure 12 Figure 2 is a schematic diagram of a KM survival curve for a validation set: based on the risk score, the KM survival curve of the high-risk group and the low-risk group is drawn; upper: the horizontal coordinate is the survival time of each patient, the vertical coordinate is the survival probability, red is the high-risk group, green is the low-risk group; lower: the horizontal coordinate is the survival time of each patient, the vertical coordinate is the high-risk group and the low-risk group;

[0042] Figure 13 Figure 3 is a schematic diagram of a risk factor correlation graph for a validation set: A: risk score curve, horizontal coordinate for each patient, vertical coordinate for risk score, red for high-risk group, blue for low-risk group; B: shows the survival status of patients in the training set, horizontal coordinate for each patient, vertical coordinate for survival time, red for survival, blue for death; the dotted line represents the median risk score, and the patients are divided into low-risk group and high-risk group; C: expression distribution of prognostic genes, the redder the color, the higher the expression value, and vice versa, the horizontal coordinate is the patient, and the vertical coordinate is each gene;

[0043] Figure 14 Figure 4 is a schematic diagram of a ROC curve for a validation set; the horizontal X-axis is the false positive rate, the closer the X-axis is to zero, the higher the accuracy; the vertical Y-axis is called the true positive rate, the larger the Y-axis, the true positive and false positive can be constructed by evaluating different thresholds, a curve is formed, which extends from the lower left to the upper right, and bends to the upper left. A classifier with no discrimination between positive and negative classes will form a diagonal line, at false positive rate 0 and true positive rate 0 and false positive rate 1 and true positive rate. Different colors represent different times; AUC is the area under the curve, the larger, the higher the accuracy;

[0044] Figure 15 Figure 5 is a schematic diagram of risk scores for different categories of different clinicopathological features; the horizontal coordinate is the category of different clinicopathological features, and the vertical coordinate is the risk score, p < 0.001:***, p < 0.01:**, p ≤ 0.05: *;

[0045] Figure 16 Figure 6 is a schematic diagram of KM curves for different categories of different clinicopathological features;

[0046] Figure 17Forest plot for regression analysis; A: forest plot for univariate regression analysis; B: forest plot for multivariate regression analysis. The first column is the relevant feature, the second column is the statistical P value of the corresponding feature, the third column is the numerical change of hazard ratio (HR) and 95% confidence interval, where the hazard ratio is used to estimate the risk change of the outcome event due to the existence of a certain factor, the fourth column is the graphical representation of HR and its confidence interval, the vertical line in the middle of the figure is the null line, i.e. HR = 1, indicating that there is no statistical association between the studied factor and the outcome; if the horizontal line falls to the left of the null vertical line, it is considered that the studied factor is beneficial to the occurrence of the outcome; if the horizontal line falls to the right of the null line, it is considered that the studied factor is not beneficial to the occurrence of the outcome. Red represents HR greater than 1, and green represents HR less than 1;

[0047] Figure 18 Nomogram for integrated risk score; the nomogram includes three parts, the first part is the variable name in the prediction model, such as risk score, period and gender in the figure, each variable corresponding to the line segment is marked with a scale, representing the value range of the variable, and the length of the line segment reflects the contribution of the factor to the outcome event; the second part is the score, including the single score, i.e. points in the figure, representing the single score corresponding to each variable at different values; the total score, i.e. Total Point, representing the total score of the sum of the single scores corresponding to all variables after taking values; the third part is the prediction probability, i.e. 1-year-Survival, 2-year-Survival and 3-year-Survival;

[0048] Figure 19 Calibration plot for different survival periods; the x-axis represents the probability predicted by the nomogram, and the y-axis represents the actual probability. Perfect prediction corresponds to the gray dashed line. Different colors represent different years, and the bias is corrected by Bootstrapping (1000 times repetition), which represents the observed performance of the nomogram;

[0049] Figure 20 CXCL5-GSEA enrichment plot; A: GO; B: KEGG; the figure mainly shows part of the enrichment items, and the figure is divided into three parts, the uppermost part shows the enrichment score (ES) value calculation process, and an ES value is calculated for each gene from left to right, and connected into a line. There is a particularly obvious peak value on the left side, which is the ES value of the gene set phenotype. The middle part of the figure represents a gene in the gene set and its ranking position in the gene list. The lowermost part shows the matrix of the association between genes and phenotypes, and the first 5 with the smallest P value are selected;

[0050] Figure 21Figure 6A and 6B show the distribution of different immune cell infiltrations between high and low risk groups; A: the horizontal axis represents different cells, the vertical axis represents infiltration scores, and the * between different groups represents significance, p < 0.0001: ****, p < 0.001: ***, p < 0.01: **, p < 0.05: *; B: the horizontal axis represents different samples, the vertical axis represents different immune cells, and the color represents the infiltration condition. The redder the color, the higher the proportion, and vice versa;

[0051] Figure 22 Figure 7 shows the correlation heat map between different immune cells; the vertical axis represents different immune cells, the horizontal axis represents the correlation value, the color represents the P value, the redder the color, the smaller the P value, and vice versa. The size of the point represents the absolute value of the correlation. The larger the point, the stronger the correlation;

[0052] Figure 23 Figure 8 shows the expression distribution of different immune checkpoint expression values between high and low risk groups; the horizontal axis represents different immune checkpoints, and the vertical axis represents the expression value of the immune checkpoint. The * between different groups represents significance, p < 0.0001: ****, p < 0.001: ***, p < 0.01: **, p < 0.05: *;

[0053] Figure 24 Figure 9 shows the correlation between immune checkpoints and risk scores that meet the correlation conditions; the outermost circle is the immune checkpoint and the risk score, and the two are connected by a line, representing the correlation between the two. The one shown is in the case where P < 0.05 and the absolute value of the correlation is greater than 0.3;

[0054] Figure 25 Figure 10 shows the distribution of TIDE between high and low risk groups; the horizontal axis represents different groups, and the vertical axis represents TIDE. The * between different groups represents significance, p < 0.0001: ****, p < 0.001: ***, p < 0.01: **, p < 0.05: *;

[0055] Figure 26 Figure 11 shows the immune therapy response; left: the number of immune therapy responders in high and low risk groups. The horizontal axis represents different high and low risk groups, and the vertical axis represents the number. Right: the distribution of risk scores in different immune therapy responders;

[0056] Figure 27Figure 1 shows the distribution of IC50 values for drugs in high-risk groups; the horizontal axis is different high and low risk groups, the vertical axis is the log10 (estimated IC50 value distribution), the middle is the size of P value, the smaller the P value, the more asterisks, the significance between different groups is represented, p < 0.0001: ****, p < 0.001: ***, p < 0.01: **, p ≤ 0.05: *;

[0057] Figure 28 Figure 2 shows the expression levels of prognostic genes in control and LIHC samples; A: training set; B: validation set; the abscissa is different genes, and the ordinate is the expression level of the gene

[0058] Figure 29 Figure 3 shows the overall flowchart of the present application. DETAILED DESCRIPTION

[0059] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0060] Embodiment 1

[0061] The method for constructing a hepatocellular carcinoma prognosis model based on immunogenic cell death-related genes (IRGs) described in this embodiment comprises:

[0062] Molecular subtype typing: based on the TCGA-LIHC database of HCC samples, unsupervised clustering was performed on 34 IRGs (source: PMID: 27057433) using ConsensusClusterPlus package, the best clustering number K = 2 was determined according to the cumulative distribution function (CDF), and the samples were divided into subtype 1 and subtype 2;

[0063] Differential gene screening: DEseq2 package was used to screen the differentially expressed genes DEGs1 between HCC tissues and normal tissues (threshold: |log2FC|>1.5, adjust.P<0.05), and 2166 genes were obtained;

[0064] Screening of differentially expressed genes DEGs2 between subtype 1 and subtype 2 (threshold: |log2FC|>1.5 and adj.P<0.05), 112 genes were obtained;

[0065] The intersection of DEGs1 and DEGs2 was taken to obtain 72 IRGs-related differentially expressed genes (IRGs-DEGs);

[0066] Prognostic gene identification: univariate Cox regression analysis was performed on 72 IRGs-DEGs (R package survival, threshold P<0.01); LASSO regression analysis was performed on significant genes (R package glmnet, λ min =0.00972); prognostic genes were screened;

[0067] Risk model construction: a risk score formula was established:

[0068] Riskscore=0.03210648×EXP(KRT20)+ 0.02462363×EXP(MMP12) + 0.02566009×EXP(AGR2)+0.02882947×EXP(CXCL5) - 0.02414141×EXP(CYP3A4)+0.06227324×EXP(MAGEA8);

[0069] Risk grouping: the patients were divided into high-risk group and low-risk group with the median value of risk score of the training set (TCGA-LIHC) 0.1427 as the threshold.

[0070] The training set was 371 HCC samples and 50 normal samples in TCGA-LIHC database, and 342 samples with survival time>30 days were reserved for survival analysis;

[0071] The verification set was 240 HCC samples in ICGC LIRI-JP database, and 237 samples were reserved for survival analysis.

[0072] In this embodiment, the prognostic genes include: KRT20, MMP12, AGR2, CXCL5, CYP3A4, MAGEA8;

[0073] In this embodiment, the immunogenic cell death-related genes IRGs include ENTPD1, NT5E, CALR, HMGB1, HSP90AA1, ATG5, BAX, CASP8, PDIA3, EIF2AK3, PIK3CA, CXCR3, IFNA1, IFNB1, IL10, IL6, TNF, CASP1, IL1R1, IL1B, NLRP3, P2RX7, LY96, MYD88, TLR4, CD4, CD8B, CD8A, FOXP3, IFNG, IFNGR1, IL17A, IL17RA, PRF1, a total of 34 genes;

[0074] On the other hand, the present invention proposes a combination of biomarkers for prognostic assessment of hepatocellular carcinoma, including the following 6 genes: KRT20, MMP12, AGR2, CXCL5, CYP3A4, and MAGEA8.

[0075] In this embodiment, the biomarker combination calculates prognostic risk by detecting gene expression levels in liver tissue samples and substituting them into a risk scoring formula:

[0076] RiskScore=0.03210648×EXP(KRT20)+ 0.02462363×EXP(MMP12) + 0.02566009×EXP(AGR2)+0.02882947×EXP(CXCL5) - 0.02414141×EXP(CYP3A4)+0.06227324×EXP(MAGEA8);

[0077] A RiskScore > 0.1427 indicates that the patient belongs to the high-risk group (significantly reduced overall survival);

[0078] A RiskScore ≤ 0.1427 indicates that the patient belongs to the low-risk group (increased response rate to immunotherapy).

[0079] On the other hand, the present invention proposes a kit for detecting combinations of biomarkers as described above, comprising: primer pairs or probes for specifically amplifying the genes KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8; and reagent components thereof for qRT-PCR detection.

[0080] Example 2

[0081] Consistent clustering analysis

[0082] First, the expression values ​​of IRGs in normal and LIHC samples from the TCGA-LIHC database are presented, such as... Figure 1 As shown, compared with cancer samples, ENTPD1, ​​CALR, HSP90AA1, BAX, CASP8, PDIA3, etc. are expressed at lower levels in the disease, while MYD88, TLR4, PRF1, IL1B, IL10, etc. are expressed at higher levels in the disease.

[0083] Secondly, disease samples in the TCGA-LIHC cohort were genotyped using Consensus Clustering. The data were then divided into different clusters. While there were significant differences in molecular patterns between samples in different clusters, the molecular patterns of samples within each cluster were quite similar, thus achieving the goal of molecular genotyping for a large sample cohort.

[0084] First, unsupervised clustering analysis was performed using the ConsensusClusterPlus package with all genes in the training set, and the samples were genotyped based on the CDF curve. The results are as follows: Figure 2 As shown, all samples were clustered into two classes. KM curves were also plotted, and the results were analyzed assuming P < 0.05. Figure 3 The study showed that different subtypes were associated with prognosis, and found that subtype 2 had a worse prognosis than subtype 1.

[0085] Identification of differentially expressed mRNAs in LIHC

[0086] Normal and disease

[0087] Differentially expressed genes (DEGs) are genes whose expression varies significantly at the RNA level under different environmental stresses, time, and spatial conditions. Differential expression analysis aims to identify differentially expressed genes among different sample groups.

[0088] The DEseq2 package was used to compare the differences in mRNA expression levels between the normal group and the disease group. The differentially expressed gene screening criteria were: |log2FC|>1.5, adjust.P<0.05. Volcano plots were generated based on the results. Compared with the control group, a total of 2166 differentially expressed mRNAs (DEGs2) were obtained through differential expression analysis of TCGA-LIHC; among them, 1702 were upregulated and 495 were downregulated. The expression heatmap and volcano plots of the differentially expressed mRNAs were then displayed using the R package pheatmap. Figure 4 As shown.

[0089] Different subtypes

[0090] The DEseq2 package was used to compare the differences in mRNA expression levels among different subtypes. The differentially expressed gene screening criteria were: |log2FC|>1.5, adjust.P<0.05. Compared with the control group, a total of 112 differentially expressed mRNAs (DEGs2) were identified through differential expression analysis, of which 51 were upregulated and 61 were downregulated. The expression heatmap and volcano plot of the differentially expressed mRNAs were displayed using the R package pheatmap. Figure 5 As shown.

[0091] Different subtypes

[0092] Construct a venn graph from the two sets of DEGs (DEGs1 and DEGs2). Figure 6 The intersection was then taken, and 72 intersection genes were obtained for further analysis.

[0093] Construction of prognostic model

[0094] Regression analysis

[0095] Prognosis refers to the estimation of the possible outcome after treatment based on the current condition of the patient. The purpose of studying prognosis is to understand the rules of disease development, find methods to break the prognosis information early, create and use effective treatment methods, master the braking power of diagnosis and treatment, intervene in adverse natural prognosis, improve adverse treatment prognosis, improve medical level, and improve quality of life. Prognosis model is mainly used to inform individuals about the future development of their disease (or the risk of developing a specific disease), and to guide doctors and patients in making joint decisions in further treatment.

[0096] First, single factor Cox regression was performed on the intersection genes identified by the training set TCGA. The main purpose was to perform Cox regression analysis on individual genes to see if they were significantly related to sample survival, thereby screening genes. As shown in Figure 7 Based on P<0.01 (P value reference: PMID: 35840928), 6 genes were selected, including KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8.

[0097] Second, LASSO regression analysis was performed to obtain a more refined model by constructing a penalty function, which compresses some coefficients while setting some coefficients to zero. Therefore, it retains the advantages of subset shrinkage and is a biased estimate for handling multicollinearity data. It can select variables while estimating parameters and better solve the problem of multicollinearity in regression analysis. The model reaches the optimal when lambda is the smallest. As shown in Figure 8 The smallest partial likelihood bias (lamda.min=0.00972) was used to obtain the best lambda value, and a total of 6 genes were obtained, including KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8.

[0098] Evaluation of risk score model

[0099] Next, a risk score model was constructed, and its risk score formula: Riskscore = 0.03210648 x EXP(KRT20) + 0.02462363 x EXP(MMP12) + 0.02566009 x EXP(AGR2) + 0.02882947 x EXP(CXCL5) - 0.02414141 x EXP(CYP3A4) + 0.06227324 x EXP(MAGEA8) was used to calculate the risk scores of all samples in the data set. According to the median (0.1427), the patients were divided into a high-risk group and a low-risk group. KM survival analysis was performed on the high and low risk groups, and it was shown from Figure 9 that the survival rate of patients in the low risk score group was significantly higher than that in the high risk group. And it was shown from Figure 10 that the more the number of survival and death patients in the training set, the higher the risk score in Figure 10 A.

[0100] At the same time, the ROC curve (Receiver Operating Characteristic Curve) was used to evaluate the prognosis model, and the area under the curve was called AUC (Area Under Curve) to predict the accuracy. The larger the area under the curve, the higher the prediction accuracy. The areas under the ROC curves of the prognosis signature at 1 year, 2 years and 3 years were 0.704, 0.656 and 0.661 respectively Figure 11 , showing a certain degree of prediction ability.

[0101] Validation of the prognosis model

[0102] In order to verify the prognosis model of the screened genes as a prognosis signature, the validation set LIRI-JP was used to evaluate the results of the training set, and the risk score was calculated using the coefficients of the above genes. According to the median (-0.0578), the patients were divided into a high-risk group and a risk group. In the case of P < 0.05, the KM curve showed that the survival rate of patients in the low risk score group was significantly higher than that in the high risk group Figure 12 . It was shown from Figure 13 that the more the number of survival and death patients in the validation set, the higher the risk score. In the validation set, the areas under the ROC curves of the prognosis signature at 1 year, 2 years and 3 years were 0.764, 0.689 and 0.683 respectively Figure 14 , showing a certain degree of prediction ability.

[0103] Correlation analysis between risk score and clinical characteristics

[0104] To obtain the correlation between risk score and different clinicopathological features (T (T1, T2, T3, T4), N (N0, N1, N2), M (M0, M1), tumor_stage.diagnoses (Stage I, Stage II, Stage III, Stage IV), Grade and gender, etc.) in the training set.

[0105] To show the relationship between clinical features and risk groups, correlation test was used, Wilcoxon test for two groups (comparison between two groups, non-parametric test), Kruskal-Wallis for multiple groups (comparison between multiple groups, non-parametric test), threshold: P < 0.05, as shown in Figure 15 , the results showed that the risk score was significantly different between different clinical feature groups (tumor_stage.diagnoses, status, Grade, T and N). At the same time, the significant clinical features were refined, the risk score of each subgroup was analyzed, and the KM survival curve was drawn, and the results of the refined clinical features were shown as Figure 16 .

[0106] Independent prognostic analysis

[0107] First, in order to evaluate the ability of the risk model as an independent prognostic factor of the disease, the risk score and clinical features (M, N, T, tumor stage, etc.) were included in the single factor and multiple factor Cox regression analysis to screen independent prognostic factors, and the screening threshold was P < 0.05, as shown in Figure 17 , the risk score can be used as an independent prognostic factor and can be used as a high-risk factor.

[0108] Second, in order to establish a more reliable clinical practice prediction method, according to the results of single factor Cox regression and multiple factor Cox regression, a composite nomogram was constructed based on risk score and grade to predict the survival rate of patients at different times Figure 18 (C index = 0.701). The nomogram (Alignment Diagram) is based on multiple factor regression analysis, integrates the prediction indicators, and uses a scaled line segment to draw according to a certain proportion, so as to express the mutual relationship between each variable in the prediction model. Different survival periods show that the prediction results of the gene prognosis nomogram are consistent with the actual results to a certain extent Figure 19 .

[0109] Functional enrichment analysis

[0110] Gene Set Enrichment Analysis (GSEA) is a computational method to determine whether a set of a priori defined genes shows statistically significant and consistent differences between two biological states.

[0111] To understand the biological function and involved signaling pathways of prognostic genes, the correlation coefficient between all genes and the expression of prognostic genes was calculated, and the correlation coefficient was sorted for GSEA enrichment analysis. The screening condition was |NES| > 1, p < 0.05 and FDR < 0.25

[0112] KRT20 was mainly enriched into 19 KEGG entries and 885 GO entries, such as Cell cycle, Spliceosome, Primary bile acid biosynthesis, DNA replication, Lysosome and the like KEGG entries and mitotic sister chromatid segregation, mitotic sister chromatid segregation and the like GO entries; MMP12 was mainly enriched into 68 KEGG entries and 1774 GO entries, such as Leishmania infection, Hematopoietic cell lineage and the like KEGG entries and mitotic sister chromatid segregation, cell cycle checkpoint and the like GO entries; MAGEA8 was mainly enriched into 13 KEGG entries and 855 GO entries, such as Base excision repair, Cell cycle and the like KEGG entries and DNA damage checkpoint, G2 / M transition of mitotic cell cycle and the like GO entries; AGR2 was mainly enriched into 32 KEGG entries and 681 GO entries, such as Glycine, serine and threonine metabolism, Tryptophan metabolism and the like KEGG entries and steroid hydroxylase activity, drug metabolic process and the like GO entries; CXCL5 was mainly enriched into 62 KEGG entries and 1673 GO entries, such as Cytokine-cytokine receptor interaction, Butanoate metabolism and the like KEGG entries and adaptive immune response, immune response-activating cell surface receptor signaling pathway and the like GO entries.

[0113] Immune analysis between high and low risk groups

[0114] The infiltration of 28 kinds of immune cells in the training set high and low risk groups was analyzed using the ssGSEA algorithm, and the results were as follows Figure 21As shown, 19 kinds of cells such as Activated CD4 T cell and Effector memeory CD4 T cell were found to be significantly different between high and low risk groups, and the infiltration was higher in the high risk group.

[0115] Secondly, Spearman was used to calculate the correlation between risk score and differential immune cells, and the results were as follows Figure 22 As shown, in the case of P<0.05 and the absolute value of correlation greater than 0.3, there were 7 kinds of immune cells that were significantly positively correlated with risk score, which were Effector memeory CD4 T cell, Plasmacytoid dendritic cell, MDSC, Activated dendritic cell, Central memory CD4 T cell, Natural killer T cell and Activated CD4 T cell.

[0116] As shown Figure 23 , the expression of (CD274 (PD-L1), CTLA-4, LAG-3, LGALS9 (GAL9), HAVCR2 (TIM-3), PDCD1 (PD-1), PDCD1LG2 (PD-1LG2), TIGHT and B7H3 (reference: PMID: 32887635; PMID: 34093817) in the training set was analyzed in high and low risk groups, and it was found that all immune checkpoints were significantly different between high and low risk groups. The correlation between the differentially expressed immune checkpoints in high and low risk groups and the risk score was analyzed, as shown in Figure 24 , in the case of P<0.05 and the absolute value of correlation greater than 0.3, there were 4 pairs of genes and risk score, which were CTLA4, HAVCR2, LGALS9 and PDCD1, and there were 49 pairs of immune checkpoints that met the conditions.

[0117] Finally, the sensitivity of high and low risk groups in the training set to immunotherapy was evaluated using TIDE algorithm score, and it was found that TIDE score was significantly different between high and low risk groups ( Figure 25 ), combined with reference PMID: 32607061, by two-sided Fisher's exact test, when P was less than 0.05, it was found that the number of patients with immunotherapy response in high and low risk groups was significantly different, and the risk score was significantly distributed between respond and non-respond ( Figure 26 ).

[0118] Drug sensitivity analysis between high and low risk groups

[0119] The 289 drugs of LIHC were obtained through the Cancer Drug Sensitivity Genomics (GDSC) (https: / / www.cancerRxgene.org / ) website, and the drug sensitivity data and standardized expression data of cell lines were obtained, wherein the expression data contained 286 drugs, reference: PMID: 35198564, and oncoPredict was used to calculate the IC50 of common chemotherapy and molecular targeted drugs in each LIHC sample in the training set TCGA. The IC50 difference between the high and low risk groups was compared by Wilcoxon test, and there were 153 drugs with significant differences between the high and low risk groups at P less than 0.05, including Sorafenib, Entospletinib, etc., as shown in Figure 27 .

[0120] Expression verification of prognostic genes The expression levels of prognostic genes in TCGA-LIHC tissue samples and normal tissue samples were calculated, as shown in Figure 28 , and at P less than 0.05, 5 prognostic genes KRT20, MMP12, CXCL5, CYP3A4 and MAGEA8 were significantly expressed between the control and disease samples, and the expression trends of all genes were consistent, wherein KRT20, MAGEA8 and MMP12 were expressed higher in the disease group, and AGR2, CXCL5 and CYP3A4 were expressed higher in the normal group.

[0122] In summary, the present application establishes a prognostic model based on (Immunogenic cell death related genes, IRGs) IRGs, and verifies the relationship between the model and the immune microenvironment, and the model for predicting the prognosis of patients with hepatocellular carcinoma (Hepatocellular carcinoma, HCC).

[0123] Firstly, consensus clustering analysis was applied to TCGA-LIHC, and LIHC samples were divided into different clusters based on IRGs. Differential expression analysis was used to screen out DEGs between LIHC and normal samples, and DEGs between clusters in TCGA-LIHC, and the intersection of them was obtained as IRGs-DEGs. Then, univariate Cox and LASSO analysis were performed on IRGs-DEGs to identify prognosis-related genes, and HCC samples were divided into high / low risk groups according to risk score. Cox proportional hazards model was used for prognosis analysis. Finally, the role of prognosis model in liver cancer immune microenvironment and chemotherapy was analyzed. Liver cancer samples could be divided into two clusters, and the survival status of cluster 2 was worse. 72 IRGs were obtained by taking the intersection of 2197 DEGs and 112 DEGs in TCGA-LIHC. Six prognosis-related genes were finally determined by univariate Cox and LASSO analysis, which were KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8. A liver cancer prognosis model related to immune, nitrogen metabolism, adaptive immune response and other metabolic pathways was constructed. In addition, TIDE (Tumor Immune Dysfunction and Exclusion) score was significantly lower in low-risk group, suggesting that patients in low-risk group were more likely to receive immunotherapy. High-risk group was sensitive to camptothecin, sorafenib, etc., while low-risk group was sensitive to veliparib, darafinib. Through bioinformatics analysis, a liver cancer prognosis model was constructed, which was helpful for the related research of liver cancer prognosis and treatment.

[0124] The preferred embodiments of the application disclosed above are only used to help explain the application. The preferred embodiments do not describe all the details, nor limit the application to the specific embodiments described. Obviously, many modifications and variations can be made according to the content of the present application. The present application selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present application, so that those skilled in the art can well understand and utilize the present application. The present application is limited only by the claims and their full scope and equivalents.

Claims

1. A method for constructing a prognostic model for hepatocellular carcinoma based on immunogenic cell death-related genes, characterized in that, Includes the following steps: S1: Obtain gene expression data of hepatocellular carcinoma tissue samples and normal tissues from the TCGA-LIHC database as the training set, and obtain gene expression data of hepatocellular carcinoma tissue samples from the LIRI-JP project in the ICGC database as the validation set. S2: Consensus clustering analysis was performed on the training set samples based on immunogenic cell death-related genes. The optimal number of clusters was determined according to the cumulative distribution function. The samples were classified as subtype 1 and subtype 2. S3: Use the DEseq2 package to screen for differentially expressed genes DEGs1 between HCC tissues and normal tissues, with a threshold of |log2FC|>1.5 and adjust.P<0.05; screen for differentially expressed genes DEGs2 between subtype 1 and subtype 2, with a threshold of |log2FC|>1.5 and adj.P<0.05; take the intersection of DEGs1 and DEGs2 to obtain IRGs-related differentially expressed genes. S4: Perform univariate Cox regression analysis on the IRGs-DEGs obtained in step S3; screen out significant genes and perform LASSO regression analysis to determine prognostic genes; S5: Establish a risk scoring formula: Riskscore=0.03210648×EXP(KRT20)+ 0.02462363×EXP(MMP12) + 0.02566009×EXP(AGR2)+0.02882947×EXP(CXCL5) - 0.02414141×EXP(CYP3A4)+0.06227324×EXP(MAGEA8); S6: Using the median risk score in the training set as a threshold, patients were divided into high-risk and low-risk groups. In the validation set, patients were grouped according to the median risk score, and the survival difference between the high-risk and low-risk groups was validated using Kaplan-Meier survival curves.

2. The method for constructing a hepatocellular carcinoma prognostic model based on immunogenic cell death-related genes as described in claim 1, characterized in that: The prognostic genes include: KRT20, MMP12, AGR2, CXCL5, CYP3A4, and MAGEA8.

3. The method for constructing a hepatocellular carcinoma prognostic model based on immunogenic cell death-related genes as described in claim 1, characterized in that: The immunogenic cell death-related genes (IRGs) mentioned in step S2 include ENTPD1, ​​NT5E, CALR, HMGB1, HSP90AA1, ATG5, BAX, CASP8, PDIA3, EIF2AK3, PIK3CA, CXCR3, IFNA1, IFNB1, IL10, IL6, TNF, CASP1, IL1R1, IL1B, NLRP3, P2RX7, LY96, MYD88, TLR4, CD4, CD8B, CD8A, FOXP3, IFNG, IFNGR1, IL17A, IL17RA, and PRF1, totaling 34 genes.

4. A combination of biomarkers for prognostic assessment of hepatocellular carcinoma, characterized in that: It includes the following 6 genes: KRT20, MMP12, AGR2, CXCL5, CYP3A4, and MAGEA8.

5. The biomarker combination for prognostic assessment of hepatocellular carcinoma as described in claim 4, characterized in that: The biomarker combination calculates prognostic risk by detecting gene expression levels in liver tissue samples and then substituting these levels into a risk scoring formula. RiskScore=0.03210648×EXP(KRT20)+ 0.02462363×EXP(MMP12) + 0.02566009×EXP(AGR2)+0.02882947×EXP(CXCL5) - 0.02414141×EXP(CYP3A4)+0.06227324×EXP(MAGEA8); A RiskScore > 0.1427 indicates that the patient belongs to the high-risk group; A RiskScore ≤ 0.1427 indicates that the patient belongs to the low-risk group.

6. A kit for detecting the combination of biomarkers as described in claim 4, characterized in that: include: Primer pairs or probes for specific amplification of the KRT20, MMP12, AGR2, CXCL5, CYP3A4 and MAGEA8 genes; And its reagent components used for qRT-PCR detection.

Images