Machine learning-based esophageal cancer prognosis risk analysis method, system and medium

Abstract

The application discloses an esophageal cancer prognosis risk analysis method and system based on machine learning and a medium, relates to the fields of artificial intelligence technology and bioinformatics technology, and comprises the following steps: step one, collecting patient multi-modal data and standardizing the processing, and constructing a stable feature set; step two, constructing multiple machine learning models, modeling based on the stable feature set to realize high and low risk group prediction of each machine learning model, so as to determine the machine learning model and the optimal feature set; step three, calculating the SHAP interaction value between each feature in the optimal feature set, drawing an interaction value curve based on the SHAP interaction value, obtaining a TopA interaction feature pair, and taking the interaction feature pair as a newly constructed feature to be included in the original feature set in step one, repeating the feature screening process in step one and step two, so as to obtain an optimal machine learning model for predicting high and low risk groups of patients; and the specificity of the prognosis risk analysis method for detecting high and low risk groups of patients with locally advanced esophageal cancer is greatly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of artificial intelligence technology and bioinformatics technology, in particular to an esophageal cancer prognosis risk analysis method and system based on machine learning and a medium. BACKGROUND

[0002] Esophageal cancer is a high-incidence malignant tumor, and the treatment of locally advanced patients is difficult and the prognosis is poor. Precise risk stratification and prognosis evaluation are key problems that need to be solved in clinical practice. At present, the methods relying on TNM staging, clinical pathological features and serum markers, among which TNM staging is the most commonly used system for staging tumors internationally, and is also the standard method for staging malignant tumors in clinical practice. The above methods have limitations: TNM staging is difficult to reflect the biological behavior of tumors and individual differences, clinical indicators are highly subjective, and single markers have insufficient sensitivity and specificity.

[0003] Although imageomics can extract quantitative features from CT, existing researches are mostly limited to single modality data, and the feature screening is not systematic and the model stability is insufficient. Clinical and transcriptomic data contain important information, but lack of deep mining of multi-modal correlation. Therefore, a method is needed to integrate multi-modal data, optimize features and models, and reveal macro-micro correlation, so as to realize accurate prediction and analysis. SUMMARY

[0004] Based on the technical problems existing in the background art, the present application provides an esophageal cancer prognosis risk analysis method and system based on machine learning and a medium, which greatly improves the specificity of detection of high-risk and low-risk groups of locally advanced esophageal cancer patients.

[0005] The esophageal cancer prognosis risk analysis method based on machine learning provided by the present application comprises:

[0006] Step one, collecting multi-modal data of patients and standardizing the data, extracting imageomics features from the standardized data, and determining a stable feature set after feature screening from an original feature set composed of imageomics features and clinical features, wherein the multi-modal data includes CT images and clinical features;

[0007] Step two, constructing multiple machine learning models, modeling based on the stable feature set after screening to realize high-risk group and low-risk group prediction of each machine learning model, selecting the machine learning model with the highest comprehensive evaluation score and the best feature number, and determining the optimal feature set;

[0008] Step three, calculate the SHAP interaction value between each feature in the optimal feature set and average it after absolute value processing to construct a feature interaction strength matrix, according to which an interaction value curve is drawn, and TopA interaction feature pairs are obtained, the two features in each interaction feature pair are normalized and multiplied, and then the square root is taken as a new constructed feature to enter the original feature set in step one, and the feature screening process in step one and step two are repeated, so as to obtain the optimal machine learning model for predicting the high-risk group and the low-risk group of patients, and for esophageal cancer prognosis risk analysis, wherein A is an integer.

[0009] Further, the patient multi-modal data is collected and standardized, and the standardization of the CT image is specifically:

[0010] The CT image of the patient before radiotherapy is collected, and the CT image has been outlined by a doctor.

[0011] Spatial standardization: B-spline interpolation is used to resample the CT image to a uniform pixel spacing; for the clinical target volume mask, nearest neighbor interpolation is used for synchronous resampling to ensure that the mask and the CT image match in spatial dimension;

[0012] Gray scale standardization: The gray scale value of the CT image is normalized by linear transformation to scale it to a set standard range;

[0013] Gray scale discretization: Set Q gray scale intervals to discretize the continuous gray scale value into discrete intervals;

[0014] Image transformation preprocessing: Based on the CT image after gray scale discretization, multi-modal image representation is generated by wavelet transform, local binary pattern, square transform, exponential transform, gradient transform, square root transform and logarithmic transform method, which is used as the standardized data.

[0015] Further, the imageomic features are extracted from the standardized data, specifically:

[0016] A multi-process parallel processing framework is used to construct a feature extractor based on the pyradiomics library, all feature categories are enabled, and the standardized data are subjected to feature extraction to obtain imageomic features.

[0017] Further, the original feature set composed of imageomic features and clinical features is subjected to a feature screening process to determine the screened stable feature set, wherein the feature screening process is specifically:

[0018] (a1) Remove low-variation features with a variance less than a set variance;

