Computer-aided differentiates benign from malignant IPMN and MCN with a novel feature selection algorithm

1.   Introduction

In recent years, due to a higher morbidity and the malignant potential of intraductal papillary mucinous neoplasm (IPMN) and mucinous cystic neoplasm (MCN), these two pancreatic cystic neoplasms (PCNs) have drawn great attention of clinical researchers ^[1,2,3]. As there is a tremendous difference in the prognosis between benign and malignant cases, recognizing malignant IPMN/MCN preoperatively is very important for tailoring an optimal strategy of treatment for a patient. Patients with benign tumors can receive a function-preserving surgery or even can escape surgery and stick to a close surveillance plan, while patients with malignant tumors may undergo a more radical surgery. However, even for experienced clinicians, it's currently always challenging to distinguish benign lesions from malignant ones, because they show similar imaging appearance and clinical indices and the potential difference cannot be easily recognized by human in routine clinical practice.

Multi-detector row computed tomography (MDCT) is recommended for the initial diagnosis and assessment of IPMN and MCN ^[4,5]. And the incidence of pancreatic cancer raised dramatically, owing to the wide application of MDCT examination ^[6]. In clinical practice, there are a lot of high-risk imaging findings associated with cancerization. According to the international consensus guidelines published by Tanaka M et al. ^[7,8], the diffuse dilation of main pancreatic duct (MPD) > 5 mm is an important signal of malignant IPMN and MCN, while the large cyst > 3 cm, the thickened cystic wall, and the papillary mural nodules are also prominent features of malignant tumors. In addition to radiologic approaches, studies also have shown that certain laboratory results and clinical symptoms can also indicate the malignant transformation of pancreatic neoplasms, e.g., the raising of serum carbohydrate antigen (CA) 19-9, CA12-5 or carcinoembryonic antigen (CEA), the increasing of fasting plasma glucose, and even abdominal pain ^[9]. However, due to the similar imaging appearances and the subjective clinical diagnosis, there is a solid need from clinicians for an assistant system to help them make a more accurate and objective diagnosis on benign or malignant IPMN/MCN. As a result, some studies have emerged on computer-aided diagnosis (CAD).

Radiomics-based CAD method constructs an objective image evaluation system and the whole process is highly automated, time-saving and effective. At present, there are some CAD methods and radiomics researches related to pancreatic disease. Jayasree et al. ^[10] predicted cancer risk in branch duct (BD)-IPMN through the extraction of 135 features, then applied wilcoxon rank-sum test (WRST) based p-value feature selection and obtained the AUC of 0.77. Park et al. ^[11] extracted 431 3D CT radiomics features and used minimum-redundancy maximum-relevancy (mRMR) algorithm to reduce feature dimension, then realized the differential diagnosis of autoimmune pancreatitis and pancreatic ductal adenocarcinoma. Zhang et al. ^[12] extracted 251 radiomics features from 2D and 3D PET/CT images, and used support vector machine recursive feature elimination (SVM-RFE) to select effective features, while support vector machine (SVM) was choose as classifier to discriminate autoimmune pancreatitis (AIP) and pancreatic ductal adenocarcinoma (PDAC). Wei et al. ^[13] extracted 409 MDCT radiomics features and used bootstrapping repeated least absolute shrinkage selection operator (LASSO) regression to select effective features. Finally, they distinguished between serous cystic neoplasm (SCN) and other subtypes of PCN, which obtained the AUC of 0.84.

Throughout the abovementioned CAD methods, there are some urgent problems to be solved. First, the above feature extraction modules did not form a complete and comprehensive feature system optimized for PCN. Second, most of the feature selection methods only selecting once, which leads to the instability and contingency of the results. Third, the CAD methods have not made many attempts focusing on the task of the differentiation between benign and malignant of both IPMN and MCN. Therefore, a robust and improved radiomics CAD system is highly needed to solve these problems.

In this paper, an enhanced CAD system was proposed, which combines the radiomics features with the clinical indices and forms a complete, comprehensive and PCN-optimized feature system. In addition, a fashion and high-performance feature selection algorithm named Bootstrapping repeated LASSO with Random selections (BLR) was used in our CAD system. Finally, the features selected by BLR were sent to SVM to classify benign and malignant of IPMN and MCN.

