Reliable gene mutation prediction in clear cell renal cell carcinoma through multi-classifier multi-objective radiogenomics model

2.1.1. Patients

2.1.2. CT image features

CT images from UTSW were acquired by GE LightSpeed VCT (GE Healthcare, Waukesha, WI) or TOSHIBA Aquilion ONE (Canon Medical Systems USA, Tustin, CA). CT image size was 512  ×  512 with a pixel size of 0.7–0.9 mm, and slice thickness was 3 or 5 mm.

The primary tumor contour was delineated by a radiation oncologist with 4 years of experience and reviewed by a radiation oncologist with nine years of experience. Contrast enhanced CT images acquired during the corticomedullary phase were used in all 57 cases for image analysis. A region of interest (ROI) was drawn along the outer contour of the mass using the Velocity 3.2.0 software excluding adjacent tissues (e.g. renal parenchyma, peri-renal fat) (figure 1).

Zoom In Zoom Out Reset image size

We resampled all images of the same slice thickness at 5 mm. We only considered primary tumors and defined 43 quantitative image features describing tumor characteristics, including 13 geometry features, nine intensity features, and 21 texture features (table 2). Geometry features describing tumor shape and size were calculated according to the actual pixel size. Features Size_X, Size_Y and Size_Z (table 2) describe the tumor size along the X, Y, and Z axes of the digital imaging and communications in medicine (DICOM) coordinate system (figure 2(a)). The shape and location of the tumors differed from patient to patient. To intuitively describe the size of the tumor, a principal component analysis (PCA) was applied to the tumor contour points to transform the data into a new coordinate system. Size_P1 is the maximum three-dimensional diameter of the tumor, measured as the largest pairwise Euclidean distance between the voxels on the surface of the tumor volume; size_P2 and size_P3 are the tumor size along the directions orthogonal to the direction of maximum size (figure 2(b)).

Table 2. Quantitative CT image features.

Geometry features	Intensity features	Texture features
Volume	Minimum	Auto correlation
Size_X	Maximum	Contrast
Size_Y	Mean	Correlation
Size_Z	Median	Cluster prominence
Size_P1	Sum	Cluster shade
Size_P2	Variance	Dissimilarity
Size_P3	Standard deviation	Energy
Roundness	Skewness	Entropy
Surface area	Kurtosis	Homogeneity
Compactness_1		Maximum probability
Compactness_2		Variance
Spherical disproportion		Sum average
Surface to volume ratio		Sum variance
		Sum entropy
		Difference variance
		Difference entropy
		Information measure of correlation_1
		Information measure of correlation_2
		Inverse difference
		Inverse difference normalized
		Inverse difference moment normalized

Zoom In Zoom Out Reset image size

Intensity features are first-order statistics that describe the distribution of the voxel intensities within the tumor on the CT image through commonly used basic metrics (Aerts et al 2014). Intensity features provide information related to the gray-level distribution of the image; however, they do not provide any information on the relative position of the various gray levels over the image. Therefore, we included textural features that either describe the patterns or the spatial distribution of voxel intensities, which were calculated from the gray level co-occurrence matrix (GLCM) (Haralick et al 1973). Texture matrix representation requires the voxel intensity values within the volume of interest to be discretized. In this work, voxel intensities were resampled into 64 equally spaced bins using a bin-width of 25 hounsfield units (HU). The detailed methodology of extracting intensity and texture features was previously described by Aerts et al (2014) (supplementary document). The calculated features of the 57 patients were normalized using the Z-scores method (Cheadle et al 2003).

2.2.1. Evidential reasoning based classifier fusion

The evidential reasoning (ER) approach was used for fusing the individual classifier probability output (Yang and Xu 2002, 2013). Our study is a binary classification problem (mutation or non-mutation). Assuming there are $M$ classifiers, for a test sample, the output probability of each classifier is denoted by ${{P}_{i}}=\left\{P_{i}^{1},P_{i}^{2} \right\},~i=1,\cdots, M$, which satisfies:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{i}^{1}+P_{i}^{2}=1.\nonumber \end{align} \tag{ 1 }$$