[0019] (a2) performing independent sample t test on characteristics of the high-risk group and the low-risk group, retaining characteristics with p value less than a preset p threshold, and screening a first candidate characteristic set with statistically significant differences between groups, wherein the high-risk group is a patient with survival time less than 6 months, and the low-risk group is a patient with survival time greater than 36 months, and the survival time is the difference between the death time and the end date of radiotherapy;

[0020] (a3) calculating Spearman rank correlation coefficients between the first candidate characteristic set, constructing a correlation matrix, removing characteristics with absolute value of correlation coefficient greater than a preset correlation threshold, and obtaining a second candidate characteristic set;

[0021] (a4) quantifying the correlation between the second candidate characteristic set and the target variable and the redundancy between characteristics based on the maximum correlation minimum redundancy principle, screening characteristics with high label correlation and low inter-feature redundancy, and obtaining a third candidate characteristic set;

[0022] (a5) screening the third candidate characteristic set using L1 regularization, optimizing the regularization parameter through cross-validation, retaining non-zero coefficient characteristics, and obtaining a fourth candidate characteristic;

[0023] (a6) repeating (a1) to (a6) to calculate the frequency of each characteristic in the screening by setting the number of resampling, and obtaining a stable feature set after screening.

[0024] Further, step two, constructing multiple machine learning models, modeling based on the stable feature set after screening to realize high-risk group and low-risk group prediction of each machine learning model, selecting the best feature number and the machine learning model with the highest comprehensive evaluation score, and determining the optimal feature set, specifically:

[0025] Each machine learning model is: random forest, extreme random tree, gradient boosting tree, gradient boosting decision tree, support vector machine, and ensemble model:

[0026] Optimizing the hyperparameters of the random forest through grid search;

[0027] Using the parameter-optimized random forest to score the importance of each feature in the stable feature set and arranging them in descending order to obtain an ordered feature set;

[0028] Starting from the first feature in the ordered feature set, traversing the subsequent features and constructing subsets in turn, obtaining a number of subsets equal to the total number of features in the ordered feature set, and each subset includes a single feature or consecutive features;

[0029] Each machine learning model traverses the feature quantity range of the corresponding subset to find the balance point of "feature quantity - model performance", and obtains a machine learning model that can stably predict the high-risk group and the low-risk group of patients and the optimal feature set.

[0030] Further, it is also used for patient prognosis evaluation, in particular:

[0031] The optimal feature set is combined with transcriptome data to mine gene expression-imageomics feature pairs, based on which the patient prognosis is analyzed, and the prognosis-related biomarkers are found.

[0032] Further, the optimal feature set is combined with transcriptome data to mine gene expression-imageomics feature pairs, in particular:

[0033] Transcriptome data of tumor tissues of patients are collected to obtain a hub gene expression dataset;

[0034] The optimal feature set is preprocessed to remove non-feature columns, and patient identifiers and imageomics feature columns are retained;

[0035] The hub gene expression dataset is combined with the preprocessed optimal feature set with the patient identifier as the association key, to form a combined dataset containing patient identifiers, gene expression data and imageomics features;

[0036] The gene expression columns and imageomics feature columns in the combined dataset are extracted as input variables for correlation analysis, a correlation calculation function is defined, each gene expression-imageomics feature pair combination is calculated, and gene expression-imageomics feature pairs meeting the conditions are selected according to the correlation threshold and the statistical significance threshold.

[0037] Further, at least three correlation calculation functions are defined and arranged side by side to calculate each gene expression-imageomics feature pair combination.

[0038] A computer system comprising a memory, a processor and a computer program stored on the memory, wherein the processor executes the computer program to implement the above method.

[0039] A computer readable storage medium, wherein a plurality of classification programs are stored on the computer readable storage medium, and the plurality of classification programs are used to be called and executed by a processor to implement the above method.

[0040] The esophageal cancer prognosis risk analysis method, system and medium based on machine learning provided by the application have the advantages that: based on multi-modal and multi-dimensional modeling analysis, the one-sidedness problem of single-index analysis is solved, the specificity of detection of the high-risk group and the low-risk group of locally advanced esophageal cancer patients is greatly improved through comparison of multiple machine learning models, and the prognosis of locally advanced esophageal cancer patients is better explained through fusion of multi-dimensional data analysis. BRIEF DESCRIPTION OF DRAWINGS

[0041] Figure 1 This is a schematic diagram of the process of the present invention;

[0042] Figure 2 A schematic diagram of the Top 15 features retained after Bootstrap resampling feature filtering;

[0043] Figure 3 The PDP image with the maximum probability is obtained by multiplying the inverse variance of the gray-level co-occurrence matrix of the wavelet-LLL layer by the gray-level co-occurrence matrix of the wavelet-HHL layer.

[0044] Figure 4 A schematic diagram illustrating the optimal number of features and the ROC of the optimal model;

