Network-based multi-class classifier to identify optimized gene networks for acute leukemia cell line classification

TL;DR

This paper introduces a new computational method to classify acute leukemia cell lines using optimized gene networks, improving both accuracy and biological relevance.

Contribution

The novel GRN-multiClassifier method jointly optimizes gene network estimation and multi-class classification for cell line categorization.

Findings

GRN-multiClassifier outperforms existing methods in classifying acute leukemia cell lines into distinct categories.

The identified gene interactions align with known acute leukemia markers and biological pathways.

Key gene interactions like HBA1&HBB and ACTB inhibition are highlighted as potentially significant for disease mechanisms.

Abstract

Unraveling the genetic regulatory networks that underlie diseases is essential for comprehending the intricate mechanisms of these conditions. While various computational strategies were developed, the approaches in the existing studies concerning network-based prediction and classification are based on the pre-estimated gene networks. However, the gene network that is pre-estimated fails to yield biologically meaningful explanations for classifying cell lines into particular clinical states. The reason for this limitation is the lack of inclusion of any information about the clinical status of cell lines during the process of network estimation. To achieve effective cell line classification and ensure the biological validity of the cell lines classification, we develop a computational strategy referred to as GRN-multiClassifier for network-based multi-class classification. The GRN-multiClassifier estimates gene network in a manner that simultaneously minimizes both the network estimation error and the negative log-likelihood function of multinomial logistic regression. That is, our strategy estimates optimized gene network to enable the multi-class classification of cell lines into specific clinical conditions. Monte Carlo simulations demonstrate the efficacy of the GRN-multiClassifier. We applied our strategy to network-based classification of acute leukemia cell lines into three distinct categories of acute leukemia. Our strategy shows outstanding performance in the classification of acute leukemia cell lines. The results for the acute leukemia marker identification are strongly supported by existing literature. The implications of our findings suggest that potential pathways involving the inhibition of ACTB and the molecular interactions between “HBA1&HBB,” “HBB&HBA1,” “IGKV1-5&IGHV4-31,” “IGHV4-31&IGKV1-5,” “HLA-DRA&CD74” and “ACTB&ACTB” could offer significant insights into the underlying mechanism of acute leukemia.

Figures

Overall framework of GRN-multiClassifier.

Leukemia cell line classification results,where GNmC: GRN-multiClassifier, prNW: multi-class classification model grounded in a pre-estimated network, Classification results based on expression levels by using LA: lasso, ELA: elastic net, KSVM: kernelized support vector machine, RF: random forest. The AML cell lines classification is performed by 50, 100, 200, 300, 400, 500 genes having largest variances of expression levels.

Heatmap of the identified leukemia markers in the B-cell acute lymphoblastic leukemia, T-cell acute lymphoblastic leukemia, and acute myeloid leukemia cell lines.

Molecular interplays of the identified acute leukemia markers in B-cell acute lymphoblastic leukemia (BALL), T-cell acute lymphoblastic leukemia (TALL), acute myeloid leukemia cell lines (AML) and the estimated network without consideration of classification (Ordinary network).

GO term pathway analysis of the genes in the identified crucial common edges for three type of acute leukemia cell lines classification: enriched GO term with p.value<0.01 and -log(p.value).

ALL-MLL-AML classification results.

Lung cancer subtype classification results.

Equations

Full text

Introduction

Gene expression level analysis alone is insufficient to comprehend the intricate mechanisms driving certain disease states, as diseases involved in complex mechanism arise from disturbances in the specific operations of molecular networks rather than anomalies in individual genes [1]. In the field of biomedical research, one of the most impactful technique is heterogeneous gene regulatory network analysis and the study of gene networks has gained significant research attention in various fields of research. The effectiveness of network-based analysis has been validated in various research works, e.g., cancer prediction, drug combinations identification, and protein-protein interaction [2–4]. In particular, various investigations have delved into network-based prediction and classification, and their effectiveness has been verified [5, 6]. Although various computational methods have been developed and employed to predict and/or classify the cancer-related status of cell lines, previous studies performed the gene network analysis based on the pre-estimated gene networks. In other words, the existing studies estimated gene networks in advance based on gaussian graphical modeling, Bayesian network, etc. [7, 8], and then used the pre-estimated gene network to classification, prediction, uncovering disease-related mechanisms. However, the analysis conducted using pre-estimated gene networks lacks the capacity to provide biologically reliable results, as it does not incorporate crucial biological mechanisms, such as those related to cancer, into the network estimation process.

In order to achieve biologically reliable results of classification and gene network analysis, we develop computational approach termed Gene Regulatory Network-based multi-class Classifier (GRN-multiClassifier), designed for gene network-based multi-class classification, following the methodology proposed by Park et al. [9]. The objective function of the GRN-multiClassifier is based on the negative log-likelihood function for the multinomial logistic regression model and error term of gene network estimation. A crucial aspect to highlight is that our approach conducts gene network estimation and multi-class classification simultaneously. In our strategy, the input of multi-class classification model is the gene network that is estimated to minimize both the error in network estimation and the negative log-likelihood function, simultaneously. That is, the gene network estimation is concurrently performed with cell line classification by minimizing not only network estimation error but also classification error. It implies that the estimated gene networks are iteratively updated to optimize multi-class classification. The optimized gene networks may have crucial information to characterize each class and it can be useful tool to identify markers of specific phenotypes. Thus, we can efficiently perform for biologically reliable interpretation of the results arising from gene network estimation and the classification of cell lines. On the other hand, pre-estimated networks are estimated by minimizing only the network estimation error without consideration of classification. Thus, the pre-estimated networks cannot describe the optimized molecular interplays to classify the cell lines.