Where $P_{i}^{1}$ is output probability of gene mutation and $P_{i}^{2}$ is output probability of non-mutation. Assume the relative weight of each classifier as $\boldsymbol{w}=\{{{w}_{1}},{{w}_{2}},\cdots ,{{w}_{M}}\}$, which satisfies the following constraint:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{i=1}{\overset{M}{\mathop \sum }}\,{{w}_{i}}=1,\quad 0\leqslant {{w}_{i}}\leqslant 1.\nonumber \end{align} \tag{ 2 }$$

The final output probabilities $P_{fin}^{\,j},j=1,2$ are obtained by classifier fusion through the ER approach (Yang and Xu 2002, 2013):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{fin}^{\,j}={\rm ER}\left( P_{i}^{\,j},{{w}_{i}} \right),\quad i=1,\cdots, M, ~ j=1,2,\nonumber \end{align} \tag{ 3 }$$

where ER represents the ER analytic algorithm (Wang et al 2006), which is calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{fin}^{\,j}=\frac{\mu \times \left[ \mathop{\prod }_{i=1}^{M}\left( {{w}_{i}}P_{i}^{\,j}+1-{{\omega }_{i}} \right)-\mathop{\prod }_{i=1}^{M}\left( 1-{{w}_{i}} \right) \right]}{1-\mu \times \left[ \mathop{\prod }_{i=1}^{M}(1-{{w}_{i}}) \right]},\quad j=1,2,\nonumber \end{align} \tag{ 4 }$$

where $\mu$ is calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu ={{\left[ \underset{j=1}{\overset{2}{\mathop \sum }}\,\prod\nolimits_{i=1}^{M}{\left( {{w}_{i}}P_{i}^{\,j}+1-{{w}_{i}} \right)}-\prod\nolimits_{i=1}^{M}{\left( 1-{{w}_{i}} \right)} \right]}^{-1}}.\nonumber \end{align} \tag{ 5 }$$

For a binary classification problem, if $P_{fin}^{1}>P_{fin}^{2}$, the test sample belongs to class 1; if $P_{fin}^{1}<P_{fin}^{2}$, the test sample belongs to class 2; if $P_{fin}^{1}=P_{fin}^{2}$, the test sample belong to either class. In this study, we sought to predict either the presence or absence of gene mutation. Therefore, if $P_{fin}^{1}>P_{fin}^{2}$, we considered the test sample had mutation; if $P_{fin}^{1}\leqslant P_{fin}^{2}$, we considered the test sample had non-mutation.

2.2.2. Reliable outcome prediction based on output probability similarity

To obtain more reliable predictive results, reliable outcome prediction (RCP) is proposed and defined as maximize the similarity between predicted output probability and true label vector. For example, we assume two models that predict the VHL gene mutation with the label vector [1, 0]. Model A has prediction probabilities (0.8, 0.2) (the probability of mutation is 0.8, the probability of non-mutation is 0.2, and the threshold is 0.5), and model B has prediction probabilities (0.55, 0.45). As 0.8 is closer to 1 than 0.55, the result of model A is more reliable than that of model B. In other words, the similarity between the probability output of model A and the label vector is higher than those observed for model B, which means model A is more reliable than model B in this prediction.

In RCP, the aim is to maximize the similarity between predicted output probability and true label vector T while training the single classifier model and weights. For a training sample, its label vector is denoted by T  =  $\left[ {{T}_{1}},{{T}_{2}} \right].$ $T$ is a binary vector, T $=\left[ 1,0 \right]~$ (mutation) or T $=\left[ 0,1 \right]$ (non-mutation). Assuming that the predictive model has q parameters denoted by R $=\{{{R}_{1}},{{R}_{2}},\cdots ,{{R}_{q}}\}$, the objective function is expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{~ }{\mathop{f=\underset{\boldsymbol{w},\boldsymbol{R}}{\mathop{{\rm max}}}\,}}\,\underset{k=1}{\overset{K}{\mathop \sum }}\,{\rm sim}\left( {{P}^{k}},{{T}^{k}} \right),\nonumber \end{align} \tag{ 6 }$$