[0045] Figure 5 This is a topological map showing the correlation between gene expression and radiomics features. Detailed Implementation

[0046] The technical solution of the present invention will now be described in detail through specific embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0047] like Figures 1 to 5 As shown, the machine learning-based prognostic risk analysis method for esophageal cancer proposed in this invention includes:

[0048] Step 1: Collect multimodal patient data and standardize it. Extract features from the standardized data to obtain radiomics features. Determine the stable feature set after screening by the original feature set composed of radiomics features and clinical features through a feature screening process. The multimodal data includes CT images and clinical features.

[0049] Step 2: Construct multiple machine learning models, and based on the selected stable feature set, model and implement high-risk and low-risk group predictions for each machine learning model. Select the machine learning model with the best number of features and the highest comprehensive evaluation score, and determine the optimal feature set.

[0050] Step 3: Calculate the SHAP interaction values ​​between each feature in the optimal feature set, and construct a feature interaction strength matrix after averaging and absolute processing. Plot the interaction value curve to obtain the Top A interaction feature pairs. Perform hierarchical normalization on the two features in each interaction feature pair, multiply them, and then take the square root. This result is used as a newly constructed feature and incorporated into the original feature set in Step 1. Repeat the feature selection process in Step 1 and Step 2 to obtain the optimal machine learning model, which is used to predict the high-risk and low-risk groups of patients for esophageal cancer prognostic risk analysis, where A is an integer.

[0051] The embodiment is based on multi-modal and multi-dimensional modeling analysis, solves the one-sidedness problem of single index analysis, greatly improves the specificity of detecting high-risk group and low-risk group of locally advanced esophageal cancer patients through comparison of multiple machine learning models, and has better interpretability for the prognosis of locally advanced esophageal cancer patients through fusion of multi-dimensional data analysis.

[0052] The embodiment mainly realizes risk stratification and risk assessment of prognosis of locally advanced esophageal cancer, wherein (1) the machine learning model obtained in step three is used for risk stratification, (2) the optimal feature set is analyzed in combination with transcriptome data to mine gene expression-imaging omics feature pairs, so as to analyze the prognosis of patients by mining the relationship between image and gene, and find out the prognosis-related biomarkers, thereby performing risk assessment, which will be described in detail below.

[0053] I. Risk stratification of prognosis of locally advanced esophageal cancer;

[0054] First stage: constructing a high-risk group and low-risk group prediction model of locally advanced esophageal cancer patients based on multi-modal data;

[0055] Specifically, it includes:

[0056] S1, data acquisition:

[0057] The positioning CT of the esophageal cancer patient before radiotherapy is collected, and the esophageal cancer patient receives enhanced CT scanning in a free breathing state. The scanning is completed on the Brilliance CT BigBore device produced by Philips Medical Systems Company in Cleveland, USA, and the CT simulation positioning technology is adopted. Preferably, the scanning parameter setting is set to 120 kV, the effective millisecond range is 300-350, the beam column width is 16x1.5 mm, the image resolution is 512x512, the pitch range is 0.9-1.0, and the rotation time range is 0.7-0.8 seconds. The thickness of each scanning layer is 2-5 millimeters, and the enhanced chest CT and arterial phase CT are integrated. Then, the CT image is imported into the CT simulation workstation to accurately delineate the region of interest.

[0058] On the positioning CT, the three-dimensional volume of interest of GTV and CTV is manually delineated layer by layer on the Elekta Monaco planning system V3.8.0 by a doctor with five years of experience in esophageal imaging (for example, a senior radiologist). Then, the images are reviewed by a radio-oncology doctor specializing in upper gastrointestinal malignancies. If there is any opinion difference, it needs to be agreed through discussion.

[0059] GTV (Gross Tumor Volume) and CTV (Clinical Target Volume) are core terms used to define the treatment area in radiotherapy. GTV refers to the tumor volume visible through imaging (such as CT, MRI), including the primary lesion and metastatic lymph nodes, and is the "visible target" of radiotherapy, usually corresponding to the part of the tumor with the highest density. CTV is an extension of GTV to cover possible subclinical lesions (such as microscopic infiltration, edema, etc.), ensuring that all potential malignant tissues are treated. The range is usually larger than GTV, and the specific extension amount is determined according to the type and location of the tumor (e.g., CTV for lung cancer is usually extended by 5-10 mm).

[0060] Clinical data collection relies on hospital clinical big data system to collect patient age, BMI, clinical stage, medical history, blood test indicators, tumor marker indicators; experienced clinicians follow up on patient survival, survival is the difference between death time and radiotherapy end date; record patient survival time (month), whether they have died (die: 1 indicates death, 0 indicates survival).

[0061] According to the survival period of patients after radiotherapy, patients with a survival period of less than 6 months are included in the high-risk group (label 1), and patients with a survival period of more than 36 months are included in the low-risk group (label 0);

[0062] Standardized processing of CT images strictly follows the international imaging biomarker standardization initiative (IBSI) specifications, which includes:

