An improved MOPSO approach with adaptive strategy for identifying biomarkers from gene expression dataset

1.   Introduction

Cancer is recognized as one of the serious challenges that endanger human health due to its high incidence and mortality. Because the patient does not show obvious symptoms, cancer is not easy to be detected in the early stage, which result in a poor prognosis. Therefore, it is urgent to use scientific and effective methods to identify biomarkers of cancer for early diagnosis and treatment. With the rapid progress of sequencing technology, research in the field of cancer treatment from the gene level has become a hot topic.

Through microarray technology, genes associated with disease can be found in a short time and used as biomarkers for early diagnosis, but due to the existence of some redundant features and inherent noise in microarray data, the selection of feature genes still face great challenges ^[1,2]. Some researchers think that the generation and development of cancer are related to some specific genes, and have proposed several methods to find genes related to different stages of cancer development ^[3,4], however, due to the lack of detailed biological processes, the results obtained by gene expression-based methods cannot be proved to be completely correct ^[5]. Hence, pathway-based methods for feature gene extraction were proposed to obtain biomarkers from microarray data, which helps to better understand the differences between phenotypes ^[6,7]. Some scholars have used particle swarm optimization to infer pathway activity and achieved good results, a popular pathway activity index t-score was chosen as the objective function in ^[8]. ^[9] proposed an approach based on binary particle swarm optimization, in this approach, partial genes are chosen automatically from microarray gene data for inferring pathway activity and the remaining genes are excluded. Reference ^[10] developed a multi-objective particle swarm optimization technology that using protein interaction scores and two improved indicators of pathway activity as objective functions to infer pathway activity, all genes are involved in pathway activity inference with different weights in this method.

Multi-objective particle swarm optimization algorithm has become the main research direction of multi-objective optimization due to the advantages of high efficiency and speed. However, due to the lack of selection pressures, feature solutions distribute around the Pareto front, which will result in the lack of diversity of optimal solution set ^[11]. Some researchers use the space decomposition strategies to maintain diversity, which divide the target space into multiple regions, such as grid division method ^[12] and angle division method ^[13]. some authors put forward the idea of population decomposition, which divide the entire population into multiple subgroups, and the subgroups guide the population update in parallel according to their respective leaders ^[14,15]. But fixed selection pressures provided by decomposition approaches unable achieve ideal results when dealing with complex Pareto fronts problems ^[16]. Reference ^[17] developed an adaptive feature selection approach, which adjust the selection pressure of archive and particles by adjusting the parameters in the decomposition method.

Motivated by the above researches, in this article, a pathway-based multi-objective particle swarm optimization algorithm with adaptive strategy for feature gene extraction is proposed. Two weighted pathway activity indicators are used as objective functions, and find feature solutions by optimizing them simultaneously. In addition, in order to improve the diversity of feature solutions, the adaptive penalty boundary intersection decomposition method is introduced in the process of optimization.

The remainder of this paper is organized as follows. Section 2 provides related methods contained a general description of multi-objective optimization (MOO), particle swarm optimization, and the PBI decomposition approach. The proposed adaptive feature selection method is elaborated in Section 3. In Section 4, the detailed introduction and the preprocessing of datasets are given, and the experimental results of the proposed IMOPSO-PBI approach are presented and compared with 6 existing methods. Finally, concluding remarks are given in Section 5.

2.   Related methods

2.1.   MOO

In a multi-objective optimization problem, multiple objective functions need to be optimized simultaneously. The optimization problem is not only used in the real-time design of production scheduling, urban transportation, network communication and other systems, but also involves intelligent planning problems such as engineering design, data mining, and capital budgeting ^[18]. The mathematical definition of multi-objective optimization is as Eq (1).

where $x$ is a vector in an n-dimensional decision space, and ${f}_{i}\left(x\right)$ is the objective function of the multi-objective optimization problem, $m$ is the scale of the objective function.

