An Immune-Related Gene Pairs Signature for Predicting Survival in Glioblastoma

Background: Glioblastoma (GBM) is the frequently occurring and most aggressive form of brain tumors. In the study, we constructed an immune-related gene pairs (IRGPs) signature to predict overall survival (OS) in patients with GBM.

Methods: We established IRGPs with immune-related gene (IRG) matrix from The Cancer Genome Atlas (TCGA) database (Training cohort). After screened by the univariate regression analysis and least absolute shrinkage and selection operator (LASSO) regression analysis, IRGPs were subjected to the multivariable Cox regression to develop an IRGP signature. Then, the predicting accuracy of the signature was assessed with the area under the receiver operating characteristic curve (AUC) and validated the result using the Chinese Glioma Genome Atlas (CGGA) database (Validation cohorts 1 and 2).

Results: A 10-IRGP signature was established for predicting the OS of patients with GBM. The AUC for predicting 1-, 3-, and 5-year OS in Training cohort was 0.801, 0.901, and 0.964, respectively, in line with the AUC of Validation cohorts 1 and 2 [Validation cohort 1 (1 year: 0.763; 3 years: 0.786; and 5 years: 0.884); Validation cohort 2 (1 year: 0.745; 3 years: 0.989; and 5 years: 0.987)]. Moreover, survival analysis in three cohorts suggested that patients with low-risk GBM had better clinical outcomes than patients with high-risk GBM. The univariate and multivariable Cox regression demonstrated that the IRGPs signature was an independent prognostic factor.

Conclusions: We developed a novel IRGPs signature for predicting OS in patients with GBM.

Introduction

Glioma is the most common and the most aggressive cancer of the central nervous system (1). Grade IV glioma: glioblastoma (GBM) is the major subtype of glioma, accounting for more than 48.6% of glioma cases in the USA (2). Although advanced therapies, like surgery, chemotherapy, targeted therapy, and so on, had been applied for GBM, the 5-year relative survival probability of GBM is only limited to 6.8% (3), and the median lifespan is only 19.2 months (4). Due to these discouraging survival data, the development of more accurate tumor-specific biomarkers of GBM for further constructing novel diagnostic signature and guiding clinical treatment is still a priority in GBM management.

There is emerging evidence showing that the immune system plays an important role in the tumorigenesis, development, and metastasis of cancers (5). Immune cells could identify and eradicate tumor cells by recognizing tumor-associated antigen, a key factor for immune surveillance (6). However, tumor cells could act on the immune system to attenuate the effects of immune surveillance and to promote self-development (7). Immunotherapy, a novel treatment method for tumors based on this theory, has been used for the therapy of many cancers (8). In recent years, two-phase III trials of immunotherapy (NCT02017717 and NCT02617589) have been initiated and raised new hopes for the GBM therapy. Nevertheless, owing to the blood–brain barrier, the central nervous system becomes an immunologically privileged situation, and immune escape challenges the effectiveness of immunotherapy in GBM (9, 10). Therefore, the combination management to achieve more efficacious therapeutic benefits is required for complex immune evasion strategies of GBM.

Abnormal gene expression is a common phenomenon in malignant tumors, which may promote the progression of tumors (11). Recently, omics technology and bioinformatics analysis have created an opportunity to identify novel molecular biomarkers and understand potential mechanisms in many tumors. For example, Zhao et al. (12) identified a six-gene risk signature for the outcome prediction of GBM. Several studies have reported immune-related gene pairs (IRGPs) signature in ovarian cancer (13), liver cancer (14), and colorectal cancer (15) and showed well-predictive accuracy. However, up to now, there is no study reporting an IRGP signature in GBM.

