EgoNet identifies differential ego-modules and pathways related to prednisolone resistance in childhood acute lymphoblastic leukemia

ABSTRACT

Purpose: To extract feature ego-modules and pathways in childhood acute lymphoblastic leukemia (ALL) resistant to prednisolone treatment, and further to explore the mechanisms behind prednisolone resistance.

Materials and methods: EgoNet algorithm was used to identify candidate ego-network modules, mainly via constructing differential co-expression network (DCN); selecting ego genes; collecting ego-network modules; refining candidate modules. Afterwards, statistical significance was calculated for these candidate modules. Biological functions of differential ego-network modules were identified using Reactome database. To verify this proposed method can lead to truly positive findings in clinical settings, support vector machine (SVM) was utilized to compute the AUC values for each significant pathway using 3-fold cross-validation method. To predict the reliability of our findings, another established method (attract) was used to identify the differential attractor modules using the same microarray profile.

Results: After eliminating the modules with classification accuracy < 0.9 and node number < 15, only ego-network module 30 was eligible. After significance calculation, module 30 was significant. Module 30 was enriched in APC/C-mediated degradation of cell cycle proteins. The AUC for the significant pathway of APC/C-mediated degradation of cell cycle proteins was 0.915. Although the attract method obtained more modules, the module identified by our proposed method owned more gene nodes, and had more classification ability (AUC = 0.915).

Conclusion: One differential ego-network module identified in childhood ALL resistance to prednisolone based on DCN and EgoNet, might be helpful to reveal the mechanisms underlying prednisolone resistance in childhood ALL.

KEYWORDS:

• Childhood acute lymphoblastic leukemia
• EgoNet
• prednisolone resistance
• module
• pathway

1. Background

Acute lymphoblastic leukemia (ALL), as an acute form of leukemia, or cancer of the white blood cells, is characterized by the accumulation and overproduction of cancerous, immature white blood cells [1]. Of note, ALL is most common childhood cancer, with a peak incidence at 2–5 years old [1]. Treatment of childhood ALL includes the application of several classes of chemotherapeutic drugs, such as glucocorticoids and anthracyclines. Prednisolone, as a glucocorticoid, can diffuse passively into the cell. In the treatment of pediatric ALL, prednisolone is used as the spearhead agent, because of its ability to inhibit cell cycle progression and induce apoptosis in ALL cells [2,3]. Over the past four decades, enormous improvements have been achieved in ALL treatment with successful disease-free survival rates of about 80%. Although the significant successful survival rate, the remaining 20% of ALL patients still show treatment failure during or after completion of contemporary treatment protocols for newly diagnosed pediatric ALL [4]. In addition, several surviving patients frequently present toxicity [5]. Notably, relapsed ALL patients often obtain prednisolone resistance [6]. Significantly, resistance towards prednisolone is closely related with worse prognosis in pediatric ALL. Thus, therapeutic intervention to conquer prednisolone resistance is urgent. Nevertheless, up to now, the molecular mechanisms behind prednisolone resistance in childhood ALL is still poorly defined.

Recently, several studies have documented that gene expression profiles can identify treatment-specific changes in ALL cells [7,8]. Using gene expression profiles, Holleman et al. [9] have identified several biomarkers predictive of outcome in ALL patients, and these biomarkers are associated with the resistance or sensitivity to chemotherapeutic agents, for example, prednisolone and daunorubicin. The gene profile of GSE19143 offered by Stam et al. [10] have identified differentially expressed genes (DEGs) from ALL samples either resistant or sensitive to prednisolone. Using microarray data of GSE19143, Sarang et al. [11] have investigated the altered expression of autophagy-related genes which might contribute to glucocorticoid-resistant phenotype in glucocorticoid-sensitive and glucocorticoid-resistant ALL cells. However, these studies only focus on the individual genes. In the current analysis, complementally, EgoNet strategy [12] was utilized to trace the dysfunctional ego-network modules associated with the resistance of prednisolone.

Thus, in our study, with the goal of investigating the gene expression pattern that might result in prednisolone resistance in more detail, EgoNet was used to extract feature ego-modules and pathways in childhood ALL resistant to prednisolone treatment. Specifically, ego-network modules were derived the origin from ego genes in the network based on snowball sampling, and pathway enrichment analysis of differential modules was also carried out. Candidate modules and pathways identified by this novel method might provide the groundwork for a therapy for childhood ALL.

2. Materials and methods

2.1. Data pre-processing

Gene expression profile, recruited from the GPL96 platform of Affymetrix Human Genome U133A Array, for 52 children diagnosed with ALL either sensitive (27 patients with LC50 values <0.1 lg/ml prednisolone) or resistant (25 patients with LC50 values >150 lg/ml prednisolone) to in vitro prednisolone treatment were downloaded from Gene Expression Omnibus repository (accession no. GSE19143).