[0063] Spatial standardization: B-spline interpolation method is used to resample CT images to a uniform pixel spacing (e.g., 1mm x 1mm x 1mm), ensuring spatial scale consistency of different patient images; for clinical target volume (CTV) mask, nearest neighbor interpolation method is used for synchronous resampling to ensure that the mask and CT image are strictly matched in spatial dimension;

[0064] Gray scale standardization: The gray scale value of the CT image is normalized by linear transformation to scale it to a set standard range (e.g., 0-100), eliminating gray scale offset caused by different scanning devices and parameters, and enhancing feature comparability;

[0065] Gray scale discretization: Set Q (e.g., 25) gray scale intervals to discretize continuous gray scale values into discrete intervals, reduce image noise interference, and stabilize texture feature calculation, where Q is an integer, and its specific value is determined according to needs;

[0066] Image transformation preprocessing: On the basis of the CT image after gray discretization, 8 kinds of transformation methods such as wavelet transformation, local binary pattern (LBP2D / LBP3D), square transformation, exponential transformation, gradient transformation, square root transformation and logarithmic transformation are used to generate multi-modal image representation, maximize the mining of tumor image information, and serve as standardized data after processing;

[0067] S2, the original feature set composed of imaging features and clinical features is determined by the feature screening process to determine the stable feature set after screening:

[0068] A multi-process parallel processing framework is used to construct a feature extractor based on the pyradiomics library, and all feature categories (including first-order statistical features, shape features, gray level co-occurrence matrix, gray level run length matrix, etc.) are enabled. pyradiomics is an open source Python library specifically designed to extract high-throughput quantitative radiomics features from medical images (such as MRI, CT, PET, etc.).

[0069] The feature screening process is a progressive screening of "significance test → redundancy filtering → correlation optimization → sparsity compression → stability verification", and the specific steps include (a1) to (a6):

[0070] (a1) The first step removes low-variation features with a variance less than a set variance (e.g. 1e-6) (to avoid interference from non-discriminatory features), and retains the initial feature set with actual information value;

[0071] (a2) The second step performs an independent sample t-test (Welch test, without assuming equal variance) on the features of the high-risk group and the low-risk group (labels 1 / 0), and retains features with a p-value less than a preset p-threshold (e.g. threshold 0.05), screening out candidate features with statistically significant differences between groups, and excluding noise features with no significant difference between groups, as the first candidate feature set;

[0072] (a3) The third step calculates the Spearman rank correlation coefficient between the first candidate feature set, and constructs a correlation matrix; remove features with an absolute value of the correlation coefficient greater than a preset correlation threshold (threshold e.g. 0.9), eliminate highly collinear redundant features, retain a subset with higher independence between features, and obtain the second candidate feature set;

[0073] (a4) The fourth step is based on the principle of maximum relevance minimum redundancy (mRMR), and the mutual information (MIC) is used to quantify the correlation between the second candidate feature set and the target variable and the redundancy between the features. The features with high correlation with the label and low redundancy between the features are selected (Formula 1), and the third candidate feature set is obtained. The number of retained features is 1 / 10 of the sample size (at most 19), which balances the amount of feature information and the calculation efficiency;

[0074] Formula 1: , wherein S represents a feature subset in the second candidate set, , and , respectively represent the average correlation between the feature subset S and the target variable (maximization term) and the average redundancy between the features in the feature subset S (minimization term). The third candidate feature set is obtained by using a greedy strategy in the calculation process.

[0075] (a5) The fifth step further screens using L1 regularization (Lasso) (Formula 2). The regularization parameter (alpha) is optimized by 10-fold cross-validation, and the features with non-zero coefficients are retained. The sparsity of the L1 penalty is used to maximize the compression of the feature space while ensuring the prediction ability, and finally the core feature set with high discriminability is obtained.

[0076] Formula 2:

[0077] , wherein is the regularization parameter, the second term is the L1 penalty term, is the coefficient vector of the Lasso model (sparse regression model), is the number of samples (patients), is the index of the sample (patient) (from 1- ), is the label value corresponding to the th sample (patient), is the intercept term of the Lasso model, is the index of the feature, is the number of features in the third candidate feature set, is the coefficient vector corresponding to the th feature.

[0078] This embodiment is based on the mathematical property that "L1 penalty will shrink some coefficients to 0" (sparse solution property in convex optimization), and the features with non-zero coefficients are retained.

[0079] (a6) By setting the number of resampling times, steps (a1) to (a6) are repeated to obtain the frequency of each feature in the screening, and the stable feature set after screening is obtained.

[0080] ​As an implementation, the screening process of (a1) to (a6) is repeated by 500 times of Bootstrap resampling (sampling with replacement), and the frequency of each feature appearing in the screening is calculated, as shown in Figure 2 to ensure the stability and robustness of the screening results.

