Introduction
============

Leukemia is a prevalent hematologic malignancy and one of the most common causes of cancer deaths in the developed countries[@B1], [@B2]. The overall incidence of leukemia is 14 per 100000 people in the United States in 2015 and is projected to continue rising. Based on the origin, leukemia can be classified into myeloid leukemia or lymphoid leukemia, which can be subdivided into acute or chronic according to the degree of cellular differentiation[@B3], [@B4].

Many of the symptoms of leukemia are non-specific and vague, which could not be diagnosed by conventional blood tests and bone marrow examination[@B5], [@B6]. Plenty of efforts have been devoted to investigate the molecular alterations in leukemogenesis. Next generation sequencing of human genomes and exomes has revealed somatic mutations, aberrantly expressed genes, microRNAs and DNA methylations with putative roles in leukemia[@B7]-[@B9]. However, most of the individual molecules suffer from low reproducibility and high false-positive rates. Few of them have been translated to the clinic for diagnostic application.

It is well recognized that cancer is a complex disease caused not by the malfunction of single molecules but their collective behavior in the network [@B10]-[@B15]. Therefore, network biomarkers are considered to better characterize leukemia than individual molecules and have recently attracted much attention. A number of protein interaction sub-networks have been proposed for early diagnosis, prognosis and efficacy prediction of cancers[@B16]-[@B19].