In recent years, evolutionary algorithms have integrated biological information into meta-heuristic algorithms. With its unique update mechanism, many breakthrough research results have been achieved in the fields of combinatorial optimization and numerical optimization ^[19]. Multi-objective evolutionary algorithms include: multi-objective particle swarm algorithm ^[20], multi-objective bee colony algorithm ^[21], multi-objective ant colony algorithm ^[22], multi-objective immune algorithm ^[23], multi-objective differential algorithm ^[24], etc.

2.2.   Particle swarm optimization

Particle swarm optimization is an evolutionary computing technology whose basic idea is to find the optimal solution through the cooperation and information sharing among individuals in the group. The individuals in the population are abstracted as particles, and the particles are affected by the combined effects of themselves and the state of the population at each iteration ^[25]. The velocity of the $i$ th particle at the next moment is determined by the current velocity, the personal best position $\left(pbest\right)$ and the global best position $\left(gbest\right)$, and the particle moves from the current position to the new position at the updated velocity. The updating process of particle velocity and position is shown in Eqs (2) and (3).

where $\omega$ represent inertia weight, ${r}_{1}$, ${r}_{2}$ are the random numbers between 0 and 1, ${c}_{1}$, ${c}_{2}$ are learning factors generally set to 2.

2.3.   PBI approach

Decomposition based multi-objective evolutionary algorithm decomposed the target into a set of scalar optimization subproblems firstly, and then optimized these subproblems simultaneously. In the existing decomposition methods, weighted sum (WS) is limited to dealing with convex Pareto front problems, for non-convex Pareto frontiers, optimal solutions cannot be obtained completely. Tchebycheff approach (TCH) overcome that problem, but for a continuous multi-objective optimization problem, its aggregation function is not smooth, and the resulting Pareto front is not smooth too ^[26]. Compared with the above two decomposition methods, PBI approach can obtain more evenly distributed solutions ^[27]. The mathematical expression of PBI decomposition approach is as Eq (4).

where 𝜆 is the predefined weight vector, and 𝜃 is the penalty parameter, and ${d}_{1}$, ${d}_{2}$ are shown in Eq (5).

As shown in Figure 1, ${z}^{\mathrm{*}}$ is the ideal reference value, the red line is the Pareto frontier (PF), the point 𝑃 is the projection of $f\left(x\right)$ on the line $L$, and the intersection of PF and the line $L$ is the optimal solution under the decomposition weight vector λ. ${d}_{1}$ represent the distance between the point 𝑃 and the ideal point ${z}^{\mathrm{*}}$, and ${d}_{2}$ represent the distance between $f\left(x\right)$ and the line $L$. It can be seen that ${d}_{1}$ determines the degree of convergence of the population. The smaller ${d}_{1}$ is, the closer its corresponding solution $x$ is to the ideal PF, and ${d}_{2}$ with a penalty parameter determines the diversity of the population. Figure 2 shows the contours of $g$ for three different penalty parameters, as can be seen from the figure, the smaller the penalty parameter 𝜃, the larger the area of the update area, the more favorable the population is to approach the ideal PF. The larger the penalty parameter 𝜃, the smaller the area of the update area, and the closer the updateable solution is to the weight vector, which is more conducive to the diversity of the population. That is to say, the penalty parameter 𝜃 plays an important role in the performance of penalty boundary interaction methods.

Many scholars set the parameter to a fixed value of 5, however, it is difficult to select the solution with better convergence and diversity with a fixed penalty value. In ^[28], the authors introduced a subproblem-based penalty mechanism and an adaptive penalty mechanism to analyze the influence of different penalty values on the feature solution set, but this method did not consider the evolution state of each subproblem. ^[29] proposed a method to adjust the penalty value by the angle between the solution and the weight vector, but this increases the time complexity of the algorithm. In ^[17], an adaptive penalty mechanism based on PBI parameters is proposed to reduce the time complexity of calculating penalty values while maintaining diversity.

3.   The proposed IMOPSO-PBI algorithm

This section mainly describes the proposed pathway-based adaptive feature gene selection algorithm. Firstly, the overall flow of the proposed algorithm is summarized and the computational complexity is described. For more comprehensibility, details about each component of IMOPSO-PBI method are described in the following sections. The process of inferring pathway activity and particle coding is briefly introduced; two objective functions are described, and the adaptive updating mechanism of penalty value based on penalty boundary interaction is introduced.

