Multicellular gene network analysis identifies a macrophage-related gene signature predictive of therapeutic response and prognosis of gliomas
Journal of Translational Medicine volume 17, Article number: 159 (2019)
The tumor-associated microenvironment plays important roles in tumor progression and drug resistance. However, systematic investigations of macrophage–tumor cell interactions to identify novel macrophage-related gene signatures in gliomas for predicting patient prognoses and responses to targeted therapies are lacking.
We developed a multicellular gene network approach to investigating the prognostic role of macrophage–tumor cell interactions in tumor progression and drug resistance in gliomas. Multicellular gene networks connecting macrophages and tumor cells were constructed from re-grouped drug-sensitive and drug-resistant samples of RNA-seq data in mice gliomas treated with BLZ945 (a CSF1R inhibitor). Subsequently, a differential network-based COX regression model was built to identify the risk signature using a cohort of 310 glioma samples from the Chinese Glioma Genome Atlas database. A large independent validation set of 690 glioma samples from The Cancer Genome Atlas database was used to test the prognostic significance and accuracy of the gene signature in predicting prognosis and targeted therapeutic response of glioma patients.
A macrophage-related gene signature was developed consisting of twelve genes (ANPEP, DPP4, PRRG1, GPNMB, TMEM26, PXDN, CDH6, SCN3A, SEMA6B, CCDC37, FANCA, NETO2), which was tested in the independent validation set to examine its prognostic significance and accuracy. The generation of 1000 random gene signatures by a bootstrapping scheme justified the non-random nature of the macrophage-related gene signature. Moreover, the discovered gene signature was verified to be predictive of the sensitivity or resistance of glioma patients to molecularly targeted therapeutics and outperformed other existing gene signatures. Additionally, the macrophage-related gene signature was an independent and the strongest prognostic factor when adjusted for clinicopathologic risk factors and other existing gene signatures.
The multicellular gene network approach developed herein indicates profound roles of the macrophage-mediated tumor microenvironment in the progression and drug resistance of gliomas. The identified macrophage-related gene signature has good prognostic value for predicting resistance to targeted therapeutics and survival of glioma patients, implying that combining current targeted therapies with new macrophage-targeted therapy may be beneficial for the long-term treatment outcomes of glioma patients.
Gliomas are among the most malignant cancers, commonly inducing profound and progressive disability and causing death in most cases . Although many cancer cell-targeted therapeutic agents have been developed, intrinsic or acquired resistance to such therapies often emerges during long-term treatment . Therefore, the preexisting or newly developed tolerance of cancer cells to molecularly targeted therapeutic drugs is a main cause of the eventual failure of most existing targeted therapies.
Several forms of cancer cell-intrinsic mechanisms of drug resistance have been revealed, including genetic/epigenetic mechanisms , posttranslational mechanisms [4,5,6], cellular mechanisms [7, 8], and metabolic mechanisms . Recently, an increasing number of experiments indicated that the tumor microenvironment may play important roles in cancer progression and drug resistance . Therefore, uncovering the intercellular interactions between cancer cells and microenvironmental cells is crucial for identifying effective biomarkers for predicting drug resistance and cancer progression, as well as understanding the mechanisms of acquired resistance and prioritizing potential drug targets. Importantly, therapies targeting the tumor microenvironment have been proposed as a promising approach for treating cancers, including gliomas [11, 12]. Several macrophage-targeted therapies for gliomas have been developed , further highlighting the importance of tumor–microenvironment interactions in determining glioma outcome.
Systems biology approaches are powerful for quantitatively studying various forms of drug resistance at multiple scales, which can provide insights into underlying mechanisms and experimentally testable predictions . Computational methods have also been used to identify prognostic biomarkers and predict drug resistance. However, conventional methods for glioma biomarker identification often focus on individual cell types or molecules rather than on the interactions between different cell types and molecules. Exploring the interactions between tumor cells and microenvironmental cells, such as macrophages, from a network view can more insightfully help understand the mechanisms that underlie cancer progression and drug resistance and discover more robust and accurate biomarkers .
In this study, we developed a multicellular gene network-based approach to investigating intercellular gene associations between tumor cells (TCs) and tumor-associated macrophages (TAMs) and to identify biomarkers of prognosis and drug resistance in gliomas. We used RNA-seq data from mice bearing gliomas treated with BLZ945 to construct drug-sensitive and drug-resistant multicellular gene networks. Based on the differential network, a macrophage-related gene signature was discovered using a dataset of glioma patients from the Chinese Glioma Genome Atlas (CGGA) database. An independent dataset from The Cancer Genome Atlas (TCGA) database was used for validation. Time-dependent receiver operator characteristics (ROC) analysis was used to assess the prognostic significance and accuracy of the gene signature for predicting survival of glioma patients. Moreover, the macrophage-related gene signature was found to be predictive of the targeted therapy outcome in glioma patients and outperformed existing gene signatures including conventional EGFR signature [15, 16] and an immune-related gene signature , which were verified using the validation dataset. Moreover, we showed that the network perturbation analysis improved the model-building process beyond the conventional methods of identifying signature genes based on differential expression.
Preclinical RNA-seq data analysis
RNA-seq data (FPKM) of TCs and TAMs in mice gliomas were downloaded from the Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE69104. The preclinical experiment  investigated the role of TAMs in glioma immunotherapy with BLZ945, a CSF1R inhibitor, where all TC-TAM paired samples were divided into 3 groups, i.e., Vehicle (Veh, 5 samples), Endpoint (EP, 6 samples, i.e., drug-sensitive), and Rebound (Reb, 4 samples, i.e., drug-resistant) tumors (Additional file 1: Table S1).
As described in detail in Additional file 1: Text S1, the expression level of each gene in the drug-resistant and drug-sensitive samples were normalized to their mean values in the vehicle samples. We selected significantly differentially expressed genes (DEGs) between the drug-resistant samples (Reb) and the drug-sensitive samples (Ep) for both TCs and TAMs by calculating the fold changes (FCs) and false discovery rate (FDR)-adjusted p values using t test. A gene with |FC| > 1.5 and adjusted p value less than 0.05 was considered as a DEG, resulting in a list of 1141 DEGs.
Multicellular gene network construction and differential network analysis
The top 50 DEGs with the largest absolute fold change value in each type of cells were selected for gene network construction. We built multicellular gene association networks between TAMs and TCs to investigate the changes between these cell types during the emergence of drug resistance from a network perspective. We classified the above DEGs into TC-specific and TAM-specific DEGs and computed the Pearson correlation coefficients (PCCs) for each pair of these DEGs. An edge is added to the TC-TAM gene network if the corresponding |PCC| > 0.95 and p value < 0.05. In this way, we built two-types of edges, including ‘intracellular edges’ within TCs or TAMs and ‘intercellular edges’ between TCs and TAMs, to construct the multicellular gene network.
Based on the work of one of our authors , we developed a network perturbation analysis technique for multicellular gene networks. As described in detail in Text S1, first we used the reference samples (N = 6) in the Ep group to construct a sensitive network (Fig. 1a) for pairs of correlated DEGs in TCs and TAMs (|PCC| > 0.95 and p value < 0.05). We then added each single drug-resistant sample (Rebi, i = 1, 2, 3, 4) to the reference samples in the Ep group to construct 4 sets of sample-specific perturbed samples, respectively (Additional file 1: Table S1). We used each set of perturbed samples (N = 7) to construct perturbation networks (Fig. 1b). The addition of each single Reb sample is the cause of differences between the sensitive and perturbation networks. If the gene expression profile in the added sample, Rebi, was similar to that in the Ep samples, the perturbation of the PCC was negligible. However, if some gene expression levels were remarkably different between the single Reb sample and the Ep samples, significant changes in the PCCs of certain gene pairs were induced upon the addition of Rebi to the reference samples.
Based on such network perturbation analysis (Fig. 1c), the robust significantly different PCCs (ΔPCC, see Text S1) of gene pairs were chosen to construct a differential network. Specifically, we selected significantly different edges, represented by ΔPCC, by setting a threshold of |ΔPCC| > 0.05 for each perturbation network versus the sensitive network. If a gene-pair was differentially correlated in at least 3 perturbation networks compared to the sensitive network, this edge was selected as a robust differential edge. We defined three types of edges in the differential network: correlation-gained edges (ΔPCC > 0), correlation-lost edges (ΔPCC < 0), and correlation-invariant edges (ΔPCC = 0). The differential network reflected the changes in both intracellular and intercellular edges in the multicellular gene networks of mice with distinct therapeutic responses to the CSF1R inhibition.
Functional enrichment analysis
The functional enrichment of genes in the differential network was analyzed using the newest version of Metascape (http://metascape.org), a gene annotation and enrichment analysis database that integrates ontology sources, including the KEGG Pathway, GO Biological Processes, Reactome Gene Sets, Canonical Pathways and CORUM databases. The whole genome was used as the enrichment background. Terms of pathways or processes with a p value < 0.01 (accumulative hypergeometric test), a minimum of 3 genes, and an enrichment factor > 1.5 were collected and grouped into clusters based on their membership similarities.
Prognostic signature identification and validation
The differential network might capture robust topological differences between the CSF1R inhibitor-sensitive network and inhibitor-resistant networks and reflect potential changes in the gene interactions across TCs and TAMs during the acquisition and development of drug resistance. It is therefore reasonable to speculate that the genes in the differential network play critical roles in promoting tumor growth, even under drug pressure. We hypothesized that the expression levels of genes in the differential network are associated with the survival outcomes of glioma patients. As such, we developed a differential network-based signature identification method to select prognostic biomarkers for glioma patients. We collected the clinical information and RNA-seq gene expression data from glioma patients in the CGGA database (http://www.cgga.org.cn/) and TCGA database (https://cancergenome.nih.gov/). The names of genes in the differential network were mapped to those of genes in homo sapiens, resulting in a list of 29 candidate genes. By matching both patient sample IDs and gene names from the clinical information and the gene expression data, a learning set of 310 samples from CGGA dataset and an independent validation set of 690 samples from TCGA dataset were created for prognostic signature identification, validation and further analysis. The details of the sample information were listed in Additional file 2: Table S2.
We used a multivariate COX proportional hazards (PH) model  and LASSO regression method  to select prognostic genes from the differential gene association network (see details in Text S1). A tenfold cross-validation was performed to select the optimal values of the tuning parameter for minimizing the mean cross-validation error.
The regression coefficients at the optimal tuning parameter were computed as risk coefficients, which were used to formularize a risk signature. The same risk signature was used to compute the risk scores (RSs) for patients in the independent validation set. The patients in each dataset were classified into a high-risk group and a low-risk group according to the optimal cutoff value using the ROC method. The Kaplan–Meier (K–M) curves for patients in the high- and low-risk groups were analyzed, and the statistical significance of the difference was assessed using the two-sided log-rank test. Time-dependent ROC analysis  was further conducted to evaluate the prognostic accuracy of the above RSs with respect to the 3- and 5-year survival predictions of patients in both the learning and independent validation sets.
Prediction of response to targeted therapeutics
The survival and gene expression data from glioma patients who received targeted therapies in the independent set were extracted to evaluate the predictive effectiveness of the macrophage-related gene signature for classifying patients into drug-sensitive and drug-resistant groups. We used the observed 3- or 5-year survival status to substitute for the latent drug-resistance status in these patients. Specifically, the 3- or 5-year survival status (alive or dead) was defined as the outcome (sensitive or resistant) of targeted drug treatment. The above risk signature was used to classify each patient into a sensitive group (i.e., low-risk group) or a resistant group (i.e., high-risk group) according to the optimal cutoff value of the RS using the ROC method. To further compare the powers of different gene signatures to predict responses to targeted therapy, we calculated the area under the curves (AUCs) of the time-dependent ROCs to assess their accuracies.
Comparison with other methods and other related signatures
We compared the macrophage-related gene signature produced from the above multicellular gene network perturbation method with other commonly used methods, including LASSO Cox regression model  and correlation network-based biomarker identification method without network perturbation. First, a LASSO Cox signature was trained from the full list of 1141 DEGs based on the CGGA set. In addition, a weighted correlation network (without perturbation) was constructed by calculating the pairwise Person correlation coefficients (PCCs) across all pairs of DEGs within both TCs and TAMs for the full set of Veh, Reb and Ep samples (30 samples). Based on the ranking of node strength score (defined as the sum of the absolute values of correlation coefficients of the node with other nodes), we selected the top ranked 12 DEGs as a gene set for defining a prognostic signature using the CGGA dataset.
Moreover, to assess whether the macrophage-related gene signature was independently correlated with the prognosis of glioma patients, we conducted univariate and multivariate COX regression analyses of clinicopathological factors and available gene signatures. Clinicopathological information, including age, gender and grade, was available for glioma patients in both the learning and validation sets. We also included the following gene signatures in the multivariate COX regression analyses: signature 1—the macrophage-related gene signature newly proposed in this study; signature 2—an EGFR gene signature studied by many groups [15, 16]; and an immune-related gene signature for predicting the prognosis of glioma patients, i.e., signature 3—the Cheng et al. signature .
The above risk factors were further extracted for construction of a combined signature using the LASSO COX model [19, 20]. As a result, we defined the combined signature as follows: CS = (0.008974621 × Age) + (1.617859481 × Grade) + (0.940077644 × Signature_1) + (0.006408624 × Signature_3). Here, Grade = 1 for lower grade glioma, and Grade = 2 for high grade glioma.
We tested the robustness of the macrophage-related gene signature obtained using the multicellular network perturbation analysis in comparison with LASSO Cox signature and correlation network-based signature using a bootstrapping approach. We generated 100 random datasets by randomly sampling 50% of the samples from the TCGA datasets (i.e., validation set). The AUC values of ROC with respect to overall survival, 3-year survival and 5-year survival were computed. Wilcoxon rank sum test (one-tailed) p values were computed to assess the significance of the difference between the probability distributions of AUC values of competing signature.
Construction and analysis of multicellular gene networks
The heatmaps of DEGs in TAMs (Additional file 1: Fig S1A) and TCs (Additional file 1: Fig S1B) across all samples demonstrated that the gene expression profile of the Reb samples was significantly different from that of Ep samples. The constructed sensitive network and sample-specific perturbation networks along with their topological attributes are shown in Additional file 1: Fig S2 and Fig S3–S6, respectively. A significant difference in the number of nodes and edges was observed between the sensitive and perturbation networks (Additional file 1: Fig S7). Other network topological attributes, including the network diameter, network centralization, and average number of neighbors, of the perturbation networks (Additional file 1: Fig S3–S6) were mainly higher than those of the sensitive network (Additional file 1: Fig S2). These results suggest that the perturbation networks had more complexity and that network rewiring might emerge during the acquisition and development of resistance to CSF1R inhibition.
A robust differential network (Fig. 2a) containing 42 nodes and 30 edges was thus constructed. The genes in the differential network are listed in Additional file 1: Table S3. Interestingly, all edges in the differential network were correlation-gained edges, indicating that the addition of Reb samples improved the correlation coefficients, which is in accordance with the above results (Additional file 1: Fig S2–S7).
Functional enrichment analysis of genes in the differential network revealed that some pathways, including the chemokine-mediated signaling and receptor-type tyrosine-protein phosphatase pathways as well as the cell-cell adhesion via plasma membrane adhesion molecule pathway, were significantly enriched during the acquired resistance to CSF1R inhibition treatment (Fig. 2b). This result suggests that cell–cell interactions between TCs and TAMs might contribute to the resistance of glioma to CSF1R inhibition, which was supported by the previous experimental results [10, 23,24,25].
Identification of the macrophage-related gene signature
We selected prognostic genes from the differential network using the LASSO regression method (see “Methods” section and Text S1). Figure 3a shows the selection of optimal tuning parameter (λ) of LASSO regression based on tenfold cross-validation of the learning set. Figure 3b shows the LASSO coefficient profiles of the 29 genes in the differential network. Each curve corresponds to evolution of the coefficient of each gene with respect to the change of the tuning parameters during the LASSO regression. Figure 3c lists the 12 genes selected under the optimal tuning parameter of LASSO. Five genes were located in macrophages (MФ) (i.e., ANPEP, DPP4, PRRG1, GPNMB, TMEM26), and 7 genes were located in TCs (i.e., PXDN, CDH6, SCN3A, SEMA6B, CCDC37, FANCA, NETO2). The coefficients of each gene in the COX PH model and the corresponding hazard ratios (HRs) were also listed and used to define a macrophage-related gene signature for predicting the prognosis of glioma patients. Accordingly, we formulated the following RS for each patient based on the expression levels of the selected genes: RS = (0.001695826 × ANPEP) + (0.001351164 × DPP4) + (0.828492221 × PRRG1) + (0.002693736 × GPNMB) + (0.572250065 × TMEM26) + (0.011112329 × PXDN) + (0.000861924 × CDH6) − (0.877296902 × SCN3A) − (0.042307865 × SEMA6B) + (0.019673956 × CCDC37) + (1.184362541 × FANCA) + (0.101032334 × NETO2). Furthermore, the univariate COX regression analysis (Additional file 1: Table S4) demonstrated that each signature gene was significantly associated with the survival of glioma patients.
Validation of prognostic significance and accuracy of the macrophage-related gene signature
We then investigated the prognostic significance and accuracy of the macrophage-related gene signature in both the learning and independent validation sets. Figure 4a shows the K–M curves for high-risk (blue) and low-risk (red) glioma patients in the learning set; significant differences were assessed with the log-rank test (p value less than 0.0001). The high-risk group of glioma patients tended to have a shorter survival time than the low-risk group. Figure 4b shows the prognostic accuracy of the risk signature evaluated by the AUCs of the time-dependent ROCs with respect to the 3- and 5-year survival rates of glioma patients in the learning set (AUC at 3 year: 0.92; AUC at 5 years: 0.891).
Moreover, we tested the prognostic significance of the macrophage-related gene signature using the independent validation set (Fig. 4c). The K–M curves demonstrated a significant difference between the survival rates of high- and low-risk patients (log-rank test p value less than 0.0001). Figure 4d shows the prognostic accuracy of the macrophage-related gene signature validated by the AUCs of the time-dependent ROCs with respect to 3- and 5-year survival rates of glioma patients in the validation set (AUC at 3 years: 0.818; AUC at 5 years: 0.836). These results demonstrated good prognostic value of the macrophage-related gene signature in glioma patients.
Furthermore, we tested the significance of prognostic accuracy of the macrophage-related gene signature evaluated using a bootstrapping approach. We randomly selected 12-gene sets from the whole transcriptome for 1000 times and compared their prognostic accuracy with that of the macrophage-related gene signature. CGGA dataset was used for training and TCGA dataset for validation. Figure 5a shows ROC curve of the macrophage-related gene signature (red) against ROC curves of bootstrapped 12-gene signatures (blue, only 100 curves were randomly shown). Figure 5b shows the probability distributions of the AUC values of 1000 sets of random 12-gene signatures. The p value (0.001) from the permutation test justified the statistical significance of the prognostic accuracy of the macrophage-related gene signature.
For comparison, we also generated gene signatures using LASSO Cox regression model and correlation network-based method (see Method section). Additional file 1: Fig S8 shows the prognostic accuracies of LASSO Cox signature and correlation network-based signature with respect to 3-year and 5-year survival prediction. Although these two signatures performed well on the learning set, their predictive accuracies on the test sets were less than the macrophage-related gene signature from the network perturbation analysis. We further employed a bootstrapping approach to test and compare the robustness of these gene signatures derived from different methods. Figure 6 shows the AUC values of ROCs of these signatures with respect to overall survival, 3-year survival and 5-year survival. The macrophage-related gene signature exhibited better accuracy and robustness than the LASSO Cox signature and the correlation-network-based signature.
The macrophage-related gene signature is predictive of the responses of gliomas to targeted therapeutics
Considering the important roles of macrophages in glioma progression and drug resistance [10,11,12, 26], a subset of TCGA patients who received targeted therapy were used to test the predictive power of the above macrophage-related gene signature for predicting drug sensitivity or resistance in glioma patients. The drug-resistance status of these patients was latent, we thus used the associated survival profiles as surrogate. The 3- or 5-year survival status (alive or dead) of each patient after the targeted therapy was used to evaluate the treatment outcome of the molecularly targeted therapeutics (sensitive or resistant). Each patient was predicted to be drug-sensitive (i.e., low-risk) or drug-resistant (i.e., high-risk) according to the optimal cutoff value of the RS using the time-dependent ROC method evaluated at 3-year or 5-year. Figure 7a, b shows the K–M survival curves of the predicted drug-sensitive patients (blue) and drug-resistant patients (red) as evaluated by 3-year survival ROC (Fig. 7a) or 5-year survival ROC (Fig. 7b), verifying distinct survival profiles between the two predicted groups of patients (log-rank test p values less than 0.0001).
To further assess the accuracy of different gene signatures for predicting the sensitivity or resistance to the targeted therapies, we compared 3 risk signatures that were designed for glioma patients: signature 1—the macrophage-related gene signature newly proposed in this study; signature 2—a conventional EGFR gene signature studied by many groups [15, 16]; signature 3—an immune-related gene signature, i.e., the Cheng et al. signature (FOXO3, IL6, IL10, ZBTB16, CCL18, AIMP1, FCGR2B, and MMP9) ; Signature 4—a gene signature based on IGF1/IGF1R-mediated pathways between macrophages and glioma cells (Quail et al. ), including several gene families in these pathways (PIK3R1, PIK3R2, PIK3R3, PIK3R4, PIK3R5, PIK3R6, PIK3AP1, AKT1, AKT2, AKT3, IGF1, IGF1R, IL4, IL4R, NFATC1, NFATC2, NFATC3, NFATC4, NFAT5, STAT6, CSF1 and CSF1R). We calculated the AUCs of the ROCs of these signatures to quantitatively assess and compare the accuracies of different gene signatures. Figure 7c shows the AUCs of the ROCs of these signatures for predicting drug sensitivity as evaluated by the 3-year survival outcomes (AUC of signature 1: 0.8295; AUC of signature 2: 0.6475; AUC of signature 3: 0.7155; AUC of signature 4: 0.5802). Figure 7d shows the AUCs of the ROCs of these signatures for predicting drug sensitivity as evaluated by the 5-year survival outcomes (AUC of signature 1: 0.7973; AUC of signature 2: 0.6462; AUC of signature 3: 0.6718; AUC of signature 4: 0.5877). Small p values computed using a bootstrap method  further demonstrated a superior predictive power of the macrophage-related gene signature compared with that of other signatures.
The macrophage-related gene signature is an independent prognostic signature
The univariate and multivariate COX regression analyses for both the learning and validation datasets revealed that the macrophage-related gene signature retained prognostic significance for glioma patients when adjusted for both clinicopathologic risk factors (age, gender, grade) and other existing gene signatures (EGFR gene signature, Cheng et al. signature and IGF1/IGF1R pathways signature) (Table 1). These results indicated that the macrophage-related gene signature was independently correlated with the overall and 5-year survival of glioma patients.
To further explore the prognostic value of the macrophage-related gene signature in stratified cohorts, patients were first classified by 2 important clinicopathological factors, age and grade, that significantly correlated with the prognosis of glioma patients (Table 1). Figure 8a–d shows the prognostic significance of the macrophage-related gene signature in different glioma cohorts stratified by age (Fig. 8a, b) or grade (Fig. 8c, d). In all of these cohorts, patients were classified as high- versus low-risk groups using cut point from ROCs, and the high-risk patients had a significantly shorter overall survival than low-risk patients. Subsequently, patients who received pharmaceutical therapy and radiotherapy were utilized to validate the prognostic significance of the macrophage-related gene signature. Figure 8e–h demonstrates that the macrophage-related gene signature retained prognostic significance for glioma patients treated with or without pharmaceutical therapy (Fig. 8e, f) and radiotherapy (Fig. 8g, h). These results indicated that the macrophage-related gene signature could accurately identify patients with an unfavorable prognosis regardless of their clinicopathological and treatment characteristics.
We further examined whether combining the macrophage-related gene signature with clinicopathological risk factors and other existing gene signatures could significantly improve the prognostic accuracy of the risk signature. The time-dependent ROC curves (Fig. 9) compared the prognostic accuracy by age, grade, signature 1 (i.e., macrophage-related gene signature), signature 3 (i.e., Cheng et al. signature) and the combined signature (see the definition in the Methods section). The AUCs of ROC curves in both the learning dataset (Fig. 9a, b) and the validation dataset (Fig. 9c, d) showed that the combined signature outperformed other risk signatures except the macrophage-related gene signature for predicting the 3- and 5-year survival rates of glioma patients. Noticeably, the AUC of the ROC of the macrophage-related gene signature was rather close to that of the combined signature in the learning dataset (Fig. 9a, b), or even higher in the validation dataset (Fig. 9c, d). These results indicated that the macrophage-related gene signature is an independent prognostic signature and possessed convincingly strong and robust prognostic power in comparison to the clinicopathological factors and other existing risk signatures.
Several studies have reported that the tumor microenvironment plays important roles in glioma progression and therapeutic response [10,11,12, 28,29,30,31]. TAMs, a type of immune cell, account for the majority of nonneoplastic cells in the glioma microenvironment . However, limited studies have investigated macrophage-based molecular biomarkers for glioma prognostication and therapeutic response prediction. In this study, we first identified and validated a macrophage-related gene signature that has a strong prognostic significance and accuracy for glioma patients.
The conventional node biomarker method identifies biomarkers using a set of single genes or molecules and does not consider interactions between these genes [32,33,34,35,36,37,38]. In recent years, a network biomarker method has been proposed to identify biomarkers of cancers or drug resistance using a network-based analysis approach, which takes into account gene interactions and uses aggregates of genes/proteins in networks for prediction [39,40,41]. However, these previous network-based approaches mainly focused on intracellular gene networks within TCs and did not explicitly or systematically consider intercellular interactions, particularly tumor–microenvironment interactions. In this study, we built multicellular gene networks connecting TAMs and TCs to identify gene signatures for predicting prognosis and drug resistance in glioma patients, which provided systems biology insights into tumor cell–microenvironment interactions in cancer progression and therapeutic resistance. Furthermore, a robust differential gene network was constructed using the network perturbation analysis method, which is suitable for small sample sizes. The resulting differential gene pairs (Fig. 2a) not only reflected the intracellular gene coexpression patterns within TCs and TAMs but also revealed the intercellularly correlated genes between these cell types during the acquisition of drug resistance.
We compared the macrophage-related gene signature with other existing signatures, and our signature showed better accuracy for predicting the survival outcomes of glioma patients who received molecularly targeted therapies (Fig. 7). We speculate that the superiority of our signature lies in the fact that it considered the interaction between TCs and the microenvironment during tumor progression. Although we considered only TAMs, which constitute a fraction of microenvironmental cell types, the identified TAM-TC gene signature outperformed the existing signatures that focused on TCs or immune genes only, which is an intriguing finding for investigating the prognostic significance of microenvironment-mediated gene signatures. The multicellular gene network approach developed herein provides a promising strategy to systematically integrate genes in both TCs and microenvironmental cells and thus identify more robust and accurate signatures or biomarkers.
In addition, we assessed the association between the macrophage-related gene signature and frequent genetic and genomic alterations in gliomas by exploring other genomic characteristics that are available within the TCGA or CGGA databases. Figure 10 shows the distribution of the risk scores evaluated by macrophage-related gene signature in patients from the TCGA set stratified by IDH1 mutation status (Fig. 10a) or methylator CIMP status (Fig. 10b), and in patients stratified by EGFR mutation status (Fig. 10c) or PTEN mutation status (Fig. 10d) with the information available from the CGGA set. The signature value was different between cases stratified by IDH1 mutation status, CIMP status, EGFR mutation status and PTEN status. We also examined the distribution of the risk scores in patients stratified by molecular subtypes with the information available from the TCGA set (Fig. 10e). Significant differences of risk values were observed between several pairs of different subtypes (classical vs. neural, classical vs. proneural, mesenchymal vs. neural, and mesenchymal vs. proneural). These results suggested the associations of the macrophage-related gene signature with genetic and genomic alterations in gliomas and molecular subtypes. In future studies, we will develop an ensemble model that integrates multiple gene signatures designed for different subgroups of glioma patients stratified by molecular subtypes or major genomic mutations in IDH1, CIMP, EGFR and PTEN in glioma cohorts to further improve the predictive accuracy.
The established macrophage-related gene signature includes both protective genes (SCN3A, SEMA6B) and risk-increasing genes (PXDN, CDH6, ANPEP, CCDC37, DPP4, GPNMB, FANCA, NETO2, PRRG1, TMEM26) for predicting the prognosis of glioma patients (Fig. 3). This signature could be regarded as a macrophage-related protective and risky pattern of TC–TAM interactions in gliomas, which is consistent with the balance between the antitumorigenic M1 phenotype and protumorigenic M2 phenotype of TAMs in gliomas.
Additional file 1: Table S5 lists the experimental evidences for functional roles of the five macrophage-related genes in cancer progression and/or drug resistance. For instance, ANPEP was reportedly involved in cell migration and tumor metastasis [42,43,44,45,46,47,48]. ANPEP and TMEM26 were associated with drug response [49, 50]. Remarkably, the important role of DPP4 in anti-tumor immunity  has been revealed in many cancers including glioblastomas . The inhibition of DPP4 was found to improve both naturally occurring tumor immunity and immunotherapy by enhancing lymphocyte trafficking . Our results may therefore implicate new macrophage-targeting treatment strategies to improve clinical outcomes. The genes in the macrophage-related gene signature could be employed as stand-alone targets or in combination with the existing targeted therapies, attributing to their prognostic significance and association with drug resistance. For example, GPNMB has been found to be highly upregulated in human glioma-associated microglia/macrophages, which are the predominant source of GPNMB transcripts . High GPNMB expression was found to be associated with poor prognosis in human glioblastoma. Furthermore, GPNMB has also been indicated as a potential molecular therapeutic target in patients with glioblastoma .
The advantages of our study originate from the use of a multicellular network-based gene screening approach, large population databases for learning and validation and robust risk signature identification methods. In future studies, we will investigate the functions and mechanisms of the 12 genes alone and in combination to verify their clinical applicability. The capability of the signature identified herein for predicting drug resistance should be further validated by prospective studies in the future. We will also leverage single cell RNA-seq data to construct multilayer networks that connect tumor microenvironmental cells to tumor cells , which was anticipated to improve the biological interpretation and predictive accuracy of the biomarkers for predicting prognosis and therapeutic response.
In summary, we developed a multicellular gene network approach to investigate the role of TC–TAM interactions in the progression and therapeutic responses of gliomas. The identified macrophage-related gene signature showed good accuracy for predicting the prognosis and targeted therapeutic responses of glioma patients. The multicellular gene network-based signature provided mechanistic insights into microenvironment-mediated drug resistance and implied that combining current targeted therapies with macrophage-targeted therapy might improve the long-term treatment outcomes of glioma patients.
Availability of data and materials
The datasets including the RNA-seq data in mice (GSE69104), clinical information and the gene expression data of glioma patients analyzed during the current study are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69104), TCGA database (https://cancergenome.nih.gov/) and CGGA database (http://www.cgga.org.cn/). The source R scripts for analysis in this study were available at https://github.com/dongbusun/MCGN-MS.
tumor-associated macrophages and microglia
differentially expressed gene
Pearson correlation coefficient
area under curve
receiver operating characteristic curve
Yiu G, He Z. Glial inhibition of CNS axon regeneration. Nat Rev Neurosci. 2006;7(8):617–27.
Brown C. Targeted therapy: an elusive cancer target. Nature. 2016;537(7620):S106.
Robert B, Edward C, Luca M, Wilhelm-Benartzi CS, Jane B. Poised epigenetic states and acquired drug resistance in cancer. Nat Rev Cancer. 2014;14(11):747–53.
Lee H-J, Zhuang G, Cao Y, Du P, Kim H-J, Settleman J. Drug resistance via feedback activation of Stat3 in oncogene-addicted cancer cells. Cancer Cell. 2014;26(2):207–21.
Wagle N, Van Allen EM, Treacy DJ, Frederick DT, Cooper ZA, Taylor-Weiner A, Rosenberg M, Goetz EM, Sullivan RJ, Farlow DN. MAP kinase pathway alterations in BRAF-mutant melanoma patients with acquired resistance to combined RAF/MEK inhibition. Cancer Discov. 2014;4(1):61–8.
Pazarentzos E, Bivona T. Adaptive stress signaling in targeted cancer therapy resistance. Oncogene. 2015;34:5599.
Dean M, Fojo T, Bates S. Tumour stem cells and drug resistance. Nat Rev Cancer. 2005;5(2):275–84.
Cojoc M, Mäbert K, Muders MH, Dubrovska A: A role for cancer stem cells in therapy resistance: cellular and molecular mechanisms. In: Seminars in cancer biology: 2015, Elsevier; 2015: 16–27
Pasipanodya JG, Srivastava S, Gumbo T. Meta-analysis of clinical studies supports the pharmacokinetic variability hypothesis for acquired drug resistance and failure of antituberculosis therapy. Clin Infect Dis Off Publ Infect Dis Soc Am. 2012;55(2):169–77.
Quail DF, Bowman RL, Akkari L, Quick ML, Schuhmacher AJ, Huse JT, Holland EC, Sutton JC, Joyce JA. The tumor microenvironment underlies acquired resistance to CSF-1R inhibition in gliomas. Science. 2016;352(6288):aad3018.
Junttila MR, de Sauvage FJ. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature. 2013;501(7467):346–54.
Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med. 2013;19(11):1423–37.
Sun X, Hu B. Mathematical modeling and computational prediction of cancer drug resistance. Brief Bioinform. 2017;19:1382–99.
Wang E, Zou J, Zaman N, Beitel LK, Trifiro M, Paliouras M. Cancer systems biology in the genome sequencing era: part 2, evolutionary dynamics of tumor clonal networks and drug resistance. Semin Cancer Biol. 2013;23(4):286–92.
Etienne MC, Formento JL, Lebrunfrenay C, Gioanni J, Chatel M, Paquis P, Bernard C, Courdi A, Bensadoun RJ, Pignol JP. Epidermal growth factor receptor and labeling index are independent prognostic factors in glial tumor outcome. Clin Cancer Res. 1998;4(10):2383.
Li J, Liang R, Song C, Xiang Y, Liu Y. Prognostic significance of epidermal growth factor receptor expression in glioma patients. Oncotargets Therapy. 2018;11:731–42.
Cheng W, Ren X, Zhang C, Cai J, Liu Y, Han S, Wu A. Bioinformatic profiling identifies an immune-related risk signature for glioblastoma. Neurology. 2016;86(24):2226–34.
Liu X, Wang Y, Ji H, Aihara K, Chen L. Personalized characterization of diseases using sample-specific networks. Nucleic Acids Res. 2016;44(22):e164.
Cox DR, Cox DR, Oakes D. Analysis of survival data. New York: Chapman and Hall; 1984.
Tibshirani R. The lasso method for variable selection in the cox model. Stat Med. 1997;16(4):385–95.
Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2015;56(2):337–44.
Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J Stat Softw. 2011;39(05):1–13.
Merzak A, Koocheckpour S, Pilkington GJ. CD44 mediates human glioma cell adhesion and invasion in vitro. Can Res. 1994;54(15):3988.
Navis AC, Eijnden MVD, Schepens JTG, Huijsduijnen RHV, Wesseling P, Hendriks WJAJ. Protein tyrosine phosphatases in glioma biology. Acta Neuropathol. 2010;119(2):157.
Ohnishi T, Izumoto S, Arita N, Hiraga S, Taki T, Hayakawa T. Expression and biological functions of L1 cell adhesion molecule in malignant glioma cells. Berlin: Springer; 1996.
Hambardzumyan D, Gutmann DH, Kettenmann H. The role of microglia and macrophages in glioma maintenance and progression. Nat Neurosci. 2015;19(1):20.
Xavier R, Natacha T, Alexandre H, Natalia T, Frédérique L, Jean-Charles S, Markus M. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. 2011;12(1):77.
Zheng Y, Bao J, Zhao Q, Zhou T, Sun X. A spatio-temporal model of macrophage-mediated drug resistance in glioma immunotherapy. Mol Cancer Ther. 2018;17(4):814–24.
Liang WZ, Yongjiang; Zhang, Ji; Sun, Xiaoqiang: Multiscale modeling reveals angiogenesis-induced drug resistance in brain tumors and predicts a synergistic drug combination targeting EGFR and VEGFR pathways. BMC Bioinform 2019, https://doi.org/10.1186/s12859-019-2737-1
Sun X. Multi-scale agent-based brain cancer modeling and prediction of TKI treatment response: incorporating EGFR signaling pathway and angiogenesis. BMC Bioinform. 2012;13(1):218.
Sun X, Bao J, Shao Y. Mathematical modeling of therapy-induced cancer drug resistance: connecting cancer mechanisms to population survival rates. Sci Rep. 2016;6:22498.
Aksoy BA, Demir E, Babur Ö, Wang W, Jing X, Schultz N, Sander C. Prediction of individualized therapeutic vulnerabilities in cancer from genomic profiles. Bioinformatics. 2014;30(14):2051–9.
Alexandrov LB, Nikzainal S, Wedge DC, Aparicio SAJR, Behjati S, Biankin AV, Bignell GR, Bolli N, Borg A, Børresendale AL. Signatures of mutational processes in human cancer. Nature. 2013;500(7463):415–21.
Shukla S, Patric IRP, Thinagararjan S, Srinivasan S, Mondal B, Hegde AS, Chandramouli BA, Santosh V, Arivazhagan A, Somasundaram K. A DNA methylation prognostic signature of glioblastoma: identification of NPTX2-PTEN-NF-κB nexus. Can Res. 2013;73(22):6563–73.
Chan E, Prado DE, Weidhaas JB. Cancer microRNAs: from subtype profiling to predictors of response to therapy. Trends Mol Med. 2011;17(5):235.
Curtis C. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature. 2012;486(7403):346–52.
Grasso CS, Wu YM, Robinson DR, Cao X, Dhanasekaran SM, Khan AP, Quist MJ, Jing X, Lonigro RJ, Brenner JC. The mutational landscape of lethal castration-resistant prostate cancer. Nature. 2012;487(7406):239–43.
Zhou M, Liu Z, Zhao Y, Ding Y, Liu H, Xi Y, Xiong W, Li G, Lu J, Fodstad O. MicroRNA-125b confers the resistance of breast cancer cells to paclitaxel through suppression of pro-apoptotic Bcl-2 antagonist killer 1 (Bak1) expression. J Biol Chem. 2010;285(28):21496.
Gao S, Tibiche C, Zou J, Zaman N, Trifiro M, O’Connormccourt M, Wang E. Identification and construction of combinatory cancer hallmark-based gene signature sets to predict recurrence and chemotherapy benefit in stage II colorectal cancer. Jama Oncol. 2015;2(1):1–9.
Li J, Lenferink AE, Deng Y, Collins C, Cui Q, Purisima EO, O’Connor-McCourt MD, Wang E. Identification of high-quality cancer prognostic markers and metastasis network modules. Nat Commun. 2010;1:34.
McGee SR, Tibiche C, Trifiro M, Wang E. Network analysis reveals a signaling regulatory loop in PIK3CA -mutated breast predicting survival outcome. Genom Proteom Bioinform. 2017;15(2):121–9.
Hee DS, Kwang-Pyo L, Dongjun J, Chang-Jin K, Kyung-Sook C, Young KJ, Bum-Chan P, Sup PS, Seon-Young K, Ki-Sun K. GPR171 expression enhances proliferation and metastasis of lung cancer cells. Oncotarget. 2016;7(7):7856–65.
Wang T, Han S, Wu Z, Han Z, Yan W, Liu T, Wei H, Song D, Zhou W, Yang X. XCR43 promotes cell growth and migration and is correlated with bone metastasis in non-small cell lung cancer. Biochem Biophys Res Commun. 2015;464(2):635–41.
Kim M, Rooper L, Xie J, Rayahin J, Burdette JE, Kajdacsy-Balla AA, Barbolina MV. The lymphotactin receptor is expressed in epithelial ovarian carcinoma and contributes to cell migration and proliferation. Mol Cancer Res. 2012;10(11):1419.
Eleonora D, Roberto R, Liliana GR, Barbu EM, Hitomi H, John LS, St. Molldrem JJ, Angelo C, Sidman RL, Wadih A. CD13-positive bone marrow-derived myeloid cells promote angiogenesis, tumor growth, and metastasis. Proc Natl Acad Sci USA. 2013;110(51):20717–22.
Rabindranath B, Chih-Yung C, Ming-Chin Y, Jei-Ming P, Chung-Ru H, Chih-Yun H, Hsiao-Ling H, Ho UY, Shi-Ming L, Yu-Jr L. Functional genomics identified a novel protein tyrosine phosphatase receptor type F-mediated growth inhibition in hepatocarcinogenesis. Hepatology. 2014;59(6):2238–50.
Du WW, Ling F, Minhui L, Xiangling Y, Yaoyun L, Chun P, Wei Q, O’Malley YQ, Askeland RW, Sugg SL. MicroRNA miR-24 enhances tumor invasion and metastasis by targeting PTPN9 and PTPRF to promote EGF signaling. J Cell Sci. 2013;126(6):1440–53.
Scrima M, Marco CD, Vita FD, Fabiani F, Franco R, Pirozzi G, Rocco G, Malanga D, Viglietto G. The nonreceptor-type tyrosine phosphatase PTPN13 is a tumor suppressor gene in non-small cell lung cancer. Am J Pathol. 2012;180(3):1202–14.
Azimi A, Tuominen R, Costa Svedman F, Caramuta S, Pernemalm M, Frostvik Stolt M, Kanter L, Kharaziha P, Lehtiö J, Hertzman Johansson C, et al. Silencing FLI or targeting CD13/ANPEP lead to dephosphorylation of EPHA2, a mediator of BRAF inhibitor resistance, and induce growth arrest or apoptosis in melanoma cells. Cell Death Dis. 2017;8:3029.
Nass N, Dittmer A, Hellwig V, Lange T, Mirjam BJ, Leyh B, Ignatov A, Weiβenborn C, Kirkegaard T, Lykkesfeldt AE. Expression of transmembrane protein 26 (TMEM26) in breast cancer and its association with drug response. Oncotarget. 2016;7(25):38408–26.
Jang JH, Baerts L, Waumans Y, Meester ID, Yamada Y, Limani P, Gil-Bazo I, Weder W, Jungraithmayr W. Suppression of lung metastases by the CD26/DPP4 inhibitor Vildagliptin in mice. Clin Exp Metastasis. 2015;32(7):677–87.
Matrasova I, Busek P, Balaziova E, Sedo A. Heterogeneity of molecular forms of dipeptidyl peptidase-IV and fibroblast activation protein in human glioblastomas. Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub. 2017;161(3):252–60.
Barreira da Silva R, Laird ME, Yatim N, Fiette L, Ingersoll MA, Albert ML. Dipeptidylpeptidase 4 inhibition enhances lymphocyte trafficking, improving both naturally occurring tumor immunity and immunotherapy. Nat Immunol. 2015;16:850.
Szulzewsky F, Pelz A, Feng X, Synowitz M, Markovic D, Langmann T, Holtman IR, Wang X, Eggen BJ, Boddeke HW. Glioma-associated microglia/macrophages display an expression profile different from M1 and M2 polarization and highly express Gpnmb and Spp1. PLoS ONE. 2015;10(2):e0116644.
Kuan CT, Wakiya K, Dowell JM, Herndon NJ, Reardon DA, Graner MW, Riggins GJ, Wikstrand CJ, Bigner DD. Glycoprotein nonmetastatic melanoma protein B, a potential molecular therapeutic target in patients with glioblastoma multiforme. Clin Cancer Res Off J Am Assoc Cancer Res. 2006;12(7 Pt 1):1970.
Zhang J, Guan M, Wang Q, Zhang J, Zhou T, Sun X. Single-cell transcriptome-based multilayer network biomarker for predicting prognosis and therapeutic response of gliomas. Brief Bioinform. 2019. https://doi.org/10.1093/bib/bbz040.
We would like to acknowledge Prof. Xun Lan at School of Medicine, Tsinghua University, Prof. Jun Chen at Zhongshan School of Medicine, Sun Yat-sen University, and Dr. Lei Huang at Houston Methodist Hospital Research Institute for valuable discussions.
X.S. was funded by the National Natural Science Foundation of China (11871070, 61503419), the Guangdong Nature Science Foundation (2016A030313234, 2014A030310355) and the Opening Project of Guangdong Province Key Laboratory of Computational Science at the Sun Yat-Sen University (2018003). Y.S. was funded by the NIH/NCI Grant 5P30CA016087. X.D.Z was supported by research Grants SRG2016-00083-FHS, FHS-CRDA-029-002-2017 and MYRG2018-00071-FHS at the University of Macau. The funding body played no role in the design of the study; collection, analysis, and interpretation of data; or writing of the manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Detailed description of the methods. Table S1. Scheme of sample grouping. Table S3. List of genes in the differential network. Table S4. Univariate Cox regression analysis for the association of each signature gene with the overall survival of glioma patients. Table S5. Experimental evidences for functional roles of the eight macrophage-related genes in cancer progression and/or drug resistance. Figure S1. Differential gene expression profiles of TAMs and TCs between the Endpoint and Rebound groups. Figure S2. TC-TAM sensitive network and its topological attributes. Figure S3. EpUReb1 sample-specific TC-TAM network and its topological attributes. Figure S4. EpUReb2 sample-specific TC-TAM network and its topological attributes. Figure S5. EpUReb3 sample-specific TC-TAM network and its topological attributes. Figure S6. EpUReb4 sample-specific TC-TAM network and its topological attributes. Figure S7. Network topological analysis of the sensitive network and 4 sample-specific resistant networks. Figure S8. Prognostic accuracies of LASSO Cox signature and correlation network-based signature. Supplementary references.
Supplementary tables for the details of patient sample information from CGGA and TCGA.
About this article
Cite this article
Sun, X., Liu, X., Xia, M. et al. Multicellular gene network analysis identifies a macrophage-related gene signature predictive of therapeutic response and prognosis of gliomas. J Transl Med 17, 159 (2019). https://doi.org/10.1186/s12967-019-1908-1