The workflow of the proposed CAD system is shown in Figure 1. Our major contributions can be summarized as follows:

1) The feature extraction modules in the existing CAD system are not systematic, nor optimized for PCN. Besides, it lacks of some important clinical symptoms or laboratory indicators. Therefore, a comprehensively, completely, and PCN-optimized feature system has been built including 436 radiomics feature and 9 clinical indices to describe the image characteristics of PCNs automatically and accurately.

2) A novel algorithm named BLR was proposed, which increases the stability of feature selection through a large number of bootstrapping repetitions while reduces the bias and overfitting, and maximizes the diagnostic potential of features by means of random selections.

3) A robust and precise preoperative CAD system was established aimed to classify the benign and malignant of IPMN/MCN, which forms a reliable and objective image evaluation and diagnostic method.

2.   Materials and methods

2.1.   Dataset and segmentation

This study is a retrospective study. The 107 patients were enrolled at Department of Pancreatic Surgery, Huashan Hospital of Fudan University, Shanghai, China, from December 2007 to August 2016. There are 73 IPMN and 34 MCN cases, including 71 benign and 36 malignant. All patients have signed informed consent, and our study has been approved by the ethics committee of the Huashan Hospital. The characteristics of enrolled patients are shown in Table 1.

All patients underwent preoperative abdominal enhance MDCT scans with the thickness of 1.5 mm and already received surgery, the postoperative pathological examination result was taken as the ground-truth and final diagnosis of the patient. Due to the higher image quality of portal venous phase, the single 2D image with the largest tumor cross section in portal venous phase of abdominal enhance MDCT sequence was chosen as the input. The portal venous phase is the images collected in 60–65 s after intravenous injection of nonionic-iodinated contrast agent with the reagents concentration of 370 mg I/mL and the injection rate of 4 mL/s. Preoperative data used in this study of these patients is complete and available.

The size of all images was 512 × 512 and the tumors were outlined by experienced radiologists manually. The 90 cases were randomly selected as the training cohort with 5-fold cross-validation, while the remaining 17 cases formed the independent testing cohort. The data processing was implemented in Matlab R2018b (Mathworks, Inc, Natick, Massachusetts).

2.2.   High-throughput feature design and extraction

A PCN-optimized feature system was designed according to the clinical guidelines published by Sahani et al. ^[14] in 2013, which practically covers all the features included in the clinical guidelines.

Specifically, 436 radiomics features were designed and extracted, including 21 structure features, 16 intensity features, 67 texture features, and 332 wavelet features. Besides, 9 clinical indices were added to our feature system to promote the effects of radiomics-based CAD scheme and make it closer to the real clinical diagnosis model. The extracted features are summarized in Table 2, while the detailed features are shown in Appendix. The correspondences between the extracted features and the clinical guidelines features are shown in Table 3.

Structure features are used for describe the shape, edge, size and internal structure of the tumor ^[13,15]. Shape features comprehensively reflect the shape information of tumors, such as tumor size, shape, edge roughness and so on. The internal structure features focus on PCN such as internal cyst, calcification, wall thickness and central scar.

Intensity features reflect the intensity and histogram information of the tumors.

Texture features calculate the texture information of tumor region, including Inner regional comparison features and Other texture features.

2.2.1.   Inner regional comparison features

To better describe the gray level and texture distribution disparity between the different inner regions of the tumor, 6 texture features were designed, named inner regional comparison features. The greatest advantage of inner regional comparison features is that it can calculate the inner-regional disparity of the tumor with any shape, in the meanwhile, without the serious interference of the other tissues outside of the mass or compared areas. Figure 2 reflects the process of obtaining the internal regions of tumor. Figure 2(a) is the original image with the size of R × C. The R and C are both 512 in this study; (b) is the mask of the tumor; (c) is the tumor region; (d) and (e) represent ${ROI1}_{ij}$ and ${ROI2}_{ij}$, which are the upper left corner and the lower right corner of (c) respectively and with the same size of $R-i$ × $C-j$; (f) shows the intersection of ${ROI1}_{ij}$ and ${ROI2}_{ij}$; (g) is the mask of intersection region; (h) and (i) are ${Region1}_{ij}$ and ${Region2}_{ij}$, which are the point multiplication of (g) and (d) or (e), respectively. In simple terms, ${Region1}_{ij}$ and ${Region2}_{ij}$ are the different inner regions of the tumor to be compared, all of the other subpictures are the intermediate processes to obtain ${Region1}_{ij}$ and ${Region2}_{ij}$. Notice that the size of ${Region1}_{ij}$ and ${Region2}_{ij}$ vary with the value of $i$ and $j$.