Through Monte Carlo simulations, we show the effectiveness of the GRN-multiClassifier in both multi-class classification and gene network estimation. We applied the GRN-multiClassifier to the acute leukemia gene expression dataset comprising 72 cell lines and 7129 genes [10], facilitating gene network-based classification of acute leukemia cell lines. Our strategy shows effective performance in the classification of cell lines associated with acute leukemia. We then identify the markers relevant to acute leukemia as well as their associated molecular interplays and the validity of the identified markers is confirmed through existing literature. While biological knowledge was not incorporated into the identfied markers and their networks, our data-driven strategy yields biologically reliable outcomes for the identification of acute leukemia markers. The results of our study suggest that the increased expression of ACTB and the molecular interplays involving “HBA1&HBB," “HBB&HBA1," “IGKV1-5&IGHV4-31," “IGHV4-31&IGKV1-5," “HLA-DRA&CD74" and “ACTB&ACTB" might contribute to the progression of acute leukemia. In light of this, we suggest that the suppression of ACTB and the associated molecular interactions could provide vital insights into understanding potential strategies for both preventing and treating acute leukemia.

The remainder of this paper is organized as follows: In the Methods section, we outline the computational methods for gene network estimation and as well as the novel strategy we have developed for network-based multi-class classification. Subsequently, we introduce the numerical resolution process for the GRN-multiClassifier. In the section of Monte Carlo simulation, we demonstrate the results derived from the executed simulation studies. The findings of gene network-based acute leukemia cell line classification are detailed in the section titled Gene Regulatory Network-based Acute Leukemia Subtypes Classification. Concluding remarks and insights are presented in the Discussion section.

Methods

Gene regulatory network estimation

Suppose $[eqn]$ is an $[eqn]$ data matrix describing the expression levels of p regulator genes that control the $[eqn]$ target gene transcription $[eqn]$ , $[eqn]$ .

The gene regulatory network can be represented by the following linear regression model,

$$[eqn]$$

where $[eqn]$ is the regression coefficient that represents the effect of p regulator genes $[eqn]$ on $[eqn]$ target gene yij and $[eqn]$ is a random error vector for the $[eqn]$ target gene. The regularization methods have been often used to estimate the gene networks,

$$[eqn]$$

where $[eqn]$ is the regularization penalty, e.g., ridge [11], lasso [12], elastic net [13], etc.,

$$[eqn]$$

and $[eqn]$ are the regularization parameters of $[eqn]$ .

The gene networks are crucial to understand complex cancer-related mechanisms, because the mechanisms are involved in dysregulation and dysfunction of molecular networks. Although numerous computational methodologies have been formulated and applied to gene network analysis, the methods cannot provide biologically reliable interpretation of the results. This limitation arises from the fact that the process of estimating gene networks did not consider cancer-associated mechanisms, i.e., the networks are pre-determined, and these pre-estimated gene networks are subsequently applied for tasks like classification, prediction, and the exploration of disease-related mechanisms. In order to identify molecular interplays that are optimized to explain cancer-related status of cell lines (e.g., clones), we developed a novel computational approach named Gene Regulatory Network based multi-class Classifier (GRN-multiClassifier).

Preliminaries

Let Zi is categorical response variable having G>2 levels to describe phenotypes of $[eqn]$ cell line and $[eqn]$ is expression of p gene in $[eqn]$ cell line. Suppose we have n independent observations $[eqn]$ , where Zi is categorical response variable having G>2 levels. The multinomial logistic regression model for multi-class classification assumes that

$$[eqn]$$

where $[eqn]$ . We suppose the $[eqn]$ element in $[eqn]$ row of $[eqn]$ equals one when the $[eqn]$ cell line belongs to phenotype g, and zero otherwise, i.e., $[eqn]$ represents the multinormial trial for subject i with zig = 1 when the response is in phenotype g and zig = 0 otherwise. In other words, $[eqn]$ serves to phenotypes of cell lines and is a random variable with multinormial distribution, such that $[eqn]$ , having the following probability mass function,

$$[eqn]$$

where $[eqn]$ . The log-likelihood has the follows form

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ . Although maximum likelihood method has been used to estimate the multinormial logistic regression model, the method frequently yields unstable estimation results with significant variation. This instability is particularly prominent when multicollinearity exists among predictor variables or when dealing with datasets of high dimensionality [14]. To address the issue, we consider the following penalized log likelihood methods,

$$[eqn]$$

Most of the previous studies in the domain of multi-class classification for cell lines were based on the gene expression levels as the input for the multinomial logistic regression model, i.e., the input of the multinormial logistic regression model was expression levels of genes $[eqn]$ . However, a single gene-based analysis is insufficient in generating biologically reliable results and meaningful interpretations for the classification of cell lines. This is due to the fact that the complex mechanisms of diseases are intricately linked to molecular networks, rather than being solely reliant on the perturbation of individual genes.

Gene regulatory network based multi-class classifier: GRN-multiClssifier

We consider the network-based multi-class classification and develop a computational strategy that performs network estimation and classification, simultaneously. For the categorical response variable Zi describing phenotypes of cell lines, the network-based multi-class classification is performed by the expression levels of genes $[eqn]$ and their regulatory networks $[eqn]$ as follows,

$$[eqn]$$ $$[eqn]$$

In the network-based multi-class classification, the log-likelihood function takes on the following form

$$[eqn]$$

The network-based multi-class classification in (9) is based on the pre-estimated gene network $[eqn]$ . Thus, we cannot effectively interpret the classification results based on the network, because any information of the status of cell line was incorporated in to the network estimation process.

To derive a gene regulatory network that is optimized for multi-class classification, we consider the following strategy that involves the simultaneous estimation of gene networks and cell line classification, following the approach presented by Park et al. [9],

$$[eqn]$$

Our approach aims to estimate the gene regulatory network $[eqn]$ in a manner that minimizes not only the errors associated with gene network estimation but also the negative log-likelihood pertaining to multi-class classification. That is, the resultant estimated network $[eqn]$ represents the optimized gene regulatory network tailored for cell line classification. Consequently, we are able to carry out biologically meaningful interpretations for both the outcomes of gene network estimation and the multi-class classification of cell lines.

In order to address the limitations of the maximum likelihood method, we consider penalized multinormial logistic regression based on the regularization approaches. Moreover, we integrate the following insights from network biology into the statistical model to attain outcomes that are more biologically reliable results.