For quality control, we first preprocess the gene expression data of GSE19143 prior to analysis. In detail, the annotation data in TXT format was employed to map the association between probes and genes. Next, RMA [13] was used to implement background adjustment, and quantiles were employed to carry our normalization [14]. Then, PM/MM correction was implemented through Microarray Suite (MAS) [15], following by summarization of the expression measures via median polish. Eventually, the probe sets were mapped to gene symbols. After that, an expression profile including 12 442 genes was obtained.

2.2. Protein–protein interactions data

Protein–protein interactions (PPIs) are very important for understanding the biological processes. Thus, compiling PPI networks provides a valuable framework for a better understanding of the protein function. In the current study, the original PPI network involved in 787 896 interactions (16 730 genes) was downloaded from STRING database [16]. Then, the common interactions between the original PPI network and microarray data were extracted. Finally, a new PPI network was obtained, which included 41 244 interactions and 6722 genes.

2.3. Identification of ego-network modules

In order to identify the candidate modules, EgoNet method was applied. This approach included the following four key steps: construction of differential co-expression network (DCN); selection of ego genes; collecting ego-network modules; and refinement of candidate modules.

2.3.1. Construction of DCN

First, the absolute value of Pearson correlation coefficient (PCC) was computed between the interactions of gene i and gene j. The interactions whose absolute value of PCC were higher than the pre-defined threshold δ (δ = 0.8) were picked out to build the co-expression network.

Second, we used one-side t-test to calculate p values of each gene expression between the two groups. Based on the p values, we assigned the weight values to the edges in the binary co-expression network.

2.3.2. Selection of ego genes

In this section, we used a function g(s) to calculate the importance of node s in the DCN. The function g(s) was shown as follows:

In which, g(s) was the importance of node s in the DCN, N_k(s) was on behalf of the neighbors of gene s; was the degree normalized weighted adjacency matrix; and D stood for diagonal matrix. For each gene, g(s) = z-score was obtained, and the genes were sorted according to z-score. The top 5% ranked genes with degree >1 were identified as the ego genes.

2.3.3. Collecting ego-network modules

This step was carried out with the help of ego genes. Taking each ego gene as the initial node, the ego-network module stage exhaustively sought all genes x whose inclusion could cause the maximum increase of the classification accuracy till there was no increase. Regarding a given ego gene x ∈ X, we treated it as an ego-network module C = {x}. For each vertex y in its neighborhood in the network, N(x) was defined as U_iN_i(x), where N_i(x) was the neighbor set as the candidate for ego-network module C. For each y ∈ N(x), a new ego-network module C' = C{y}. The classification accuracy increase between C' and C was calculated based on the formula △S (C', C) = S(C)-S(C').

