TMB high group and immunityhigh group have a better overall survival of OC.
To explore if a high load of somatic mutation is associated with enhanced tumor immune cell infiltration and result in prolonged overall survival, datasets comprising RNAseq from 379 TCGA tumors and mutation annotation files from 436 TCGA tumors (272 which overlapped with RNAseq data) were obtained from the TCGA Data Portal (https://portal.gdc.cancer.gov/) in March 2019. We evaluated all OC cases with complete gene expression data and clinical information in TCGA ( n = 417, except for the blank values). Based on the somatic mutation data across all the OC cases in this study, calculated TMB scores range from 0.0263 to 6.5260. Interestingly, the TMB scores of OC cases at Grade 3 and Grade 4 were significantly higher than those of Grade 1 and Grade 2 OC cases (Fig. 2A, p = 0.048). In addition, a multivariate Cox regression analysis of clinical factors of OC patients includes age, grade, and tumor residual. The results show that high TMB was associated with a good prognosis (HR = 0.843, p = 0.026), while the age of patients and tumor residual was correlated with a poor prognosis. TMB can be used as an independent prognostic indicator for OC patients (Fig. 2B). Then based on the ssGSEA scores of OC patients, we identified three immune clusters by unsupervised clustering, 193 cases were defined as immunityhigh, while 149 cases were considered as immunitylow by comprehensive infiltration level of all 29 immune subsets (Fig. 2C). Moreover, the PD-L1 expression of the immunityhigh group is significantly higher than that of immunitylow group (Fig. 2C). Furthermore, the tumor purity in immunityhigh group was significantly lower than that in immunitylow group, which verified the reliability of our classification (Fig. 2D).
We further accessed whether there is any association between the TMB or immune infiltration profile and the overall survival. Kaplan-Meier survival curves showed that the overall survival time of TMBhigh or immunityhigh group was significantly longer than TMBlow or immunitylow group (Fig. 2E, 2.58 years vs. 1.92 years, p = 0.0295 in the log-rank test. Figure 2F, 2.50 years vs. 2.45 years, p = 0.00082 in log-rank test). We then downloaded and analyzed the expression of reference genes representing 22 immune cell subgroups from CIBERSORT (see Materials and Methods), and further evaluated the abundance of 22 different immune cell subsets in each of the TMB subgroups or the immune subgroups. The results showed that there were significant differences in the abundance of memory B cells and NK cells between the TMBhigh and TMBlow subgroups, while there were substantial differences in the abundance of 7 immune cell subgroups between the immunityhigh and immunitylow subgroups (Figure S1A and S1B).
Functional enrichment analysis of up-regulated differentially expressed genes (Up-DEGs).
Tumor mutation burden (TMB) or tumor immune infiltration have been demonstrated as a potential biomarker to predict responses to immunotherapy in several types of solid tumors[14–16]. Considering the close correlation between TMB, or the immune cell infiltration and the immune cell infiltration. We comprehensively compared the upregulated DEGs in the TMBhigh and immunityhigh group to looking for common differential genes. We reasoned that the genes enriched in both these two subgroups would be involved in important biological functions. We crossed the top 200 genes enriched in the TMBhigh group with the 427 genes significantly up-regulated in the immunityhigh group through the Venn diagram, and finally got 14 Up-DEGs (Fig. 3A). The expression level of 14 Up-DEGs in each case within the TMB group and the immune group were shown in Figure S2A. Next, we downloaded and analyzed the expression of 14 Up-DEGs in immune cells (including T cells and B cells) and ovarian cancer cells on the CCLE website (https://portals.broadinstitute.org/ccle/about). IGLL1, CXCL13, FCRLA, MS4A1, and PLA2G2D are expressed higher in immune cells than in OV cells (Figure S2B). Interestingly, we also found CCL11, IL27, and CXCL13 in the signature genes of the “chemokine receptor (CCR)” in Table S1. CXCL13 was also found in the signature genes of “follicular helper T cells (Tfh)”. PLA2G2D was found in the signature genes of “CD8+ T cells”. MS4A1 found in the signature genes of “tumor-infiltrating lymphocytes (TIL)”. Finally, we conducted unsupervised hierarchical clustering of KEGG pathways related to 14 Up-DEGs and found that most of these KEGG pathways were immune-related, such as “Cytokine cytokine receptor interaction”, “Chemokine signaling pathway”, as well as “Primary immunodeficiency” (Fig. 3B). Moreover, to predict the potential function of the DEGs, we performed the functional enrichment analysis of the 14 Up-DEGs. GO terms and KEGG pathway analysis also demonstrate that the Up-DEGs function in the immune and inflammatory response, cytokine activities, and T cell proliferation (Fig. 3C and 3D). These results indicated that the gene sets shared by both TMBhigh and immunehigh groups are involved in active immune functions.
Protein-protein interactions (PPI) of Up-DEGs
To understand the interrelationship between the 14 Up-DEGs, we used the Search Tool for the Retrieval of Interacting Gene (STRING) tool to establish a PPI network (Fig. 4A). By counting the number of connected nodes, it shows that CXCL13 6 interconnected proteins, which are the most key links with other members of the module in the network. Furthermore, CXCL13, several key genes associated with T cell proliferation, inflammation, and immune response were located at the center of the module, such as PLA2G2D and FCRLA (Fig. 4B). These results indicated that CXCL13, FCRLA, MS4A1, PLA2G2D, and ADAMDEC1 interact more with other molecules, suggesting that these genes are more likely to be key genes in regulating tumor restricting functions in OCs.
Identification of a gene set correlated with better survival.
Both TMBhigh group and immunityhigh group were correlated with better survival. However, the detection of either TMB or tumor immune profile takes much time and resources, which limits their prevalence in the clinic. Based on our previous analysis, we reasoned that there could be potential prognosis predicting biomarkers among the 14 shared Up-DEGs defined above that could be applied to predict the survival. To this end, we analyzed the potential role of Up-DEGs in the overall survival time of OC patients. Among 14 Up-DEGs, four genes, including CXCL13, FCRLA, MS4A1, and PLA2G2D were validated to be significantly relevant to good overall survival and prognosis outcomes (log-rank test, p < 0.05) (Fig. 5A-D). Next, we downloaded from the CPTAC database and analyzed the difference in protein expression of FCRLA and MS4A1 in OC patients with different stages. The results showed that the protein expression of FCRLA and MS4A1as the stage of tumors increase (Figure S3A, B). Furthermore, We analyzed the prognosis of CXCL13, FCRLA, MS4A1, and PLA2G2D on the Human Protein Atlas (HPA, https://www.proteinatlas.org) website. The results showed that the high expression of CXCL13, FCRLA, MS4A1, and PLA2G2D were correlated with a good prognosis of OC (Figure S3C-F). To validate whether the four genes have an equal prognostic value in other OC studies, we downloaded and analyzed the expression datasets for different cohort of the four genes from the GEO database, including GSE30161 (n = 58), GSE9891 (n = 285), GSE63885 (n = 101), GSE26712 (n = 195), GSE15622 (n = 69), GSE19829 (n = 70), GSE18520 (n = 63), GSE26193 (n = 107), GSE27651 (n = 49). Meta-analysis further indicated that OC cases with high expression of these four genes showed a significant survival benefit (Fig. 5E-H). These results indicated that the transcription levels of CXCL13, FCRLA, MS4A1, and PLA2G2D could be applied to predict the prognosis of OC patients.
CXCL13, FCRLA, MS4A1, and PLA2G2D correlated with immune status in OC.
Tumor purity is usually considered when obtaining immunotherapy expression markers from transcriptome data. Some studies have shown that tumor samples with low purity tend to have more immune cells and a higher mutation load. These patients with low purity responded were better to immunotherapy. Therefore, to further validate if the four signature genes can predict tumor immune infiltration, we used TIMER (Tumor Immune Estimation Resource, https://cistrome.shinyapps.io/timer/), an online database, to analyze the infiltration of different immune cells in tumor tissues based on available RNA-seq data. Our results demonstrated that CXCL13, FCRLA, MS4A1, and PLA2G2D have a negative correlation with tumor purity in OC (Fig. 6A-D). Furthermore, CXCL13 expression showed a very weak relationship with B cell infiltration level in OC (Figs. 6A), while FCRLA, MS4A1, PLA2G2D have no significant correlations with B cell infiltration level in OC (Fig. 6B-D). However, there were moderate to strong positive relationships between the expression levels of CXCL13, FCRLA, MS4A1, PLA2G2D, and infiltration level of CD4+ T cells, as well as have a prominent positive correlation between expression level and infiltration level of CD8+ T cells in OC, especially CXCL13, FCRLA, MS4A1, and PLA2G2D (Figs. 6A-D). More information can be found in Table S2. In conclusion, CXCL13, FCRLA, MS4A1, and PLA2G2D might have multiple and closely function in immune infiltration.
To identify the pathways that these four genes might be involved in gene set enrichment analysis (GSEA) was performed in the published TCGA ovarian cancer database (n = 304). The results indicated that all of four genes were related to immune cell signaling pathways. They were all associated with the B cell receptor, chemokine, cytokine-cytokine receptor interaction, natural killer cell-mediated cytotoxicity, and T cell receptor signaling pathways (Fig. 6E-H). Interestingly, CXCL13, FCRLA, and PLA2G2D were also associated with primary immunodeficiency (Fig. 6E-H). MS4A1 was associated with autoimmune thyroid disease (Fig. 6G). All of these findings revealed that CXCL13, FCRLA, MS4A1, and PLA2G2D’s expression was closely related to immune function. More information can be found in Table S3. These results indicated that CXCL13, FCRLA, MS4A1, and PLA2G2D were associated with immune pathways.
Altogether, all the results indicated that CXCL13, FCRLA, MS4A1, and PLA2G2D were associated with immune effector cell infiltrations and involved in immune pathways.
CXCL13, FCRLA, MS4A1, and PLA2G2D correlated with immunotherapy.
Immune infiltration and active inflammations are prerequisites of anti-tumor immune reaction and better response to immunotherapy. We have shown that CXCL13, FCRLA, PLA2G2D, and MS4A1 are up-regulated in the TMBhigh group and immunityhigh group, involved in immune functions, and are associated with better survival. We then asked if we could use the transcription levels of these four genes to predict the response to immunotherapy. Since the data of OC response to immunotherapy are not available, to solve this problem, we explored an online database—TIDE, then analyzed the data collected from the melanoma patients treated with anti-PD-1 or anti-CTLA4 alone from Gide et al. in the TIDE website. Previous reports showed that some tumors have a high level of infiltration by cytotoxic T cells[24, 25]. Then the CTL levels were used in the TIDE website to evaluate the tumor immune infiltration of these four Up-DEGs. As shown in Fig. 7A-D, our results showed that the CXCL13, FCRLA, MS4A1, and PLA2G2D were positively correlated with the level of CTL infiltration in the melanoma patients treated with anti-PD-L1 or anti-CTLA4 (Fig. 7A-D). Moreover, we also found that the high expression levels of CXCL13, FCRLA, MS4A1, and PLA2G2D were positively associated with the better response of melanoma patients treated with anti-PD1 and anti-CTLA4 (Fig. 7E-H). These results indicated that CXCL13, FCRLA, MS4A1, and PLA2G2D could be applied to predict responses to immunotherapies.