In the last two decades, genetic analyses identified point mutation of the gene: isocitrate dehydrogenase (IDH) was an early and common event in gliomagenesis, especially in Grade II and Grade III gliomas (16, 17). In recent years, GBM was defined as a diffuse astrocytic glioma without IDH-mutation (18, 19). Moreover, IDH-mutation in patients with glioma predicts better overall survival (OS) and influences the infiltration of immune cells and immunity (20–22), which may bias the analysis in the study. Hence, in the present study, we did not enroll the patients with GBM with IDH-mutation. Herein, we collected the gene expression matrix from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov) and immune-related gene (IRG) set from the ImmPort database (https://www.immport.org/home) to establish IRGPs and an IRGP signature in GBM with IDH-wildtype. Meanwhile, the signature was validated with data from the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn). Then, we explored the relationship between the signature and clinical characteristics. Finally, we proved the predictive effectiveness and accuracy of this signature by comparing the area under the curve (AUC) of this signature with that of clinical indexes.

Materials and Methods

Data Source and Pre-processing

The gene expression profile and relevant clinical characteristics were downloaded from the databases of TCGA and CGGA. Patients with follow-up time <1 mouth or lacking complete OS would be removed. Three independent cohorts including 514 cases were gathered in the present study [160 cases in Training cohort (TCGA-GBM); 260 cases in Validation cohort 1 (CGGA-693); and 94 cases in Validation cohort 2 (CGGA-301)]. In addition, 472 IRGs were collected from the ImmPort database (https://www.immport.org/home). All raw data from CGGA were transformed with log₂.

Establishment of an IRGP Signature

As described in the previous study (23), we performed the pairwise comparison in all samples. Score 1 represented that the expression level of first IRG was higher than that of second IRG in a specific IRGP, and score 0 represented the opposite. Additionally, an IRGP would be deserted if the score of which was 0 or 1 in over 80% of samples. Then, in the Training cohort, the univariate regression analysis with p-value <0.001 as the threshold was used to preliminary filtrate IRGPs. Later, the least absolute shrinkage and selection operator (LASSO) regression analysis with iteration = 1,000 and 10-fold cross-validations to prevent overfitting were used to screen up the remaining IRGPs. Ultimately, with p-value <0.05 as the cut-off criterion, we applied multivariate Cox regression analysis to determine top OS-related IRGPs, which would be used to develop an IRGP signature and immune risk-score formula for calculating risk score. Based on the risk cut-off score: the median value, patients were classified into two subgroups (high- and low-risk groups).

Evaluation and Validation of the IRGPs Signature

In the same way, the risk score of each case in Validation cohort 1 and Validation cohort 2 was also calculated. The receiver operating characteristic curve (ROC) was carried out to evaluate the predictive capacity of this IRGPs signature by calculating the AUC.

Difference of Tumor-Infiltrating Immune Cells Between Two Subgroups

As described previously (24), the CIBERSORT algorithm was utilized to assess the landscape of 22 tumor-infiltrating immune cells (TIICs) in the tumor microenvironment. Cases with p < 0.05 were selected for further research. The distribution of 22 TIICs between subgroups was also pictured.

Gene Set Enrichment Analysis

In the Training cohort, the gene set enrichment analysis (GSEA) was used to assess the potential biological pathways related to this IRGPs signature. False discovery rate (FDR) <0.01 was chosen as the cut-off criterion.

Statistical Analysis

Table 1 presents all clinical and pathological information of patients with GBM in three cohorts. For categorical data, the differences among different groups were compared with the chi-square test, whereas, for measurement data, the differences were compared with the t-test. The univariate and multivariate Cox regression analyses were used to determine independent prognosis-related factors. Survival difference was compared with the logrank test and shown with the Kaplan–Meier method. All statistical analyses were performed using R 3.6.3 (https://www.r-project.org).

TABLE 1

Table 1. The baseline characteristics of patients with glioblastoma (GBM) in this study.

Results

Establishment of the IRGPs Signature

A total of 19,916 IRGPs with score of 0 or 1 in <80% were paired. With the univariate regression analysis, 48 IRGPs were filtered in the Training cohort. Then, the LASSO regression analysis screened out 24 IRGPs for further analysis (Figure 1A). Finally, 10 IRGPs were identified and used to develop an IRGP signature (Table 2 and Figure 1B). Risk score = (−0.4987 × Score_BMP2|IL13RA2) + (−0.4115 × Score_BMP2|NOV) + (−0.3985 × Score_BMP2|OXTR) + (0.8071 × Score_CHGB|FGFR2) + (0.6000 × Score_CHGB|GDF11) + (0.4410 × Score_FCGR2B|PRKCB) + (−0.3854 × Score_{FGF12|SLC11A1}) + (−0.7438 × Score_IL10RA|OSMR) + (0.6128 × Score_LEFTY2|MSTN) + (0.8191 × Score_{SLC11A1|TRIM22}). Based on the cut-off risk score (1.089), patients were divided into two subgroups (high- and low-risk groups). The distribution of risk score and survival status of the IRGPs signature in Training cohort were shown in Supplementary Figure 1A. A nomogram combining risk score and clinical characteristics was also constructed for the prediction of 1-, 3-, and 5-year survival probability (Figure 1C).

FIGURE 1

Figure 1. Construction of immune-related gene pairs (IRGPs) signature. (A) “Leave-one-out-cross-validation” for the parameter selection in the least absolute shrinkage and selection operator (LASSO) regression. (B) The forest map of multivariate Cox regression analysis. (C) The prognostic nomogram for the prediction of 1-, 3- and 5-year overall survival (OS) in glioblastoma (GBM). IRGPs, immune-related gene pairs. *p < 0.05, **p < 0.01, ***p < 0.001.

TABLE 2

Table 2. Information on the 10 immune-related gene pairs (IRGPs).

Evaluation and Validation of the IRGPs Signature

We calculated the AUC of this IRGPs signature in the Training cohort to assess the predictive accuracy and effectiveness of the IRGPs signature. The AUC of IRGPs signature in the Training cohort for 1-, 3-, and 5-year OS was 0.801, 0.901, and 0.964, respectively (Figure 2A), which were markedly higher than the AUCs of clinical indexes (Supplementary Figures 2A–C). Also, we calculated the immune risk score of each patient in the Validation cohort 1 and Validation cohort 2 and divided them into high- and low-risk groups according to the cut-off risk score. The distribution of risk score and survival status of the IRGPs signature in the Validation cohort 1 and Validation cohort 2 were shown in Supplementary Figures 1B,C. Then, we validated the predictive accuracy of the IRGPs signature in the Validation cohort 1 and Validation cohort 2. The AUC in Validation cohort 1 reached 0.763 at 1-year, 0.786 at 3-year, and 0.884 at 5-year OS (Figure 2B), and in Validation cohort 2 reached 0.745 at 1-year, 0.989 at 3-year, and 0.987 at 5-year OS (Figure 2C). Similar to the results in the Training cohort, the AUCs of the IRGPs signature at 1-, 3-, and 5-year OS in Validation cohorts 1 and 2 were obviously superior to the AUCs of clinical indexes (Supplementary Figures 2D–I).

FIGURE 2

Figure 2. Evaluation of immune-related gene pairs (IRGPs) signature. (A) The area under the receiver operating characteristic (ROC) curve (AUC) for 1-, 3-, and 5-year overall survival (OS) of patients with glioblastoma (GBM) in the Training cohort. (B) In Validation cohort 1. (C) In Validation cohort 2.

Here, we constructed a signature with 10 IRGPs, which included 16 IRGs (BMP2, CHGB, FCGR2B, FGF12, IL10RA, LEFTY2, SLC11A1, IL13RA2, NOV, OXTR, GDF11, PRKCB, FGFR2, OSMR, MSTN, and TRIM22). Among 16 genes, the expressions of 9 genes (CHGB, FCGR2B, IL10RA, SLC11A1, IL13RA2, NOV, OXTR, OSMR, and LEFTY2) in the high-risk group was significantly upregulated, and the expressions of two genes (FGFR2 and MSTN) were opposite (Supplementary Figure 3).

The IRGPs Signature Is an Important Prognosis-Related Factor

The Kaplan–Meier plot in the Training cohort indicated patients with GBM with high-risk scores exhibited worse clinical outcomes than patients with low-risk scores (Figure 3A). Similar results were observed in Validation cohorts 1 and 2 (Figures 3B,C). Moreover, the stratified analysis showed the prognosis of patients with high-risk scores was poorer than that of patients with low-risk scores (Figure 4).

FIGURE 3

Figure 3. Survival difference between high- and low-risk groups: (A) In Training cohort. (B) In Validation cohort 1. (C) In Validation cohort 2.

FIGURE 4

Figure 4. Survival difference between high- and low-risk patients in different subgroup. (A) Age <50. (B) Age >50. (C) Female. (D) Male. (E) Primary tumor. (F) Recurrent tumor. (G) KPS ≤70. (H) KPS >70. (I) Radiotherapy: No. (J) Radiotherapy: Yes. (K) Chemotherapy: No. (L) Chemotherapy: Yes. (M) Mesenchymal. (N) Non-mesenchymal. KPS, Karnofsky.

For comparing the IRGPs signature with clinical characteristics, we assessed the importance of the IRGPs signature for prognosis of patients with GBM. The univariable Cox regression analysis showed the IRGPs signature was a meaningful factor influencing the prognosis of patients with GBM [Training cohort: HR = 4.036, 95% CI (2.668, 6.104), p < 0.001, Figure 5A; Validation cohort 1: HR = 3.023, 95% CI (2.227, 4.103), p < 0.001, Figure 5C; Validation cohort 2: HR = 19.023, 95% CI (5.750, 62.936), p < 0.001, Figure 5E]. In addition, the multivariable Cox regression analysis indicated that the IRGPs signature was an independent prognosis-related factor [Training cohort: HR = 4.536, 95% CI (2.733, 7.528), p < 0.001, Figure 5B; Validation cohort 1: HR = 3.180, 95% CI (2.288, 4.421), p < 0.001, Figure 5D; Validation cohort 2: HR = 21.331, 95% CI (6.248, 72.826), p < 0.001, Figure 5F].

FIGURE 5

Figure 5. The immune-related gene pairs (IRGPs) signature is an independent prognostic factor of overall survival (OS) of patients with glioblastoma (GBM). (A,B) The result of univariable and multivariable Cox regression analysis in the Training cohort. (C,D) In Validation cohort 1. (E,F) In Validation cohort 2.

Association Between the IRGPs Signature and Clinical Indexes

Clinicopathological data, including age, gender, tumor molecular subtype, KPS (Karnofsky), recurrent, radiotherapy, and chemotherapy, were collected. As shown in Figure 6A, the distribution of age (p = 0.047) and tumor subtype (p = 0.036) between high- and low-risk groups was significantly different. In addition, the risk scores of patients aged ≥50 were significantly increased (p = 0.034, Figure 6B). An analogous phenomenon was observed in patients with mesenchymal GBM (p = 0.004, Figure 6F). Compared with patients without chemotherapy, patients treated with chemotherapy had lower risk scores (Figure 6H). However, there was no difference in the risk scores between female patients and male patients (p = 0.790, Figure 6C), patients with KPS ≤70 and KPS >70 (p = 0.120, Figure 6D), primary tumor and recurrent tumor (p = 0.600, Figure 6E), and, patients treated with radiotherapy and without radiotherapy (p = 0.180, Figure 6G).

FIGURE 6

Figure 6. Correlation between risk score and clinical characteristics. (A) Heat map for the distribution of clinicopathological features between high- and low-risk groups. (B) The difference of risk scores among patients aged <50 and aged ≥50. (C) Between female and male. (D) Between patients with KPS ≤70 and >70. (E) Between primary tumor and recurrent tumor. (F) Between different molecular subtypes. (G) Between patients with radiotherapy and without radiotherapy. (H) Between patients with chemotherapy and without chemotherapy. *p < 0.05, ***p < 0.001.

Relationship Between the IRGPs Signature and TIICs

The CIBERSORT method was applied to quantify the proportions of 22 TIICs. A total of 248 samples (116 low-risk patients and 132 high-risk patients) met the cut-off value: p < 0.05. The highest quantity of immune cell among 22 TIICs was macrophages (M0, M1, and M2) (49.73%), followed by monocytes (12.53%) and resting memory CD4 T cells (10.45%) (Figures 7A,B). Compared with high-risk patients, the fractions of macrophages M1 (p = 0.014, Figure 7D) and resting memory CD4 T cells (p = 0.033, Figure 7E) in low-risk patients were significantly elevated, whereas, the fractions of macrophages M0 (p < 0.001, Figure 7C) and regulatory T cells (Tregs) (p = 0.021, Figure 7F) were decreased.

FIGURE 7

Figure 7. Relationship between immune-related gene pairs (IRGPs) and tumor-infiltrating immune cells (TIICs). (A) The proportions of 22 TIICs in 248 samples were quantified by CIBERSORT. (B) The difference of 22 TIICs between high- and low-risk groups. (C) The difference of macrophages M0. (D) The difference of macrophages M1. (E) The difference of resting CD4 memory T cells. (F) The difference of regulatory T cells (Tregs) between high- and low-risk groups. *p < 0.05, ***p < 0.001.

Then, we investigated the distributions of TIICs between score:1 and score:0 of 10 IRGPs. As shown in Supplementary Figure 4, we found that the proportion of macrophages M0 in IL10RA|OSMR:0, LEFTY2|MSTN:1, and SLC11A1|TRIM22:1 was increased. The fractions of macrophages M1 and resting memory CD4 T cells in FCGR2B|PRKCB:0, LEFTY2|MSTN:0, and SLC11A1|TRIM22:0 were also increased. Moreover, the content of Tregs in FCGR2B|PRKCB:1, LEFTY2|MSTN:1, SLC11A1|TRIM22:1, BMP2|NOV:0, and FCGR2B|PRKCB:0 was upregulated. In addition, the difference of the distributions of several TIICs between score:1 and score:0 of 10 IRGPs was observed (Supplementary Table 1).

Gene Set Enrichment Analysis

The results of GSEA illustrated that a total of 22 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway-related gene sets were enriched. Among them, a few pathways were closely correlated with the immune system, such as “cytokine–cytokine receptor interaction,” “chemokine signaling pathway,” “Janus Kinase (JAK)-signal transducer and activator of transcription (STAT) signaling pathway,” “complement and coagulation cascades,” “Nucleotide-binding oligomerization domain (NOD)-like receptor signaling pathway,” and “intestinal immune network for IgA production” (Figure 8).

FIGURE 8

Figure 8. Gene set enrichment analysis (GSEA) between high and low immune risk groups. (A) 22 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway-related gene sets. (B) Cytokine–cytokine receptor interaction. (C) Chemokine signaling pathway. (D) Complement and coagulation cascades. (E) Intestinal immune network for IgA production. (F) JAK–STAT signaling pathway. (G) NOD-like receptor signaling pathway.

Discussion

Glioblastoma is the frequently occurring and most aggressive malignancy in the central nervous system (25). The 5-year survival rate of patients with GBM after diagnosis is still limited to 6.8% (3). It is important to seek the tumor-specific markers and build risk stratification for the assessment of prognosis in GBM, which may facilitate the development of novel strategies for the diagnosis and therapy of GBM.

The immune system is a considerable influencing factor for the progression of cancers (5). Immunotherapy as a novel treatment method has revolutionized cancer therapy (8), which also provides new hope for the treatment of GBM. In this study, we downloaded the IRG expression matrix from the TCGA database to establish IRGPs and develop an IRGP signature for the prediction of clinical outcomes of patients with GBM. Based on the cut-off value, patients were segmented into two subgroups. The survival analysis showed that patients with high-risk scores had a more unfavorable prognosis. In addition, the univariable and multivariate Cox regression analysis in three cohorts demonstrated the IRGPs signature was an independent prognostic factor in patients with GBM. However, the HR of univariable and multivariate Cox regression analysis in the Validation cohort 2 was 19.023 and 21.331, respectively, substantially higher than that in the Training cohort and Validation cohort 1. It may be caused by the patients with low-proportion low-risk scores in the Validation cohort 2 (high-risk score:low-risk score = 71: 23).

Later, we estimated the predicted capacity of this IRGPs signature. The AUCs of the IRGPs signature for predicting OS probability in the Training cohort reached 0.800 at 1 year, 0.901 at 3 years, and 0.948 at 5 years, which were significantly superior to the AUCs based on the clinical characteristics, such as age, gender, chemotherapy, recurrent, and so on (Supplementary Figures 2A–C). The AUCs for 1-, 3-, and 5-year OS in the Validation cohorts 1 and 2 were similar to that in the Training cohort and were also obviously higher than that of clinical indexes (Supplementary Figures 2D–I). All data demonstrated that the IRGPs signature was suitable to predict 1-, 3- and 5-year OS in GBM.

Recently, several studies explored prognostic biomarkers and evaluated the prognosis of GBM with immune-related features. For example, Zhang et al. (26) used a random forest algorithm to screen the tumor microenvironment (TME)-related genes and developed a TME score with five genes using the multivariate Cox regression analysis. The Kaplan–Meier plot showed a significant survival difference between patients with high- and low-TME scores. A study utilized the univariable Cox regression analysis and LASSO regression analysis to establish an immune-based prognostic scoring model for predicting survival probability in GBM with the AUC = 0.657 at 1 year, 0.667 at 3 years, and 0.667 at 5 years (27). Another study constructed a Foxp3-related immune prognostic signature in patients with GBM using the LASSO regression analysis (28). The AUC for predicting 1- and 3-year OS was 0.633 and 0.695, respectively. In 2018, Zhou et al. (29) performed the survival analysis to select prognostic immune-related-lncRNA and then carried out the multivariate Cox regression analysis to develop a six immune-related lncRNA prognostic risk model. The AUCs of the model for predicting 5-year OS was 0.842. All of those studies illustrated the well-predictive ability of immune-related features. Different from previous studies, herein, we constructed IRGPs in GBM and applied the univariable Cox regression analysis and LASSO regression analysis to preliminary filtrate IRGPs. Finally, we executed the multivariate Cox regression analysis to further screen and develop IRGPs signature. Up to now, this is the first time using IRGPs, a novel immune-related feature, to build a signature for predicting OS probability in GBM. The AUCs of the signature were obviously superior to that in the study of Qin et al. (27) and Guo et al. (28) and slightly preferable to the AUC in the study of Zhou et al. (29), showing the IRGPs signature is robust and reliable.

The IRGPs signature consisted of 10 gene pairs, including 16 IRGs (BMP2, CHGB, FCGR2B, FGF12, IL10RA, LEFTY2, SLC11A1, IL13RA2, NOV, OXTR, GDF11, PRKCB, FGFR2, OSMR, MSTN, and TRIM22). Between low- and high-risk groups, the expression of 11 genes (CHGB, FCGR2B, IL10RA, SLC11A1, IL13RA2, NOV, OXTR, OSMR, LEFTY2, FGFR2, and MSTN) was with statistical difference (Supplementary Figure 3). Fibroblast growth factor receptor 2 (FGFR2) and fibroblast growth factor 12 (FGF12) are the members of the fibroblast growth factor family, which is involved in the activation of the Kirsten Rat Sarcoma (RAS)-Mitogen-Activated Protein Kinase (MAPK) and the Threonine Kinase (AKT) and the PI3K–AKT pathway, and the cell proliferation, differentiation, and apoptosis of cancers (30). The studies demonstrated that they could influence tumor development via regulating macrophages, enhance PD-1 and PD-L1 expression, and promote inflammation (31–34). Fc fragment of IgG receptor IIb (FCGR2B) encoded a receptor for the Fc region of immunoglobulin gamma complexes. It is involved in several regulatory functions, such as modulating B cells producing antibodies and the phagocytosis of immune complexes (35). Moreover, it could bring about the downmodulation of cell activation triggered through regulating antigen receptors on B cells and T cells and upregulation of IgG responses, which contribute to humoral immune responses (35, 36). Protein kinase C beta (PRKCB) is a kind of serine-specific protein kinases, and it is involved in various cellular signaling processes, such as the activation of B cells and cell proliferation (37). Moreover, it also could mediate the differentiation of dendritic cells (37). Tripartite motif-containing 22 (TRIM22) and solute carrier family 11 member 1 (SLC11A1) participate in host resistance to the pathogen and interferon-mediated antiviral effects (38–40). SLC11A1 also participates in the activation of macrophages to promote macrophage immune effector functions, such as inhibiting immune inflammation and altering the processing and presentation of antigens by enhancing the autoimmune T cell (39, 40). TRIM22 could regulate macrophage autophagy (38). Interleukin 10 receptor subunit alpha (IL10RA), interleukin 13 receptor subunit alpha 2 (IL13RA2), and oncostatin M receptor (OSMR) are cytokine receptors and could mediate the immunosuppressive signal and the synthesis of proinflammatory cytokines. Previous studies showed IL10RA was highly related to the clinical outcomes of colorectal cancer and melanoma (41, 42). IL13RA2 is overexpressed in over 75% of patients with GBM and related to prognosis and participates in the downregulation of immune response mediated by helper T cells (43, 44). Immunization with IL13RA2 DNA vaccine could inhibit tumor growth and induces tumor immunity (45). OSMR could transduce signaling induced by oncostatin M, IL6, and IL31 and regulated immune cell-mediated metabolism (46). Basic research showed that OSMR contributed to the adjusting of extracellular matrix process and local immune response in GBM, which was highly associated with the development of GBM (47). Chromogranin B (CHGB), oxytocin receptor (OXTR), and cellular communication network factor 3 (CCN3) also belong to cytokine and cytokine receptor, which play a vital role in the homeostasis of innate and adaptive immunity and remarkably affect anti-tumor immunity through activating cytotoxic T cells and cancer immunotherapy (48, 49). Bone morphogenetic protein 2 (BMP2), myostatin (MSTN), growth differentiation factor 11 (GDF11), and left–right determination factor 2 (LEFTY2) belong to the transforming growth factor-beta superfamily, which has pleiotropic effects on cell proliferation, differentiation, invasion, metastasis, and apoptosis of cancer cells. It also influences inflammation, neurogenesis, and immune response via regulating immune cell differentiation and the activity of immune components (50–52).

Nowadays, emerging studies show treatment targeting the tumor microenvironment is an inspiring way to overcome therapeutic escape issues (53). Herein, we depicted the landscape of immune cells in the tumor microenvironment of GBM and found macrophages was the highest proportion immune cell, which was in line with previous findings (53, 54). Furthermore, we identified that the proportion of macrophages M1 in the low-risk group was significantly increased, whereas, macrophages M0 was opposite. After induced by lipopolysaccharide and interferon-γ stimulation, macrophages M0 may polarize into M1 that could promote GBM cell proliferation, invasion, and metastasis (55–57). Additionally, in the high-risk group, the fractions of Tregs were also significantly increased. The accumulation of Tregs in the tumor microenvironment functions as immunosuppression via inhibiting proliferation and degranulation of other immune cells (58). Increased proportion of Tregs was highly associated with tumor recurrence and decrease of survival and efficacy of immunotherapeutic strategies in GBM (59).

In this study, the IRGPs signature demonstrated an excellent predictive capability for patients with GBM. However, there are still several limitations that need to be improved. (1) As all cases were obtained from open databases, the possibility of selection bias cannot be excluded. (2) The IRGPs signature was developed based on RNA-seq data and microarray expression. It is costly and time-consuming. (3) Due to all samples were collected from the public database, the potential selection bias could not be eliminated, and some clinical information, such as treatments, molecular subtype, and so on, were missing in some samples, which may result in information bias. Meanwhile, the study was a retrospective study. Hence, validation of prospective samples, especially samples with complete clinical information, was still needed. (4) In the Training cohort, we did not obtain the clinical information: the status of IDH-mutation, which may present in enrolled patients with GBM, especially the patients with recurrent tumor. (5) No experimental studies were performed to testify the finding in this study. Thus, further basic experiments are warranted to validate the finding of this research.

Conclusions

In conclusion, in this study, we constructed a 10-IRGPs signature for predicting 1-, 3-, and 5-year OS in GBM. The IRGPs signature would be an effective and convenient means to predict the clinical outcome of GBM and provide a novel way for risk assessment.

Data Availability Statement

Publicly available datasets were analyzed in this study, this data can be found at: The Cancer Genome Atlas (https://portal.gdc.cancer.gov/).

Author Contributions

SW: data collection, data analysis, and writing. XX: study design, data analysis, and revised the manuscript. All authors reviewed the manuscript.

Funding

This study was supported by National Natural Science Foundation of China (No. 81770739), the Clinical Big Data System Construction Project of Central South University (No. 46), and Hunan Key Laboratory of Liver viscercal manifestation (No. 420010087).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.564960/full#supplementary-material

Supplementary Figure 1. Distribution of risk score and survival status. (A) Training cohort, (B) Validation cohort 1, (C) Validation cohort 2.

Supplementary Figure 2. The area under the curves (AUCs) of clinical indexes for predicting 1-, 3-, and 5-year OS. (A–C) in the Training cohort. (D–F) in Validation cohort 1. (G–I) in Validation cohort 2. AUC, the area under the receiver operating characteristic curve; OS, overall survival.

Supplementary Figure 3. The different expression of 16 immune-related genes (IRGs) between high- and low-risk groups. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Supplementary Figure 4. Comparison of 22 immune cells between 10 IRGPs scores:1 and 0. (A) The distributions of Macrophages M0 between score:1 and score:0 of IL10RA/OSMR, (B) LEFTY2/MSTN, (C) SLC11A1|TRIM22. (D) The distributions of Macrophages M1 between score:1 and score:0 of FCGR2B|PRKCB, (E) LEFTY2|MSTN, (F) SLC11A1|TRIM22. (G) The distributions of resting memory CD4 T cells between score:1 and score:0 of BMP2|NOV, (H) FCGR2B|PRKCB, (I) LEFTY2|MSTN, (J) SLC11A1|TRIM22. (K) The distributions of Tregs cells between score:1 and score:0 of BMP2|NOV (L) FCGR2B|PRKCB, (M) LEFTY2|MSTN, (N) SLC11A1|TRIM22, (O) FCGR2B|PRKCB.

Supplementary Table 1. P-values of Wilcoxon test for the distribution of 22 immune cells between 10 IRGPs score:1 and 0.

References

1. Ferlay J, Colombet M, Soerjomataram I, Mathers C, Parkin DM, Piñeros M, et al. Estimating the global cancer incidence and mortality in 2018: GLOBOCAN sources and methods. Int J Cancer. (2019) 144:1941–53 doi: 10.1002/ijc.31937

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ostrom QT, Patil N, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2013-2017. Neuro Oncol. (2020) 22(Suppl. 2):iv1–iv96. doi: 10.1093/neuonc/noaa200

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2012-2016. Neuro Oncol. (2019) 21(Suppl. 5):v1–v100. doi: 10.1093/neuonc/noz150

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Tesileanu CMS, Dirven L, Wijnenga MMJ, Koekkoek JAF, Vincent AJPE, Dubbink HJ, et al. Survival of diffuse astrocytic glioma, IDH1/2 wildtype, with molecular features of glioblastoma, WHO grade IV: a confirmation of the cIMPACT-NOW criteria. Neuro Oncol. (2020) 22:515–23. doi: 10.1093/neuonc/noz200

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gerada C, Ryan KM. Autophagy, the innate immune response and cancer. Mol Oncol. (2020) 14:1913–29. doi: 10.1002/1878-0261.12774

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Brown NF, Carter TJ, Ottaviani D, Mulholland P. Harnessing the immune system in glioblastoma. Br J Cancer. (2018) 119:1171–81. doi: 10.1038/s41416-018-0258-8

CrossRef Full Text | Google Scholar

7. Corthay A. Does the immune system naturally protect against cancer? Front Immunol. (2014) 5:197. doi: 10.3389/fimmu.2014.00197

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lakshmi Narendra B, Eshvendar Reddy K, Shantikumar S, Ramakrishna S. Immune system: a double-edged sword in cancer. Inflamm Res. (2013) 62:823–34. doi: 10.1007/s00011-013-0645-9

PubMed Abstract | CrossRef Full Text | Google Scholar

9. McGranahan T, Therkelsen KE, Ahmad S, Nagpal S. Current state of immunotherapy for glioblastoma. Curr Treat Options Oncol. (2019) 20:24. doi: 10.1007/s11864-019-0619-4

CrossRef Full Text | Google Scholar

10. De Felice F, Pranno N, Marampon F, Musio D, Salducci M, Polimeni A, Tombolini V. Immune check-point in glioblastoma multiforme. Crit Rev Oncol Hematol. (2019) 138:60–9. doi: 10.1016/j.critrevonc.2019.03.019

CrossRef Full Text | Google Scholar

11. Wang S, Yu ZH, Chai KQ. Identification of EGFR as a novel key gene in clear cell renal cell carcinoma (ccRCC) through bioinformatics analysis and meta-analysis. Biomed Res Int. (2019) 2019:6480865. doi: 10.1155/2019/6480865

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhao J, Wang L, Hu G, Wei B. A 6-Gene risk signature predicts survival of glioblastoma multiforme. Biomed Res Int. (2019) 2019:1649423. doi: 10.1155/2019/1649423

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhang L, Zhu P, Tong Y, Wang Y, Ma H, Xia X, et al. An immune-related gene pairs signature predicts overall survival in serous ovarian carcinoma. Onco Targets Ther. (2019) 12:7005–14. doi: 10.2147/OTT.S200191

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Sun XY, Yu SZ, Zhang HP, Li J, Guo WZ, Zhang SJ. A signature of 33 immune-related gene pairs predicts clinical outcome in hepatocellular carcinoma. Cancer Med. (2020) 9:2868–78. doi: 10.1002/cam4.2921

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Wu J, Zhao Y, Zhang J, Wu Q, Wang W. Development and validation of an immune-related gene pairs signature in colorectal cancer. Oncoimmunology. (2019) 8:1596715. doi: 10.1080/2162402X.2019.1596715

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Watanabe T, Nobusawa S, Kleihues P, Ohgaki H. IDH1 mutations are early events in the development of astrocytomas and oligodendrogliomas. Am. J. Pathol. (2009) 174:1149–53. doi: 10.2353/ajpath.2009.080958

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hartmann C, Meyer J, Balss J, Capper D, Mueller W, Christians A, et al. Type and frequency of IDH1 and IDH2 mutations are related to astrocytic and oligodendroglial differentiation and age: a study of 1,010 diffuse gliomas. Acta Neuropathol. (2009) 118:469–74. doi: 10.1007/s00401-009-0561-9

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Weller M, van den Bent M, Preusser M, Le Rhun E, Tonn JC, Minniti G, et al. EANO guidelines on the diagnosis and treatment of diffuse gliomas of adulthood. Nat Rev Clin Oncol. (2020) 18:170–86. doi: 10.1038/s41571-020-00447-z

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Brat DJ, Aldape K, Colman H, Figrarella-Branger D, Fuller GN, Giannini C, et al. cIMPACT-NOW update 5: recommended grading criteria and terminologies for IDH-mutant astrocytomas. Acta Neuropathol. (2020) 139:603–8. doi: 10.1007/s00401-020-02127-9

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Cancer Genome Atlas Research Network, Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, et al. Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N Engl J Med. (2015) 372:2481–98. doi: 10.1056/NEJMoa1402121

PubMed Abstract | CrossRef Full Text

21. Lapointe S, Perry A, Butowski NA. Primary brain tumours in adults. Lancet. (2018) 392:432–46. doi: 10.1016/S0140-6736(18)30990-5

CrossRef Full Text | Google Scholar

22. Klopfenstein Q, Truntzer C, Vincent J, Ghiringhelli F. Cell lines and immune classification of glioblastoma define patient's prognosis. Br J Cancer. (2019) 120:806–14. doi: 10.1038/s41416-019-0404-y

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li L, Ouyang Y, Wang W, Hou D, Zhu Y. The landscape and prognostic value of tumor-infiltrating immune cells in gastric cancer. PeerJ. (2019) 7:e7993. doi: 10.7717/peerj.7993

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Alexander BM1, Cloughesy TF. Adult glioblastoma. J Clin Oncol. (2017) 35:2402–9. doi: 10.1200/JCO.2017.73.0119

CrossRef Full Text | Google Scholar

26. Zhang J, Xiao X, Zhang X, Hua W. Tumor microenvironment characterization in glioblastoma identifies prognostic and immunotherapeutically relevant gene signatures. J Mol Neurosci. (2020) 70:738–50. doi: 10.1007/s12031-020-01484-0

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Qin Z, Zhang X, Chen Z, Liu N. Establishment and validation of an immune-based prognostic score model in glioblastoma. Int Immunopharmacol. (2020) 85:106636. doi: 10.1016/j.intimp.2020.106636

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Guo XY, Zhang GH, Wang ZN, Duan H, Xie T, Liang L, et al. A novel Foxp3-related immune prognostic signature for glioblastoma multiforme based on immunogenomic profiling. Aging (Albany NY). (2021) 13:3501–17. doi: 10.18632/aging.202282

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhou M, Zhang Z, Zhao H, Bao S, Cheng L, Sun J. An immune-related six-lncRNA signature to improve prognosis prediction of glioblastoma Multiforme. Mol Neurobiol. (2018) 55:3684–97. doi: 10.1007/s12035-017-0572-9

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang C, Fu L, Fu J, Hu L, Yang H, Rong TH, et al. Fibroblast growth factor receptor 2-positive fibroblasts provide a suitable microenvironment for tumor development and progression in esophageal carcinoma. Clin Cancer Res. (2009) 15:4017–27. doi: 10.1158/1078-0432.CCR-08-2824

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bohrer LR, Schwertfeger KL. Macrophages promote fibroblast growth factor receptor-driven tumor cell migration and invasion in a CXCR2-dependent manner. Mol Cancer Res. (2012) 10:1294–305. doi: 10.1158/1541-7786.MCR-12-0275

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Li P, Huang T, Zou Q, Liu D, Wang Y, Tan X. FGFR2 promotes expression of PD-L1 in colorectal cancer via the JAK/STAT3 signaling pathway. J Immunol. (2019) 202:3065–75. doi: 10.4049/jimmunol.1801199

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Huang Y, Wang F, Li H, Xu S, Xu W, Pan X, et al. Inhibition of fibroblast growth factor receptor by AZD4547 protects against inflammation in septic mice. Inflammation. (2019) 42:1957–67. doi: 10.1007/s10753-019-01056-4

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Nakayama F, Müller K, Hagiwara A, Ridi R, Akashi M, Meineke V. Involvement of intracellular expression of FGF12 in radiation-induced apoptosis in mast cells. J Radiat Res. (2008) 49:491–501. doi: 10.1269/jrr.08021

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Su K, Wu J, Edberg JC, Li X, Ferguson P, Cooper GS, et al. A promoter haplotype of the immunoreceptor tyrosine-based inhibitory motif-bearing FcgammaRIIb alters receptor expression and associates with autoimmunity. I. Regulatory FCGR2B polymorphisms and their association with systemic lupus erythematosus. J Immunol. (2004) 172:7186–91. doi: 10.4049/jimmunol.172.11.7186

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Xiu Y, Nakamura K, Abe M, Li N, Wen XS, Jiang Y, et al. Transcriptional regulation of Fcgr2b gene by polymorphic promoter region and its contribution to humoral immune responses. J Immunol. (2002) 169:4340–6. doi: 10.4049/jimmunol.169.8.4340

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Farren MR, Carlson LM, Lee KP. Tumor-mediated inhibition of dendritic cell differentiation is mediated by down regulation of protein kinase C beta II expression. Immunol Res. (2010) 46:165–76. doi: 10.1007/s12026-009-8118-5

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lou J, Wang Y, Zheng X, Qiu W. TRIM22 regulates macrophage autophagy and enhances Mycobacterium tuberculosis clearance by targeting the nuclear factor-multiplicity κB/beclin 1 pathway. J Cell Biochem. (2018) 119:8971–80. doi: 10.1002/jcb.27153

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Dai YD, Marrero IG, Gros P, Zaghouani H, Wicker LS, Sercarz EE. Slc11a1 enhances the autoimmune diabetogenic T-cell response by altering processing and presentation of pancreatic islet antigens. Diabetes. (2009) 58:156–64. doi: 10.2337/db07-1608

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Nairz M, Fritsche G, Crouch ML, Barton HC, Fang FC, Weiss G. Slc11a1 limits intracellular growth of Salmonella enterica sv. Typhimurium by promoting macrophage immune effector functions and impairing bacterial iron acquisition. Cell Microbiol. (2009) 11:1365–81. doi: 10.1111/j.1462-5822.2009.01337.x

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zadka Ł, Kulus MJ, Kurnol K, Piotrowska A, Glatzel-Plucińska N, Jurek T, et al. The expression of IL10RA in colorectal cancer and its correlation with the proliferation index and the clinical stage of the disease. Cytokine. (2018) 110:116–25. doi: 10.1016/j.cyto.2018.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Cheng S, Li Z, Zhang W, Sun Z, Fan Z, Luo J, et al. Identification of IL10RA by weighted correlation network analysis and in vitro validation of its association with prognosis of metastatic melanoma. Front Cell Dev Biol. (2021) 8:630790. doi: 10.3389/fcell.2020.630790

CrossRef Full Text | Google Scholar

43. Zeng J, Zhang J, Yang YZ, Wang F, Jiang H, Chen HD, et al. IL13RA2 is overexpressed in malignant gliomas and related to clinical outcome of patients. Am J Transl Res. (2020) 12:4702–14.

PubMed Abstract | Google Scholar

44. Sattiraju A, Solingapuram Sai KK, Xuan A, Pandya DN, Almaguel FG, Wadas TJ, et al. IL13RA2 targeted alpha particle therapy against glioblastomas. Oncotarget. (2017) 8:42997–3007. doi: 10.18632/oncotarget.17792

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Nakashima H, Fujisawa T, Husain SR, Puri RK. Interleukin-13 receptor α2 DNA prime boost vaccine induces tumor immunity in murine tumor models. J Transl Med. (2010) 8:116. doi: 10.1186/1479-5876-8-116

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Luo P, Wang PX, Li ZZ, Zhang XJ, Jiang X, Gong J, et al. Hepatic oncostatin m receptor β regulates obesity-induced steatosis and insulin resistance. Am J Pathol. (2016) 186:1278–92. doi: 10.1016/j.ajpath.2015.12.028

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sharanek A, Burban A, Laaper M, Heckel E, Joyal JS, Soleimani VD, et al. OSMR controls glioma stem cell respiration and confers resistance of glioblastoma to ionizing radiation. Nat Commun. (2020) 11:4116. doi: 10.1038/s41467-020-17885-z

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Estrellas KM, Chung L, Cheu LA, Sadtler K, Majumdar S, Mula J, et al. Biological scaffold-mediated delivery of myostatin inhibitor promotes a regenerative immune response in an animal model of Duchenne muscular dystrophy. J Biol Chem. (2018) 293:15594–605. doi: 10.1074/jbc.RA118.004417

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Chikuma S, Kanamori M, Mise-Omata S, Yoshimura A. Suppressors of cytokine signaling: potential immune checkpoint molecules for cancer immunotherapy. Cancer Sci. (2017) 108:574–80. doi: 10.1111/cas.13194

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Zhao H, Wei J, Sun J. Roles of TGF-β signaling pathway in tumor microenvirionment and cancer therapy. Int Immunopharmacol. (2020) 89(Pt B):107101. doi: 10.1016/j.intimp.2020.107101

CrossRef Full Text | Google Scholar

51. Huynh LK, Hipolito CJ, Ten Dijke PA. Perspective on the development of TGF-β Inhibitors for cancer treatment. Biomolecules. (2019) 9:743. doi: 10.3390/biom9110743

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Bai X, Yi M, Jiao Y, Chu Q, Wu K. Blocking TGF-β signaling to enhance the efficacy of immune checkpoint inhibitor. Onco Targets Ther. (2019) 12:9527–38. doi: 10.2147/OTT.S224013

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Quail DF, Joyce JA. The microenvironmental landscape of brain tumors. Cancer Cell. (2017) 31:326–41. doi: 10.1016/j.ccell.2017.02.009

CrossRef Full Text | Google Scholar

54. Engler JR, Robinson AE, Smirnov I, Hodgson JG, Berger MS, Gupta N, et al. Increased microglia/macrophage gene expression in a subset of adult and pediatric astrocytomas. PLoS ONE. (2012) 7:e43339. doi: 10.1371/journal.pone.0043339

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Martinez FO, Gordon S, Locati M, Mantovani A. Transcriptional profiling of the human monocyte-to-macrophage differentiation and polarization: new molecules and patterns of gene expression. J Immunol. (2006) 177:7303–11. doi: 10.4049/jimmunol.177.10.7303

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Mantovani A, Sozzani S, Locati M, Allavena P, Sica A. Macrophage polarization: tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol. (2002) 23:549–55. doi: 10.1016/S1471-4906(02)02302-5

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Poon CC, Sarkar S, Yong VW, Kelly JJP. Glioblastoma-associated microglia and macrophages: targets for therapies to improve prognosis. Brain. (2017) 140:1548–60. doi: 10.1093/brain/aww355

PubMed Abstract | CrossRef Full Text | Google Scholar

58. DiDomenico J, Lamano JB, Oyon D, Li Y, Veliceasa D, Kaur G, et al. The immune checkpoint protein PD-L1 induces and maintains regulatory T cells in glioblastoma. Oncoimmunology. (2018) 7:e1448329. doi: 10.1080/2162402X.2018.1448329

PubMed Abstract | CrossRef Full Text | Google Scholar

59. See AP, Parker JJ, Waziri A. The role of regulatory T cells and microglia in glioblastoma-associated immunosuppression. J Neurooncol. (2015) 123:405–12. doi: 10.1007/s11060-015-1849-3

PubMed Abstract | CrossRef Full Text | Google Scholar