[0081] S3, The modeling process adopts multi-model comparison and system optimization strategy, combined with data preprocessing, feature importance sorting and cross-validation, to realize the precise prediction of high-risk group and low-risk group of patients with locally advanced esophageal cancer. The specific steps are as follows: (b1) to (b6):

[0082] (b1) Load the image features and clinical features, select the preset core feature subset (including radio features such as wavelet_LLL_glcm_InverseVariance and clinical features such as Lymphocyte Percentage), remove non-feature columns such as patient ID and invalid data rows, and ensure the validity of the input data;

[0083] (b2) Feature standardization: use the standardization tool StandardScaler to scale the features to a distribution with mean 0 and standard deviation 1, to eliminate the influence of dimension difference on the model;

[0084] (b3) Class imbalance processing: balance the sample ratio of high-risk group and low-risk group (label 0 / 1) by synthesizing minority class samples through SMOTE algorithm, to avoid model bias towards majority class, where SMOTE (Synthetic Minority Over-sampling Technique) is an oversampling algorithm for solving data class imbalance problem, which balances the dataset by synthesizing minority class samples, so as to improve the recognition ability of the classification model for minority class;

[0085] (b4) Data set division: divide the resampled data into training set (for model training) and test set (for preliminary performance evaluation) according to the ratio of 8:2, to ensure the randomness and stability of data division (the random seed is fixed as 42);

[0086] (b5) Construct 6 machine learning models;

[0087] Construct 6 classifiers, including random forest, extreme random tree, gradient boosting tree (XGBoost), gradient boosting decision tree (LightGBM), support vector machine (SVM) and integrated model (VotingClassifier, which uses soft voting to fuse the probability prediction results of random forest and extreme random tree);

[0088] The random forest is optimized by grid search (5-fold cross-validation) for hyperparameters (total number of decision trees n_estimators takes 50 / 100 / 200, maximum depth of decision tree max_depth takes None / 10 / 20).

[0089] (b6) Calculate feature importance based on the random forest method (quantify the contribution of tree node splitting);

[0090] The random forest after parameter optimization is used to score the importance of each feature in the stable feature set and arrange them in descending order, from which the ordered feature set is obtained.

[0091] (b7) Starting with the first feature in the ordered feature set, traverse the subsequent features and build subsets in turn, obtaining a number of subsets equal to the total number of features in the ordered feature set, each subset including a single feature or consecutive features;

[0092] It should be noted that, for example, there are 5 features in the ordered feature set, then the first feature is taken as the starting point to build subsets, a total of 5 subsets are built, which are: {1}, {1, 2}, {1, 2, 3}, {1, 2, 3, 4}, {1, 2, 3, 4, 5}.

[0093] (b8) Each machine learning model finds the balance point of "feature number-model performance" by traversing the feature number range of the corresponding subset, from which the machine learning model that can stably predict the high-risk group and low-risk group of patients and the optimal feature set are obtained.

[0094] It should be noted that, in the traversal of the feature number range of the subset, the random forest is specifically the random forest after parameter optimization.

[0095] Example 1: Evaluation of 6 machine learning models in the first stage;

[0096] Multi-dimensional indicators are used to evaluate the model, including AUC (area under the ROC curve, reflecting the discrimination ability), precision (precision), recall (recall), and F1 score (all using macro average, suitable for unbalanced data). The best model is selected by the comprehensive score (mean) of the four indicators to find the optimal machine learning model.

[0097] Load the dataset containing the imageomics features, survival time (month), and outcome event (die, 1 indicates the occurrence of the event, and 0 indicates the deletion).

[0098] Quantify the association between imageomics features and patient survival prognosis through Kaplan-Meier survival analysis and log-rank test, where Kaplan-Meier is a widely used non-parametric statistical method for survival analysis, mainly used to estimate the survival function and compare the survival rates between different groups.

[0099] The quantile (25%, 50%, 75%) of the feature value sets the grouping threshold, divides the patients into "high value group" (feature value > threshold) and "low value group" (feature value < threshold), and realizes the dichotomous grouping of samples;

[0100] When determining the threshold, ensure that the sample proportion of "high value group" and "low value group" is close to 1:1;

[0101] In the optimal feature set found, 12 features were significantly associated with patient survival, with p values <0.05, further supporting the robustness and stability of the feature set.

[0102] It was found that the image fusion clinical data analysis predicted the high-risk group and low-risk group of patients, and the AUC was about 7% higher than that of a single image model, and about 13% higher than that of a single clinical model.

[0103] Second stage: Based on the optimal feature set found in the first stage, analyze the "and" relationship between each feature in the optimal feature set, including:

[0104] Based on the optimal machine learning model obtained by training in the first stage, the SHAP interaction value (shap_interaction_values) between image features in the optimal feature set is calculated using shap.TreeExplainer (shap.TreeExplainer is an explainability analysis tool in SHAP library for tree models), see formula 3;

[0105] Formula 3:

[0106] ;