In this study, we proposed a framework (Figure [1](#F1){ref-type="fig"}) that integrates protein-protein interaction (PPI) data and microarray-based gene expression profiles to construct network biomarkers for accurate prediction of leukemia. The network biomarkers prove to be effective in distinguishing leukemia from normal samples.

Materials and Methods
=====================

Data collection
---------------

We used two different types of datasets, protein-protein interaction data and disease annotation of the protein-coding genes to reconstruct the leukemia-specific PPI network. PPI data was extracted from the Protein Interaction Network Analysis (PINA) v2.0 platform [@B20]. PINA is a unified database of protein-protein interaction

that collects 14454 genes and 108470 interactions from six manually curated public databases (listed in Table [1](#T1){ref-type="table"}). The leukemia-associated genes were extracted from the commercial knowledge database Metacore^TM,^ which is developed by GeneGo.

The public gene expression data were downloaded from the Gene Expression Omnibus (GEO) database. All the gene expression data were obtained using Affymetrix Human Genome arrays. The samples in each GEO datasets are divided into three categories: Leukemia (including AML, CLL, T-PLL and B-CLL), others and Normal. The others samples are filtered out in this study since they are not associated with leukemia. Detailed information for GEO datasets is summarized in Table [2](#T2){ref-type="table"}. The six groups of expression datasets were analyzed to get statistics values. Additional three sets of expression datasets were used for further verification (Table [3](#T3){ref-type="table"}).

Reconstruction of leukemia-specific PPI network
-----------------------------------------------

Human leukemia-specific protein-protein interaction network was first downloaded from PINA and then refined with the 1495 leukemia-associated gene from GeneGo. Only the interactions formed between leukemia-associated genes were selected to form a leukemia-specific PPI network.

Integration with gene expressing profiles
-----------------------------------------

The statistical analysis was invoked through the limma (Linear Models for Microarray Data) R package [@B36] and the affy(Methods for Affymetrix Oligonucleotide Arrays) R package in R software platform [@B37]. Student t-test was used to identify the significant difference level (P-value) of each considerable gene in each dataset. To enhance the accuracy, we applied the empirical Bayesian statistical method to moderate the standard errors and then utilized the method proposed by Benjamini et al. to adjust the multi-testing [@B38], and got the adjusted P-values simultaneously.

To integrate the gene expression data and leukemia-specific PPI network, the adjusted P-value of each gene was mapped onto its corresponding gene in the leukemia-specific PPI network to obtain a dataset-specific weighted PPI network, with adjusted P-value as weight factor.

Active module subtraction
-------------------------

In general, the network integration analysis was performed in 3 steps. At the first step, we converted the adjusted P-values to Z score through using the inverse normal cumulative distribution. Higher Z score indicates more important role in leukemogenesis. Given the Z score, we performed a greedy search to identify the modules with a locally maximal Z score. The candidate modules were seeded with a single gene and then a neighbor within a distanced=3 from the seed were iteratively added. If the neighbor added to the Z score, it was incorporated into the module. The search terminated when no addition increased the Z score over the improvement rate r. The parameter r was set as 0.05 to avoid over fitting. At last the top 10 modules with the highest Z-score identified from each run were merged and iteratively searched for 3-5 times, until the module reached the optimal size of 70-80 nodes. We used jActiveModules [@B39] to select active modules from the weighted PPI network since it is a fashionable method for this kind of investigation. jAM is a plug-in of Cytoscape which evaluated module activity with Z score.

Network-based biomarkers construction
-------------------------------------

At last, as 6 optimized modules include 290 genes in total, which are too large and loosely interconnected for further analysis, we carried out the overlapping analysis to find out the number of enriched genes shared by each optimized modules. We overlapped the six modules and selected the genes shared by at least two networks to construct the final network-based biomarker.

Pathway enrichment analysis
---------------------------

We performed pathway enrichment analysis in Ingenuity Pathway Analysis (IPA) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [@B40]to provide functional insight into the identified network marker. The statistical significance of the enrichment was calculated using hypergeometric test and adjusted by FDR method (P-value \< 0.05).

Statistical significance assessment
-----------------------------------

Hypergeometric test was used to test whether the network biomarkers are significantly enriched with leukemia-related genes. Known mutation genes related to cancers were obtained from the Catalogue of Somatic Mutations in Cancer (COSMIC), which is a cancer gene census [@B41]. 213 of the COSMIC genes are found in common with the GeneGO database. An empirical P-value was calculated to evaluate the statistical significance. P-value was obtained according to the following equation:

Where, N represents the number of genes in the leukemia-specific PPI network; M is the number of known leukemia related genes in COSMIC; n denotes the number of genes in the final network biomarkers; k represents the known leukemia related genes in the final network biomarkers.

Performance evaluation
----------------------

We employed the receiver-operating characteristic (ROC) analysis to evaluate the prediction performance of the network biomarkers in distinguishing leukemia samples from the normal controls. The epicalc R package (<http://CRAN.R-project.org/package=epicalc>) was used to produce the ROC curves. A 5-fold cross validation was performed on three gene expression dataset listed in Table [3](#T3){ref-type="table"}. Normal samples were set as 0 and cancer samples were set as 1. The classification performance was represented as the area under curve (AUC). We also provided sensitivity, specificity and accuracy for the network biomarkers.

Results and Discussion
======================

Sub-network involved in leukemogenesis
--------------------------------------

The leukemia-specific PPI network was reconstructed by integrating PPI from PINA and 1495 leukemia-associated genes from GeneGo. As a result, the leukemia-specific PPI network consists of 4136 interactions among 978 genes.

As is described in Methods section, gene expression profiles of 6 independent GEO datasets were overlaid to the reconstructed leukemia-specific PPI network and 6 correspondent dataset-specific sub-networks were obtained (marked orderly as PPI_GSE9476_raw, PPI_GSE6691_raw, PPI_GSE9476_raw, PPI_GSE22529_raw, PPI_GSE23293_raw and PPI_GSE26725_raw in Figure [1](#F1){ref-type="fig"}). Greedy search was performed for 6 sub-networks respectively. After 3 iterations for GSE5788, GSE9476 and GSE22529, 4 iterations for GSE23293 and GSE26725, 5 iterations for GSE6691, finally we obtained 6 active modules with a locally maximal Z score by jActiveModules (marked orderly as PPI_GSE9476_TR, PPI_GSE6691\_ TR, PPI_GSE9476_TR, PPI_GSE22529_TR, PPI_GSE23293_TR and PPI_GSE26725_TR in Figure [1](#F1){ref-type="fig"}). The number of nodes and edges in each module is summarized in Table [4](#T4){ref-type="table"}.

After overlap analysis, a total of 97 genes along with their interactions were incorporated into the final network-based biomarkers, as illustrated in Figure [2](#F2){ref-type="fig"}. Genes with previous evidence in leukemia are marked yellow.

Functional analysis of candidate network biomarkers
---------------------------------------------------

The network biomarkers were most enriched for molecular mechanisms of cancer (IPA) and pathways in cancer (KEGG). Leukemia-specific pathways such as Chronic Myeloid Leukemia (KEGG) and Acute Myeloid Leukemia Signaling (both IPA and KEGG) were also enriched and showed high statistical significance. It indicates that genes in the biomarker network are closely associated with the development of different types of leukemia. Besides, in He\'s study, P13K/AKT Signaling (IPA) was also proved to be involved in chronic myeloid leukemia[@B42]. Irwin et al. found that ErbB inhibitors played important roles in Philadelphia chromosome-positive acute lymphoblastic leukemia (Ph(+)ALL) and ErbB signaling (KEGG) was a complementary molecular target in Ph(+)ALL [@B43]. The top-ranked pathways in both IPA and KEGG displayed apparent correlation between leukemia and the network biomarkers, which implied the potential accuracy of our result. Figure [3](#F3){ref-type="fig"} shows the top 10 most significantly enriched IPA and KEGG pathway respectively.

We used the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [@B44] for the Gene ontology (GO) annotation in three domains: molecular function, biological process, and cellular component. The top 10 most significantly enriched items for each domain are shown in Figure [4](#F4){ref-type="fig"}. These results indicate that genes in the network are closely associated with the biological processes in the development of different types of leukemia, such as cell death [@B45] and apoptosis [@B46]. This indicated the accuracy of the predicted network biomarkers to a certain extent.

Network biomarkers are significantly associated with leukemia
-------------------------------------------------------------

We further investigated whether the genes in the network biomarkers were randomly obtained. The statistical significance was checked using hypergeometric test and a significant p-value of 0.008987933 was obtained. This indicates that the candidate network biomarkers are enriched with known leukemia-related genes and could not be obtained randomly.

As illustrated in Figure [5](#F5){ref-type="fig"}(A), the blue circle represents the 978 genes in the leukemia-specific PPI network; the red circle includes the 522 known leukemia-related genes in COSMIC. The leukemia-specific PPI contains 195 known leukemia-related genes in COSMIC. The purple circle represents 97 genes in final network biomarkers, among which 29 genes belong to the known leukemia-related category.

Sub-network marker with higher classification accuracy
------------------------------------------------------

To evaluate the performance of network biomarkers in classifying leukemia and normal gene expression profiles, we used three independent gene expression datasets listed in Table [3](#T3){ref-type="table"} as tested datasets to produce the ROC curves. We compared the network biomarker with three reported gene biomarkers: CD38[@B47], BCL2 [@B48] and IGFBP7 [@B49]. The reasons we chose these three markers for comparison are as follows, 1) these biomarkers are all well-studied and all of them have been validated by clinical experiments. 2) The marker CD38 is a member of our network whereas the remaining two are not. We included two others for fair evaluating the performance of our network biomarker. Figure [5](#F5){ref-type="fig"} shows the ROC curves for network biomarkers and 3 known biomarkers. Network-based biomarker has higher AUC than any of the single markers which means network-based biomarker could more effectively discriminate the leukemia from the normal controls. The sensitivity, specificity and accuracy of each dataset are given in Table [5](#T5){ref-type="table"}.

It is worth noting that for network biomarkers from GSE8835 has a relatively lower AUC than the other two datasets. This may be caused by the difference of platform and method between GSE8835 and the training datasets. Anyhow, the accuracy of the network biomarkers is still higher than other three single biomarkers.

The result indicates that the putative network biomarkers could diagnose leukemia samples more accurately and could be used as putative biomarker to aid in early diagnosis of leukemia.

Conclusions
===========

In conclusion, we developed a network approach for molecular investigation and diagnosis of leukemia. The constructed network biomarkers not only achieve higher accuracy rate of diagnosis compared to known single biomarkers but also provide systematic insights into the leukemogenesis process. We noticed that we only considered the combination of genes (or the nodes) in the network for the prediction of leukemia. The interactions among genes can also provide valuable biological signatures for diagnosis of diseases. We will take the edge-variation in the network into the account for the further improving of the leukemia prediction.

Availability of Data and Materials
==================================

PINA: <http://cbg.garvan.unsw.edu.au/pina/interactome.stat.do.>

GeneGo: <https://portal.genego.com/cgi/data_manager.cgi.>

GSE9476: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9476.>

GSE6691: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6691.>

GSE5788: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE5578.>

GSE22529: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22529.>

GSE26725: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26725.>

GSE23293: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE23293.>

GSE8835: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE8835.>

GSE24739: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE24739.>

GSE39411: <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39411.>

COSMIC: <http://cancer.sanger.ac.uk/cosmic/download>

IPA: <http://www.ingenuity.com/products/login>

DAVID: <https://david.ncifcrf.gov/.>

BS, XY and JC were supported by National Natural Science Foundation of China (NSFC) (grant no. 31670851, 31470821, 31400712, 91530320).

PPI

:   Protein-protein interaction.

PINA

:   Protein Interaction Network Analysis.

GEO

:   Gene Expression Omnibus.

IPA

:   Ingenuity Pathway Analysis.

KEGG

:   Kyoto Encyclopedia of Genes and Genomes.

COSMIC

:   Catalogue of Somatic Mutations in Cancer.

ROC

:   Receiver-operating characteristic.

AUC

:   area under curve.

DAVID

:   Database for Annotation, Visualization, and Integrated Discovery.

GO

:   Gene ontology.

![The flowchart of network biomarkers identification for leukemia diagnosis.](jcav08p0278g001){#F1}

![**The final network-based biomarker for leukemia.** The known cancer related genes in final network are marked yellow.](jcav08p0278g002){#F2}

![**IPA and KEGG pathway enrichment analysis for network biomarkers.**The top 10 most significantly enriched IPA and KEGG pathway are shown in panel (A) and (B) respectively.](jcav08p0278g003){#F3}

![**Gene ontology annotation for the network biomarkers.** The network biomarkers identified by our method were annotated with DAVID tools at three levels of gene ontology: Molecular Function, Biological Process, and Cellular Component. The top 10 most significantly enriched items for each level are shown.](jcav08p0278g004){#F4}

![**Validation of the network biomarkers.**(A) Distribution of the leukemia-associated genes in the network. (B-D) ROC curves obtained with the network biomarkers tested by three gene expression datasets. Panel (B), (C) and (D) represent respectively the results of the gene expression datasets in series of GSE8835, GSE24739 and GSE39411.](jcav08p0278g005){#F5}

###### 

Source databases of PINA.

  Original database   Version         Ref.     Link
  ------------------- --------------- -------- -----------------------------------------------
  IntAc               Oct 4,2012      \[21\]   <http://www.ebi.ac.uk/intact/>
  BioGRID             3.1.93          \[22\]   <http://thebiogrid.org/>
  MINT                Dec 21,2010     \[23\]   <http://mint.bio.uniroma2.it/mint/Welcome.do>
  DIP                 June 14,2010    \[24\]   <http://dip.doe-mbi.ucla.edu/dip/Stat.cgi>
  HPRD                April 13,2010   \[25\]   <http://www.hprd.org/download>
  MIPS/Mpact          Oct 1,2008      \[26\]   <http://mips.helmholtz-muenchen.de/>

###### 

Leukemia-associated gene expression datasets used for analysis.

  Series     Platform   No. Samples   Leukemia   Others   Normal   Ref.             
  ---------- ---------- ------------- ---------- -------- -------- ------ ---- ---- --------
  GSE9476    GPL96      64            26                                       38   \[27\]
  GSE6691    GPL96      56                       11                       32   13   \[28\]
  GSE5788    GPL96      14                                6                    8    \[29\]
  GSE22529   GPL96      52                       41                            11   \[30\]
  GSE26725   GPL570     17                                         12          5    \[31\]
  GSE23293   GPL570     41                       7                        18   16   \[32\]

###### 

Leukemia-associated gene expression datasets used for validation.

  Series     Platform   No. of samples   CML   CLL   Normal   Ref.
  ---------- ---------- ---------------- ----- ----- -------- --------
  GSE8835    GPL96      66                     42    24       \[33\]
  GSE24739   GPL570     24               16          8        \[34\]
  GSE39411   GPL570     152                    104   48       \[35\]

###### 

Detailed information of the active modules.

          PPI_GSE5788_TR   PPI_GSE6691\_ TR   PPI_GSE9476\_ TR   PPI_GSE22529\_ TR   PPI_GSE23293\_ TR   PPI_GSE23293\_ TR
  ------- ---------------- ------------------ ------------------ ------------------- ------------------- -------------------
  Nodes   77               71                 75                 73                  71                  75
  Edges   205              186                166                193                 126                 188

###### 

Detailed information of ROC curves.

  Series               Biomarker   Sensitivity   Precision   Specificity   Accuracy   AUC
  -------------------- ----------- ------------- ----------- ------------- ---------- -------
  GSE8835              CD38        0.913         0.700       0.212         0.658      0.629
  BCL2                 0.885       0.650         0.166       0.623         0.604      
  IGFBP7               0.965       0.665         0.148       0.668         0.584      
  Network biomarkers   0.851       0.686         0.316       0.657         0.698      
  GSE24739             CD38        0.893         0.662       0.088         0.625      0.431
  BCL2                 0.943       0.654         0.004       0.630         0.237      
  IGFBP7               0.938       0.652         0.001       0.625         0.197      
  Network biomarkers   0.874       0.986         0.976       0.908         0.966      
  GSE39411             CD38        0.999         0.725       0.177         0.740      0.981
  BCL2                 0.915       0.931         0.853       0.895         0.966      
  IGFBP7               0.886       0.797         0.513       0.768         0.880      
  Network biomarkers   0.996       0.999         0.998       0.997         0.999      

[^1]: Conflict of Interest: The authors declare no conflict of interest.