where $K$ represents the number of training samples and sim is the similarity measure. Since the above problem can be considered as the similarity of probability distribution and the dice coefficient (Cha 2007) is effective for measuring similarity, it is used here:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm sim}\left( {{P}^{k}},{{T}^{k}} \right)=\frac{2\sum\limits_{j=1}^{2}{P_{j}^{k}T_{j}^{k}}}{\sum\limits_{j=1}^{2}{{{\left( P_{j}^{k} \right)}^{2}}}+\sum\limits_{j=1}^{2}{{{\left( T_{j}^{k} \right)}^{2}}}}.\nonumber \end{align} \tag{ 7 }$$

Since a single objective may not be a good measure when the training dataset is imbalanced, we consider sensitivity and specificity simultaneously as a better solution as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{f}_{sen}}=\frac{{\rm TP}}{{\rm TP}+{\rm FN }~ },\quad {{f}_{spe}}=\frac{{\rm TN}}{{\rm TN}+{\rm FP }~ },\nonumber \end{align} \tag{ 8 }$$

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives, and FN is the number of false negatives. In our previous study, ${{f}_{sen}}$ and ${{f}_{spe}}$ were both considered as objective functions (Zhou et al 2017).

However, ${{f}_{sen}}$ and ${{f}_{spe}}$ are label based measures, while we sought to maximize the similarity of probability output and the true label vector. Therefore, we defined the new similarity-based sensitivity and specificity denoted by ${{f}_{sim\_sen}}$ and ${{f}_{sim\_spe}}$, respectively. Assume that $\left\{P_{tp}^{1},P_{tp}^{2},~\cdots, P_{tp}^{TP} \right\}$ represents the probability output of true positives and the corresponding true label vector is $\left\{T_{tp}^{1},T_{tp}^{2},~\cdots, T_{tp}^{TP} \right\}$. The similarity of true positives ${\rm T}{{{\rm P}}_{sim}}$ is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm T}{{{\rm P}}_{sim}}=\underset{k=1}{\overset{{\rm TP}}{\mathop \sum }}\,{\rm sim}\left( P_{tp}^{k},T_{tp}^{k} \right)=\underset{k=1}{\overset{{\rm TP}}{\mathop \sum }}\,\frac{2\sum\limits_{j=1}^{2}{P_{tp,j}^{k}T_{tp,j}^{k}}}{\sum\limits_{j=1}^{2}{{{\left( P_{tp,j}^{k} \right)}^{2}}}+\sum\limits_{j=1}^{2}{{{\left( T_{tp,j}^{k} \right)}^{2}}}}.\nonumber \end{align} \tag{ 9 }$$

In gene mutation prediction, j  =  1 represents mutation, j  =  2 represents non-mutation, and ${{P}^{k}}$ is the mutation probability. Also, $P_{tp,1}^{k}={{P}^{k}}$, $P_{tp,2}^{k}=1-{{P}^{k}}$, $T_{tp,1}^{k}=1,$ and $T_{tp,2}^{k}=0$. Therefore, equation (9) can be simplified as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm T}{{{\rm P}}_{sim}}=\underset{k=1}{\overset{{\rm TP}}{\mathop \sum }}\,\frac{{{P}^{k}}}{{{\left( {{P}^{k}} \right)}^{2}}-{{P}^{k}}+1}.\nonumber \end{align} \tag{ 10 }$$

