Mining of prognosis-related genes in cervical squamous cell carcinoma immune microenvironment

Database sources and pre-processing

The overall flow diagram of present study was summarized in Fig. 1. The RNA-seq counts data, SNP data, and clinical follow-up information were downloaded from the TCGA database. The Reads PerKilobase Million (FPKM) data of RNA-Seq were transformed into Transcripts PerKilobase Million (TPM) expression profiles. In consistent with previous studies, 13 metagenes (shown in ImmuneScore.genes.ids.txt, Supplemental File) corresponded to various immunocyte types, reflecting the different immune functions.

Computational methods of multiple immune scores and result determination

The scores of each sample in the 13 types of metagenes were calculated based on the log2-transformed expression of each gene member in the immune metagene (shown in immune.meta.score.txt, Supplemental File) (Safonov et al., 2017). TIMER (https://cistrome.shinyapps.io/timer/) database (immune.immu.score.txt, Supplemental File) was utilized to calculate the scores of each sample in the immunocyte infiltration (six categories in total) (Li et al., 2016). Moreover, the ImmuneScore (that represents the infiltration of immune cells in tumor tissue), StromalScore (that captures the presence of stroma in tumor tissue) and ESTIMATEScore (that infers tumor purity) of each sample (immune.est.score.txt, Supplemental File) were calculated by ESTIMATE function of R software package (Yoshihara et al., 2013). Detailed description of above 3 scores could be found on the website below: https://bioinformatics.mdanderson.org/public-software/estimate/. Finally, R software package MCPcounter was utilized for the calculation of the abundances of ten immune-related cell (eight categories of immune cells, endothelial cells and fibroblasts) populations in the tumor microenvironment (Supplemental File).

Survival analysis

Patients were divided into several groups according to each specific parameter (including ImmuneScore, StromalScore, ESTIMATEScore and gene expression level). Afterwards, the association between the gene expression level (or level of ImmuneScore, StromalScore, ESTIMATEScore) and overall survival was analyzed by univariate Cox regression model.

The construction of immune scores-related gene modules through WGCNA

To begin with, transcripts with over 75% TPM of >1 and median absolute deviation (MAD) of >median were chosen from the expression profile data of all the obtained samples. Hierarchical clustering for cluster analysis of the samples was also adopted. Subsequently, samples with a distance of over 80,000 were taken as the outlier samples for screening. Moreover, the distance between any two transcripts was calculated by Pearson correlation coefficient, the establishment of the distance between any two transcripts was performed by the R software package Weighted Gene Co-expression Network Analysis (WGCNA) (Xia et al., 2019), and the soft threshold was set as eight for the screening of the co-expression modules. The co-expression network has been suggested to conform to the scale-free network. In other words, the logarithm of node with the connectivity of k (log(k)) should be negatively correlated with the logarithm of the occurrence probability of the specific node (log(P(k))), and the correlation coefficient should be >0.85. Proper β value was selected in order to ensure the network as a scale-free network. The expression matrix was subsequently transformed into the adjacent matrix, and the latter was further transformed into the topological matrix for gene clustering based on Topological Overlap Matrix (TOM) utilizing the average-linkage hierarchical clustering method in accordance with the mixed dynamic shear tree standard. In addition, the gene number of each gene network module was set at least 30. The dynamic shear method was used to determine the gene module, followed by calculation of the eigengene value of each module in succession. Afterwards, clustering analysis was performed on the modules, in which, modules close to each other were merged into a new module, with re-set appropriate height, deepSplit and minModuleSize values. Finally, the association of the acquired gene modules with ImmuneScore, StromalScore and ESTIMATEScore were separately calculated, in order to explore the gene modules with high correlation for further research.

Establishment of the gene interaction network and functional analysis

Genes were mapped into the String database (Szklarczyk et al., 2019). The gene-gene interactions were acquired at the score threshold of >0.4, followed by visualization using Cytoscape software (Shannon et al., 2003). Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis was performed by utilizing the clusterprofile R package (Xu et al., 2019) to examine the signaling pathways affected by these genes.

The immune-related scores based on ESTIMATE algorithm is the most suitable immune scoring method for CSCC

To be specific, we retrieved CSCC samples from the TCGA database and analyzed their scores in 23 types of scoring systems, including 13 types of metagenes scores, six types of immunocyte infiltration scores, three types of immune-related scores according to ESTIMATE algorithm (ImmuneScore, StromalScore and ESTIMATEScore) and 10 types of abundances of immune-related cell. In addition, Spearman’s correlation coefficient was used to calculate the correlations among these scoring systems (shown in Fig. 2). As shown in Fig. 2A, the average correlation between different types of immune-related scores was greater than 0.4. Among which, ImmuneScore (R = 0.59), Co_inhibition (R = 0.59) and LCK (R = 0.62) had the highest relevance with other immune scores. These findings showed that there were fine consistency and comparability between different immune scoring systems. The clustering heat maps of various types of scoring systems were shown in Fig. 2B, suggesting the great correlation among the scoring systems MHC1, MHC2, Monocytic lineage, Dendritic, Macrophages, ESTIMATEScore, ImmuneScore, Tfh, LCK, Co_stimulation, Co_inhibition, Mete_ImmuneScore, Neutrophil and STAT1. We further investigated the average correlation among immune scores according to four different algorithms. As shown in Fig. 2C, the immune-related scores calculated by the ESTIMATE algorithm harbored the highest average correlation with the other three algorithms, which is greater than 0.52 on average. These findings implicated that the immune-related scores based on ESTIMATE algorithm were the most representative immune scoring methods for CSCC.

It is widely accepted that HPV infection has a significant association with the occurrence and progression of CSCC (Ding et al., 2019). Therefore, we separately analyzed the ImmuneScore, StromalScore and ESTIMATEScore distribution among CSCC patients with or without HPV infection. As shown in Figs. 3A–3C, the three immune-related scores in CSCC with HPV infection were significantly higher than those without HPV infection. It should be noted that ImmuneScore was most significantly correlated with the infection status of HPV (p < 0.05).

Subsequently, in order to investigate the association between the above three immune-related scores and prognosis, samples were sorted based on the median of scores of all samples. And then, prognostic difference was analyzed by Kaplan–Meier method (Zou, O’Malley & Mauri, 2007). As a result, the prognosis of samples in different groups was significantly different (shown in Fig. 4). And the five-year survival rate of samples with high ImmuneScore and ESTIMATEScore were significantly superior in comparison with those with low scores, suggesting that the three immune-related scores on the basis of ESTIMATE algorithm could be accepted as promising novel prognostic markers for CSCC.

A large number of somatic mutations of HLA genes have been reported in CSCC, strongly indicating that loss of function due to HLA mutations is tightly correlated with the immune escape of cancer cells (Xiao et al., 2013). It is of great significance for us to analyze the changes of HLA gene sequence in tumor patients. In addition, the mutation of TP53, a tumor suppressor gene, can induce unlimited proliferation and apoptosis resistance of tumor cells (Laprano et al., 2014; Li et al., 2015). Next, we focused on analyzing the associations of three immune-related scores with mutations of HLA and TP53. To this end, we extracted the mutation data of HLA-A, HLA-B, HLA-C and TP53 from the mutect-processed SNP database and then calculated the three immune-related scores based on ESTIMATE algorithm in HLA-A, HLA-B, HLA-C and TP53 mutation and non-mutation groups. As shown in Fig. 5, there was higher level of ImmuneScore in HLA-A and HLA-B mutation groups compared with wild-type groups, while there was also higher level of ESTIMATEScore in HLA-B mutation groups but lower level in TP53 mutation groups comparison with that in wild-type groups.

In summary, we demonstrated that the immune-related scores on the basis of ESTIMATE algorithm were the most proper immune scoring method for CSCC. Additionally, the co-expressed genes with remarkable correlation with these three immune-related scores might be considered as the representative genes in CSCC immune microenvironment, which could be further validated as potential prognostic markers and novel therapeutic targets of CSCC.

Screening of the representative genes in the immune scores-related gene modules

In this section, clustering analysis was first conducted through hierarchical clustering. As shown in Fig. 6A, a total of 296 samples were finally screened out among all the outlier samples, which had a distance of larger than 80,000. Subsequently, the weight co-expression network was constructed by WGCNA with β=8 to guarantee the scale-free network (Figs. 6B and 6C). Afterwards, dynamic shear method (Dong & Horvath, 2007) was utilized to determine the gene modules, and clustering analysis was performed on these modules. Additionally, modules with close distance were further merged into the new module, having height, deepSplit and minModuleSize set to 0.25, 2 and 30, respectively. Finally, a total of 30 modules were acquired (Fig. 6D). Of note, the grey module indicated gene sets that could not be clustered into other modules. The transcripts of each module were counted and displayed in Table 1. In total, 6,679 transcripts were allocated to 29 co-expression modules. The correlations of the eigenvectors of these 30 modules with ImmuneScore, StromalScore and ESTIMATEScore were subsequently calculated, respectively. As shown in Fig. 6E, the yellow module obviously harbored extremely high association with these three immune-related scores based on ESTIMATE algorithm containing 422 genes.

The gene functions in the yellow module were subsequently analyzed. Meanwhile, KEGG and GO enrichment analysis was also conducted using the clusterProfiler of R software package, with flase discovery rate (FDR) set as <0.05. The detailed enrichment results were shown as supporting information file (yellow enrich.txt). As a result, the genes in the yellow module were enriched into 50 KEGG pathways, 670 GO biological processes (BP), 85 GO cellular components (CC) and 74 molecular functions (MF). The most significant top 20 KEGG pathways and GO terms were shown in Fig. 7. The enriched pathways mainly included Th1 and Th2 cell differentiation, cytokine-cytokine receptor interaction and so on. And the enriched biological processes primarily included T cell activation, leukocyte cell–cell adhesion and so on. The enriched cell components mainly included MHC class II protein complex and T cell receptor complex, and so on. The enriched molecular functions mainly included cytokine receptor activity and MHC class II receptor activity, and the rest. Intriguingly, these enriched pathways and GO term have previously been reported to have close association with CSCC and its immune microenvironment (Roca et al., 2019; Wang et al., 2017; Yasmeen et al., 2010; Zehbe et al., 2005).

Finally, to further mine the immune scores-related genes, the weight co-expression relationship between genes in the yellow modules was calculated, with the weight threshold greater than 0.2. Cytoscape software was used for derivation and visualization of the co-expression network of these genes (as shown in Fig. 8A). Afterwards, we further analyzed the topological properties of the network, which contained 244 nodes and 4,083 edges, indicating that genes with greater association with modules had more close correlation with other genes in the network. As shown in Fig. 8B, the degree distribution of the network was further analyzed, suggesting that the degree of the majority of nodes was extremely small, while the degree of a few nodes was rather large, which was consistent with the characteristics of biological network. The correlation between the gene and the module was further calculated. As shown in Fig. 8C, the correlation between most genes and the module was over 0.6, suggesting a high expression similarity between the genes in the module. Moreover, a total of 26 genes (Table 2 and lst.genes.txt as Supplemental File) with a correlation over 0.9 and a degree over 50 in the network were selected, with seven members of LCK Metagenes, and one member of Co_inhibition Metagenes. Thus, 18 new representative immune microenvironment-related genes were finally screened.

Function analysis of 18 novel representative immune microenvironment-related genes in CSCC patients

Firstly, to further analyze the functions of these 18 novel representative immune microenvironment-related genes, the R software package clusterProfiler was utilized for KEGG and GO enrichment analysis, with the significance FDR set at <0.05. The detailed results were summarized in lst enrich.txt (Supplemental File). In brief, these 18 genes were enriched into 11 KEGG pathways, 202 GO biological processes, 8 GO cell components, 19 molecular functions. The most significant 20 KEGG pathways and GO terms were shown in Fig. 9, the majority of which were involved in the proliferation, growth and differentiation of T cells. Intriguingly, LAPTM5, EVI2A and MS4A6A were not enriched in any signaling pathways and GO term, indicating that the functions of these three genes remained completely unclear, which is the focus of our further studies.

Secondly, to further investigate the potential roles of the 18 novel representative immune microenvironment-related genes in clinical practice, the R package corrgram was utilized for the calculation of the association between these genes and immune checkpoints (PDCD1, CD274, PDCD1LG2, CTLA4, CD86, CD80, CD276, VTCN1). As shown in Fig. 10, apart from CD276 and VTCN1, the other 6 immune checkpoints were significantly related to these 18 genes, with an average correlation coefficient over 0.5, which indicated that these immune microenvironment-related genes might be promising targets for immunotherapy.

Finally, the prognostic significance of 18 novel representative immune microenvironment-related genes was assessed. According to the median of gene expression, samples were categorized into high and low expression groups. And then the differences of prognosis between these groups were analyzed. As shown in Fig. 11, high expression of 13 genes were significantly associated with better overall survival according to the threshold of p < 0.05, suggesting that these genes might be closely associated with patient prognosis.

Validation of the correlations of 18 immune microenvironment-related genes with ImmuneScore for CSCC patients by using external dataset

External database was used for further validation of the correlations of 18 immune microenvironment-related genes with the immune-related scores according to ESTIMATE algorithm for CSCC patients. Standardized expression matrix was downloaded and extracted from an independent dataset GSE44001 (Lee et al., 2013) from Gene Expression Omnibus (GEO). R packages hgu133plus2.db was utilized to map a probe for gene to extract the expression profiles of these 18 genes, followed by the calculation of the ImmuneScore for each sample using R software package ESTIMATE. Subsequently, the Pearson correlation was calculated between expression of these genes and the level of ImmuneScore for every CSCC sample in this dataset. As shown in Fig. 12, apart from CCR5 (P = 0.867, R = 0.01), the other 17 genes were significantly associated with the ImmuneScore, which was consistent with our previous findings.