A prognostic analysis method for non-small cell lung cancer based on the computed tomography radiomics

The data used in the experiments were derived from two data sets, which are lung 1 data set (called data set 1 below) of TCIA NSCLC-Radiomics project⁵ and lung 3 data set (called data set 2 below) of TCIA NSCLC Radiogenomics project⁶. The data inclusion-exclusion criteria is as follows: INCLUDE single tumor (In order to control the variables and ensure the influence of a single lesion on the prognosis of patients, we only selected the data of patients with one primary lesion) and EXCLUDE patients who had a cut-off survival time of less than three years and were still alive at that time (Such patients cannot be given a proper label so that they do not meet our experimental requirements).

Based on the above requirements, we screened a total of 173 eligible patient data Among them, 124 patients from data set 1 are used for model training, and 49 patients from data set 2 are used for model testing. Besides, the database also included many clinical information such as age, gender and TNM stage. The full information of the eligible data is shown in table 1.

The clinical experiments indicated that the chance of recurrence will be greatly reduced when lung cancer patients have survived for more than 3 years. If the patient's disease-free survival (DFS) was more than 5 years, it was considered to be probably recovered. These time points were important for patient treatment design and review options in clinical practice (Zeng et al 2018). Therefore, according to the postoperative survival time of patients in the data set, the survival period of 3 years was selected as grouping threshold clinically. The experimental data from data set 1 was divided into 'long survival group' (LSP) and 'short survival group' (SSP), which correspond to survival time more than 1100 d (from 1141 to 1926 d with the median time of 1479 d) and less than 900 d (from 10 to 896 d with the median time of 418 d) respectively. After grouping, lung tumors were segmented and the features were extracted and optimized at first. Then, various classifiers were applied to classify the feature optimized data, and the prediction model of NSCLC prognosis survival was constructed based on CT radiomics. At the same time, the survival-related radiomics prognostic factors were obtained. Finally, the effectiveness of radiomics prognostic factors was verified and evaluated for performance of the classification model. The process of experiments is shown in figure 1.

Zoom In Zoom Out Reset image size

Tumor segmentation is the basis for tumor radiomics feature extraction. In order to extract accurately the radiomics features from tumor regions, the tumor regions should be firstly segmented accurately from CT image sequence to prevent interference of surrounding background (Shakibapour et al 2019).

For the requirements of segmentation, a semi-automatic segmentation method was used in this paper. Firstly, tumor area sequence images were selected that contain the entire tumor through the interactive medical imaging software RadiAntViewer. Secondly, a bicubic interpolation algorithm along the axial direction was applied to generate isotropic scans in order to solve the problem of anisotropic 3D spatial resolution CT scanning. Then, the semi-automatic segmentation method was used to segment the tumor, and different segmentation schemes were adopted for tumors with different types: (1) For the solitary tumors, the gray threshold segmentation algorithm and three-dimensional region growing algorithm were used to segment; (2) For the pleural tumors, the edge of lung area was repaired based on concave hulls algorithm (Soltaninejad et al 2016) and chain code algorithm (Sun and Wang 2016). On the basis of a complete repair of lung boundary, the same segmentation method as solitary pulmonary tumor segmentation method was used; (3) For the vascular tumors, the above-described solitary pulmonary tumor segmentation method is used for coarse segmentation at first. Moreover, the method of Geodesic distance histogram (Sun et al 2014) and fuzzy C-means clustering were adopted to achieve fine segmentation, and the segmentation outcome was the intersection of the above two segments. The final results were selected from the average of several times of segmentations, which are shown in figure 2.

Zoom In Zoom Out Reset image size

The extraction of radiomics features is an important method for calculating and quantifying tumor image information. In order to reveal the correlation between CT radiomics features and NSCLC prognosis, quantitative indicators of tumors must be extracted, as many as possible, to avoid the loss of information and increase the accuracy of predictions (Way et al 2006). In this paper, a total of 258 radiomics features based on gray, shape, texture and other information were extracted (Han et al 2015, Dhara et al 2016, de Carvalho Filho et al 2017, Ferreira et al 2017). In addition to traditional gray features, shape features, texture features, we also used mathematical modeling methods to accurately quantify medical signs commonly used in clinical practice. The features extracted by the experiment are shown in table 2.