Similarly, we define the similarity of true negatives ${\rm T}{{{\rm N}}_{sim}}$, false positives ${\rm F}{{{\rm P}}_{sim}}$, and false negatives ${\rm F}{{{\rm N}}_{sim}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm T}{{{\rm N}}_{sim}}=\underset{k=1}{\overset{{\rm TN}}{\mathop \sum }}\,\frac{1-{{P}^{k}}}{{{\left( {{P}^{k}} \right)}^{2}}-{{P}^{k}}+1},\nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm F}{{{\rm P}}_{sim}}=\underset{k=1}{\overset{{\rm FP}}{\mathop \sum }}\,\frac{1-{{p}^{k}}}{{{\left( {{P}^{k}} \right)}^{2}}-{{P}^{k}}+1},\nonumber \end{align} \tag{ 12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm F}{{{\rm N}}_{sim}}=\underset{k=1}{\overset{{\rm FN}}{\mathop \sum }}\,\frac{{{p}^{k}}}{{{\left( {{p}^{k}} \right)}^{2}}-{{p}^{k}}+1}.\nonumber \end{align} \tag{ 13 }$$

Then, ${{f}_{sim\_sen}}$, ${{f}_{sim\_spe}}$ are calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{f}_{sim\_sen}}=\frac{{\rm T}{{{\rm P}}_{sim}}}{{\rm T}{{{\rm P}}_{sim}}+{\rm F}{{{\rm N}}_{sim}}~},\ {{f}_{sim\_spe}}=\frac{{\rm T}{{{\rm N}}_{sim}}}{{\rm T}{{{\rm N}}_{sim}}+{\rm F}{{{\rm P}}_{sim}}~}.\nonumber \end{align} \tag{ 14 }$$

Our aim was to maximize the two similarity-based objective functions simultaneously as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{f}_{sim}}=\underset{\boldsymbol{w},\boldsymbol{R}}{\mathop{\max }}\,\left( {{f}_{sim\_sen}},~{{f}_{sim\_spe}} \right).\nonumber \end{align} \tag{ 15 }$$

Once training is finished, the Pareto-optimal solution set is generated, and the best model parameters and weights are selected based on the clinical needs. In the following subsection, we describe a new algorithm that was developed to solve the similarity-based multi-objective optimization problem.

Multi-objective evolutionary algorithms (MOEA) have demonstrated the superior performance for multi-objective optimization (Deb 2001). Based on MOEA, we have proposed an iterative multi-objective immune algorithm (IMIA), which adopts the traditional sensitivity and specificity as the optimized objective functions (Zhou et al 2017). Based on IMIA, we propose a new SMO algorithm. The major difference between IMIA and SMO is that the reliability is measured through similarity between output probability and label vector during the training process. Additionally, weighting coefficients needs to be optimized in SMO. For conciseness, we just gave a brief description of SMO and focused on the difference between SMO and IMIA. For a full and detailed algorithm, please refer to our previous paper (Zhou et al 2017).

As IMIA, SMO consists of the seven steps: initialization, cloning, mutation, deletion, solution updating, termination and best solution selection. In initialization, model parameters R and weights w were both initialized. We generated the initial solution set $D\left( t \right)=\left\{{{d}_{1}},\cdots ,{{d}_{H}} \right\}\left(t=0 \right)~$ randomly, ${{d}_{i}}(i=1,\cdots, H)$ is a particular solution, H is the number of solutions, and $t$ is the number of generation. In cloning step, there is a big difference. In SMO, new similarity-based proportional cloning operation was proposed, where the solution with higher similarity was reproduced multiple times. Specifically, the clonal time ${\rm CL}{{{\rm T}}_{i}}$ for each solution is calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm CL}{{{\rm T}}_{i}}={{n}_{c}}\times \frac{{\rm sim}\left( {{d}_{i}} \right)}{\sum\limits_{i=1}^{H}{{\rm sim}\left( {{d}_{i}} \right)}},\nonumber \end{align} \tag{ 16 }$$

where ${{n}_{c}}$ is the expected value of the clonal solution set and $\left\lceil {} \right\rceil$ is the ceiling operator. The similarity measure for solution ${{d}_{i}}$ denoted by ${\rm sim}\left( {{d}_{i}} \right)$ is calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm sim}\left( {{d}_{i}} \right)=\underset{k=1}{\overset{K}{\mathop \sum }}\,\frac{2\sum\limits_{j=1}^{2}{P_{j}^{k}T_{j}^{k}}}{\sum\limits_{j=1}^{2}{{{\left( P_{j}^{k} \right)}^{2}}}+\sum\limits_{j=1}^{2}{{{\left( T_{j}^{k} \right)}^{2}}}},\nonumber \end{align} \tag{ 17 }$$