3.1.   Framework of the proposed IMOPSO-PBI algorithm

The flowchart of our proposed IMOPSO-PBI method is shown in Figure 3, and the whole process of IMOPSO-PBI method is shown in Table 1. Firstly, the original datasets are preprocessed, the top 1000 feature genes are selected and the corresponding pathway has been obtained through the David database. Then, the pathway is encoded according to the encoding strategy introduced in Section 3.2, and an initial population containing $S$ particles has been obtained. Next, the population is initialized, the position of each particle in the population is randomly initialized between 0 and 1, the target value of each particle is calculated through the objective functions ${f}_{1}$ and ${f}_{2}$, and the speed and position of the particle are determined by Eqs (1) and (2) to update. The leader archive is updated through Algorithm 1, and the penalty value is also updated during the archive update process, and the crowding distance between particles is calculated. Finally, the solution in the lead archive is the desired result.

3.2.   Inferring pathway activity and encoding strategy

At first, the top 1000 genes taken from the preprocessed data (described in Section 4.1) are put into the DAVID Bioinformatics Resources ^[30], then we can obtain the corresponding KEGG pathway. Because each pathway contains several genes, so the pathway activity can be calculated by Eq (6).

Here, M is the number of genes in the pathway, and N is the number of samples of each gene.

After calculating the pathway activity, we encode particles based on the pathway. As shown in Figure 4, particle ${P}_{k}$ contains three pathways, ${p}_{1} = \left\{{g}_{1}^{1}、{g}_{2}^{1}、{g}_{3}^{1}\right\}$, ${p}_{2} = \left\{{g}_{1}^{2}、{g}_{2}^{2}、{g}_{3}^{2}\right\}$, ${p}_{3} = \left\{{g}_{1}^{3}、{g}_{2}^{3}\right\}$, ${g}_{i}^{j}$ represents the $i$ th gene on the $j$ th pathway. The numerical value on each cell represents the degree of relevance of that gene in inferring pathway activity.

3.3.   Objective functions

In this paper, weighted t-score and weighted z-score are selected as the objective function of the multi-objective optimization algorithm. In general, the t-score is used to measure the ability to differentiate the cumulative expression of the constituent genes of a given pathway ^[31]. The mathematical expression of the t-score is as follows.

where ${p}_{i}$ indicates the pathway activity level of a given pathway, $\mu \left(x\right)$ and $\mu \left(y\right)$ represent the mean of pathway activity for classes ${C}_{1}$ and ${C}_{2}$ respectively, $\sigma \left(x\right)$ and $\sigma \left(y\right)$ represent the standard deviation of pathway activity for classes ${C}_{1}$ and ${C}_{2}$ respectively. ${s}_{1}$ and ${s}_{2}$ represent the number of samples in the two classes. In this paper, the concept of weighted t-score is introduced, which assumes that all genes in the pathway participate in the inference of pathway activity with a certain weight. The formula of weighted t-score is defined as Eq (8).

where ${\mu }_{\omega }\left(x\right) = \frac{{\sum }_{i = 1}^{{s}_{1}}~~~\left({x}_{i}^{k}*{\omega }_{i}\right)}{{\sum }_{i = 1}^{{s}_{1}}~~~{\omega }_{i}}$, ${\mu }_{\omega }\left(y\right) = \frac{{\sum }_{i = 1}^{{s}_{2}}~~~\left({y}_{i}^{k}*{\omega }_{i}\right)}{{\sum }_{i = 1}^{{s}_{2}}~~~{\omega }_{i}}$ represent the weighted mean of pathway activity for classes ${C}_{1}$ and ${C}_{2}$ respectively, and the weighted standard deviation of pathway activity for classes ${C}_{1}$ and ${C}_{2}$ are described as

Suppose a particle consists of n pathways, then the weighted t-score of the particle should be expressed as Eq (9).

The higher the weighted t-score is, the greater the differentiation ability is. Therefore, the weighted t-score of the particle should be maximized. The proposed algorithm is designed to be a minimization problem, so the first objective function can be expressed as Eq (10).