So, the calculation method and formula of each feature are as follows:

1) Mean of inner-regional dissimilarity ($Dis{s}_{m}$):

where the $sum$ means to sum all pixels of an image; N is the number of pixels of foreground in (g); M and K are the positive integers, which can be regarded as the maximum offset when using Figure 2(c) to get (d) and (e). Their values should be appropriate, not too large or too small. In this study, M and K are both equal to 5. The dissimilarity of different inner regions of the tumor is proportional to the value of this feature.

2) Standard deviation of inner-regional dissimilarity ($Dis{s}_{std}$):

This feature reflects the stability of inner-regional dissimilarity.

3) Mean of inner-regional contrast (${Con}_{m}$):

where the $\mathrm{*}$ represents point multiplication of two images; $BW$ and $Tumor$ represent the Figure 2(b), (c), respectively. $min$/$max$ are the minimum and maximum value of the tumor region in (c); $Tumo{r}_{n}$ is the normalized image of the tumor region in (c). ${ROI1\_n}_{ij}$ and ${ROI2\_n}_{ij}$ are the upper left and the lower right corner of $Tumo{r}_{n}$ respectively, and with the same size of $R-i$ × $C-j$. This feature reflects the contrast inside the tumor.

4) Standard deviation of inner-regional contrast (${Con}_{std}$):

This feature reflects the stability of inner-regional contrast.

5) Mean of inner-regional covariance (${Cova}_{m}$):

P is the number of pixels with the value of 1 in the Figure 2(b); $\alpha$ represents the Figure 2(g), which is the mask of intersection region of ${ROI1\_n}_{ij}$ and ${ROI2\_n}_{ij}$. This feature calculates the covariance of the tumor's different inner-regions.

6) Standard deviation of inner-regional covariance (${Cova}_{std}$):

This feature reflects the stability of inner-regional covariance.

2.2.2.   Other texture features

ROI-based features calculate the relevant texture characteristics in the bounding-box of the tumor ^[13,16]. Besides, gray-level co-occurrence matrix (GLCM) ^[17], gray-level run-length matrix (GLRLM) ^[18], gray-level size zone matrix (GLSZM) ^[19], and neighborhood gray-tone difference matrix (NGTDM) ^[20] are common high-order matrices to describe image texture.

In order to get the deeper intensity and texture information of the image, wavelet transform was applied to all two-dimensional images and four components were obtained, which are Low pass/Low pass (LL), Low pass/High pass (LH), High pass/Low pass (HL), and High pass/High pass (HH). The 83 intensity and texture features were extracted from four component respectively, and obtained 332 wavelet features in total.

Clinical indices record the patients' demographic information, preoperative laboratory indicators, and clinical symptoms comprehensively. The quantization methods of each feature are described below:

1) Sex: 0 and 1 represent male and female, respectively.

2) Age: 1–4 represent age range from 1 to 19, 20 to 39, 40 to 59 and above 60.

3) Preoperative fasting plasma glucose (mmol/L): range from 4.1 to 11.3.

4) Serum cancer antigen 199 (CA199, u/mL): range from 0.6 to 999.

5) Serum cancer antigen 125 (CA125, u/mL): range from 2.0 to 116.9.

6) Carcinoembryonic antigen (CEA, ug/mL): range from 0.2 to 6.1.

7) tumor size (cm): range from 0.2 to 15.

8) Location: 1, 1.5, 2, 2.5, 3, 3.5, 4, in which 1, 2, 3, 4 represent the head, neck, body, tail of pancreas, respectively. While 1.5 represents between the head and the neck of the pancreas, 2.5, 3.5 are similar.

9) Clinical symptoms: 0–6 represent asymptomatic, abdominal pain, abdominal distension, weight loss, jaundice, mass, and pancreatitis, respectively.

To sum up, a complete, comprehensive and PCN-optimized feature system with 445 high-throughput features has been established in this study.

2.3.   BLR feature selection algorithm