where $K$ is the number of training samples and $T_{j}^{k}$ is the label vector.

The mutation and deletion in SMO are the same as those in IMIA. In this paragraph, 'mutation' refers to the operation performed on the cloned solution set, not the gene mutation. The mutation probability threshold MP is determined empirically and an operation probability ${\rm R}{{{\rm P}}_{i}}$ is generated randomly. If ${\rm R}{{{\rm P}}_{i}}>{\rm MP}$, a mutation operation is performed in which a new solution $d_{i}^{m}~$ was generated randomly and replace the original solution ${{d}_{i}}$. Then, a newly generated mutated solution set $M\left( t \right)$ and solution set $D\left( t \right)$ constitute the new solution set denoted by $F\left( t \right)$. Same solutions in $F\left( t \right)~$ were removed and a new solution set ${\rm DF}\left( t \right)$ is generated.

In solution updating step of SMO, the aim is to select $H$ solutions from ${\rm DF}\left( t \right)$ to maintain the population size. For each solution, we can obtain the similarities based on the probability outputs of all the training samples, according to equation (17). ${{f}_{sim\_sen}}$ and ${{f}_{sim\_spe}}$ can also be calculated according to equation (14). Using the MOEA (Deb 2001, Deb et al 2002), we selected H solutions. Unlike most traditional MOEAs, the solution in ${\rm DF}\left( t \right)$ is sorted according to the similarity of each solution. Then, the new solution set ${\rm UD}\left( t \right)$ is generated.

When t reaches the maximal number of generations ${{G}_{max}}$, the algorithm terminates. The best solution is selected from the Pareto-optimal solution set ${\rm UD}\left( {{G}_{max}} \right)$ according to clinical needs. We selected the best solution according to the similarity-based sensitivity, specificity, and AUC. First, the thresholds T_{sim_sen} and T_{sim_spe} are determined for similarity-based sensitivity and specificity according to clinical needs. Second, for each solution ${{d}_{i}}$ in ${\rm UD}\left( {{G}_{max}} \right)$, we calculate its similarity-based sensitivity $f_{sim\_sen}^{i}$ and specificity $f_{sim\_spe}^{i}$. If $f_{sim\_sen}^{i}>~{{T}_{sim\_sen}}$ and $f_{sim\_spe}^{i}>~{{T}_{sim\_spe}}$, ${{d}_{i}}$ is selected as a candidate solution. Third, we select the solution with the highest AUC from the candidate solutions as our final solution.

The training process mainly consists of three stages: feature calculation, feature selection, and predictive model construction. To achieve optimal performance for each classifier, we adopted our multi-objective feature selection method (Zhou et al 2017). After selecting the features for each classifier, model parameters R_i and weights ${{w}_{i}}\left( i=1,2,\ldots M \right)$ were trained. The workflow is illustrated in figure 3.

Zoom In Zoom Out Reset image size

The testing process consists of three stages (figure 4). For a test sample, first, the features for each classifier are selected; second, each classifier outputs a probability $P_{i}^{\,j}\left(\,j=1,2;~i=1,2,\ldots M \right)$; third, the final mutation probability $P_{fin}^{\,j}$ is obtained by combining all $P_{i}^{\,j}$ and ${{w}_{i}}$ using the ER approach. Then the label can be determined.

Zoom In Zoom Out Reset image size

We used six different classifiers in the MCMO model, including support vector machine (SVM) (Keerthi and Lin 2003), logistic regression (LR) (Freedman 2009), discriminant analysis (DA) (Hastie and Tibshirani 1996), decision tree (DT) (Breiman 2001), K-nearest-neighbor (KNN) (Keller et al 2012), and naive Bayesian (NB) (Goldszmidt 1997). Since SVM has two model parameters and other classifiers use default parameters, we train eight parameters $\left\{{{R}_{SVM\_1}},{{R}_{SVM\_2}},~{{w}_{1}},\ldots ~{{w}_{6}}~ \right\}$ for the predictive model.