[0107] Among them, is all image features in the optimal feature set, is a feature subset that does not contain , is the optimal machine learning model, represents the output of only retaining features in the subset , are feature indexes, is a feature in the feature subset , is the total number of features in the feature subset , , , represent adding feature , adding feature and The predicted value at that time Features in the sample With features The interaction value.

[0108] Average and absolute value of SHAP interaction. (Formula 5), ​​and remove diagonal elements (only retain the interaction effects between features off-diagonal, and exclude the main effects of a single feature), and construct the feature interaction intensity matrix;

[0109] Formula 4: ;

[0110] in, This refers to the sample size (i.e., the number of patients). For the first Features in a sample (patient) and The interaction value.

[0111] Create an interaction heatmap: Use color gradients to display the interaction intensity of all feature pairs, intuitively presenting the global interaction pattern;

[0112] Plot the interaction value curve: Mark significant inflection points and define the threshold boundary of significant interactions. That is, the interaction threshold is derived from the interaction value curve. The selected threshold is the value when the curve has the largest rate of change - the inflection point of the curve. That is, the threshold boundary is determined to be 0.0068. Retain the top A interaction feature pairs, for example, obtain the top 5 interaction feature pairs, where A is an integer.

[0113] The five interaction feature pairs include:

[0114] Wavelet-LLL layer_grayscale dependency matrix_dependence variance (wavelet-LLL_gldm_DependenceVariance) × Lymphocyte Percentage: capturing the synergistic effect of "tumor macrostructural heterogeneity" and "body's immune status";

[0115] Lymphocyte Percentage × Wavelet-LLL Layer Gray-Level Dependence Matrix Large Dependence High Gray-Level Emphasis: Focusing on the combined effects of "immune status" and "range of highly active tumor regions";

[0116] Wavelet-LLL_gldm_LargeDependenceHighGrayLevelEmphasis x Wavelet-HHL_glcm_MaximumProbability: capture the synergistic effect of the tumor "macro-structure activity" and "local typical pattern", reflect the "global-local" structural consistency of the tumor;

[0117] Exponential_glrlm_RunEntropy x Wavelet-HHL_glcm_MaximumProbability: focus on the joint effect of the tumor "local texture disorder" and "typical structure pattern", reflect the relationship between the "chaos degree" and "dominant pattern" of the tumor microstructure;

[0118] Wavelet_LLL_glcm_InverseVariance x Wavelet_HHL_glcm_MaximumProbability: reflect the synergistic effect of the tumor "macro-heterogeneity" and "local typical pattern".

[0119] Generate dependence plot: for the Top5 interaction feature pairs, show the trend of the SHAP value of one feature with the value of another feature, visualize the influence law of the joint action of two features on the prediction result, a total of 5 dependence plots, 5 dependence plots correspond to the Top5 interaction feature pairs, the drawing is to visualize the interaction between the interaction features, one of the dependence plots is as shown in Figure 3

[0120] PDP plot of Wavelet_LLL_glcm_InverseVariance x Wavelet_HHL_glcm_MaximumProbability;

[0121] Take the 5 interaction feature pairs screened out, multiply the two features in the feature pair after hierarchical normalization, then take the square root, and then take the square root as a newly constructed feature into the first stage original feature set.

[0122] After the same feature screening and modeling process in the first stage;

[0123] ​Wavelet - LLL layer_Grey Co-matrix_Inverse Variance x Wavelet - HHL layer_Grey Co-matrix_Maximum Probability is the only reserved interaction feature pair in Bootstrap resampling (sampling with replacement) feature selection and ranks second in importance in the model (calculated by random forest algorithm).

[0124] Wavelet - LLL layer_Grey Co-matrix_Inverse Variance describes the "uniformity of global contour";

[0125] Wavelet - HHL layer_Grey Co-matrix_Maximum Probability describes the "dominant mode of local details";

[0126] When the two features interact, the model can focus on "regularity of global shape" and "typical mode of local details" at the same time, capturing collaborative information that single features cannot cover.

[0127] As shown in the results, after adding the interaction feature, the model performance AUC improved by about 2%, and Precision and Recall improved by about 3%, as shown in Figure 4 It also further improves the interpretability of the prediction of the high-risk group and the low-risk group of patients.

[0128] II. Risk assessment of prognosis of locally advanced esophageal cancer

[0129] Third stage: joint analysis of the optimal feature set and transcriptome data to mine gene expression-imageomics feature pairs, analyze patient prognosis, and find prognosis-related biomarkers. The specific process includes:

[0130] Collecting patient tumor tissue transcriptome data

[0131] Using betweenness centrality (BC betweenness centrality), closeness centrality (CC closeness centrality), degree centrality (DC degree centrality), and maximal clique centrality (MCC maximal clique centrality) to obtain Top20 pivot genes.

[0132] Load the pivot gene expression dataset (gene_df) and the optimal feature set, respectively.

[0133] Preprocess the optimal feature set, remove non-feature columns (such as the label column label), and retain patient identification (patient_id) and imageomics feature columns.