LASSO regression has been widely used in the feature extraction because of its efficiency and high performance ^[21,22]. It can shrink the coefficients of the redundant features to zero. But selecting features only once can be extremely occasional and cause overfitting. Therefore, in our previous works ^[13,16], LASSO regression was combined with a large number of bootstrapping repetitions to increase the reliability of feature selection meanwhile prevent overfitting. To be exact, all of the input features were ranked according to their occurrence number in the multiple bootstrapping repeated LASSO regression, and the features with top 10% reproducibility were selected as the final feature subset.

However, it is far from sufficient if only combines features with top reproducibility mechanically, because it ignores the potential ability of features with non-top but high reproducibility and will reduces the classification accuracy. Therefore, in this study, to obtain better performances, an extremely huge number of random selections were implemented among the features with high reproducibility, then chose the best feature combination according to the comprehensive classification effect. This algorithm named BLR.

The steps of proposed BLR algorithm are as follows:

1) Use the LASSO regression model on the training cohort with 5-fold cross-validation to select the effective features, and bootstrapping 300 times.

2) Sort the features according to the reproducibility in above LASSO selections, and the top M features were selected to form the feature subset1.

3) Select N (N < M) features from the feature subset1 randomly K times (K is extremely large) to obtain K feature combinations, feed them to SVM (use 5-fold cross-validation, besides, bootstrapping repeat 100 times and take the average) to get their classification ability.

4) The threshold method was used for select L (L < < K) feature combinations with excellent classification performance from K feature combinations obtained in step 3).

5) Select the best one from L feature combinations obtained in step 4) manually according to the comprehensive classification performance, and get the corresponding final feature subset.

BLR increases the reliability of LASSO regression through a lot of bootstrapping repetitions. Besides, it avoids simple feature combination, which picks up the features with the non-top reproducibility and maximizes their diagnostic potential through a large number of random extractions. Meanwhile, it explores various possibilities of feature combinations. Therefore, the BLR algorithm has a great probability to obtain better classification effect.

2.4.   Diagnosis of benign and malignant

The final feature subset selected by BLR was sent to SVM with linear kernel to complete the differentiation between benign and malignant of IPMN and MCN. All selected features were normalized to [−1, +1].

Some important details of SVM were configured as follows: The maximal number of optimization iterations was 1,000,000. In order to solve the slight sample imbalance in the training cohort, we doubled the cost of the misjudged malignant tumors in the training stage. In the K-fold training set, the classification results were averaged as the final classification scores.

3.   Experiments

In this section, a series of experiments were conducted to illustrate the performance of our method.

3.1.   Compared methods

Our proposed method was compared with some state-of-the-art methods in four aspects: the comprehensive feature sets, the feature selection methods, the performance of various classifiers, and the overall pancreas related CAD methods.

As for the feature sets, the effect of inclusion and exclusion of clinical indices on the proposed CAD system was contrasted.

As for the feature selection methods, our BLR algorithm was compared with three most commonly feature selection algorithms, which were p-value, relief and logistic regression, respectively. In addition, the influences of two important factors of BLR were observed, which are bootstrapping repetitions and random selections. For a fair comparison, the feature sets fed into various selection methods were the same. The condition of each method is as follows:

1) P-value algorithm: The WRST based p-value feature selection method with p < 0.01 was applied.

2) Relief and logistic regression algorithm: The features with top 20 weights were selected.

3) LASSO: The lasso fit was constructed by 10-fold cross-validation, and the features corresponding to the minimum cross-validated mean squared error (MSE) were combined as final feature subset.

4) BL: The effect of bootstrapping repetitions on the LASSO model was verified and repeated this step 300 times. Then the features with top 10% reproducibility were selected as final feature subset.

5) BLR: N features were selected randomly from the top M features in BL. Repeated this step K times and obtained K feature combinations, fed them to SVM. Threshold method was used for selecting L combinations with excellent classification effect, then the best combination was selected manually as the final feature subset. The hyperparameters like N, M, K, and thresholds were optimized according to the experimental conditions.

As for the classifiers, the effect of different classifiers was compared on the 5-fold cross-validation cohort. The image preprocessing was exactly the same. The proposed feature extraction module without clinical indices was used. And the feature selection algorithm was the baseline LASSO regression. The selected features were sent to the compared classifiers, including Decision Tree (DT), K-Nearest Neighbor (KNN), Back Propagation Neural Network (BPNN), Naive Bayes (NB), and SVM.

