A hybrid feature selection algorithm and its application in bioinformatics

The mRMR algorithm

The mRMR (Peng, Long & Ding, 2005) algorithm, which uses mutual information to assess the relevance between features, has been used in bioinformatics (Ding & Peng, 2005; Li et al., 2012; Mundra & Rajapakse, 2010). Mutual information is widely used to analyze the correlation between two variables, and it can be expressed as Eq. (1). (1)${\displaystyle I\left(X,Y\right)=\sum _{x\in X}\sum _{y\in Y}p\left(x,y\right)log\frac{p\left(x,y\right)}{p\left(x\right)p\left(y\right)}.}$

In the Eq. (1), P represents the probability and X, Y represent the feature vector or class vector. The relevance V and redundancy W of the mRMR can be expressed using Eqs. (2) and (3). (2)${\displaystyle V=\frac{1}{\left|S\right|}\sum _{x_i\in S}I\left(y;x_i\right).}$ (3)${\displaystyle W=\frac{1}{{\left|S\right|}^2}\sum _{x_i,x_j\in S}I\left(x_i;x_j\right).}$

In the Eqs. (2) and (3), y is the target variable, S is candidate feature set and x_i, x_j are arbitrary variables of S. To calculate the final score of relevance, the MIQ can be used, as shown in Eq. (4). (4)${\displaystyle MIQ:argmax\left(\frac{V}{W}\right).}$

Maximal information coefficient

As a measure of dependence for two-variable relationships in a large dataset, MIC has been widely used in various fields, including global health, gene expression, human gut microbiota and identify novel relationships due to its ability to capture a wide range of functional and non-functional associations. The definition of MIC is: (5)${\displaystyle MIC\left(x,y|D\right)=ma{x}_{i\ast j<B\left(n\right)}\left\{\frac{I^{\ast }\left(x,y,D,i,j\right)}{logmin\left(i,j\right)}\right\}.}$

In the Eq. (5), x and y are the pairs(x, y) of the dataset and I*(x, y, D, i, j) denotes the maximum mutual information of D—_G with the i-by-j grid, where the default B(n) = n^0.6.

Particle swarm optimization algorithm

The particle swarm optimization (PSO) algorithm is a heuristic search algorithm that originated from studies of bird predation behavior (Eberhart & Kennedy, 1995). The first step of PSO is to initialize a group of particles. It then iterates until it finds the best solution. The particles update themselves in each iteration by tracking two extreme values. The first is the particle’s determination of the individual extreme value Pbest. The other is called Gbest, which is determined by the entire particle swarm. After determining the Gbest and Pbest values, the particle updates its speed and position based on Eqs. (6) and (7). (6)${\displaystyle V_{k+1}=\omega V_k+c_1{r}_1\left(Pbest-X_k\right)+c_2{r}_2\left(Gbest-X_k\right).}$ (7)${\displaystyle X_{k+1}=X_k+V_k.}$

In the above equations, k represents the number of iterations; V_k and X_k represent the particle’s current velocity and position, respectively; r₁, r₂ are random values between [0, 1]; c₁, c₂ are the learning factors; ω is the inertia weight, which is used to control the influence of the last iteration’s speed on the current speed. A smaller and larger ω can strengthen the PSO algorithm’s local or global search ability, respectively.

The hybrid algorithm for feature selection

In this study, we proposed the MMPSO hybrid method. First, the dataset needed to be preprocessed. On the one hand, the aim of preprocessing is to remove some features that contain a large quantity of noisy data, such as features that contain many zeros. On the other hand, if the proportion of samples is clearly unbalanced, it is necessary to balance the samples. Here, we employed random oversampling technology to address this issue. A random over-sampler randomly copies and repeats the minority class samples, eventually resulting in the minority and majority classes having the same number. The next step was to rank the features. MIC was used to measure the correlation of two features, resulting in a more accurate ranking result of features based on the mRMR framework (Cao et al., 2021). Considering the numerous features and the complexity of MIC, we used the multithreading method in paper (Tang et al., 2014) to speed up the calculation. After performing the mRMR based on MIC method, we obtained the ranking features and use the top K features as the input of the next step to reduce the computational load for the PSO. The K features were used in the third step to initialize the particle swam and calculate the fitness of each particle. For the wrapper feature selection algorithm, only the classification accuracy is used as the fitness function to guide the feature selection process, which will lead to a larger scale of the selected feature subset (Liu et al., 2011). Therefore, some studies combine the classification accuracy and the number of selected feature subsets to form a fitness function (Xue, Zhang & Browne, 2012). Here, the fitness we defined is shown in Eq. (8). V_error was the error, which is measured by a classification method of k-nearest neighbor (KNN) (Chen et al., 2021a). N_selected and N_all were the numbers of selected features and the entire features, respectively. α and β were parameters whose sum is 1. The larger α is, the more features will be chosen; otherwise, fewer features will be chosen. When the specified number of iterations is reached, the PSO program terminates, and the final selected features will be available. (8)${\displaystyle cost=\alpha V_{error}+\beta \frac{N_{selected}}{N_{all}}.}$

The validation method