Genes linked in the networks may have analogous biological functions.The hub genes, involved in interactions with a multitude of other genes, play key roles in governing the function and expression of multiple target genes. Disruption or dysfunction in these genes can undermine the intricate balance of gene networks, resulting in profound subsequent influences on cellular functions and disease phenotypes [15]

The knowledge of network biology can be incorporated by using the network-constrained regularization [16]. The estimated gene network from the second term in (10) can be represented by a weighted graph $[eqn]$ , where $[eqn]$ is the set of vertices corresponding to p genes and $[eqn]$ is the set of edges (i.e., pair (i,j), where $[eqn]$ and $[eqn]$ . $[eqn]$ is the edge weight. The normalized Laplacian matrix $[eqn]$ for the graph is given as [16].,

$$[eqn]$$

where di is the degree of each gene, which is given as $[eqn]$ .

In our method, the estimated effect of regulators on their target genes is presented within the matrix $[eqn]$ , where each row and column of $[eqn]$ corresponds to the index of a regulator and a target gene, respectively. We compute the weight of edges $[eqn]$ based on the effect of the $[eqn]$ gene to $[eqn]$ gene (i.e., $[eqn]$ ) and the $[eqn]$ gene to $[eqn]$ gene (i.e., $[eqn]$ ) as follows,

$$[eqn]$$

We calculate the Laplacian matrix using the edge weights $[eqn]$ , and subsequently, we integrate the estimated network into the penalized multinomial logistic regression model utilizing the Laplacian matrix $[eqn]$ .

We then proposed the following Gene Regulatory Network based multi-class Classifier (GRN-multiClassifier),

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ with $[eqn]$ [17]. The proposed GRN-multiClassifier can estimate gene regulatory network that is optimized to multi-class classification. As a result, we are enabled to carry out biologically credible interpretation of the classification outcomes, grounded in the gene network. Furthermore, our method encourages similarity in gene coefficients for genes that share common edges. This is facilitated by giving considerable weight to the distinction between coefficients of genes that share a substantial number of edges. Consequently, our approach has the ability to simultaneously identify connected genes within the network. Furthermore, our approach enforces a relatively slight penalty on hub genes making it possible for these hub genes allowing hub genes to be readily recognized as key features for the classification of cancer-related statuses. In summary, our model integrates network biology knowledge, enabling robust gene network analysis and accurate multi-class classification of cell lines with a strong biological foundation. We expect that the proposed GRN-multiClassifier will be a useful tool for identifying crucial molecular interplays to characterize diseases related phenotypes. Figure 1 shows overall framework of the proposed GRN-multiClassifier.

Overall framework of GRN-multiClassifier.

Implementation

The optimization of the unknown parameters $[eqn]$ is nonlinear and the equation does not have explicit solution. The solution, $[eqn]$ in this case can be obtained by using iterative algorithm. We consider Fisher’s score method to estimate the proposed GRN-multiClassifier.

By refer to the Fisher’s scoring algorithm in Appendix section, the quadratic approximation of the objective function of the GRN-multiClassifier in (13) is given as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ .

We propose the coordinate descent algorithm to optimize $[eqn]$ and $[eqn]$ . The coordinate update to optimize $[eqn]$ given $[eqn]$ has the following form,

$$[eqn]$$

where

$$[eqn]$$

and $[eqn]$ is a soft thresholding operator with value

$$[eqn]$$

The sign adjusted Laplacian matrix is computed as $[eqn]$ . Then, the parameters $[eqn]$ are estimated by the following coordinate-wise update

$$[eqn]$$

where $[eqn]$ is $[eqn]$ element in c row of $[eqn]$ and $[eqn]$ .

Given $[eqn]$ and $[eqn]$ , coordinate $[eqn]$ is given as

$$[eqn]$$

The parameters $[eqn]$ and $[eqn]$ are cyclically updated until convergence.

Regularization parameters selection

The results of the proposed GRN-multiClassifier heavily rely on the regularization parameters $[eqn]$ and $[eqn]$ . The $[eqn]$ and $[eqn]$ are regularization parameters of $[eqn]$ and $[eqn]$ , respectively. When $[eqn]$ and $[eqn]$ ( $[eqn]$ , $[eqn]$ ), the penalty of $[eqn]$ ( $[eqn]$ ) is reduced to the lasso. On the other hand, for $[eqn]$ and $[eqn]$ ( $[eqn]$ and $[eqn]$ ) retains the good properties of both sparsity and smoothness for edge estimation in a network, i.e., $[eqn]$ , (predictor variable estimation in classification model, i.e., $[eqn]$ ).

We consider the following Bayesian Information Criterion (BIC) to select the regularization parameters,

$$[eqn]$$

where $[eqn]$ the set of indexes of the validation dataset, $[eqn]$ is a number of observations of the validation dataset, $[eqn]$ are estimated parameters based on training dataset, and df is the degree of freedom of the GRN-multiClassifier. We use the number of nonzero elements in $[eqn]$ as an estimate of the degrees of freedom in line with [9, 18, 19]. We apply our method to acute leukemia subtypes classification based on dataset from the R package golubEsets (https://jokergoo.github.io/cola_examples/Golub_leukemia/).

Monte Carlo simulations

We demonstrated the effectiveness of the GRN-multiClassifier through the utilization of Monte Carlo simulations. We assumed that each transcription factor gene (TF) regulates a set of 10 genes, and the expression levels of these transcription factor genes were generated from N(0,1). The expression levels of each of the regulated genes ( $[eqn]$ , $[eqn]$ ) by the TF ( $[eqn]$ ) were generated in accordance with the expression level of the TFs as follows,

$$[eqn]$$

where $[eqn]$ .

The response variable $[eqn]$ is generated by taking into account the regulatory effect of genes, that is, both the gene expression levels denoted as $[eqn]$ and the effect of regulators on their targets represented by the matrix $[eqn]$ .

We conducted simulations with 50 datasets generated from the following true model,

$$[eqn]$$

We set G = 3 and assume that the data were simulated based on the scenarios established as benchmarks in prior studies [9, 16].

Scenario 1:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Scenario 2:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Scenario 3:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Scenario 4:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

We consider the number of observations n = 150 consisting of training, validation and test dataset with 80%, 10%, and 10% of 150 observations, respectively. In each scenario, we consider number of TFs (i.e, T) as 10 and 20. We assess the performance of the proposed Gene Regulatory Network based multi-class Classifier (GNmC) by comparing it with a multi-class classification model grounded in a pre-estimated network (prNW). The gene network for the prNW model is obtained through the utilization of the lasso. We further conduct a comparison among classification methods based on not the network, but the expression levels of genes $[eqn]$ , where lasso (LA), elastic net (ELA), kernelized support vector machines (KSVM) and random forest (RF) are used for classification of cell lines. That is, the GNmC and prNW are the network-based classification, while LA, ELA, KSVM and RF are the expression-based classification.

We also consider the scenarios for networks with the regulators comprising both activators and inhibitors, as described below:

Scenarios 1 and 2:

$$[eqn]$$ $$[eqn]$$

Scenario 3 and 4:

$$[eqn]$$ $$[eqn]$$

In the situation for networks consisting of activators and inhibitors, the $[eqn]$ is given in same as the scenarios for networks consisting of activators only.

We compare the multi-class classification accuracy in Table 1. As indicated in Table 1, the network-based methods (i.e., GRN-multiClassifier and prNW) exhibit exceptional performance in terms of multi-class classification accuracy compared with expression levels based approaches (i.e., ELA, LA, KSVM and RF). Furthermore, it can be seen that the proposed GRN-multiClassifier provides the most effective multi-class classification results in overall. We also evaluate our method for the feature selection accuracy, encompassing true positive rates, true negative rates, and their averages, for both $[eqn]$ in the multinomial logistic regression model and $[eqn]$ in the network estimation. [Simul_FeatureS]Table 2 shows the feature selection results. From the perspective of feature selection accuracy for $[eqn]$ , the network-based methods demonstrate effective outcomes in comparison to approaches based on expression levels. Our approach demonstrates remarkable outcomes in terms of feature selection, encompassing not only multinomial logistic regression but also the selection of edges in gene regulatory network estimation. The results clearly indicate that the proposed GRN-multiClassifier delivers exceptional performance across gene network estimation, feature selection, and multi-class classification of cell lines.

Table 1: Multi-class classification accuracy, where GNmC: GRN-multiClassifier, prNW: multi-class classification model grounded in a pre-estimated network, Classification results based on expression levels by using LA: lasso, ELA: elastic net, KSVM: kernelized support vector machine, RF: random forest.

Table 2: Simulation results of networks selection in multinormial logistic regression and edge selection in network estimation, where GNmC: GRN-multiClassifier, prNW: multi-class classification model grounded in a pre-estimated network, Classification results based on expression levels by using LA: lasso, ELA: elastic net, KSVM: kernelized support vector machine, RF: random forest.

Gene regulatory network based acute leukemia subtypes classification

Acute leukemias are malignancies that originate from either the lymphoid or myeloid cell line and are characterized by rapid and uncontrolled proliferation of immature blood cells. AML is a cancerous condition that affects myeloid cells, responsible for generating certain white blood cells, whereas Acute Lymphocytic Leukemia (ALL) primarily impacts lymphocytes, a critical component of the immune system’s white blood cells. AML is a highly aggressive cancer that tends to progress rapidly and can be particularly deadly and thus uncovering the molecular interplays that play key role in AML-related mechanism is a crucial issue. We apply the proposed GRN-multiClassifier to estimate a gene network that is optimized to provide insights into the characteristics of acute leukemia cell lines. We consider the well-known acute leukemia gene expression dataset [10], which has been explored in various fields of research on multiclass cancer analysis [21, 22]. The data set consist of 72 cell lines for 7129 genes that are available in the R package golubEsets (https://jokergoo.github.io/cola_examples/Golub_leukemia/).

Diagnosing acute leukemia can be framed as a tri-classification problem, i.e., B-cell acute lymphoblastic leukemia (BALL), T-cell acute lymphoblastic leukemia (TALL), and AML [23, 24], where the dataset consisting of 38 BALL cell lines, 9 TALL cell lines, and 25 AML cell lines. In this study, we consider the classification of cell lines into the tree types of acute leukemia (g = 1: AML; g = 2:BALL; g = 3:TALL) based on the gene regulatory networks.

It is known that extremely high dimensional data situation can cause poor performance for machine learning model [20] and lead to difficulty in visualization of the results. Thus, we peroform multi-class classification based on 50, 100, 200, 300, 400, and 500 genes with the highest variance in 72 cell lines. The training and test dataset consist of 57 (80%) and 15 (20%) cell lines, respectively.

Figure 2 shows the classification results. The proposed GRN-multiClassifier shows effective results in the classification of acute leukemia cell lines in overall, even though the random forest show better classification accuracy than our method.

Leukemia cell line classification results,where GNmC: GRN-multiClassifier, prNW: multi-class classification model grounded in a pre-estimated network, Classification results based on expression levels by using LA: lasso, ELA: elastic net, KSVM: kernelized support vector machine, RF: random forest. The AML cell lines classification is performed by 50, 100, 200, 300, 400, 500 genes having largest variances of expression levels.

Our strategy shows the effective classification results (i.e., 100% accuracy) based on the classifiers with 100, 300, 500 genes. Other methods (i.e., prNW, KSVM and RF) also show the effective results (i.e., 100% accuracy) based on 300, 500, 400 genes, respectively. It implies that increasing number of features cannot always improve classification accuracies. The classification utilizing the gene network comprising 50 genes achieves a perfect accuracy of 100%, where the large number of genes do not guarantee improvement of accuracy. Thus, we consider that the gene network with 50 genes is enough to explain status of cell lines, because the network achieves the perfect classification result. We proceed to interpret the outcomes of acute leukemia cell line classification based on this gene network of 50 genes. The optimized gene networks (i.e., edges) with 50, 100, 300, 500 genes estimated by our method are given in Supplementary file.

We identify genes that correspond to the largest absolute values of $[eqn]$ , with $[eqn]$ , as crucial markers, i.e., we extract the top 5 genes for each category ( $[eqn]$ ) based on their largest absolute $[eqn]$ values. Table 3 shows the crucial markers and their leukemia related evidences. As shown in Table 3, 11 identified markers have been uncovered in previous studies and only a marker has no evidences.

Table 3: The identified markers for Leukemia cell line classification and their evidences related AML markers uncovered in previous studies.

ELANE In the study conducted by Yanli et al. [25], it was observed that heightened expression of ELANE correlated with comparatively shorter survival durations among leukemia patients. In light of this, the authors proposed that ELANE serves as an oncogene driving leukemia development. Patients harboring the ELANE mutation may experience the onset of myelodysplastic syndrome (MDS), AML, or, in rarer cases, acute lymphoblastic leukemia (ALL) [26]. Although it is vital to keep track of patients with ELANE-germline mutations, the initial leukemogenesis process in ELANE-neutropenia patients is characterized by the presence of CSF3R mutations, potentially induced by GCSF therapy and is adequate to cause myelodysplasia and acute leukemia [26]. CXCL8 Research findings indicate that the expression level of CXCL8 is positively correlated with recurrence probability in AML and CXCL8 plays significant role in promoting AML cell growth through the activation of the ERK1/2 signal pathway [28]. The CXCL8 derived from mesenchymal stromal cells promotes the survival and proliferation of AML cells by activating the PI3K/AKT pathway [29]. Significant upregulation of CCL3, CCL4, CXCL8, and IL-17A messenger RNA (mRNA) expression levels was observed in the adult T-cell leukemia/lymphoma (ATLL) groups, as reported by Soltani et al. [30]. The significant association between CXCL8 and ATLL was evident, and the upregulation of CXCL8 was associated with an increased odds ratio of ATLL. CST3 CST3 has been identified as a biologically relevant gene in the context of leukemia [31], and it has also been noted as one of the up-regulated genes in patient/donor pairs with AML [32]. Elevated expressions of CSTA, CSTB, CST3, and CST7 were found to be correlated with higher percentages of monocytes and neutrophils in the peripheral blood of patients with AML, as reported by Yuita et al. [34]. CST3 was selected in the Top 10 ranked informative genes from the leukaemia dataset [34]. CST3 and MPO were experimentally demonstrated to exhibit a correlation with acute lymphoblastic leukemia (ALL) or AML [35]. In the study by Wang et al. [35], it was shown that the gene groups CST3, MPO, and IGL are strongly associated with the development of leukemia. MPO Itonaga et al. [36] revealed that MPO is linked to myeloid lineage commitment and is indicative of a favorable prognosis in patients with AML. The presence of MPO serves as an informative marker for distinguishing a distinctive and crucial DNA methylation profile in CD34-positive AML cells. Patients with B-ALL-isoMPO exhibit a greater risk of relapse than those with B-ALL [38]. MPO plays a crucial role in determining the susceptibility of leukemia cells to parthenolide-induced apoptosis [39]. Parthenolide holds promise as a potential therapy targeting leukemic stem cells, particularly for AML cases characterized by elevated levels of MPO expression. Expression of MPO was detected in a majority of AML [40]. MPO gene expression has the potential to function as an additional marker for distinguishing various types of acute leukemias. Moreover, it could aid in the identification of leukemic cells arrested during the initial stages of the myeloid differentiation pathway. The presence of elevated MPO expression, indicative of a more differentiated cellular profile, has been associated with positive clinical outcomes in AML, as outlined by Kumari et al. [41]. HLA-DRA In specific human leukemia cases, the expression of HLA-DR antigens is indicative of cellular differentiation, as discussed in the work by Tobe et al. [42]. Dorak et al. [43] showed a highly significant association of a homozygous HLA-DR genotype in childhood ALL with a strong gender effect. Multiple sclerosis risk markers in HLA-DRA, HLA-C, and IFNG genes are associated with sex-specific childhood leukemia risk [44]. CD74 The expression of CD74 was detected in various AML cell lines and patient samples that exhibit sensitivity to cytotoxicity induced by the anti-CD74 treatment (milatuzumab), as demonstrated by Le et al. [45]. It was also demonstrated that CD74 is expressed more frequently and at higher levels on AML blasts compared to lymphocytes. The ease of targeting the phenotype and the presence of the anti-CD74 antibody milatuzumab suggest that further investigate into the role of CD74 in AML biology [40]. TCL1A Aggarwal et al. [46] demonstrated the prognostic relevance of TCL1A expression in patients with chronic lymphocytic leukemia and mantle cell lymphoma. Expression and/or rearrangement of TCL1A is a useful marker to diagnosing T-cell prolymphocytic leukemia [47]. Upregulation of PAX5, CD72, CSRP2, LOC100130458, TCL1A and EBF1 genes is associated with patients diagnosed with ALL patients in a phenotype-related signature. DEFA1 Overexpression of FLT3 and DEFA1 genes retained independent prognostic significance for B-ALL outcome [48]. The overexpressions of DEFA1-3 may be associated with an increase in malignancy during relapsed leukemia and could potentially serve as predictive markers for therapy resistance during relapse, as suggested by Te et al. [49]. Te Kronnie et al. [49] illustrated that relatively high expression levels of DEFA1-3 are linked to an unfavorable response to chemotherapy, which in turn leads to early relapse in leukemia patients. The overexpression of DEFA1-3 was found to be associated with elevated expression of both MPO and ELA2, where these latter genes are preferentially transcribed within the most immature granulocytes. There was a suggestion that the simultaneous overexpression of DEFA1-3, MPO, and ELA2 might contribute to a specific granulopoiesis signature. TRBV19/TCRB The presence of TCRB gene rearrangements was detected in 35% of precursor-B-ALL patients, encompassing both children and adult [50]. ACTB ACTB was identified as a member of set for use as control/reference for the analysis of gene expression in peripheral blood and bone marrow samples from patients with acute leukemias [51]. The genes CLUS, CERU, APOE, APOA4, APOA1, GELS, S10A9, AMBP, ACTB, CATA, and AFAM have been identified as important factors in leukemia prognosis, with the potential to act as unique biomarkers for gauging the aggressiveness of leukemia or as suppressor proteins specifically in cases of high-risk acute lymphoblastic leukemia (HR-ALL), as suggested by the study referenced as [52]. HLA-B Fernandez et al.[53] uncovered the contradictory effects of the HLA-B*40 allele in terms of genetic susceptibility to develop ALL or AML. The somatic mutation identified in the HLA-B gene of leukemic cells was responsible for the observed variations in typing and sequencing alterations in the peripheral blood sample, as stated in the study cited as [54]. Vikash et al. [55] provided the evidence of the HLA-B allele is associated with leukemia in the North Indian population.

It can be seen through literature that the crucial genes identified for the classification of acute leukemia cell lines have strong evidences as acute leukemia markers.

Figure 3 shows the heatmap of the identified markers in the BALL, TALL, and AML cell lines. The markers show different pattern in the BALL, TALL, and AML cell lines. X00274 $[eqn]$ at, M13560 $[eqn]$ s $[eqn]$ at, L19779 $[eqn]$ at show relatively high expression levels in AML and BALL cell lines than TALL cell lines. Y00787 $[eqn]$ s $[eqn]$ at, M27891 $[eqn]$ at, M19507 $[eqn]$ at, M27783 $[eqn]$ s $[eqn]$ at, M20203 $[eqn]$ s $[eqn]$ at are up-regulated in AML cell lines, while X00437 $[eqn]$ s $[eqn]$ at is up-regulated in TALL cell lines. D49824 $[eqn]$ s $[eqn]$ at, AFFX.HSAC07.X0035 $[eqn]$ M $[eqn]$ at, X00351 $[eqn]$ f $[eqn]$ at and AFFX.HSAC07.X00351 $[eqn]$ 5 $[eqn]$ at show high expression levels in all acute leukemia cell line. This suggests that the ACTB gene (Accession ID: AFFX.HSAC07.X00351 $[eqn]$ M $[eqn]$ at, X00351 $[eqn]$ f $[eqn]$ at, and AFF X.HSAC07.X00351 $[eqn]$ 5 $[eqn]$ at) holds the potential to be regarded as a critical marker for acute leukemia. Furthermore, the outcomes imply that inhibiting ACTB could potentially offer valuable insights into comprehending the mechanisms underlying acute leukemia.

Heatmap of the identified leukemia markers in the B-cell acute lymphoblastic leukemia, T-cell acute lymphoblastic leukemia, and acute myeloid leukemia cell lines.

The proposed GRN-multiClassifier incorporates the L1 penalty of $[eqn]$ into multinomial logistic regression, and thus provides sparse estimation result for $[eqn]$ . It implies that we can identify crucial markers for the classification of acute leukemia cell lines. We present the networks of the identified crucial markers, i.e., the networks consist of the selected target genes ( $[eqn]$ ) and their regulator genes. We estimate gene networks for AML, BALL, and TALL cell lines separately using the genes within the extracted network. These networks are constructed using the specific cell lines corresponding to AML, BALL, and TALL. To effectively visualize, we extract edges with an absolute value exceeding 0.5. Figure 4 shows the crucial molecular interactions for acute leukemia cell line classification in AML, BALL, TALL cell lines and the estimated gene network without consideration of classification (ordinary network), where the ordinary network describes not optimized molecular interplays for leukemia cell line classification but the original network without consideration of the classification. As shown in Figure 4, the networks within AML cell lines display sparse molecular interactions, a dense gene network is observed in TALL cell lines, where numbers of edges are 18, 26, and 28 in the networks of AML, BALL and TALL cell lines, respectively. The ordinary network shows the relatively sparse network, i.e., the network consists of 18 edges. The AML and BALL cell lines exhibit analogous molecular interactions, specifically sharing 11 common edges in their respective networks. Conversely, the TALL cell lines feature distinct gene networks when compared to the networks of AML and BALL cell lines. In detail, the networks of AML and TALL cell lines share 6 common edges, while the networks of BALL and TALL cell lines share 7 common edges. The ordinary network relatively larger number of common edges with the networks of AML and BALL cell lines, while the network of TALL cell lines shows differentially regulated gene network with those of the ordinary network. The results imply that the ordinary network describes the general molecular interplays for the network of AML, BALL and TALL cell lines, where 12 edges among the 18 edges of the ordinary network are existed in the networks of AML and BALL cell lines. On the other hand, the optimized gene networks for leukemia cell line classification show relatively distinguishing molecular interplays for each cell lines. Especially, the T-cell acute lymphoblastic leukemia may possess distinct molecular characteristics that set it apart from AML and BALL.

Molecular interplays of the identified acute leukemia markers in B-cell acute lymphoblastic leukemia (BALL), T-cell acute lymphoblastic leukemia (TALL), acute myeloid leukemia cell lines (AML) and the estimated network without consideration of classification (Ordinary network).

For the genes in the optimized networks for the leukemia cell line classification (i.e., genes in the networks of AML, BALL and TALL cell lines in Figure 4), we perform differentially expression genes analysis. Table 4 show the results of analysis of variance (ANOVA). Among the 33 genes in the optimized networks, 17 genes (51.2%) show significant difference (Pvalue<0.05) of expression levels in AML, BALL and TALL cell lines, while only 23.2% of the non-selected genes show significantly different expression levels. The results imply that our strategy for the optimized gene networks analysis for cell line classification can identify crucial genes to characterize status of cell lines in the viewpoint of expression levels. It can be considered that the molecular interplays of the remain 16 genes corresponding Pvalue larger than 0.05 are crucial features for classification, even though the expression levels of the genes are not differently expressed between the AML, BALL and TALL cell lines. Thus, the genes can be considered as novel founding (i.e., candidate markers) for leukemia cell line classification,which cannot be revealed by widely used single gene-based analysis (e.g., differentially expressed gene analysis).

Table 4: Differentially expressed genes analysis

We extract the common edges from the three gene networks for AML, BALL and TALL cell lines in Table 5. To comprehend the biological processes implicated in the shared gene regulatory system of acute leukemia cell lines, we conducted a Gene Ontology (GO) term pathway analysis of the genes present in the common edges. Figure 5 illustrates the enriched pathways, with the p-value represented as -log(p-value). As shown in Figure 5, “blood microparticle" is the most enriched pathway of the genes in common edges for AML, BALL and TALL cell lines. The immune response related pathways (“immune response" and “adaptive immune response") are also enriched for the common markers of acute leukemia cell lines. Furthermore, lumenal side of endoplasmic reticulum membrane-related pathways (“integral component of lumenal side of endoplasmic reticulum membrane" and “lumenal side of endoplasmic reticulum membrane") are also identified as enriched pathways for the common markers.

GO term pathway analysis of the genes in the identified crucial common edges for three type of acute leukemia cell lines classification: enriched GO term with p.value<0.01 and -log(p.value).

Table 5: Identified common edges in gene networks for AML, BALL and TALL cell lines shown in Figure 3. The common edges are existing in three gene networks estimated by BALL, TALL, AML cell lines.

Table 6 show the leukemia related evidences for the identified GO terms of the crucial common edges. The uncovered biological pathways for common genes in the networks were identified as crucial biological functions to AML mechanism by the Gene Ontology (GO) term analysis of differentially expression genes analysis.

Table 6: Evidences for the Leukemia of the identified biological pathways.

Blood microparticle The GO term Blood microparticle was revealed as an enriched pathway related differentially expressed genes (DEGs) related leukemia cells as follow. A comparison of differential mRNAs and miRNAs between the pAML and control groups demonstrated that pediatric AML DEGs are predominantly enriched in the Blood microparticle and immunoglobulin complex pathways [56]. The biological pathway Blood microparticle was identified as relevant to the DEGs in high-white and low-white count B-cells as well [57]. According to Huang et al. [58], cellular component enrichment analysis revealed that the DEGs distinguishing TP53 mutation from wild-type AML patients were significantly associated with the Blood microparticle pathway. The Blood microparticle pathway was also identified as being enriched with DEGs related to chronic lymphocytic leukemia [59]. Nitric-oxide synthase binding According to Brandao et al. [60], AML patients exhibited elevated levels of nitric oxide synthase expression relative to controls. Immune response Fu et al. [61] discovered that upregulated leukemia-promoting genes were significantly enriched in biological processes such as cell activation involved in Immune response, cytokine production, and leukocyte migration. The differentially expressed probe sets that exhibit sex-specific variations in AML patients are significantly enriched in pathways related to extracellular space, immune response, protein binding [62]. Yang et al. [63] showed that the DEGs distinguishing high-risk from low-risk AML patients were predominantly associated with immune response signaling pathways. Deepak et al. [64] revealed that the network of DEGs in AML compared to normal samples was enriched in immune response-related processes. Adaptive immune response The GO term Adaptive immune response was highlighted as an enriched pathway for the differentially expressed long non-coding RNAs in AML [65]. Iumenal side of endoplasmic reticulum membrane Hu et al. [66] identified the GO term Iumenal side of endoplasmic reticulum membrane as one of the top ten pathways associated with aberrant methylation DEGs in AML. Antigen processing and presentation of exogenous peptide antigen via MHC class II It was revealed that the enriched pathway for DEGs in common myeloid progenitor cells was the GO term Antigen processing and presentation of exogenous peptide antigen via MHC class II [67].

It can be seen through the literature that the identified biological pathways for genes linked the extracted the common edges in BALL, TALL and AML networks are key biological functions and may provide crucial clue for understanding AML mechanism.

From our results, we suggest that the suppression of ACTB could hold the key to comprehending the progression of acute leukemia. We also suggest that targeting the molecular interactions among the genes extracted from the shared edges, as displayed in Table 5, could offer pivotal insights for the prevention and treatment of acute leukemia.

We also compare the classifiers of various methods. The classifiers of GNmC, LA, KSVM and RF with 100 genes show 93% classification accuracies. Thus, we compare the classifiers of the methods for 100 genes. Table 7 shows the classified cell lines by the classifiers of methods.

Table 7: Comparison of classifiers with 100 genes.

As shown in Table 7, the methods that provide same accuracy (i.e., 93%.) construct same classifier. Although the methods show same classifiers for cell line classification, we can expect that our method provides more interpretable results for the classification, as our strategy provides information of not only crucial genes but also molecular interplays.

Evaluations for multi-class cancer classification

We also evaluate our method for multi-class cancer classifications.

Mixed lineage leukemia dataset.

We consider another popular acute leukemia data, called MLL dataset, for 12582 gene expression values for 72 peripheral blood or bone marrow samples consisting of 24 ALL, 20 MLL (mixed lineage leukemia) and 28 AML samples [68]. The MLL Leukemia dataset has been widely explored in research on multi-class cancer analysis and mult-classi classification [69, 70]. For each 24 ALL, 20 MLL and 28 AML samples, we generated the training and test dataset consist of 80% and 20% samples, respectively. Similar to the BALL-TALL-AML classification, we perform ALL-MLL-AML classification based on 50, 100, 200, 300, 400, and 500 genes with the highest variance in 72 samples. Figure 6 shows the classification results of ALL, MLL and AML samples. As shown in Figure 6, our method and RF show the outstanding results for the ALL-MLL-AML classification in overall. The proposed method (i.e., GNmC), KSVM and RF show the most effect classification accuracy (i.e., 100%, 93%, 100% and 100%) based on the models with 50 genes. On the other hand, the expression levels based classification by LA and ELA shows the perfect classification accuracy based on the models with 200 genes. Although some methods shows effective accuracies by classifier with 300, 400, 500 genes (i.e., GNmC and RF), the results implies that the larger number of features does not always provide the highest classification accuracy.

ALL-MLL-AML classification results.

Lung cancer subtype classification.

Our strategy is also applied to classification of not the leukemia cell lines but lung cancer subtypes. We used the publicly available CCLE expression dataset consisting of mRNA expression levels of 19,221 genes in 1,406 cell lines from the DepMap database (https://depmap.org/portal/). From 1,406 cell lines, we extract 206 cell lines that indicate “primary disease" as “Lung Cancer" (see Table 8). The training and test dataset are randomly selected from 80% and 20% of each type of cell lines, respectively. The classification accuracies are also evaluated by classifiers with 50, 100, 200, 300, 400, and 500 genes. The results of lung cancer subtype classification are given in Figure 7. In lung cancer subtype classification, almost methods (i.e., GNmC, LA, KSVM and RF) shows effective classification performances based on not the largest genes, while classifiers with the largest number of genes (i.e., 500) based on prNW and ELA shows the effective results. Our strategy also shows the most effective results based on the classifier with 50 genes. The results also implies that increasing number of features cannot always improve classification accuracies. We can expect through the results of various cancer classification that the proposed strategies will be a useful tool for multi-class classification of disease subtype.

Lung cancer subtype classification results.

Table 8: Subtypes of Lung cancer cell lines of CCLE expression levels dataset

Discussion

The mechanisms involved in disease are related to perturbations in complex molecular networks, rather than in a single gene, thus gene networks are crucial to understand complex mechanisms of disease [1]. The single gene-based analysis cannot provide comprehensive understanding of the disease mechanism. It implies that the expression-based classification cannot describe crucial molecular interplays to understand status of cell lines and/or cancer subtype. The crucialness of the network-based approaches has been demonstrated by many previous studies. Kim et al. [71] proposed a deep learning strategy for drug response prediction based on protein-protein interaction and demonstrated that the network-based prediction can reveal the subnetworks of genes that contribute to the drug response. Rapaport et al. [72] also proposed a classification methodology by incorporating the knowledge of the gene network a priori. Their strategy was applied to dataset of transcriptional response of irradiated and non-irradiated yeast colonies, and provided the accurate and interpretable discriminative model that may lead to new biological insights. Mi et al. [73] developed a disease classification model and demonstrated that the disease classification models based on gene network enable us to look at diseases in the viewpoint of commonalities in etiology and pathology.

In this study, we have introduced a novel computational strategy for gene network analysis and multi-class classification, i.e., method for gene regulatory network-based multi-class classification. The proposed GRN-multiClassifier aims to estimate the gene regulatory network while simultaneously minimizing both the error in gene network estimation and the negative log-likelihood linked to the multinomial logistic regression model. This indicates that the gene network obtained through our approach is the fine-tuned network designed for optimal multi-class classification. This enables us to achieve a biologically meaningful interpretation of gene network analysis and classification, as the estimated gene network is optimized to elucidate the cancer-related status of cell lines. We can expect through our results and the literature that the proposed network-based classifier will be a useful to biologically reliable and interpretable classification, because our method can explain the status and subtype of cancer cell lines based on the not only expression levels of genes but also their interplay.

To demonstrate the effectiveness of the proposed strategy, we conduct Monte Carlo simulations. The simulation results clearly show that the proposed GRN-multiClassifier outperforms methods relying on pre-estimated gene networks as well as methods that use gene expression levels for multi-class classification. Additionally, our strategy yields effective results in terms of feature selection for multinomial logistic regression and the selection of edges in network estimation.

We apply the GRN-multiClassifier to the classification of acute leukemia cell lines based on gene networks. Our approach demonstrates superior performance in classifying acute leukemia cell lines across all three distinct types. Through the interpretation of the optimized estimated networks for acute leukemia cell line classification, we identify essential markers and fundamental molecular interactions that play crucial roles in achieving accurate classification. Integrating information from prior research, the identified markers demonstrate substantial evidence that underlines their significance in mechanisms associated with acute leukemia. Our findings imply that ACTB may potentially play a pivotal role in the context of acute leukemia. Furthermore, the suppression of ACTB and the molecular interactions involving pairs such as “HBA1&HBB", ”HBB&HBA1", “IGKV1-5&IGHV4-31", “IGHV4-31&IGKV1-5", “HLA-DRA&CD74", and “A CTB&ACTB" might provide essential insights into unraveling the intricate mechanisms of acute leukemia. These insights could extend to implications for both preventive and therapeutic strategies.

Our strategy provides data-driven results for gene networks analysis, where the network is the optimized molecular interplays to explain status of cell lines. Although our method can estimate the optimized gene network for describing the diseases-related status of cell lines, the use of the known a-priory network can improve interpretability of the network-based classification results. That is, we consider the use of a known a-priory network related a specific disease (e.g., the networks involved in specific biological pathway in the known database) as an initial network in network-based classification, and then the known a-priory network is estimated to optimize for the specific disease status related classification. We consider the known a-priory network-based classification as one of the future works of this study and expect that the use of the known network can provide biologically interpretable and reliable results.

Appendix

The first and second derivatives of the objective function of the GRN-multiClassifier in (13) with respect to $[eqn]$ are given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is $[eqn]$ identify matrix, $[eqn]$ is an n–dimensaional vector and $[eqn]$ and $[eqn]$ are $[eqn]$ diagonal matrices defined as

$$[eqn]$$ $$[eqn]$$

Initiating from an initial value, we achieve a numerical solution by employing the subsequent update formula:

$$[eqn]$$

The update formula is referred to as Fisher’s scoring algorithm, and the $[eqn]$ estimator $[eqn]$ is updated by

$$[eqn]$$

where $[eqn]$ [14].