Abstract
Background. In a pathological examination of pancreaticoduodenectomy for pancreatic head adenocarcinoma, a resection margin without cancer cells in 1 mm is recognized as R0; a resection margin with cancer cells in 1 mm is recognized as R1. The preoperative identification of R0 and R1 is of great significance for surgical decision and prognosis. We conducted a preliminary radiomics study based on preoperative CT (computer tomography) images to evaluate a resection margin which was R0 or R1. Methods. We retrospectively analyzed 258 preoperative CT images of 86 patients (34 cases of R0 and 52 cases of R1) who were diagnosed as pancreatic head adenocarcinoma and underwent pancreaticoduodenectomy. The radiomics study consists of five stages: (i) delineate and segment regions of interest (ROIs); (ii) by solving discrete Laplacian equations with Dirichlet boundary conditions, fit the ROIs to rectangular regions; (iii) enhance the textures of the fitted ROIs combining wavelet transform and fractional differential; (iv) extract texture features from the enhanced ROIs combining wavelet transform and statistical analysis methods; and (v) reduce features using principal component analysis (PCA) and classify the resection margins using the support vector machine (SVM), and then investigate the associations between texture features and histopathological characteristics using the Mann–Whitney U-test. To reduce overfitting, the SVM classifier embedded a linear kernel and adopted the leave-one-out cross-validation. Results. It achieved an AUC (area under receiver operating characteristic curve) of 0.8614 and an accuracy of 84.88%. Setting in the Mann–Whitney U-test, two features of the run-length matrix, which are derived from diagonal sub-bands in wavelet decomposition, showed statistically significant differences between R0 and R1. Conclusions. It indicates that the radiomics study is rewarding for the aided diagnosis of R0 and R1. Texture features can potentially enhance physicians’ diagnostic ability.
1. Background
Pancreatoduodenectomy is the main treatment for pancreatic head adenocarcinoma. Knowledge of preoperative assessment of cancer resection and excision expansion will help to choose optimal therapies for patients. Thus, it is very important to evaluate the resection margin of pancreaticoduodenectomy. In a pathological examination of pancreaticoduodenectomy, a resection margin without cancer cells in 1 mm is recognized as R0; a resection margin with cancer cells in 1 mm is recognized as R1. The preoperative identification of R0 and R1 is of great significance for surgical decision and prognosis.
Intertumoral heterogeneity is generally considered as a typical finding of malignancy. It reflects variations in tumor-cell differentiation, extracellular matrix, and cellularity angiogenesis [1]. Image-based texture analysis is a noninvasive technique for quantifying tumor heterogeneity and has been widely applied to aided diagnosis, efficacy evaluation, and prognosis [2]. This is termed radiomics [3, 4]. Computer tomography (CT) is a commonly used examination for diagnosis of pancreatic head adenocarcinoma. To the best of our knowledge, there are currently no texture-based radiomics studies yet to evaluate the aided diagnosis of R0 and R1, but there are a few similar studies of pancreatic cancer on portal-venous phase CT images. In 2017, Cassinotto et al. [5] used the Laplacian of Gaussian (LoG) filter and histogram to extract texture features to evaluate pathologic tumor aggressiveness and predict disease-free survival in patients with resectable pancreatic adenocarcinoma; Eilaghi et al. [6] used the method of gray-level co-occurrence matrix (GLCM) to extract texture features to assess whether CT-derived texture features predict survival in patients undergoing resection for pancreatic ductal adenocarcinoma; and Chakraborty et al. [7] used the methods of histogram, GLCM, gray-level run-length matrix (GLRLM), and angle co-occurrence matrix (ACM), etc. to extract texture features to predict 2-year survival of pancreatic ductal adenocarcinoma (PDAC). In 2018, Canellas et al. [8] used the LoG filter and histogram to extract texture features to assess whether CT texture analysis and CT features are predictive of pancreatic neuroendocrine tumor grade based on the World Health Organization classification and to identify features related to disease progression after surgery; Qiu et al. [9] used the methods of histogram, GLCM, wavelet transform, etc. to extract texture features on nonenhanced CT images and then explored the feasibility of discriminating pancreatic cancer from normal pancreas. In 2019, Cheng et al. [10] used the LoG filter and histogram to extract texture features to determine if CT texture analysis measurements of the tumor are independently associated with progression-free survival and overall survival in patients with unresectable PDAC.
We evaluated whether an operation was performed by R0 resection or R1 resection based on the radiomics aided diagnosis on its preoperative portal-venous CT images and investigated the differences of histopathological characteristics between R0 and R1 by using statistical significance tests of texture features. This study has been approved by the Ethics Committee of West China Hospital of Sichuan University (trial registration: NCT02928081).
2. Methods
Figure 1 illustrates the framework of our radiomics study. It consists of five stages: Stage 1: obtain ROIs (regions of interest) by preprocessing patients’ CT images Stage 2: by solving discrete Laplacian equations with Dirichlet boundary conditions, fit the ROIs to rectangular regions Stage 3: combining wavelet transform and fractional differential, enhance the textures of the rectangular ROIs Stage 4: combining wavelet transform and statistical analysis methods, extract texture features from the enhanced ROIs Stage 5: reduce features using principal component analysis (PCA), perform classification using the support vector machine (SVM) (to reduce overfitting, the SVM classifier embeds a linear kernel and adopts the leave-one-out cross-validation), and then investigate the associations between texture features and histopathological characteristics using the Mann–Whitney U-test
2.1. Patients
This study retrospectively analyzed 258 preoperative CT images of 86 patients (34 cases of R0 and 52 cases of R1) who were diagnosed as pancreatic head adenocarcinoma at West China Hospital from October 2015 to March 2018. These patients underwent pancreaticoduodenectomy. The surgeries were pathologically diagnosed as R0 resection or R1 resection. Patients were screened based on NCCN guidelines for diagnostic criteria and standard surgical procedures. We selected 3 portal-venous phase CT images from each case for analysis, which are located at the top, middle, and bottom of a tumor [11].
Abdominal scan and enhanced scan were performed using 64-slice spiral CT of American GE. Collimator was set to 0.625 mm, FOV was set to 350 mm × 350 mm, tube voltage was set to 120 kV, tube current was set to 160 mAs, and layer thickness was set to 1.250 mm. In enhanced scanning, iopamide was injected via cubital veins, and flow rate was 3 ml/s, dose was 90∼100 ml, and delayed time was 25∼30 s for scanning of the portal-venous phase. A CT image was exported as an 8-bit grayscale image.
2.2. Delineation and Segmentation
The steps for delineating and segmenting are as follows: (1) choose three portal-venous phase CT images from each case, which are located at the top, middle, and bottom of a tumor; Figure 2 [11] illustrates the locations; (2) delineate resection margins around portal veins on the chosen images, and it is shown in Figure 3; to ensure authenticity of signals, the delineated resection margins exclude edges of stent and metal artifacts; and (3) segment the delineated regions to form ROIs based on a region growing segmentation method.
(a)
(b)
Two physicians with 10 years of experience in abdominal CT diagnosis delineated all resection margins. The first physician delineated the resection margins and repeated the delineations after 2 weeks to prevent observer deviations. The other physician only delineated the resection margins once to assess whether his delineations were consistent with the delineations of the first physician.
2.3. Fitting ROIs
We fitted the strip-shaped ROIs to rectangular ROIs by solving discrete Laplacian equations with Dirichlet boundary conditions. The fitting method is abbreviated as LD. The LD method has good applications in signal fitting [121314]. Discrete Laplace equation can be defined in the following equation:
Equation (1) shows that a linear equation can be established based on a 4-neighborhood of a point . The point is to be fitted. A region to be fitted is named as a mask. If the current pixel is on an edge of the mask, then at least one of its neighbors (on the Dirichlet boundary) is known. A set of linear equations can be established along the Dirichlet boundary (along edges of the mask). The pixel values to be fitted can be obtained by solving the established set of linear equations. The solving procedure is then extended into the interiors of the mask. Figure 4 shows a mask to be fitted and its boundaries.
(a)
(b)
Figure 5 illustrates two fitted examples, where the black regions of an ROI are the mask of this ROI, and the R0 ROI is better fitted, while the fitted regions of the R1 ROI are smoother. Actually, the fitted regions are too smooth to express more information. Next, we would enhance the textures of the fitted regions.
(a)
(b)
(c)
(d)
(e)
(f)
2.4. Enhancing ROIs
We designed a texture enhancement method with reference to the Grunwald–Letnikov (G-L) fractional differential definition and wavelet transform [15, 16]. The enhancement method is abbreviated as WF. It consists of 3 steps as illustrated in Figure 6. Step 1: decompose an ROI into 4 components using wavelet transform [17]: H (horizontal), V (vertical), and D (diagonal) are the high-frequency components; A (approximate) is the low-frequency component. It is 1-level decomposition. The approximate component can be decomposed again. Step 2: convolve each high-frequency component with a fractional differential operator M. Step 3: perform wavelet inverse transform based on the convolution results of Step 2 and the approximate component in the last-level decomposition.
Wavelet inverse transform will reconstruct the ROI, which is the enhanced ROI. The steps for constructing a fractional mask are as follows:(1)Discretize G-L definition: equation (2) is the -order G-L definition of on , where is a gamma function; discretize the continuous duration equally by unit interval , where ; and it is known that , and equation (3) can be derived:(2)Expand equation (3): it is known that (unit interval), and equation (4) can be derived as follows:
We constructed a fractional differential operator named M based on the expanded coefficients of equation (4). Figure 7 demonstrates the operator M. It performs fractional differential operations in eight symmetric directions in a 5 × 5 neighborhood. at the center point position is an adjustable parameter and is called the compensation parameter. In experiments, the order and the parameter can be appropriately adjusted. Figure 5 illustrates two enhancing examples using the WF method.
2.5. Texture Analysis
We used rbio2.8 for wavelet transform. The steps for feature extraction are as follows: Step 1: fit and enhance the ROIs as described in the previous sections. Step 2: perform wavelet transform on the fitted and enhanced ROIs. A decomposition of a fitted and enhanced ROI will derive 4 components; a coefficient matrix uniquely expresses a component. Step 3: convert high-frequency components to grayscale images called sub-band images.
In the coefficient matrix of a high-frequency component, elements with larger absolute values usually represent singular value points (meaning a fast and large change). First, absolute values of coefficient matrices are calculated. Then, elements of a coefficient matrix are linearly and equally discretized into a grayscale range of [0, 255] (the range of gray level) according to the minimum and maximum values of the coefficient matrix. The calculations are shown in equations (5)(7), where is the coefficient matrix and is the discretized matrix (sub-band image): Step 4: extract features from the sub-band images using the methods of histogram, co-occurrence matrix, and run-length matrix. Considering that the sizes of sub-band images are also small, the gray level is rescaled.
2.6. Feature Reduction and Classification
Reducing features can usually improve classification performance. We used principal component analysis (PCA) for feature reduction and limited the number of features to reduce overfitting. Empirically, it is appropriate that the number of features is 1/15 or 1/10 of the number of samples, and a linear classifier allows for more features.
Support vector machine (SVM) [18] is widely used due to its outstanding performance in pattern recognition problems of small sample sizes. To reduce overfitting, we used a linear kernel and used the leave-one-out cross-validation. Linear kernel-based SVM allows more features without easily overfitting. In the vast majority of cases, especially in classification problems of small sample sizes, the model evaluated in the leave-one-out cross-validation is close to the model that expected to be evaluated using a training set. Thus, evaluation results of the leave-one-out cross-validation are often considered more accurate [19].
3. Results
We performed other texture analysis methods that are frequently used in pancreatic cancer-related radiomics studies and applied the PCA-based feature reduction method and the linear SVM classification method. Table 1 shows the texture analysis methods.
Considering the size of an ROI is small, we performed 1-level wavelet decomposition and set the distances of the co-occurrence matrix to 1 and 2. Feature values of 4 directions (0, 45, 90, and 135) of a co-occurrence matrix were averaged, so was a run-length matrix. Wavelet transform should be performed on rows and columns. Before applying the WT method and WT-HCR method, we filled ROIs into valid matrices based on interpolation methods. The linear interpolation method was first applied, and then we fill the remaining missing values using the nearest interpolation method. The LD-WF method used a reverse biorthogonal wavelet and selected rbio2.8 by experiments. Figure 8 illustrates two examples of decomposing ROIs using the rbio2.8 wavelet.
(a)
(b)
A binary classification problem can use a confusion matrix to express the results. R1 is used as the positive class, and R0 is used as the negative class. Table 2 shows the experimental results. The LD-WF method achieves the best classification performance, and its accuracy and AUC are 84.88% and 0.8641, respectively, followed by the LOG-GH method and the CTM method. Although the accuracy of CTM is lower than that of LOG-GH, its AUC value is higher than LOG-GH. The ROC (receiver operation curve) and AUC (area under the ROC) are powerful indicators for measuring a binary classification model.
Figure 9 illustrates the ROC curves of these methods. The classifier based on the LD-WF method approaches the upper-left corner faster followed by the classifier based on the CTM method and the classifier based on the LOG-GH method.
4. Discussion
This study aims to conduct a preliminary radiomics exploration to evaluate whether a surgery was performed by R0 resection or R1 resection based on its surgical margin of portal-venous CT images. To eliminate the bias of possible episodes of acute or chronic pancreatitis, which can occur concomitantly with the neoplastic evolution or the pancreatic reaction after endoscopic biopsy sampling etc., all the selected patients underwent pancreaticoduodenectomy followed by pathological diagnosis, and the pathological diagnoses were used as the gold standard. Physicians delineated the resection margins around portal veins on the chosen CT slices as the initial ROIs. It is illustrated in Figure 3.
In an R0 or R1 resection margin, an ROI is an irregular strip-shaped region, and its structure contains complex internal details such as capillary distribution, cancer cell tissue, and pancreatic cell tissue. Statistical texture analysis methods are appropriate for this. Multiresolution texture analysis methods perform well in extracting detail features. However, both statistical texture analysis methods and multiresolution texture analysis methods are limited to irregular strip-shaped and small ROIs. Figure 3 shows two examples of irregular strip-shaped regions.
An image is a two-dimensional signal based on rows and columns, and two-dimensional relationships between intensity and position usually better express texture characteristics. Thus, we used the LD method to fit the ROIs to rectangular regions. Furthermore, to make the fitted regions express more information and further improve the performance of the texture analysis method, we designed the WF method to enhance the textures of the fitted ROIs. The main purpose of texture enhancement is to highlight high-frequency contour information (detailed information, that is, portions of gray levels that change relatively more varied or more quickly) while preserve low-frequency smoothing information as much as possible. Traditional enhancement methods such as histogram equalization, integer-order differentials, and frequency enhancement filters, increase contrast or highlight contours, but they lose lots of low-frequency texture information and usually sharpen contour information. In recent years, fractional differentials compensate for the drawback of greatly losing low-frequency information, making it an effective method for texture enhancement of medical images [20–22]. Thus, we consider the following 3 factors to enhance the textures: (1) wavelet transform is appropriate for detail analysis of an image, and its characteristic of perfect inverse transform enables corrections of transform coefficients to be highlighted in the reconstructed image; (2) fractional differential can enhance contours without sharpening edges; and (3) characteristics of the details usually well characterize lesions or tissues. We designed the WF method based on these 3 factors.
After fitting and enhancing the ROIs, texture analysis methods were used to extract quantitative features: texture features. Deep learning algorithms have made significant progress in image pattern recognition. However, these algorithms are limited by the problems of small sample sizes, small targets, and so on [23, 24]. Moreover, deep learning algorithms lack pertinence in quantitative analysis of ROIs. Therefore, ROI-based radiomics is still a mainstream scheme in medical image-aided diagnoses.
In histopathology, an ROI of R1 has cancer cells, some parts of its tissue are more compact, and its capillary distribution is less, while an ROI of R0 has no cancer cells within 1 mm, it only contains pancreatic tissue, and its capillary distribution is more abundant [25, 26]. However, these differences are just qualitative in details and difficult to visually observe from CT images. Multiresolution analysis methods are advantageous in local time-frequency analysis and are appropriate for deriving detail characteristics. Statistical analysis methods can usually derive representative mathematical descriptors. It can be inferred that multiresolution analysis methods and statistical analysis methods are appropriate here. Based on the stated characteristic analysis and texture analysis of the ROIs, we combined the methods of wavelet transform, histogram, GLCM, and GLRLM to extract texture features.
Radiomics uses computer methods such as computer vision and machine learning to perform digital medical image processing, which can deeply mine the heterogeneous data at levels of tissue and molecular that contained in medical images such as CT images [2, 27, 28]. CT imaging is that X-rays penetrate different media with different attenuations to form different gray levels. Thus, grayscale patterns in CT images should be able to reflect changes of body’s pathology. From histopathological analysis, an R1 resection margin contains a large number of normal pancreatic tissue and some tumor tissue, and its capillary distribution is less than an R0 resection margin; relatively, an R0 resection margin only contains normal pancreatic tissue, and its capillary distribution is more abundant. Thus, characteristics of internal details can better characterize R0 and R1. Analogous to wavelet transform, LOG-GH is also a multiscale analysis method. Both types of methods are suitable for characterizing detail characteristics. From the classification results, the multiresolution or multiscale analysis methods behave better.
In addition, it is necessary to address some issues such as the problem of irregular strip-shaped ROIs and the problem of atypical manifestations of details (macroscopically difficult to distinguish). This radiomics study used the LD-WF method to process ROIs (fitted the ROIs and enhanced textures) followed by combining wavelet transform and statistical methods to extract descriptors on the sub-band images. The experimental results indicated that it pronouncedly improved classification performance.
We expect that some texture features should be able to reflect the differences between R0 and R1. To investigate the discriminations of texture features between R0 and R1, we performed Mann–Whitney U-tests on the texture features that are extracted based on the LD-WF method. Table 3 shows the features with , which usually means that there are statistically significant differences between the two types of samples (R0 samples and R1 samples). It demonstrates that the middle and bottom ROIs present more differences on the texture features, and the diagonal sub-band image expresses more characteristic differences in detail. The values of F4 and F6 are ≤0.01, which means that there are extremely significant differences between the two types of samples.
Three features were selected based on the ascending order of values. Table 3 shows these three features in bold: F4, F6, and F9. To test the feature values, larger or smaller, right-tailed hypothesis tests based on the Wilcoxon rank sum method were performed on F4, F6, and F9, where the alternative hypothesis states that the median of R1 samples is greater than the median of R0 samples. Table 4 demonstrates the results of right-tailed hypothesis tests.
Table 4 shows that F4-values of R1 are larger than F4-values of R0 at significant level , F6-values of R1 are larger than F6-values of R0 at significant level , and F9-values of R1 are larger than F9-values of R0 at significant level . In wavelet transform, every coefficient is in charge of an oscillation in certain scale and frequency. The discussions of these three features are as follows.
As for feature F4, (1) the sub-band image expresses the component that gray level changes more and faster in diagonal direction; (2) high gray-level run emphasis (HGRE) of the run-length matrix measures the distribution of higher gray-level values, with a higher value indicating a greater concentration of high gray-level values in an image; (3) in the diagonal component, higher gray level means larger oscillation; and (4) the test result for F4 in Table 3 indicates that points with larger oscillations appear more continuously in ROIs of R1 than those of in ROIs of R0; this should be associated with the fact that the ROIs of R1 contain normal pancreatic tissue and cancer tissue, while the ROIs of R0 only contain normal pancreatic tissue.
As for feature F6, it is similar to F4. Short run high gray-level emphasis (SRHGE) is a supplement to HGRE, indicating that points with larger oscillations (fine texture) appear more continuously.
As for feature F9, (1) the meaning of the diagonal sub-band image has explained above; (2) the cubic moment of the histogram measures skewness, higher skewness means greater degree of asymmetry; and (3) the test result for F9 in Table 3 indicates that the degree of asymmetry in R1 is greater than that in R0; it should still be associated with the fact that the ROIs of R1 contain normal pancreatic tissue and cancer tissue, while the ROIs of R0 only contain normal pancreatic tissue; because R0 has only normal pancreatic tissue, the structural changes on the diagonal component are relatively more uniform and more symmetry.
This study has some limitations and deficiencies. First, it was a retrospectively study in a single institution, patients’ population and imaging methods were basically homogeneous, and selection bias may exist, making it difficult to generalize the results to other institutions. Second, ROIs were fitted to rectangles, but the pixel size of a ROI is still small. Third, sensitivity still needs to be improved. Finally, no sufficient samples for the test led to some overfitting (although the leave-one-out cross-validation is used). Next, we will collect more samples and conduct further studies using better fitting methods.
5. Conclusions
By analyzing the histopathological characteristics of R0 and R1 and considering the deficiencies that ROIs are irregular strip-shaped and small regions, we designed the LD-WF method and conducted a preliminary radiomics study based on portal-venous CT images to identify whether a surgery was conducted by R0 resection or R1 resection. The experimental results indicate that the designed method is rewarding for discriminating R0 from R1. By analyzing statistically significant differences on texture features, it elucidates that the histopathological characteristics of R0 and R1 can be represented by the texture features of preoperative CT images. It implies that texture features can potentially enhance physicians’ diagnostic abilities.
Abbreviations
| CAD: | A computer-aided diagnosis |
| CT: | Computer tomography |
| ROI: | Region of interest |
| PCA: | Principal component analysis |
| SVM: | Support vector machine |
| AUC: | Area under receiver operating characteristic curve |
| ROC: | Receiver operating characteristic |
| LoG: | Laplacian of Gaussian |
| GLCM: | Gray-level co-occurrence matrix |
| GLRLM: | Gray-level run-length matrix |
| ACM: | Angle co-occurrence matrix |
| PDAC: | Pancreatic ductal adenocarcinoma |
| LD: | A fitting method by solving discrete Laplacian equations with Dirichlet boundary conditions |
| G-L: | Grunwald–Letnikov |
| WF: | A texture enhancing method based on fractional differential and wavelet transform. |
Data Availability
The datasets during and/or analyzed during the current study are available from the corresponding author on reasonable request with the approval of the institution and trial/study investigators who contributed to the dataset.
Ethical Approval
The study in this paper has been approved by the Ethics Committee of West China Hospital of Sichuan University (the trial registration number in ClinicalTrials.gov is NCT02928081).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Bei Hui, Jia-Jun Qiu, and Jin-Heng Liu substantially contributed to conception and design; Jin-Heng Liu and Neng-Wen Ke were responsible for acquisition of data; Jia-Jun Qiu contributed to analysis and interpretation of data; and Bei Hui and Jia-Jun Qiu drafted the article and revised it critically for important intellectual content. Bei Hui was responsible for the agreement to be accountable for all aspects of the work. All the authors approved the final version to be published.
Acknowledgments
This work was supported by the Key Study Plan for State Commission of Science Technology of China under Grant no. 2018YFC0807501.