[0134] Merge the gene expression dataset and the filtered optimal feature set with patient_id as the linking key to form a merged dataset (merged_df) containing patient_id, gene expression data, and radiomics features, ensuring sample-level data alignment.

[0135] Extract the gene expression columns (excluding patient_id) and radiomics feature columns (excluding patient_id) from the merged dataset as input variables for correlation analysis.

[0136] Define correlation calculation functions to calculate the following for each gene expression-radiomics feature pair:

[0137] Pearson correlation coefficient and corresponding p-value (measuring linear correlation);

[0138] Spearman correlation coefficient and corresponding p-value (measuring monotonic correlation, suitable for non-linear relationships);

[0139] Kendall correlation coefficient and corresponding p-value (a rank-based correlation measure, more robust to outliers).

[0140] Set correlation threshold (absolute value > 0.7) and statistical significance threshold (p-value < 0.05) to filter gene expression-radiomics feature pairs that meet the conditions.

[0141] Draw a correlation heatmap of gene expression-radiomics feature pairs to visualize the correlation relationship, as shown in Figure 5 .

[0142] Find multiple pairs of gene expression-radiomics feature pairs, including:

[0143] Original - Shape Major Axis Length (ITGB2---original_shape_MajorAxisLength), gene and image positively correlated;

[0144] Wavelet - HHL Layer Gray Level Co-occurrence Matrix Maximum Probability (IL1B---wavelet-HHL_glcm_MaximumProbability), gene and image negatively correlated;

[0145] Wavelet - HHH Layer Gray Level Co-occurrence Matrix Second Order Inertia Moment (IL1B---wavelet-HHH_glcm_Imc2), gene and image negatively correlated; Three gene-image feature pairs meet three correlation calculation methods, which can be used as consensus gene expression-image feature pairs.

[0146] Draw a correlation topology map of gene expression-radiomics feature pairs, asFigure 5 ;

[0147] Based on Kaplan-Meier survival analysis and log-rank test analysis of the relationship between gene expression-imaging features and patient prognosis, it is found that:

[0148] ITGB2 low expression may be accompanied by better prognosis;

[0149] IL1B low expression is accompanied by better prognosis (consistent with clinical logic).

[0150] The above three stage results fully prove that: imaging genomics combined with clinical multi-index data can accurately realize high and low risk prediction of patients with locally advanced esophageal cancer, and the effect is better than single imaging or single clinical data; After in-depth study of the "and" relationship between image and clinical features, it is found that the interaction feature pair can explain the high and low risk stratification of patients and the prognosis of patients from a higher dimension; The joint analysis of imaging genomics and transcriptome data can explore the subtle relationship between image and gene, thereby linking the relationship between macroscopic phenomenon and microscopic nature, realizing multi-dimensional evaluation of patient prognosis. Therefore, the locally advanced esophageal cancer high risk group and low risk group prediction and prognosis analysis method based on multi-modal data of the present application, layer by layer, from macro to micro, solves the current single modality limitation, multi-modal difficulty in series, poor interpretability problem.

[0151] The above describes only the preferred specific embodiments of the present application, but the protection scope of the present application is not limited thereto, any person skilled in the art within the technical range disclosed by the present application, according to the technical scheme and inventive concept of the present application, equivalent replacement or change, should be covered within the protection scope of the present application.

Claims

1. A machine learning-based prognosis risk analysis method for esophageal cancer, characterized by, The application relates to a method for predicting esophageal cancer prognosis risk, comprising the following steps: Step 1: collecting multi-modal data of patients and performing standardization processing, extracting image features from the data after standardization processing, and determining a stable feature set after screening from an original feature set composed of the image features and clinical features through a feature screening process, wherein the multi-modal data comprises CT images and clinical features; Step 2: constructing multiple machine learning models, modeling based on the stable feature set after screening in step 1 to realize high-risk group and low-risk group prediction of each machine learning model, selecting a machine learning model with the highest comprehensive evaluation score and the best feature number, and determining an optimal feature set; Step 3: calculating SHAP interaction values between features in the optimal feature set, constructing a feature interaction intensity matrix after average absolute processing, drawing an interaction value curve diagram, obtaining TopA interaction feature pairs, multiplying two features in each interaction feature pair after hierarchical normalization, and then taking the square root to obtain a newly constructed feature, which is added to the original feature set in step 1, and the feature screening process in step 1 and step 2 are repeated, so as to obtain an optimal machine learning model for predicting high-risk group and low-risk group of patients, and esophageal cancer prognosis risk analysis is performed based on the optimal machine learning model, wherein A is an integer.