When △S (C', C) was greater than 0, which suggested the inclusion of gene y enhanced the categorizing accuracy of the ego-network module C. If so, the gene y was added to ego-network C. If there were more than one node which were contained at each stage, we randomly selected one. We would stopped this spreading stage if the classification accuracy dropped. As known, the above process of growing ego-network module was snowball sampling [17].

2.3.4. Refinement of candidate modules

In this step, we removed the ego-networks with classification accuracy smaller than 0.9. Moreover, the ego-networks with gene size <15 were also eliminated.

2.4. Estimating statistical significance of modules

Based on the null distribution of classification accuracy generated by random permutation test, we calculated the statistical significance of ego-network modules. With the goal of obtaining the null distribution, randomized permutation test on the same ego-network modules was implemented for 1000 times. Based on the null distribution, the p value of an ego-network module was calculated as the probability of the candidate ego-network module having the smaller classification accuracy value. Next, Benjamini and Hochberg method [18] was utilized to correct p values. In our study, false discovery rate (FDR) of 0.01 was set as the significance criteria.

2.5. Functional profiling

In order to gain insight into the biological processes of the given ego-network modules, the biological functions of differential ego-network modules were identified relying on Reactome pathway databases. In the current analysis, we recruited all reference pathways in human from Reactome database. Then, the common genes of each reference pathway and the new PPI network were extracted. Next, several informative pathways were selected for subsequent analysis when the pathways with less than 5 genes or more than 100 genes were removed. Afterwards, genes of differential ego-network modules were mapped to each informative pathway to identify the pathways enriched by each differential ego-network module. Statistically, significant enrichment of module genes within terms was determined by Fisher’s exact test adjusted by Benjamini–Hochberg method for multiple testing. In our paper, pathways with an FDR <0.01 were defined as the pathways enriched by a given module. Subsequently, the pathways of each differential ego-network module were sorted based on FDR values, and the pathway having the smallest FDR was considered as the significant pathway for a given ego-network module.

2.6. Validation test

Through the literature, area under the curve (AUC) is an available index which is frequently employed to assess the predictive ability of machine learners [19]. In order to verify that this proposed method can lead to truly positive findings in clinical settings, support vector machine (SVM) was utilized to compute the AUC values for each significant pathway using threefold cross-validation method. The AUC of 1.0 devotes a perfect classification ability, and the AUC value >0.7 is regarded good [20].

2.7 Comparison to other established method (attract method)

To predict the reliability of our findings, another established method (attract) was used to identify the differential attractor modules using the same microarray profile. Specific steps were described in the study provided by Wang et al. [21]. Then, pathway analysis was conducted for the genes in the attractor modules based on FDR < 0.01.

3. Results

3.1. A framework for detection of ego-network modules

Taking the PPI network and microarray profile as input, EgoNet repeatedly sought all genes with two or more neighbors in the network. An ego-network module is a part of a network which refers to a particular gene, and we call this ego gene. In addition to the ego gene, the network also includes a neighborhood containing all nodes to which the ego is connected to at a certain path length. Starting with every ego gene as initial node, the score of the level-one ego-network module was calculated based on classification accuracy. Then, it spread outward from the ego gene progressively to include more genes in the predictive model. This process quitted when the classification accuracy reduced. We defined this procedure of growing ego-network as snowball sampling. Afterwards, the significance was assessed by randomized permutation test. This algorithm mainly included four steps: DCN construction, ego genes selection, ego-network modules identification via snowball sampling, as well as refinement of ego-network modules.

3.2. Construction of DCN

In an attempt to construct the DCN, edges were picked out according to the absolute value of the PCC of the expression profiles of two genes. Totally, 605 genes (845 edges) were selected to construct the DCN when δ was greater than 0.8. The composition of the DCN was shown in Figure 1. However, this DCN was so big that the composition of the DCN was not seen clearly. Thus, we have provided the specific content in form of Excel (Supplementary Table 1), where genes and interactions were included.

Figure 1. The composition of DCN which included 245 nodes and 845 interactions. Yellow nodes stood for the ego genes.

3.3. Identification of ego-network modules

According to the z-score distribution of each node in DCN, a total of 30 ego genes were identified, as shown in Table 1. Among these ego genes, the top five were respective RPS16 (z-score = 28.853), NUSAP1 (z-score = 26.205), KIF11 (z-score = 24.456), RPL27 (z-score = 24.248), and NCAPG (z-score = 23.984).

Table 1. The 30 ego genes extracted from the DCN and their z-scores.

Ego genes	z-score	Ego genes	z-score
RPS16	28.853	U2AF2	16.961
NUSAP1	26.205	RPS24	16.676
KIF11	24.456	RPS23	16.487
RPL27	24.248	RPL9	15.785
NCAPG	23.984	SNRNP70	15.707
KIAA0101	20.788	UBA52	15.529
RAD51AP1	20.719	RPLP1	14.997
PABPN1	20.145	RPL22	14.727
RPL19	19.693	RPS6	14.399
MAD2L1	19.366	SF3A1	14.223
EIF3E	18.034	RPS27	14.197
RPS12	17.970	RPS7	13.753
RPL38	17.737	FUS	13.668
RPS10	17.572	MKI67	13.582
CCNA2	17.269	RPLP0	13.369

Next, starting with every ego gene, we progressively grew the ego-network modules by levels on the basis of the classification accuracy increase △S (C', C). Based on △S (C', C) > 0, a total of 30 ego-network modules were screened out and the mean gene size in an ego-network module was 7. After eliminating the modules with classification accuracy <0.9 and the number of nodes <15, we observed that only module 30 was remained.

3.4. Evaluating the significance of candidate ego-network modules

After random permutation test with 1000 times, we found that the ego-network module 30 was significant with FDR = 0. In order to better display the properties of this module, we presented the composition of the ego-network module 30, shown in Figure 2.

Figure 2. Composition of differential ego-network module 30 identified in the childhood ALL resistant to prednisolone treatment. Yellow node was ego gene.

3.5. Functional profiling and validation test

With the goal of better understanding of the module functions, we performed the pathway enrichment analysis based on the Reactome database to extract cellular functions specifically enriched by these differential ego-network modules. A set of pathways were selected for the given module 30 using a threshold FDR of 0.01. For this given module, we only selected the pathway with the smallest FDR as the statistically enriched pathway of a given module. Thus, we found that the genes of module 30 were enriched in the pathway of APC/C-mediated degradation of cell cycle proteins (FDR = 1.84E-09).

In order to assess the classification performance, we then used SVM to compute the AUC values for each significant pathway using threefold cross-validation method. Based on the results, we found the AUC score for the significant pathway of APC/C-mediated degradation of cell cycle proteins was 0.915, which further demonstrated that this method actually could lead to truly positive findings in clinical settings.

3.6. Comparison to other established method (attract method)

In order to prove the reliability of our results, we subsequently used attract method to identify the differential attractor modules. Based on FDR < 0.05, only two modules were identified to be differential (11 nodes in Module 1, 8 nodes in Module 2). The composition of these two attractor modules and pathways enriched by the module genes are shown in Table 2. Although the attract method obtained more modules, the module identified by our proposed method owned more gene nodes, and had more classification ability (AUC = 0.915).

Table 2. Composition of the two attractor modules and pathways enriched by the module genes.

Modules	Gene composition	Pathways
Module 1	RPN2, RPS7, RPL27A, RPL27, RPL8, EIF3D, EIF3J, SMG7, UPF1, ZNF268, POLE	Formation of a pool of free 40S subunits
Module 2	SRSF11, SNRNP70, SNRPA, CD2BP2, FUS, DHX38, GTF2F2, RNMT	mRNA splicing

4. Discussion

Although the advances have made in ALL treatment, about 20% of patients still experience a relapse and death [22]. Glucocorticoids resistance has been recognized to be a feature of poor prognosis in the treatment of childhood ALL. However, taking consideration of the importance of prednisolone in ALL treatment, it is profound that identifying discriminating biomarkers possibly associated with ALL treatment outcome, particularly those that could better explain the molecular mechanisms behind drug resistance [23,24]. Currently, as a relatively unbiased approach, genome-wide analysis of gene expression in human ALL cells is used to extract gene determinants of drug resistance [9,25]. However, many of the previous studies have focused on the identification of a set of signatures of DEGs [9,26,27]. Moreover, many of the gene signatures obtained from different researches on the same disease are frequently inconsistent [28]. Thus, to overcome this difficulty, network-based approaches are potentially more effective to evaluate the disease-related biomarkers because it addresses the organization of networks of pairwise biological interactions, instead of the individual gene [29,30]. Unfortunately, if a network is too big, a certain amount of important genes and interactions are neglected [31]. Identification of modules or sub-networks from the intricate network can avoid this problem [30,32]. Many approaches have been created to seek modules by jointly analyzing network data with gene expression profiles. For example, Vert and Kanehisa [33] utilized kernel methods to screen out the expression patterns which characterize gene sets matching pathways in a given network. Nevertheless, some hidden genes showing no difference by themselves but clustered in a module are often overlooked. Excitedly, a novel method called EgoNet [12] can identify significant modules and the hidden genes in this module whose genes are highly predictive of the disease condition.

In our study, we also used this EgoNet strategy to trace the dysfunctional ego-network modules which might be associated with the resistance of prednisolone. There were 605 nodes in DCN, and among these nodes, 30 were ego genes. Selection of differential ego-network modules was expanded on DCN and the corresponding ego genes via EgoNet algorithm. After significance calculation, we found that only module 30 was significant with FDR = 0. Pathway analysis inhibited that dysregulation of the APC/C-mediated degradation of cell cycle proteins might stand for the key in resistance or sensitivity to prednisolone. Based on the AUC results, we found the AUC score for the significant pathway of APC/C-mediated degradation of cell cycle proteins was 0.915, which further demonstrated that this method actually could lead to truly positive findings in clinical settings.

As reported, the glucocorticoids prednisolone exerts a key function in therapy protocols, since it has the ability to block cell cycle progression in ALL cells [34]. Drug-resistance genes involved in cell cycle progression have been related with the prognosis of ALL [35]. Theoretically, prednisolone can diffuse passively into the cell via binding the glucocorticoid receptor (GR). Generally, GR is defined as a cytosolic receptor, while Gametchua et al. [36] described a membrane-bound variant (mGR). The overexpression of mGR has been observed in the period of S-G2/M stage, when the cells are sensitive to the glucocorticoids-induced apoptosis [37]. MIB-1, as the marker of chemotherapy resistance, is also known as the marker of proliferation Ki-67 (MKI67). Of note, MKI 67 involved in cell cycle regulation was up-regulated in ALL [38]. In the current study, MKI67 was the ego gene in the differential ego-network module 30, and this module was enriched in the pathway of APC/C-mediated degradation of cell cycle proteins. In light of these, we infer that targeting cell-cycle-related pathway might represent a promising possibility in the treatment of prednisolone-resistant childhood ALL.

5. Conclusion

To summarize, we have successfully extracted one differential ego-network module in childhood ALL resistance to prednisolone based on DCN and EgoNet algorithm, this module might be helpful to reveal the molecular mechanisms underlying prednisolone resistance and give insight for future research of childhood ALL. Targeting cell-cycle-related pathway might represent a promising possibility in the treatment of prednisolone-resistant childhood ALL.

Disclosure statement

No potential conflict of interest was reported by the authors.