In MCMO, the population number H and the maximal generation number ${{G}_{max}}~$ were both set to 100. In the clone operation, ${{n}_{c}}$ was set to 200. In the mutation operation, the mutation probability MP was set to 0.9. The proposed MCMO predictive model was compared to three types of predictive model: (1) models with single classifiers; (2) models with different multi-objective optimization; (3) models with different fusion strategies. Because multi-objective optimization is better than the single-objective model (Zhou et al 2017), we did not compare our proposed model to models with single-objective optimization. Area under the receiver operating characteristic curve (AUC), accuracy, sensitivity, specificity, similarity-based sensitivity (Sim-sensitivity), and similarity-based specificity (Sim-specificity) were used for evaluating the model performance.

In our study, eight parameters in the predictive model need to be trained and our dataset has 57 cases. We adopted two-fold cross-validation in this work. Cross-validation is a widely used model validation technique which can test the model's ability to predict new data that were not used in training, in order to flag problems like overfitting (Kohavi 1995). One round of cross-validation involves partitioning a sample of data into complementary subsets. In our work, two-fold cross-validation was used, where for each round, half cases (training set) were selected randomly for training and the other half cases (validation set) were used for validation, then reverse. To reduce variability caused by subset partition, ten rounds of two-fold cross-validation were performed for each model, and the validation results were averaged over the ten rounds to give an estimate of the model's predictive performance (table 3). Prediction results of AUC, accuracy, sensitivity and specificity of training set were also listed in table 3. The prediction accuracy of training set is higher than that obtained from the validation set, which is considered as normal for a machine learning algorithm. On the other hand, the accuracies of the training set are not close to 1 and the predictive results on the validation set are acceptable. In the remaining of the paper, all the prediction results were from the validation set.

We obtained the classifier-specific feature subset and summarized the features selected by all (six) classifiers (table 4). Statistics using the unpaired t-test were also listed. The feature selected by all classifiers indicates that the feature is important to predict mutation. P-value smaller than 0.05 indicated significant differences between presence and absence of mutations. However, four selected features (Minimum, Contrast, Variance and Sum variance) for the PBRM1 gene and three selected features (Variance, Sum average and Sum variance) for the BAP1 gene had P-values greater than 0.05. Because we selected the optimal features according to multi-objective model, using the AUC as the figure of metric, these selected features did not necessarily have P-values smaller than 0.05. Example boxplots were plotted to show the potential of a single feature for differentiating between the presence and absence of mutations (figure 5).

Zoom In Zoom Out Reset image size

Two intensity features including Mean and Kurtosis are the most frequently selected features in VHL gene prediction. A histogram with a more elongated tail indicates smaller Kurtosis. A tumor with smaller Kurtosis is more likely to carry a VHL mutation (figures 5(a) and 6). For the PBRM1 gene prediction, nine features from intensity and texture features, were most frequently selected. A boxplot of the Mean is illustrated in figure 5(b). Six features were selected as the most prominent contributors in BAP1. Four texture features were selected as follows: Homogeneity, Variance, Sum average, and Sum variance. A lower similarity in intensity between a voxel and its neighbors led to higher Variance and Sum variance. A less uniform or more focal intensity distribution led to reduced Homogeneity. Therefore, larger Variance, and smaller Homogeneity were associated with the likelihood that a tumor carried a BAP1 mutation. Boxplot of Homogeneity is illustrated in figure 5(c).

Zoom In Zoom Out Reset image size

It is noted that one single feature may not achieve accurate predictive results. For each classifier, a feature set is necessary. For KNN classifier, a feature set consisting of 12 features (Volume, Size_P3, Minimum, Maximum, Mean, Sum, Variance, Standard deviation, Kurtosis, Cluster shade, Energy, Inverse difference) was selected to predict VHL mutation; while for LR classifier, a feature set consisting of 19 features (Volume, Size_Z, Size_P2, Roundness, Surface area, Mean, Standard deviation, Skewness, Kurtosis, Contrast, Dissimilarity, Energy, Entropy, Homogeneity, Sum entropy, Information measure of correlation_1, Information measure of correlation_2, Inverse difference, Inverse difference normalized) was selected to predict VHL mutation.