2. The analysis method according to claim 1, characterized in that, The multi-modal data of patients is collected and standardized, and the standardization processing of the CT images is specifically as follows: CT images before radiotherapy of patients are collected, and the CT images have been delineated by doctors; spatial standardization: the CT images are resampled to a uniform pixel interval by using a B-spline interpolation method; for a clinical target volume mask, a nearest neighbor interpolation method is used for synchronous resampling to ensure that the mask and the CT images match in spatial dimensions; gray scale standardization: the gray scale values of the CT images are normalized, and the gray scale values are scaled to a set standard range through linear transformation; gray scale discretization: a plurality of gray scale intervals are set, and the continuous gray scale values are discretized into discrete intervals; image transformation preprocessing: based on the CT images after gray scale discretization, a plurality of modal image representations are generated through wavelet transformation, local binary pattern, square transformation, exponential transformation, gradient transformation, square root transformation and logarithmic transformation, which are used as the data after standardization processing.

3. The analysis method of claim 1, wherein, The image features are extracted from the data after standardization processing, and the specific process is as follows: a multi-process parallel processing framework is used, a feature extractor is constructed based on a pyradiomics library, all feature categories are enabled, and the data after standardization processing are subjected to feature extraction to obtain the image features.

4. The analysis method of claim 1, wherein, The original feature set composed of the image features and the clinical features is subjected to a feature screening process to determine a stable feature set after screening, and the feature screening process is specifically as follows: (a1) removing low-variation features with a variance less than a set variance; (a2) performing independent sample t-test on the features of the high-risk group and the low-risk group, retaining features with a p value less than a preset p threshold, and screening out a first candidate feature set with statistically significant differences between groups, wherein the high-risk group is composed of patients with a survival period less than 6 months, the low-risk group is composed of patients with a survival period greater than 36 months, and the survival period is the difference between the death time and the radiotherapy end date; (a3) calculating Spearman rank correlation coefficients between the first candidate feature set, constructing a correlation matrix, removing features with absolute correlation coefficients greater than a preset correlation threshold, and obtaining a second candidate feature set; (a4) based on the principle of maximum correlation and minimum redundancy, quantifying the correlation between the second candidate feature set and the target variable and the redundancy between the features by mutual information, selecting features with high label correlation and low feature redundancy, and obtaining a third candidate feature set; (a5) using L1 regularization to screen the third candidate feature set, and optimizing the regularization parameter through cross-validation to retain non-zero coefficient features, thereby obtaining a fourth candidate feature set; (a6) by setting the number of resampling, repeating (a1) to (a6) to calculate the frequency of each feature in the screening, and obtaining a stable feature set after screening.

5. The analysis method of claim 1, wherein, Step two, constructing multiple machine learning models, modeling based on the stable feature set after screening to realize the prediction of high-risk group and low-risk group of each machine learning model, selecting the best feature number and the machine learning model with the highest comprehensive evaluation score, and determining the optimal feature set, specifically: Each machine learning model is: random forest, extreme random tree, gradient boosting tree, gradient boosting decision tree, support vector machine, and ensemble model: Optimizing the hyperparameters of the random forest through grid search; Using the parameter-optimized random forest to score the importance of each feature in the stable feature set and arranging them in descending order to obtain an ordered feature set; Starting with the first feature in the ordered feature set, traversing the subsequent features and constructing subsets in turn, obtaining a number of subsets equal to the total number of features in the ordered feature set, and each subset includes a single feature or consecutive features; Each machine learning model traverses the feature number range of the corresponding subset to find the balance point of "feature number - model performance", thereby obtaining a machine learning model that can stably predict the high-risk group and low-risk group of patients and the optimal feature set.

6. The analysis method of claim 1, wherein, It is also used for patient prognosis evaluation, specifically: Jointly analyzing the optimal feature set and transcriptome data to mine gene expression-imageomics feature pairs, thereby analyzing patient prognosis and finding prognosis-related biomarkers.

7. The analysis method according to claim 6, characterized in that, The joint analysis of the optimal feature set and the transcriptome data to mine gene expression-imageomics feature pairs is specifically: Collecting transcriptome data of patient tumor tissues to obtain a hub gene expression dataset; Pretreating the optimal feature set, removing non-feature columns, and retaining patient identifiers and imageomics feature columns; Merging the hub gene expression dataset and the pretreated optimal feature set based on the patient identifier as the association key to form a joint dataset containing patient identifiers, gene expression data, and imageomics features; Extracting the gene expression columns and imageomics feature columns in the joint dataset as input variables for correlation analysis, defining a correlation calculation function, calculating each gene expression-imageomics feature pair combination, and selecting gene expression-imageomics feature pairs that meet the conditions based on the correlation threshold and statistical significance threshold.

8. The analysis method of claim 6, wherein, Defining a correlation calculation function, and using at least three correlation calculation functions arranged side by side to calculate the correlation coefficient of each gene expression-imageomics feature pair combination.

9. A computer system comprising a memory, a processor and a computer program stored on the memory, characterized in that, The processor executes the computer program to implement the method of any one of claims 1-8.

10. A computer readable storage medium, characterized in that, The computer readable storage medium has stored thereon a plurality of classification programs for being called and executed by the processor to implement the method of any one of claims 1-8.

Images