Here, we respectively compared the classification accuracy of the MMPSO method with the results of other algorithms, including mRMR (Peng, Long & Ding, 2005), ILFS (Roffo et al., 2017), ReliefF (Liu & Motoda, 2007), Mutinffs (Zaffalon & Hutter, 2002), FSV (Bradley & Mangasarian, 1999), Fisher (Gu, Li & Han, 2012), CFS (Zeng & Cheung, 2010), UFSOL (Guo et al., 2017), to demonstrate that our method has better classification accuracy. LIBSVM (Chang & Lin, 2011) is an integrated library, which supports multi-class classification. Here, we performed classification using LIBSVM to test the accuracy with k-fold cross validation.

Summary of datasets

A total of eleven datasets were used in this study; the basic information about the datasets was shown in Table 1. The UCI Machine Learning Repository is a collection of databases, domain theories, and data generators that are used by the machine learning community for the empirical analysis of machine learning algorithms (Dua & Graff, 2017). In the beginning, ten datasets that have different numbers of features, instances and classes were downloaded from the UCI website, and they were used to evaluate the performance of our proposed method. Furthermore, we conducted a more thorough analysis to demonstrate the biological application of our method. Tremendous amount of RNA expression data has been produced by large consortium projects such as TCGA and GTEx, creating new opportunities for data mining and deeper understanding of gene functions (Tang et al., 2017). Thus, the final liver hepatocellular carcinoma (LIHC) dataset was obtained from UCSC Xean, which contains large-scale standardized public, multiomic and clinical/phenotype information (Goldman et al., 2020). The LIHC dataset used in the current study contains RNA expression data of over 60,000 genes in 531 biosamples (371 tumor samples and 160 normal samples, and the latter further containing 50 normal samples from the TCGA-LIHC cohort and 110 normal tissues from GTEx), and it is available at https://toil-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA-GTEx-TARGET-gene-exp-counts.deseq2-normalized.log2.gz.

Results of the experiment based on ten datasets

Figure 1 and Table 2 summarized the classification accuracy on the basis of the MMPSO and the compared methods. Here, we defined the threshold K was 100. When the number of original features was greater than 100, the top 100 features from the ranking result were selected as the PSO input for the MMPSO method; otherwise, all features were selected as the input. In Fig. 1, we obtained the conclusion that our method was superior to the other methods in terms of classification accuracy based on six datasets including gisette, hillvalley, isolet, madelon, scene and usps. In the other four datasets, the MMPSO method achieved the accuracies of top three rankings. Therefore, the MMPSO algorithm was outperformed other methods with respect to accuracy of classification by utilizing the same number of features.

Results of the experiment based on the biological dataset

After analyzing the first ten datasets, we employed the MMPSO method on the LIHC dataset to identify features (gene biomarkers) that can be used to distinguish the tumor group from the normal group with high accuracy. Different from the previous datasets, LIHC is an unbalanced dataset containing 371 tumor samples and 160 normal samples. Therefore, the preprocessing was needed and the number of genes was reduced from over 60,000 to 15,185 after preprocessing. When the genes were ranked by the mRMR based on MIC method, the top 100 genes were selected and input into PSO to identify the best gene signatures. Finally, we obtained a signature of 18 genes through the MMPSO algorithm, that had better classification compared to other methods. The 18 genes were ACTN1, CACHD1, ERAP2, FAM171A1, FKBP1B, HIST1H2BC, PLVAP, PRKD1, RPL7AP6, ADRA2B, DMKN, FNDC4, NPC1L1, POMC, PYROXD2, RBP1, TRIM29, and ZBED9, with the relative expression of the first nine genes significantly increasing in tumors and the last nine genes decreasing (P < 0.01 in the Wilcoxon rank sum test with continuity correction, Fig. 2). Principal component analysis (PCA) was then performed using the ‘FactoMineR’ and ‘Factoextra’ packages in R version 4.0.2 based on the expression profiles of the candidate 18 genes. As shown in Figs. 3A–3B, Dim 1 and Dim 2 were 15% and 11.9%, respectively. Figure 3C illustrated the heatmap of all samples based on the 18 gene expression profiling. The results revealed that the 18-gene signature obtained from MMPSO algorithm could effectively separate the 531 samples into two groups. Since logistic regression show that seven biomarkers of ADRA2B, ERAP2, NPC1L1, PLVAP, POMC, PYROXD2 and TRIM29 were significantly associated with Wald and P value, as shown in Table 3. We further investigated the combined diagnostic efficacy of the seven candidate genes according to the Eq. (9).

In the above Equation, PP was the functional formula for predicting the incidence of LIHC, i.e., PP-value, and the constant and coefficient were the result of logistic regression in Table 3. The results of receiver operating characteristic (ROC) curve analysis using MedCalc software based on the value of PP and classification labels of LIHC dataset was shown in the Fig. 4, which had an area under the curve (AUC) greater than 0.999 and P < 0.0001. It demonstrated that the seven genes significantly distinguished tumors from normal samples in LIHC dataset.

To explore whether the 18 genes are associated with survival time of phenotype information in LIHC dataset, the Kaplan–Meier (KM) survival curve was performed by using the “survival” and “survminer” packages in R. For each gene, the cut-off points obtaining from “survminer” package then divided gene expression values into the high (high) and the low (low) groups. We identified that higher expression levels of CACHD1, FKBP1B, PRKD1, and RPL7AP6 were associated with worse overall survival (OS) time, whereas higher expression levels of ADRA2B, PYROXD2 were associated with better OS, as shown in Fig. 5.