A large number of studies have shown that the parameters of imaging equipment affect radiomics feature values. Since the data used in this experiment is from public data sets, we cannot obtain the complete imaging protocol, but we try to consider all factors that can be obtained. In the process of feature extraction, we comprehensively considered these parameters which are easy to affect the robustness and reproducibility of radiomics features, including pixel size, layer thickness, layer spacing. We obtained this information in the Dicom file. In the process of analyzing the Dicom file information of the experimental data, we found that the device information of the CT image sequence of the data set 1 is consistent, and the type of imaging device is SIEMENS CT VA1 DUMMY. So during the training of the model, we can ensure that the CT imaging device parameters have the same impact on the imaging feature extraction for each patient. In the process of model verification using data set 2, we also considered these parameters, and tried to ensure that the radiomics feature extraction of each sample data is performed under the same operation. Finally, the feature data is normalized, and value interval of all features is scaled to the range of [0, 1].

In order to maximize the collection of information on all aspects of the tumor, a large number of radiomics features describing the characteristics of the tumor need to be extracted (Lohrmann et al 2018). However, too high dimension of features often tend to make features matching too complicated and affect operation speed. At the same time, some of the extracted features will bring irrelevant information or error information to the prediction. Therefore, in order to obtain a good performance prediction model, it is necessary to filter redundant features data and retain effective features (Lee et al 2008 Jaffar and Al Eisa 2015). So the Relief feature weighting algorithm (Hancer et al 2018) was used to retain features which are highly relevant to categories to reduce redundancy features and improve performance of the prediction model. The pseudo-code for the Relief algorithm is shown as follows.

Algorithm:. Pseudo-code for the Relief algorithm

INPUT: Sample set S, Sampling number m, Feature weight threshold τ.

OUTPUT: the vector w of estimations of the qualities of attributes and the Selected feature subset.

Begin:

1.set all weights $W\left( i \right)=0$,

2.for i  =  1 to m

Randomly select a sample R from S,

Find the k-nearest neighbors ${{H}_{j}}\left(\,j=1,2,L,k \right)$ of R from the same sample set of R,

Find k-nearest neighbors from each different class of sample set ${{M}_{j}}(C)$,

3.for i  =  1 to N

$$\begin{align} \newcommand{\e}{{\rm e}} \begin{array}{@{}llcccccccccccccccccc@{}} & W(A)=W(A)-\sum\limits_{j-1}^{k}{dif\,f(A,R,{{M}_{j}}(C))/(mk)} \\ & \quad \quad \quad \quad +\sum\nolimits_{C\notin class(R)}{[\frac{p(C)}{1-p(class(R))}\sum\nolimits_{j-1}^{k}{dif\,f(A,R,{{M}_{j}}(C))/(mk)}]} \\ \end{array}\end{align} \tag{ 1 }$$

4. Sort W and sort its corresponding features.

End

${{M}_{j}}(C)$ represents the jth nearest neighbor sample in class C.