As for the overall pancreas related CAD methods, we compared our CAD system with other pancreas related radiomics-based method, including Park et al.^[11], Zhang et al. ^[12], and Wei et al. ^[13] For Park et al., the features of volume were replaced with those of area to adapt our data dimensions. And mRMR algorithm was used to feature selection, then the random forest was used to classification. For Zhang et al., we extracted the 2D CT radiomics features mentioned in their article and used SVM-RFE for feature selection, then the selected features were sent to SVM. For Wei et al., we extracted 409 features and used their LASSO-based feature selection algorithm, then the SVM was used to diagnosis.

3.2.   Evaluation metrics

Four indices were used for evaluate the performance of established model, there are the area under the ROC curve (AUC), accuracy (ACC), sensitivity (SEN), specificity (SPE), respectively.

4.   Results and discussion

4.1.   Efficiency of different classifiers

The classification performances of 5 classifiers were compared. We extracted the proposed feature system without clinical indices and sent them to the baseline LASSO algorithm to select features. The selected features were fed to the compared classifiers, and the related experimental results are listed in Table 4.

From Table 4, in the 5-fold cross-validation cohort, the SVM obtained AUC, ACC, SEN, SPE = 0.80, 0.74, 0.68, and 0.77, respectively, which achieves the best classification effects among the compared classifiers. So, we choose SVM as the final classifier of the following experiments.

4.2.   Efficiency of feature selection algorithms

The classification efficiencies of 6 feature selection algorithms were compared in this part. The feature system without clinical indices was fed to the above 6 algorithms, and the selected feature subsets were sent to SVM to analyze their classification performance.

In a cross-validation cohort, for p-value, 24 features with p < 0.01 were fed to SVM and obtained AUC, ACC, SEN, SPE = 0.72, 0.72, 0.65, and 0.75, respectively. For relief and logistic regression, 20 features with top weights were selected as final feature subsets and obtained AUC, ACC, SEN, SPE = 0.63 versus 0.66, 0.63 versus 0.65, 0.50 versus 0.54, and 0.69 versus 0.70, respectively. For LASSO, 11 features were selected and obtained AUC, ACC, SEN, SPE = 0.80, 0.74, 0.68, and 0.77, respectively. And for BL, 13 features with top 10% reproducibility were selected and obtained AUC, ACC, SEN, SPE = 0.80, 0.74, 0.67, and 0.77, respectively. For BLR, M, N, and K were set to 23, 13, and 20,000, the thresholds were set to AUC = 0.81, ACC = 0.75, SEN = 0.70, SPE = 0.76, and selected L = 100 feature combinations. Finally, the best feature subset with 13 features was fed to SVM, obtained AUC, ACC, SEN, SPE = 0.82, 0.76, 0.71, 0.79. The above results are summarized in Table 5.

Hence, the BLR algorithm has the best performance. By a large number of bootstrapping repetitions and random selections, the BLR algorithm avoids the bias and overfitting of single LASSO regression while maximizes the potential diagnosis abilities of the features with non-top reproducibility. It should be noted that the value of N is approximately equal to the number of features selected by BL algorithm, and M should not be too large, because the excessive value of M may lead to dramatic number of combinations while consider many useless features. At the same time, the value of K should be very large to make the experiment more general.

4.3.   Analysis of the selected features

Some typical characteristics selected by BLR and their relevant information are shown in Table 6. And Figure 3 shows the boxplots of some representative features.

We found that the roundness, mean of inner-regional covariance, the run-length variance and so on are of outstanding significance for the classification of benign and malignant IPMN or MCN. Now, we will briefly observe the imaging differences between benign and malignant tumors from these important selected features.

As we can see in Figure 4, the malignant tumors are more likely to have irregular shapes, while the benign tumors tend to be rounder. This finding is in accordance with the different growth pattern of benign and malignant tumors: benign tumors have an expansible growth pattern, while malignant ones have an infiltrative pattern and tend to invade adjacent structures.