Weighted z-score is chosen as another objective function in this paper, the z-score is a measure of the distance between a data point and its overall mean, a positive z-score means that the data point is greater than the weighted mean, and a negative z-score means it is less than the weighted mean. The weighted z-score of a particle is described as Eq (11).

in the formula, ${\alpha }_{w}\left({P}_{k}\right)$, ${\mu }_{w}\left({P}_{k}\right)$, ${\sigma }_{w}\left({P}_{k}\right)$ represent the pathway activity, weighted mean and weighted standard deviation of particle ${P}_{k}$ respectively, here ${\alpha }_{w}\left({P}_{k}\right) = \frac{\sum _{i = 1}^{\mathrm{n}}~~~\alpha \left({p}_{i}\right)}{\mathrm{n}}$. The smaller the absolute value of the z-score, the closer it is to the overall mean, so the second objective function can be expressed as Eq (12).

3.4.   The mechanism for archive updating

Since the capacity of the archive is limited, it is necessary to update the archive in real time during the feature selection process. This paper updates the archive by adaptively changing the penalty value of the PBI parameter. According to Section 2.3, we can know that when the value of ${d}_{2}$ corresponding to the current weight vector is small, increasing the penalty value helps to obtain a solution closer to the direction of the current weight vector. Therefore, this section proposes a method that change the penalty value of ${d}_{2}$ adaptively. When the ${d}_{2}$ value of the solution obtained under a specific weight vector in the archive increases, the penalty value decreases exponentially. Its specific mathematical expression is as Eq (13).

where, ${d}_{i}$ is the value of ${d}_{2}$ of the corresponding weight vector at the $i$ th iteration, the outer part of the formula is set to 1 to maintain the original balance, the coefficient 𝑁 of ${d}_{i}$ is set to 30. Moreover, in order to avoid too small d value to make the convergence worse, an adaptive penalty value that varies with the number of iterations is introduced, the 𝑔𝑒𝑛𝑥 value is added to Eq (13) to ensure that the penalty value varies within the normal range, 𝑔𝑒𝑛𝑥 is expressed as Eq (14).