$dif\,f(A,{{R}_{1}},{{R}_{2}}))$ represents the difference between the sample R1 and the sample R2 in the feature A, as shown below: $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm dif\,\!f}(A,R1,R2)=\left\{\begin{array}{@{}lll@{}} \frac{\left| {{R}_{1}}[A]-{{R}_{2}}[A] \right|}{\max (A)-\min (A)}, & {\rm if}~A~ {\rm is}~{\rm continuous} \nonumber \\ 0 & {\rm if} ~ A~{\rm is}~{\rm discrete}~{\rm and}~{{R}_{1}}[A]={{R}_{2}}[A] \nonumber \\ 1, \quad \quad \quad \quad \quad & {\rm if}~A~{\rm is}~{\rm discrete}~{\rm and}~{{R}_{1}}[A]\ne {{R}_{2}}[A] \nonumber \end{array} \right..\nonumber \end{align} \tag{ 2 }$$

The feature screening was introduced into every cross-validation process to avoid errors in feature screening. According to the quantitative relationship between samples and features, the top ten features with the highest weighted value were selected for the training and testing by using the Relief feature weighting algorithm. The features sorted by weight value from high to low is shown in table 3.

In 124 case samples, the huge quantity difference between 'short survival group' and 'long survival group' results in data imbalance (sample ratio: 94/30  =  3.13). Data imbalance will cause errors in classifier training and affect classifier performance. Therefore, synthetic minority over-sampling technique (SMOTE) resampling method (Last et al 2017) was used to resample the feature data of 'long survival group' samples in the data set to balance the two types of samples.

SMOTE algorithm is a data synthesis algorithm based on interpolation to synthesize a small number of new samples. The implementation of the algorithm was shown as follows:

• (1)  
  For each sample x in a few classes, distance to all samples in a small sample set were calculated by Euclidean distance and the k-nearest neighbor should be obtained.
• (2)  
  A resampling magnification N was reset according to the sample imbalance ratio. For each minority sample x, several samples from its k nearest neighbors were randomly selected assuming the selected neighbors as ${{x}_{n}}$.
• (3)  
  For each randomly selected neighbor ${{x}_{n}}$, a new sample ${{x}_{new}}$ was constructed by the primitive sample according to equation (3).

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{x}_{new}}=x+{\rm rand}(0,1)\times (x-{{x}_{n}}).\nonumber \end{align} \tag{ 3 }$$

  Where: ${\rm rand}(0,1)$indicates a random number between 0 and 1.

By using SMOTE data resampling method, the number of 'long survival group' has increased from 30 to 90. Experimental data has been contained 94 samples of 'short survival group' and 90 samples of 'long survival group'. During the training and testing, SMOTE resampling method was also involved into verification process. Since the added 'long survival group' sample from resampling technique is not the real data, the data of this part only were used for classifier training but have not been used to test. In the process of classifier training and testing with all data sample, the classification results of only 124 real samples were obtained in the original data set. This method is widely used in the design of computer aided diagnosis (CAD) of pulmonary nodules (Aghaei et al 2016, Yan et al 2016).

The classification model of lung cancer patientsl survival time range was trained by the feature data with strong characterization ability of lung cancer regions obtained from feature selection, and model predictive performance was tested on all 124 cases by cross-validation. In order to construct a performance-optimized classification model as much as possible, seven representative classifiers with strong generalization ability and often used for small sample data set training have been selected (Messay and Rogers 2010, Cascio et al 2012, Wang et al 2013). They are decision trees (DT), discriminant analysis classifiers (DAC), logistic regression (LR), support vector machine (SVM), k-nearest neighbor (KNN), ensemble classifiers (EC), random forest (RF). In this paper, all the above biomedical data classification methods had been used to select the classification model with best performance for the experimental results, and the optimal performance model was used to predict survival of patients.

In order to explore the degree of correlation between the features obtained by the above feature selection algorithm and prognosis survival, and to select statistically significant correlation features as factors of prognosis survival analysis model, Pearson correlation analysis (Grkovski et al 2016, Bailey et al 2018) was applied to calculate the correlation between the top ten features obtained by the Relief algorithm and the prognostic survival labels. Correlation coefficient |r| and P-value were used to evaluate the 'significant degree' of correlation degree.

According to the prediction results obtained by the above performance optimal prediction model, the Kaplan-Meier survival analysis method was used to perform survival analysis on the predicted two groups of cases to verify the validity of the experiment (Wu et al 2016b).

The kaplan–meier survival analysis method used the principle of conditional probability and probability multiplication to calculate the survival rate by combining the cut-off survival time and survival status in the case information. In SPSS software, we inputed patients' cut-off survival time, patients' cut-off survival status, original grouping label and prognosis prediction results of patients into Kaplan Meier survival analysis model to draw the survival curve of NSCLC patients.