SARS-CoV-2 Pattern Provides a New Scoring System and Predicts the Prognosis and Immune Therapeutic Response in Glioma

The training set contained 667 glioma patients and 1152 normal people as a control. The RNA-Seq data and corresponding clinicopathological data of these patients were obtained from the TCGA website (http://cancergenome.nih.gov/↗, accessed on 1 September 2022). A total of 1152 normal brain tissue samples were obtained from Genotype-Tissue Expression Project (GTEX, https://www.gtexportal.org/home/index.html↗, accessed on 1 September 2022) database.

A total of 333 SARS-Cov-2-related genes were obtained from the human protein atlas (HPA) (https://www.proteinatlas.org/↗, accessed on 1 September 2022) database. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to identify pathways enriched by differentially expressed genes [17]. Then, the mRNA expression and survival data were used to further screen the above genes by univariate COX regression analysis, least absolute shrinkage, selection operator (LASSO) regression, and multivariate COX regression analysis (p < 0.05).

The regression coefficients were linearly integrated with the expression levels of seven SARS-CoV-2-related genes to establish the optimal SARS-CoV-2-related risk signature. The risk score is defined here as the SARS-CoV-2 Index (SC2I) and is calculated as follows:

Expi represents the expression value, and Coefi is the corresponding regression coefficient calculated by multivariate Cox analysis.

Age, sex, MGMT promoter methylation status, IDH mutation status, WHO grade, tumor subtype, and SC2I were combined to construct a nomogram using the R package “survival” and “rms” to predict the probability of clinical events in a personalized way. Calibration curves were adopted to evaluate the accuracy of the nomogram in predicting the 1-, 3-, and 5-year survival rates of glioma patients [18].

GO enrichment analysis and KEGG pathway analysis were used for functional enrichment analysis via GSEA and R software.

The abundance of 22 immune cells was calculated using the Cibersort algorithm [19]. The extent of 28 infiltrating immune cells was quantified using the single sample gene set enrichment analysis (ssGSEA) of the “gsva” package [20]. Spearman correlation analysis was used to obtain the relationship between risk scores and the expression levels of common immune checkpoints. The ESTIMATE Score, Immune Score, and Stromal Score were calculated using the ESTIMATE algorithm via the R package “estimate” to evaluate the role of SC2I in predicting the glioma tumor microenvironment [21]. Finally, the tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu↗, accessed on 1 September 2022) was applied to investigate the role of SC2I in predicting the response of glioma to immunotherapy.

Patient data from the CGGA database (http://www.cgga.org.cn/↗, accessed on 1 September 2022) were used for external validation. The median of the SC2I was used as the boundary to divide the sample into high- and low-risk groups. Different packages in R software were used to draw the Kaplan–Meier curve and ROC curve, and “ggplot2” software was used to draw the PCA scatter plot.

From May 2022 to August 2022, we collected clinicopathological data, including age, gender, WHO grade, IDH mutation status, MGMT promoter methylation status, and tissue samples of nine patients with glioma, from the First Affiliated Hospital of Soochow University. In this paper, we presented the results of the analysis with WHO grade and IDH mutation status as grouping indicators. This study was approved by the Medical Ethics Committee of the First Affiliated Hospital of Soochow University, and informed consent was obtained from all patients. A quantitative real-time polymerase chain reaction was performed. Tissues were stored at −80 °C, and total RNA was extracted from tissues using RNA Trizol reagent (Beyotime, Shanghai, China) following manufacturer’s instructions. The extracted total RNA was reverse-transcribed into cDNA using reverse transcription reagent, and cDNA was amplified in real-time using SYBR fluorescent quantitative PCR reagent. CT values were measured and calculated by computer software. The relative mRNA levels were normalized to β-actin as an internal control by the 2^−ΔΔCT method. Sequences of primers are shown in Supplementary Table S1.

The RNA-seq transcriptome data were preprocessed using Perl programming language (version 5.32.0, created by Larry Wall from Los Angeles, CA, USA). The R software (version 4.2.1, created by Ross Ihaka and Robert Gentlemen from Auckland, New Zealand) was applied for all statistical analyses and graph visualization. The chi-square test was used to compare categorical variables between the high- and low-risk groups. p value ≤ 0.05 was considered statistically significant.

A total of 74 differential expressed genes (DEGs) were identified by differential analysis between the normal and tumor groups (Figure 1A,B). GO enrichment analysis and KEGG pathway analysis were used to identify the biological characteristics of these DEGs. Protein targeting, protein targeting to mitochondrion, protein folding, endomembrane system organization, and cellular response to heat were identified as the most abundant processes using GO analysis (Figure 1C,E). KEGG pathway analysis suggested that DEGs were significantly enriched in protein processing in the endoplasmic reticulum, nucleocytoplasmic transport, lysosome, RNA degradation, and DNA replication (Figure 1D,F). Thirty-two genes were screened by univariate Cox regression analysis using mRNA expression and survival data from TCGA. These 32 genes were then incorporated into the least absolute shrinkage and selection operator (LASSO) regression to avoid overfitting problems, and 26 genes stood out (Figure 2B,C). Finally, multivariate Cox regression analysis was performed to screen out seven genes, which were REEP6, CEP112, LARP4B, CWC27, GOLGA2, ATP6AP1, and ERO1B (Figure 2D). Pan-cancer analysis of these seven genes mentioned above in the TCGA and GTEX showed variable expression levels in the tumor and normal groups (Supplementary Figure S1). As a result, the identification of prognostic SARS-CoV-2-related genes was completed. Among these genes, CEP112 and ATP6AP1 are risk factors for glioma with hazard ratios (HRs) > 1, while REEP6, LARP4B, CWC27, GOLGA2, and ERO1B are protective factors with HRs < 1. In addition, we explored the proteins expressed by REEP6, LARP4B, CWC27, GOLGA2, ATP6AP1 and ERO1B in glioma patients through the HPA database, and REEP6, LARP4B, GOLGA2 and ATP6AP1 showed higher staining intensity in glioma than in normal brain tissue (Supplementary Figure S2A–F). Statistical histograms are shown in Supplementary Figure S2G.

The risk score was defined as the SARS-CoV-2 Index (SC2I), which was calculated using the following formula:

SC2I = (−0.132 × REEP6 exp) + (0.895 × CEP112 exp) + (−0.622 × LARP4B exp) + (−0.347 × CWC27 exp) + (−0.184 × GOLGA2 exp) + (0.066 × ATP6AP1 exp) + (−0.473 × ERO1B exp).

Then, we analyzed the correlations between the genes from the TCGA dataset using the Pearson correlation analysis. Genes with correlation coefficients >0.3 are shown in Figure 3A, with positive correlations in blue and negative ones in red. The survival curves of each gene are shown in Supplementary Figure S3.

Samples in the TCGA dataset were divided into high-risk and low-risk groups based on the median SC2I. First, we performed a Kaplan–Meier (KM) analysis, which showed that patients in the low-risk group had a higher probability of survival than those in the high-risk group (Figure 3B). The heatmap showed that the expression of the seven modeling genes was differentially expressed in the high- and low-risk groups. Compared with the high-risk group, the expression of CEP112 was decreased in the low-risk group, while the expression level of REEP6, LARP4B, and ERO1B was relatively increased (Figure 3C). Patients in the high-risk group had a shorter survival time than those in the low-risk group (Figure 3C). Then, we performed principal component analysis (PCA) in the TCGA cohort, which showed a significant and stable distribution of differences between the high- and low-risk subgroups (Figure 3D). The time-dependent receiver operating characteristic (ROC) curve was used to evaluate the effectiveness of SC2I in predicting the prognosis of glioma patients. As shown in Figure 3E, the AUC of 1-, 3- and 5-year overall survival (OS) rates in the TCGA dataset were 0.895, 0.894, and 0.844, respectively. These results indicate that the SC2I is accurate in predicting the prognosis of patients with glioma.

A heatmap containing mRNA expression and other clinical information (subtype, WHO grade, IDH mutation status, MGMT promoter status, gender, age, and survival status) of SARS-CoV-2-related genes showed that the expression of LARP4B, GOLGA2, CWC27, REEP6, and ERO1B were positively correlated with the risk score. In contrast, CEP112 and ATP6AP1 were negatively correlated. Meanwhile, the WHO grade of gliomas was positively correlated with the SC2I, as expected (Figure 4A). SC2I was significantly higher in patients with gliomas of WHO grade IV than those of WHO grade III and II (Figure 4B). An IDH mutation status-based subgroup analysis showed that patients with IDH wild type had significantly higher SC2I than IDH mutated type (Figure 4B). Subgroup analysis based on tumor subtype suggested that the SC2I in the proneuronal type was lower than that in classical and neural type (Figure 4B). The above findings implied that SC2I could differentiate glioma patients classified based on the above key indicators. Univariate Cox and multivariate Cox of SC2I combined with clinical information showed that SC2I itself was an independent risk factor for the prognosis of glioma (Figure 4C,D). Moreover, the ROC curve of SC2I and other clinical characteristics revealed that its predictive performance (AUC = 0.959) was better than other classical indicators, including tumor subtype (AUC = 0.778), WHO grade (AUC = 0.906), IDH mutation status (AUC = 0.132), MGMT promoter status (AUC = 0.321) and age (AUC = 0.809) (Figure 4E).

By integrating SC2I, age, gender, MGMT promoter status, IDH mutation status, WHO grade, and original subtype, we established a nomogram to accurately predict the 1-, 2- and 3-year survival probability of glioma patients (Figure 5A). The calibration curves showed significant concordance between the predicted and actual rates of OS at 1, 2, and 3 years in the TCGA dataset (Figure 5B–D).

Gene set enrichment analysis (GSEA) reminds us of some key pathways. In the seven SARS-Cov-2-related genes, the up-regulated ones are enriched in cell differentiation involved in embryonic placenta development, mitotic metaphase plate congression, azurophil granule lumen, Ficolin-1-rich granule and vesicle lumen by Gene ontology enrichment analysis. Meanwhile, those downregulated genes are rich in the following pathways: neuron cell–cell adhesion, ligand-gated ion channel signaling pathway, cytoplasmic microtubule organization, negative regulation of synaptic transmission, and the neuromuscular process controlling balance (Figure 5E).

The KEGG pathway analysis was also carried out. The results revealed that the upregulated genes are enriched in systemic lupus erythematosus pathways, N-glycan biosynthesis, glutathione metabolism, pathogenic Escherichia coli infection, and leukocyte transendothelial migration. Those down-regulated genes are enriched in terpenoid backbone biosynthesis, long-term depression, taste transduction, WNT signaling pathway, and mTOR signaling pathway (Figure 5F).

We performed ssGSEA analysis to detect the correlation between the modeling genes and 28 types of tumor-infiltrating immune cells. We found that ATP6AP1 and CEP112 had an overall positive correlation with the expression of tumor-infiltrating immune cells, while the five other modeling genes showed the opposite trend (Figure 6A). Based on the immune functional enrichment analysis results, we further explored the correlation between prognostic SCRGS and glioma immune status. The results showed that the infiltration levels of CD8⁺ T cells, follicular helper T cells, M0, M1, and M2 macrophages and neutrophils in the high-risk group were higher than those in the low-risk group. In comparison, the infiltration levels of activated NK cells, monocytes, and activated mast cells were lower than those in the low-risk group. The levels of infiltration of the remaining immune cells were similar between the two groups (Figure 6B). The immune score results showed that the scores of the high-risk group were significantly higher than those of the low-risk group in 13 immune pathways, indicating that anti-cancer immunity was activated (Figure 6C).

In recent years, immune checkpoint blockers (ICB) have made breakthrough progress in the clinical treatment of a variety of malignant tumors [22,23]. Except for CD200, all immune checkpoints were significantly up-regulated in the high-risk group compared with the low-risk group (Figure 7A). CD274 (PD-L1) and PDCD1 (PD-1) are of great importance and play a crucial role in tumor immunosuppression and immunotherapy [24].

In addition, we further evaluated the association between prognostic SCRGS and the glioma microenvironment. ESTIMATE score, immune score, and stromal score in the high-risk group were higher than those in the low-risk group, suggesting a potential correlation between SCRGS and the tumor microenvironment of gliomas (Figure 7B–D). We analyzed TIDE scores in both groups to further clarify the value of prognostic SCRGS in predicting response to immunotherapy. The high-risk group had lower TIDE and exclusion scores (Figure 7E,F), while the dysfunction scores were similar between the two groups (Figure 7G). These results suggest that patients in the high-risk group are more sensitive to immunotherapy. In conclusion, SCRGS can predict immunotherapy response in patients with glioma to some extent.

To confirm the reliability of SCRGS, we performed an external validation in CGGA. Firstly, 929 glioma patients in CGGA were divided into a high-risk group and a low-risk group based on the median risk score (Figure 8B). It was found that the survival probability and time of the high-risk group were significantly lower than those of the low-risk group (Figure 8A,B). Next, ROC curve analysis was then performed to further verify the accuracy of prognostic SCRGS in predicting the prognosis of patients with glioma in CGGA. The AUC of 1, 3 and 5 years is 0.754, 0.794 and 0.795, respectively (Figure 8C). The expression of seven modeling genes in the validation cohort was also detected. It is not difficult to conclude that the expression of these genes in the CGGA cohort was consistent with that in the TCGA cohort (Figure 8B). These were all consistent with the results in the training set. Above all, these results indicate that the SARS-CoV-2-related gene signature constructed using the TCGA dataset has a good predictive effect regarding the prognosis of glioma patients.

To validate the predictive power of SCRGS in glioma samples, a qRT-PCR analysis of the collected glioma samples was performed. These samples were first divided into WHO grade II-III (4 cases) and grade IV (5 cases), and the SC2I scores were calculated for both groups. It was found that the SC2I of the WHO grade IV group was significantly higher than the corresponding score of the WHO grade II-III group (p = 0.025), which was consistent with our expectation (Figure 8D, left). Then, the tumor samples were divided into wild-type and mutant groups according to IDH mutation status. The results showed that there was no significant difference in SC2I scores between the two groups (Figure 8D, right).

Public data were from https://www.gtexportal.org/home/index.html↗ (GTEX), http://cancergenome.nih.gov/↗ (TCGA), http://www.cgga.org.cn/↗ (CGGA), https://www.proteinatlas.org/↗ (HPA), http://tide.dfci.harvard.edu↗ (TIDE), accessed on 1 September 2022.