MCMO yielded better AUC, accuracy, sensitivity, and specificity results than other single classifiers (table 5). The prediction accuracy of MCMO is 0.81, 0.78, and 0.90 for VHL, PBRM1, and BAP1 genes, respectively, with AUC  ⩾  0.86 sensitivity  ⩾  0.75 and specificity  ⩾0.80. MCMO yielded better results than other single classifiers. KNN is the best single classifier for the VHL and PBRM1 genes, with specificities of 0.66 and 0.62, respectively. MCMO can achieve specificities of 0.86 and 0.80 for VHL and PBRM1 genes, respectively. SVM and DA achieved similar results for the BAP1 gene, which are better than other single classifiers, but sensitivities were only 0.57 and 0.63, respectively; the sensitivity obtained by MCMO was 0.87. Some single classifiers achieved higher sensitivities for the VHL gene, but the corresponding specificities were poor. MCMO achieved the highest AUC and accuracy with balanced sensitivities and specificities (difference  <  0.1). Also, MCMO is a stable predictive model because its standard deviations are much smaller than those of the single classifiers.

One group of the Pareto-optimal solution set and the selected final solution in SMO is shown in figure 7. As described in section 2.3, the best solution was selected according to the similarity-based sensitivity, specificity (equation (14)), and AUC. First, thresholds T_{sim_sen} and T_{sim_spe} were determined for similarity-based sensitivity and specificity based on clinical needs. In this study, the thresholds T_{sim_sen} and T_{sim_spe} are both 0.9. The selected candidate solutions were included within the red rectangle and the selected final solution (highest AUC) was marked in red.

Zoom In Zoom Out Reset image size

We evaluated the performance of our MCMO by comparing it to the iterative multi-objective immune algorithm (IMIA), which adopts the traditional sensitivity and specificity as the optimized objective functions (Zhou et al 2017) (figure 8). The two methods were compared with the unpaired t-test at a significance level 0.05 (table 6). Results were similar based on AUC, accuracy, sensitivity and specificity (P-value  >  0.05). For VHL and PBRM1 genes, SMO achieved a little higher AUCs. For all three genes, SMO achieved significantly higher similarity scores (P-value  ⩽  0.01), indicating that these results are more reliable. For the prediction result with higher AUC, the difference of ${{f}_{sim\_sen}}$ or ${{f}_{sim\_spe}}~$ between SMO and IMIA is small. For example, for the BAP1 gene, the difference of ${{f}_{sim\_sen}}$ and that of ${{f}_{sim\_spe}}$ are 0.03 and 0.02. However, for the prediction result with lower AUC, such as for the PBRM1 gene, the difference of ${{f}_{sim\_sen}}$ is 0.13 and that of ${{f}_{sim\_spe}}$ is 0.07.

Zoom In Zoom Out Reset image size

We used the ER approach (equation (4)) for fusing the output of different classifiers. The classic weighted fusion (WF) method is used for comparison, as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P= ~ \underset{i=1}{\overset{M}{\mathop \sum }}\,{{P}_{i}}{{w}_{i}}\nonumber \end{align} \tag{ 18 }$$

where ${{P}_{i}}~$ is the individual classifier output probability and ${{w}_{i}}$ is the relative weight. SMO is used in both fusion strategies, and the comparative results are shown in figure 9. The two methods were compared with the unpaired t-test at a significance level 0.05 (table 7). For VHL and PBRM1 genes, the ER approach achieved higher AUCs (P-value  <  0.05). Also, sim-sensitivity and sim-specificity in ER are higher than WF (P-value  ⩽  0.02), which indicates that more reliable results can be obtained when using ER fusion.

Zoom In Zoom Out Reset image size