To satisfy the condition that the initial penalty value should be greater than 1, the value of ${\mathrm{\theta }}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ should be set to 1, and for the value of ${\mathrm{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, 100 is a good choice.

The process of updating archive and penalty values is shown in Table 2. Firstly, the optimal solution is selected from the former leader archive under the current weight vector. Then, remove the optimal solution from the former leader archive and add it to the updated archive. At last, calculate the ${d}_{2}$ value corresponding to the current weight vector and update the penalty value.

3.5.   Computational crowding distance

The set of non-dominated solutions distributed on the Pareto front is called the pareto optimal set, the optimal solutions can be stored in the archive. Here, the process of updating the archive and calculating the crowding distance is described. Assume that the maximum number of non-dominated solutions that the archive can accommodate at each iteration is N, suppose that by calculating the fitness in each iteration, M non-dominated solutions are obtained, if $M < N$, keep the optimal solution in the archive, else if $M > N$, calculate the crowding distance for all solutions and add solutions with higher crowding distances to the updated archive.

4.   Experiments results

4.1.   Datasets and preprocessing

Six real gene expression datasets were used in this paper, which are available from the public website (www.biolab.si/supp/bi-cancer/projections/info/).

Prostate: Gene expression measurements for samples of prostate tumors and adjacent prostate tissue not containing tumor were used to build this classification model. The number of genes in this dataset are 12, 533 and the number of samples are 102, including 50 normal tissues and 52 prostate tumor samples.

DLBCL: Diffuse large B-cell lymphomas (DLBCL) and follicular lymphomas (FL) are two B-cell lineage malignancies that have very different clinical presentations, natural histories and response to therapy. The number of genes in this dataset are 7070 and the number of samples are 77, including Diffuse large B-cell lymphoma 58 examples and Follicular lymphoma 19 examples.

GSE412(Child-ALL): The childhood ALL data set (GSE412) includes gene expression information on 110 childhood acute lymphoblastic leukemia samples. The number of genes in this dataset are 8280 and the number of samples are 110, including Diffuse large B-cell lymphoma 50 examples and Follicular lymphoma 60 examples.

GSE2535 (Chronic Myeloid leukemia Treatment): Imatinib induces complete cytogenetic response in the majority of patients with chronic myeloid leukemia (CML) in chronic phase. The number of genes in this dataset are 12, 625 and the number of samples are 28, including non-responder to imatinib treatment 12 examples and responder to imatinib treatment 16 examples.

GSE2443 (Prostate Cancer Treatment): After an initial response to androgen ablation, most prostate tumors recur, ultimately progressing to highly aggressive androgen-independent cancer. The number of genes in this dataset are 12, 627 and the number of samples are 20, including androgen dependent tumor 10 examples and androgen - independent tumor 10 examples.

GSE116959 (Lung adenocarcinoma): Based on GPL17077 platform(Agilent-039494 Sure Print G3 Human GE v28x60K Microarray 039381), the number of samples are 68, including normal Lung tissue 11 examples and LUAD 57 examples.

The above original datasets can be obtained as matrix format, whose columns are genes and rows are samples. firstly, the signal-to-noise ratio for each gene(column) needs to be calculated by Eq (15).

Compute the mean and standard deviation (S.D.) of classes 1 and 2 and put them into signal-to-noise equation, and then according to the calculated signal-to-noise ratio, the genes (column) are sorted in descending order. A high SNR indicates that these genes are distributed over a wide range of values, some genes with low signal-to-noise ratio may be considered unimportant for class labels, so we take the top 1000 genes and used them for further experiments. In order to eliminate the adverse effects caused by singular sample data, normalize the filtered data using min-max normalization method.

4.2.   Related parameter settings

In this section, the proposed algorithm is compared with some existing algorithms, to ensure the objectivity and fairness of the comparison results, the parameters of this experiment are consistent with the parameter settings in the existing algorithm. The specific experimental parameters are listed in Table 3.

4.3.   Performance metrics

Performance metrics play an important role in the process of verifying and comparing algorithm performance. Six classic performance metrics are chosen in this paper, they are sensitivity, specificity, F-score, accuracy, G-mean and AUC, and the expressed as follows respectively.

$tp$, $tn$, $fp$, $fn$ are represented true positive, true negative, false positive and false negative, respectively. In order to avoid the contingency of experimental results and improve the reliability, 10-fold cross-validation method was used in this experiment. Each dataset is divided into ten parts, of which nine parts are used for training and the remaining one is used for testing, and the results are demonstrated in Tables 4–9. In addition, the Pareto fronts of the proposed algorithm and existing algorithms on the experimental datasets are shown in Figure 5(a)–(f), respectively.

4.4.   Comparative analysis

In order to verify the performance of the proposed algorithm, some existing algorithms are selected for comparison with it. For example, a method for identifying genes using protein interaction scores ^[10], the sequential forward search method ^[5], T-test ^[32], feature selection algorithm based on feature clearness ^[33], a correlation-based feature selection method ^[34], and a multi-objective particle swarm algorithm ^[3]. Next, we will analyze the comparison results.

Table 4 shows the performance comparison of different methods for Prostate dataset, the sensitivity, specificity, F-score, accuracy, G-mean and AUC of the proposed IMOPSO-PBI algorithm are 0.96275, 0.9758, 0.96239, 0.96472, 0.9693, 0.9728 respectively, it is obvious that the performance of the proposed IMOPSO-PBI algorithm is better than that of the comparison algorithm. And as Figure 5(a) shows, when ${f}_{2}$ is the same, the non-dominant solution obtained by the proposed IMOPSO-PBI algorithm has a smaller ${f}_{1}$. Table 5 shows the performance comparison of different methods for DLBCL dataset, the sensitivity, specificity, F-score, accuracy, Gmean and AUC of the proposed IMOPSO-PBI algorithm are 0.95376, 0.98572, 0.92592, 0.97582, 0.9696, 0.9794 respectively, better than that of the comparison algorithm. And as Figure 5(b) shows, no matter for ${f}_{1}$ or ${f}_{2}$, the proposed IMOPSO-PBI algorithm has better performance. Table 6 shows the performance comparison of different methods for GSE412 (Child-All) dataset, the sensitivity, specificity, F-score, accuracy, Gmean and AUC of the proposed algorithm are 0.9236, 0.98674, 0.89742, 0.90691, 0.9546, 0.9462 respectively, the specificity of proposed IMOPSO-PBI method is slightly lower than the MOPSO-PPI algorithm, but the difference can be ignored. And as Figure 5(c) shows, the proposed IMOPSO-PBI method can get a better Pareto front.

For GSE2535 (Chronic Myeloid Leukemia Treatment) dataset, the performance comparison results are shown in Table 7, the sensitivity of the proposed IMOPSO-PBI method, MOPSO-PPI and the MOPSO algorithm are the same value as 1, the proposed IMOPSO-PBI algorithm has better performance on other metrics, and it can be seen from Figure 5(d) that most of the solutions on the Pareto front of the proposed IMOPSO-PBI method outperform other algorithms. For GSE2443 (Prostate Cancer Treatment) dataset, Table 8 describes the outcome of the proposed IMOPSO-PBI method and other techniques, the specificity of CBFS higher than the proposed IMOPSO-PBI method, apart from this, the proposed IMOPSO-PBI method has a better performance on other metrics, the proposed IMOPSO-PBI algorithm has a better Pareto front, which can be seen from the Figure 5(e). For GSE116959 (Lung adenocarcinoma) dataset, although the Pareto front of MOPSO-PPI locally is better than the proposed IMOPSO-PBI as shown in Figure 5(f), our proposed algorithm is better on the whole. From Table 9, we can know that the proposed IMOPSO-PBI algorithm has better performance on the metrics. In general, the proposed IMOPSO-PBI algorithm has better performance than the comparison algorithms on all experimental data, which proves that the IMOPSO-PBI algorithm in this paper is useful and efficient.

4.5.   Biological relevance

In order to verify the biological relevance of the selected pathway marker genes, the top 60% pathway markers are chosen and searched in the disease-gene association database (http://www.disgenet.org/). From this database we can obtain the number of PubMed citations of the disease-gene association, and a part is shown in Table 10. The first column of the Table 10 is the disease name and the second column is the corresponding gene symbol, in addition, the numbers in parentheses indicate the number of PubMed citations. From the results, we can know that the marker genes are related to particular diseases, that is to say, the selected marker genes are biologically related.

To explore the potential role of marker genes on overall survival rate, the Kaplan-Meier (KM) survival curve was presented based on the partial marker genes in this study. In Figure 6(a)–(f), the horizontal axis represents the patient's survival time and the vertical axis represents the patient's survival rate. The red and blue lines represent the high and low expression groups, respectively. Although there was no significant difference in survival rate between the two groups during the initial period, with the increase of time, the survival rate of the high expression group decreased rapidly. The low expression group had higher survival rate and longer survival time, which demonstrates that the cancer is related to some specific genes. Through identifying the biomarkers of cancers can help prediction and diagnosis cancers, which can result in a good prognosis.

5.   Conclusions

Biomarkers play an important role in the diagnosis and treatment of diseases, in this paper, a decomposition-based multi-objective optimization algorithm is proposed to screen pathway marker genes. Consider in the process of inferring pathway activity, the contribution of each gene should be different. In this research, the relevance of each gene in pathway activity inference is quantified. Moreover, in order to avoid the excessive concentration of feature solutions in most multi-objective optimization algorithms, an adaptive mechanism based on PBI decomposition is introduced. The proposed IMOPSO-PBI method is applied to six real datasets and compared with some existing algorithms, the results show that the method achieves better performance. In addition, the biological relevance of the screened marker genes was proved by biological analysis. In the future, we can optimize on this framework, such as choosing different objective functions and performance metrics, and we can also use this method on other datasets, such as HNSCC, KIRC, besides, some public available single-cell RNA-seq along with the bulk RNA-seq are also good choice.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 61803257) and Natural Science Foundation of Shanghai (No. 18ZR1417200).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.