Figure 5 shows the tumors with different H-variance, this feature reflects the intensity and histogram distribution inside the tumor, indicating that the intensity of the malignant tumor is more uneven. This might be caused by the solid components inside the cyst of the tumor. These solid nodules or components grow irregularly and are long regarded as the origin of malignant formation in the tumor.

Figure 6 reflects the image differences of the benign and malignant tumors from the perspective of texture. In the malignant tumors, the texture and density distributions are more complex and inhomogeneous, the gray level ranges are wider, and the disparities between the different inner regions of the tumor are quite distinct. This phenomenon identified reflects the fact that malignant tumors usually have a more complex composition than benign tumors. Beside cancerous cells, there are also small nourish vessels resulted from neovascularization, the characteristic fibrotic stroma of pancreatic malignancy, and so on. All of these components above have distinct difference on CT images, which make the texture of malignant tumors shows more complex on the image.

A great difference between benign and malignant tumors in shape, density, and texture are illustrated. Compared with the benign tumors, the malignant tumors usually have more irregular shape, intensity, and texture distribution. However, the thickened cystic wall, the large cyst and papillary mural nodules may lead to the above problem. It also reflects that these imaging features are very important to the classification of benign/malignant IPMN or MCN.

4.4.   Diagnostic performances

In this part, the feature system with clinical indices was fed to 6 feature selection algorithms to contrast the influence of clinical indices in the classification performance.

The parameters of p-value, relief, and logistic regression, LASSO, BL algorithm were configured same as that in section 4.1. The only difference was the final feature subsets selected by LASSO and BL contain 12 and 10 features, respectively. For the BLR algorithm, the hyperparameters were as follows: M = 21, N = 10, K = 20,000, the thresholds of classifiers were AUC, ACC, SEN, SPE = 0.82, 0.75, 0.73, and 0.76, respectively. The output parameter was L = 108, and the best subset with 10 features were selected.

The final feature subsets of the various algorithms were fed to SVM to realize the identification of benign and malignant. Table 7 compares the classification results of proposed CAD scheme using different feature selection methods with or without clinical indices, and records the number of features selected by different algorithms. The ROC curves are shown in Figure 7.

As can be seen, the BLR is the state-of-the-art method. Beyond that, the clinical indices are also very important, which improved diagnosis accuracy as well as made the CAD system closer to the real clinical diagnosis. This suggests that, from a statistical point of view, the patients with malignant tumors may have some abnormal symptoms such as abdominal pain, weight loss, and so on. However, different people may exhibit different characteristics.

4.5.   Overall comparison with other CAD methods

The overall diagnosis performances of different pancreas related CAD methods were compared. The results are compared as follows:

As can be seen from Table 8, our method achieves the state-of-the-art effects, proving our superiority. Especially for the sensitivity of diagnosis, which is usually low in clinical artificial diagnosis and other methods, our approach has made great progress. Compared to Park et al., the SEN of our method is increased from 0.52 to 0.73 in cross-validation cohort, and from 0.44 to 0.89 in independent validation cohort. The superiority of our method is due to the comprehensively, completely, and PCN-optimized feature system and the excellent feature selection algorithm.

5.   Conclusions

In conclusion, this study builds a CAD system based on radiomics, and realizes the differentiation malignant from benign of IPMN and MCN. A complete, comprehensive, and PCN-optimized feature system is formed in our research, meanwhile 9 clinical indices are added to improve the diagnostic accuracy. Besides, a novel and wonderful feature selection algorithm is proposed and achieved remarkable effects. Our CAD system forms an objective tumor evaluation method which can reduce the probability of misdiagnosis of malignant tumors, so as to provide crucial suggestions for the treatment planning and delivery for patients. In addition, the proposed CAD system provides a time- saving diagnosis scheme with low computation cost. The calculation part of our CAD system only takes less than 2 s per patient.

However, we still have a lot of room for improvement. Firstly, a small dataset may underestimate the effect and practicability of proposed CAD system. In the future work, the data set can be expanded to make the results more universal. Secondly, manual segmentation can ensure the diagnosis accuracy of proposed CAD system, but it increases the time cost, an automatic segmentation method is considered in the next step.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant 61771143, 61871135 and 81830058), NSFC-DFG Cooperation Group (GZ1456) and the Science and Technology Commission of Shanghai Municipality (Grant 18511102904 and 19511121200).

Conflict of interest

There is no conflict of interest in this work.

Appendix