The TMB scores and the immunoactivity-related gene sets are conjoined with OC’s 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 in March 2019 (https://portal.gdc.cancer.gov/). We evaluated all OC cases with complete gene expression data and clinical information in TCGA. And from mutation annotation files, we presented the Word Cloud of gene symbols identified in OC which the mutation sites > 5 based on the WordArt (https://wordart.com). We found that TP53 is the most significantly mutated gene (Figure S1A). We then selected the top 25 genes to execute some higher-level studies that we employed the DISCOVER (Discrete Independence Statistic Controlling for Observations with Varying Event Rates) [25] algorithm to examine the association between co-occurrence and mutual exclusion among these genes and we found that the pairwise combinations that the highest significance of mutual exclusivity was TTN and LRP2 (p < 0.001) (Figure S1B). In order to better visualized the mutually exclusive genes pairs, we altered 420 (96.33%) of 436 samples and in a waterfall plot shown that the most frequently mutated genes were TP53 and TTN (p < 0.05, FDR < 0.01) (Figure S1C). Based on their somatic mutation data, TMB was distributed between 0.0263 to 6.5260 and the rank order of TMB across OC cases from highest to lowest is Grade 3 & Grade 4 > Grade 1 & Grade 2, that is, the higher the OC grade, the higher the TMB.
Based on their somatic mutation data, TMB was distributed between 0.0263 to 6.5260 and the rank order of TMB across OC cases from highest to lowest is Grade 3 & Grade 4 > Grade 1 & Grade 2, that is, the higher the OC grade, the higher the TMB (Figure 2A, n=417, except for the blank values, p=0.048). And then based on the ssGSEA scores of immune cell types, we divided the immunoactivity-related gene sets into three immune infiltration clusters by unsupervised clustering, and some of which were defined as high immunity clusters (n=193) and low immunity clusters (n=149) (Figure 2B). Besides, the significant differences in tumor purity were found in the microenvironment between the high immunity clusters and the low immunity clusters. In addition, we used tumor purity in the microenvironment to validate the grouping and found that the tumor purity in the high immunity group was significantly lower than that in the low immunity group. This result proves the reliability of our classification (Figure S1D). Moreover, from this trend, it was found that there was a significant correlation between the high immunity clusters and the expression of PD-L1 (Figure 2B). To figure out the potential interrelationship between the overall survival and TMB groups and/or immunity groups, we assessed the prognosis of OC samples. In Figure 2C, Kaplan-Meier survival curves showed that the overall survival time of patients with the high TMB group was longer than the low TMB group and this was statistically significant. (2.58 years vs. 1.92 years, p=0.0295 in log-rank test). In Figure 2D, the result showed that the overall survival time of the high immunity group was higher than that of the low immunity group and this was statistically significant. (2.50 years vs. 2.45 years, p=0.00082 in log-rank test). These results indicated that the higher the TMB, the better the prognosis and the higher the OC grade. Besides, the high immunity cluster was also associated with a good prognosis of OC.
KEGG pathway and GO term analysis of DEGs associated with OC.
First, we compared the relationship between the immunity group and the immune cell infiltration in OC. Then we compared the infiltration degree of 22 different immune cells subpopulation [26] between the low immunity group and the high immunity group of cancers. We found that in the correlation between immunoactivity-related gene sets and immune infiltration, most of these immune cells have a higher infiltration degree in the high immunity group than in the low immunity group of OC (Figure S2). Next, as to reveal the correlation between the TMB, the immunoactivity-related gene sets, and global gene expression profiles, we crossed the high TMB gene expression signature with the high immune gene expression signature by Venn diagrams and finally obtained 14 DEGs (Figure 3A). Heatmaps of 14 DEGs in the up-regulated group were shown in Figure 3B. In the next work, we conducted unsupervised hierarchical clustering of KEGG pathways related to Up-DEGs (Figure 3C) and it’s shown that these KEGG pathways were able to distinguish OC from normal controls with high specificity and sensitivity. And most of these KEGG pathways are immune-related,such as “Cytokine cytokine receptor interaction”, “Chemokine signaling pathway”, as well as “Primary immunodeficiency”. Moreover, in order to summarize the potential function of the DEGs, we performed functional enrichment analysis of the 14 up-regulated genes between high TMB and high immunity group. GO terms and KEGG pathway analysis identified including immune and inflammatory response, cytokine activities and T cell proliferation in the Up-DEGs group (Figure 3D and 3E). These results indicated that the gene sets at the intersection of Higher TMB gene expression signature and higher immune gene expression signature have active immune functions.
Protein-protein interactions among genes associated with prognosis in OC
In order to better make out the interactions between the identified DEGs, we acquired PPI networks using the STRING tool. The network was consists of Up-DEGs module (Figure 4), in the Up-DEGs module, we uploaded 14 DEGs to STRING to construct PPI networks and screened nine hub genes, it indicated that 9 nodes and 10 edges when minimum required interaction score was set 1.5 and disconnected nodes were hidden in the network in STRING website (Figure 4A). We showed a connection number with other members of the module for further analysis (Figure 4). In the Up-DEGs module (Figure 4B), CXCL13 involving 6 nodes, which has the most key link with other members of the module in the network. For the Up-DEGs module (Figure 4B), several key genes associated with T cell proliferation, inflammation and immune response were located at the center of the module, including PLA2G2D, CXCL13, and FCRLA. These results indicated that PLA2G2D, MS4A1, ADAMDEC1, FCRLA, CXCL13 interact more with other molecules, suggesting that these genes are more likely to be key genes.
Correlation of expression of individual DEGs in overall survival
For the sake of probe the potential effects of individual DEGs in overall survival time, Kaplan-Meier survival curves were generated from the TCGA database. Among a total of 14 Up-DEGs, four genes were validated to be significantly relevant to a good survival and prognosis outcomes, including PLA2G2D, CXCL13, FCRLA, and MS4A1 (log-rank test, p < 0.05) (Figure 5A, 5B, 5C, 5D). To valid the prognostic of these genes, we downloaded and analyzed gene expression data from the Group on Earth Observations (GEO). Meta-analysis depicting forest plots of DEGs further indicated that OC showed a significant survival benefit with high TMB gene expression signature of the four genes and were all independent predictors of OC (Figure 5E, 5F, 5G, 5H). These results indicated that CXCL13, FCRLA, PLA2G2D, and MS4A1 were validated as prognostic biomarkers for the immunotherapy of OC.
Correlation of four genes expression with immune infiltration level in OC and pathway analysis.
Another vital aspect of this research is the correlation between DEG expression and different levels of immune infiltration in OC. We used TIMER (Tumor Immune Estimation Resource, https://cistrome.shinyapps.io/timer/), an online database, to detect the infiltration of immune cells in tumor tissues using the RNA-seq. Our results demonstrated that CXCL13, MS4A1, and PLA2G2D have a negative correlation with tumor purity in OC, especially CXCL13, MS4A1, PLA2G2D (Figure 6A, 6C, 6D). CXCL13 expression showed a very weak relationship with B cell infiltration level in OC (Figures 6A) while FCRLA, MS4A1, PLA2G2D have no significant correlations with B cell infiltration level in OC (Figure 6B, 6C, 6D). Moreover, There were moderate to strong positive relationships between the expression levels of PLA2G2D, FCRLA, MS4A1, CXCL13, and infiltration level of CD4+ T cells, and also have a prominent positive correlation between expression level and infiltration level of CD8+ T cells in OC, especially CXCL13, MS4A1, PLA2G2D, FCRLA (Figures 6A, 6B, 6C, 6D). More information can be got from Table S1.
In the next work, we conducted gene set enrichment analysis (GSEA) to analyze the four genes and found that all of four genes were related to immunity, Both CXCL13, FCRLA, MS4A1, and PLA2G2D were associated with B cell receptor, chemokine, cytokine-cytokine receptor interaction, natural killer cell-mediated cytotoxicity and T cell receptor signaling pathway (Figure 6E, 6F, 6G, 6H). Besides these, both CXCL13, FCRLA, PLA2G2D were associated with primary immunodeficiency (Figure 6E, 6F, 6H). MS4A1 was associated with autoimmune thyroid disease (Figure 6G). All of these findings revealed that CXCL13, FCRLA, MS4A1 and PLA2G2D’s expression was closely related to immunity were enriched for hallmarks of OC. More information can be got from Table S2. These results indicated that CXCL13, FCRLA, PLA2G2D, and MS4A1 were validated as prognostic biomarkers for the immunotherapy of OC. These results indicate that CXCL13, FCRLA, PLA2G2D, and MS4A1 are positively correlated with immune cell infiltration, and these genes are enriched to be associated with immune pathways by GSEA enrichment method.
Prognostic evaluation of four genes in melanoma patients treated with immunocheckpoint inhibitors.
In addition, we used an online database, TIDE(Tumor Immune Dysfunction and Exclusion, http://tide.dfci.harvard.edu/query/), to predict patient response rates to immune checkpoint inhibitors. As shown in Figure 7, our results showed that among the melanoma patients treated with anti-PD-L1 or anti-CTLA4, The CXCL13, MS4A1, FCRLA, and PLA2G2D were positively correlated with the number of CTL infiltration, the two-sided t-test p values for correlations in CXCL13, MS4A1, FCRLA, PLA2G2D and mean are1.43×10-10, 8.12×10-2, 6.6×10-10 and 2.1×10-4, respectively. And the high expression of these genes was associated with a better prognosis of patients.
The last but not least, we researched the protein expression of CXCL13, FCRLA, and MS4A1 in OC and normal ovarian tissues based on the IHC (immunohistochemistry) in The Human Protein Atlas (https://www.proteinatlas.org). We found that the staining intensity of CXCL13 and MS4A1 was slightly different in normal ovarian tissues and OC tissues, the protein expression levels in cancer were slightly lower than in adjacent cancers, but the difference was not significant. Contrary to the mRNA level of FCRLA, the level of FCRLA protein was up-regulated in OC tissues compared to the normal tissues (Figure S3). PLA2G2D could not find the relevant protein map in the database for the time being, which may indicate that the gene is a molecule worth digging and exploring in the future, and it is expected to play a greater role as a prognostic marker.