A novel tumor mutational burden-based risk model predicts prognosis and correlates with immune infiltration in ovarian cancer

Abstract

Tumor mutational burden (TMB) has been reported to determine the response to immunotherapy, thus affecting the patient’s prognosis in many cancers. However, it is unclear whether TMB or TMB-related signature could be used as prognostic indicators for ovarian cancer (OC), as its potential association with immune infiltration remains poorly understood. Therefore, this study aimed to develop a novel TMB-related risk model (TMBrisk) to predict the prognosis of OC patients on the basis of exploring TMB-related genes, and to explore the potential association between TMB/TMBrisk and immune infiltration. The mutational landscape, TMB scores, and correlations between TMB and clinical characteristics and immune infiltration were investigated in The Cancer Genome Atlas (TCGA)-OV cohort. Differentially expressed gene (DEG) analyses and weighted gene co-expression network analysis (WGCNA) were performed to derive TMB-related genes. TMBrisk was constructed by Cox regression and further validated in Gene Expression Omnibus (GEO) datasets. The mRNA and protein expression levels and biological functions of TMBrisk hub genes were verified through Gene Expression Profiling Interactive Analysis (GEPIA), GSCA Lite, the Human Protein Atlas (HPA) database, and RT-qPCR. TMBrisk-related biological phenotypes were analyzed in function enrichment and tumor immune infiltration signature. Potential therapeutic regimens were inferred utilizing the Genomics of Drug Sensitivity in Cancer (GDSC) database and connectivity map (CMap). According to our results, higher TMB was associated with better survival and higher CD8+ T cell, regulatory T cell, and NK cell infiltration. TMBrisk was developed based on CBWD1, ST7L, RFX5-AS1, C3orf38, LRFN1, LEMD1, and HMGB1. High TMBrisk was identified as a poor factor for prognosis in TCGA and GEO datasets; the high-TMBrisk group comprised more higher-grade (G2 and G3) and advanced clinical stage (stage III/IV) tumors. Meanwhile, higher TMBrisk was associated with an immunosuppressive phenotype, with less infiltration of a majority of immunocytes and less expression of several genes of the human leukocyte antigen (HLA) family. Moreover, a nomogram containing TMBrisk showed a strong predictive ability demonstrated by time-dependent ROC analysis. Overall, this novel TMB-related risk model (TMBrisk) could predict prognosis, evaluate immune infiltration, and discover new therapeutic regimens in OC, which is very promising in clinical promotion.

Introduction

Ovarian cancer (OC) is the most lethal gynecological malignancy, with 5-year survival rates below 45%. According to the World Health Organization, it is estimated that the global incidence of OC will be 225,500 cases and 140,200 patients will succumb to this disease each year, making it the seventh most common cancer and the eighth most predominant cause of cancer-related death among women (1, 2). Owing to its insidious outset, rapid development and lack of obvious symptoms, most OC patients are detected at an advanced stage (3). Although the level of diagnosis and treatment has continuously improved, OC remains a serious threat to women’s lives and a salient public concern (4). At present, the International Federation of Gynecology and Obstetrics (FIGO) stage system (3) and some common serum biomarkers, such as carbohydrate antigen 125 (CA125) (5), human epididymis protein 4 (HE4) (6), and breast cancer gene1 (BRAC1) (7), are currently used as diagnostic tools for OC. However, these markers are not proven ideal to evaluate the prognosis and curative effect of each patient precisely (8). There is an obvious need to develop novel and reliable predictive tools for accurate individual evaluation, as well as preselection of suitable treatments.

Human tumors harbor a different number of somatic mutations collectively known as tumor mutational burden (TMB) (9), which is defined as the total number of somatic coding mutations, base substitutions, and insertion–deletion errors per million bases. In recent years, emerging evidence has suggested that TMB can determine the response to immunotherapy, thus affecting the patient’s prognosis (10, 11). One of the primary explanations is that tumor types with high TMB correspond to increased degree of tumor-specific neoantigen generation and presentation, which impacts the strength of immune response (9, 12). Fortunately, OC has been considered one of the most “immunogenic tumors” (13). Immunotherapies have attracted substantial attention and shown promising potential in OC therapy (14). For example, immune checkpoint blockade (ICB) therapy, which inhibits negative regulatory immune checkpoints through various immune checkpoint inhibitors (ICIs), has shown great promise for the treatment of OC (15–17). Due to the promising future of immunotherapy in the treatment of OC and the crucial role of TMB in predicting immunotherapy effect, evaluation of TMB status may be an effective way to predict the prognosis and therapeutic benefit for patients, individually.

However, the TMB cutpoints associated with improved curative effect and survival varied markedly between cancer types, and there may not be one universal definition of low/high TMB (11). In addition, determination of TMB is related to whole exome sequencing (WES) of the selected target panel, which increases the difficulty of TMB detection. Coordination and calibration are also required to achieve optimal utility and alignment of all platforms currently in international use (18). Although some studies have established prognostic models based on TMB, the process of screening TMB-related genes can be further optimized (19, 20). Therefore, in the current study, based on the data from The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and Genotype-Tissue Expression (GTEx) databases, we screened TMB-related genes by Weighted gene co-expression network analysis (WGCNA) and developed a simplified and practical TMB-related risk model (TMBrisk) to predict survival and immune infiltration in OC patients, which may have higher clinical value. Figure 1 presents an overview of approach in this study.

Figure 1

Materials and methods

Data acquisition and processing

We collected transcriptome profiles (Workflow type: HTseq-Counts), somatic mutation data (Data type: Masked Somatic Mutation; Workflow type:VarScan2) (21), and clinical data of OC patients from Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Due to the absence of matched normal samples in the TCGA database, gene expression data for normal ovarian tissues were downloaded from GTEx portal (http://www.gtexportal.org/home/). For transcriptome profiles, Entrez ID was transformed to the corresponding official gene name with Perl language, and genes with zero expression in more than 50% samples were removed. For mutation analysis, the somatic mutation data in Mutation Annotation Format (MAF) was analyzed and visualized using the “matfools” R package, which offered a number of analysis and visualization modules commonly used in cancer genomic research (22). For clinical data, corresponding information was extracted, including age, histological grade (G1, G2, and G3), FIGO stage (stage I–IV), cancer status (tumor free/with tumor), residual disease largest nodule (no macroscopic disease, 1–10 mm, 11–20 mm, and >20 mm), anatomic neoplasm subdivision (unilateral/bilateral), venous invasion (yes/no), lymphatic invasion (yes/no), Karnofsky performance score, and survival data. Moreover, we selected three validation datasets (GSE18520, GSE26193, and GSE63885) from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) and obtained their normalized microarray gene expression data and clinical data.

TMB calculation, prognostic analysis, clinical correlation analysis, and immune and stromal scores analysis

In our study, the somatic mutation information was extracted with a Perl script (23), and the TMB score of each sample was calculated through dividing the number of variants by exon length (38 million). For prognostic analysis and clinical correlation analysis, R was utilized to merge the patients’ TMB scores with corresponding clinical information. Kaplan–Meier (K-M) analysis was conducted to compare the difference in overall survival (OS), progression-free survival (PFS), and cancer-specific survival (CSS) with the log-rank test for statistical significance. Wilcoxon rank-sum test was employed to compare clinical traits between two groups, while Kruskal–Wallis test was used among multiple groups. Immune scores, stromal scores, and tumor purity were calculated using the ESTIMATE algorithm (24), which was provided in the “estimate” R package. The correlation between TMB and immune/stromal scores was evaluated using the Pearson method.

Analysis of immune infiltration for TMB

CIBERSORT, a deconvolution algorithm, can quantify the abundance of any of 22 types of immune cells based on bulk transcriptome profiles (25). We normalized the transcriptome data of OC patients with the “limma” R package (26) and uploaded the prepared data to the CIBERSORT algorithm (R script v1.03) to evaluate the distribution of different immune cells in each sample. The Tumor Immune Estimation Resource (TIMER) database (http://timer.comp-genomics.org/) (27) includes 10,897 samples from 32 cancer types in TCGA and pre-calculates the infiltration levels of six types of immune cells. In the current study, to explore the prognostic value of immune cells, we used the “Survival” module to output tumor-infiltrating immune cell (TIIC)-related K-M plots (27). Furthermore, the “Somatic Copy Number Alterations (SCNA)” module was used to assess the relationship between different SCNAs of genes and immune infiltration. The SCNAs were defined by GISTIC 2.0, including deep deletion, arm-level deletion, diploid/normal, arm-level gain, and high amplification (28).

Differentially expressed gene screening

We utilized the “limma” R package (26) to identify DEGs between OC and normal samples with the criteria of fold change (FC) > 2 and false discovery rate (FDR) < 0.05. Since the tumor sample size (n = 379) was much larger than that of the normal sample (n = 80), we implemented a subset-based strategy to balance the sample size in order to make the screening result more objective and accurate (29). In detail, five random samplings were performed. For each random sampling, a subset of 100 patients were selected from the tumor group without repetition. DEGs were identified by analyzing differences in transcriptome profiles between each tumor subset (n = 100) and normal samples (n = 80). The intersection of DEGs obtained from five independent analyses was selected as the final DEGs between the two groups. Volcano plots of five analyses were drawn by the “ggplot2” R package. Venn diagrams were plotted by the “VennDiagram” R package to exhibit the common DEGs.

Weighted gene co-expression network construction

In the current study, a scale-free gene co-expression network was constructed using the “WGCNA” in R software to recognize co-expressed gene modules closely related to TMB (30). First, to ensure the reliability of network construction, sample clustering was carried out to distinguish and remove outlier samples. By calculating the corresponding scale independence (R²) and mean connectivity when the soft threshold ranged from 1 to 20, the optimal value of soft threshold was determined. Based on Pearson’s coefficients of pairwise gene correlations and the chosen soft threshold value, a weighted adjacency matrix was constructed and subsequently transformed into a topological overlap matrix (TOM) (31). Next, the corresponding dissimilarity matrix (1-TOM) was used to establish a hierarchically gene clustering tree. Genes in the same branch of the dendrogram were highly correlated and clustered into the same co-expression module using the dynamic tree cut method, with a minimum module size cutoff of 50 and a deepSplit value of 2. Highly similar modules were screened and merged together with the height cutoff of 0.25. In addition, the relationship between each module and clinical parameters was analyzed by Pearson correlation analysis and visualized by a heatmap. Gene modules of |correlation coefficient| > 0.2 and p value < 0.05 were considered as strong TMB-correlated modules.

Functional enrichment analysis for co-expressed gene modules

To explore the biological functions and pathways of modules, we obtained the Entrez ID of each gene via the “org.Hs.eg.db” R package and conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the “clusterProfiler,” “enrichplot,” “ggplot2”, and “Goplot” R packages (32).

Identification of TMB-related immune genes

We downloaded a list of 2,483 immune-related genes from the Immunology Database and Analysis Portal (Immport) (https://www.immport.org/shared/genelists/) (33, 34), and selected genes that overlapped with DEGs in the green module as TMB-related immune genes. The result was visualized using the “VennDiagram” R package. Then, batch survival analysis was performed via a “for cycle” R script to screen candidate genes associated with survival outcomes for further investigation.

Development of TMB-related risk model

As a preprocessing step, expression values of genes were log2 transformed after adding a pseudo-count of 1. Univariate Cox regression analyses, least absolute shrinkage and selector operation (LASSO) regression analyses, multivariate Cox regression analyses, and Akaike information criterion (AIC)-based stepwise Cox regression analyses were successively used to screen prognostic genes for constructing the risk model. LASSO is a regularization and descending dimension method that can be used in conjunction with Cox models for biomarker screening (35). The model with minimal AIC value was determined as the final model, established by multiplying mRNA levels of each gene by respective multivariate Cox regression coefficient.

Validation of expression patterns and identification of signaling pathways of hub genes

The mRNA expression patterns of the TMBrisk hub genes were verified by Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) (36). The protein expression of the hub genes between OC and normal tissues was determined using immunohistochemistry (IHC) from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/), which is a valuable database providing extensive transcriptome and proteomic data for specific human tissues and cells. GSCA Lite database (http://bioinfo.life.hust.edu.cn/web/GSCA Lite/) (37), a web-based platform for Gene Set Cancer Analysis, was employed for statistics of deletion/amplification of hetero/homozygous copy number variation (CNV) and identification of OC-related signaling pathways of hub genes. Due to the lack of RFX5-AS1 in GEPIA, HPA, and GSCA Lite, RFX5 was used as a substitute for analysis.

Development of the nomogram

We performed the nomogram with the independent prognostic factors screened by multivariate Cox analysis via the “rms” R package. To evaluate the prediction accuracy of the nomogram, the “survival” R package was used to calculate Harrell’s Concordance index (C-index) to quantify its discrimination performance. The calibration curves of survival probability for different years were plotted using the Hosmer–Lemeshow test.

Functional enrichment analysis for TMBrisk

In order to explore TMBrisk-related pathways without the restriction of DEGs, gene set enrichment analysis (GSEA) was implemented with GSEA software 4.2.2 (https://www.gsea-msigdb.org/gsea/index.jsp) (38). We used TMBrisk as the phenotype; “c5.go.mf.v7.5.1.symbols.gmt”, “c5.go.bp.v7.5.1.symbols.gmt”, and “c2.cp.kegg.v7.5.1.symbols.gmt” as the reference gene set, which were obtained from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb/) (39). The significant pathways were defined as those whose |normalized enriched score (NES)| > 1 and FDR < 0.05.

Analysis of tumor immune signatures for TMBrisk

On the one hand, we evaluated the expression of the human leukocyte antigen (HLA) gene family and the immune checkpoints (40, 41). On the other hand, the levels of infiltrating immune cells and stromal cells were calculated by seven algorithms ]TIMER (27), CIBERSORT (25), xCell (42), CIBERSORT-ABS (43), QUANTISEQ (44), MCPCOUNTER (45), and EPIC (46)]. The results are available on the TIMER database.

Prediction of treatment sensitivity and small molecule drugs

Chemotherapeutic response in OC patients was assessed utilizing the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org) (47). The 50% inhibiting concentration (IC50) value of the 138 drugs in GDSC was inferred using the pRRophetic algorithm. Possible small-molecule drugs for OC were forecasted using the Connectivity map (CMap) (https://www.broadinstitute.org/connectivity-map-cmap) (48), which was premised on DEGs between low- and high-TMBrisk groups with |FC| > 2 and FDR < 0.05.

Cell lines and cell culture

OC cells A2780, SKOV3, and OVCAR3, and normal ovarian surface epithelium cells IOSE80 were purchased from the Cell Storage Center of Wuhan University (Wuhan, China). A2780, SKOV3, and OVCAR3 cells belong to the human epithelial OC cell line, and IOSE80 cells as normal control were also from ovarian surface epithelium. A2780, IOSE80, and OVCAR3 were grown in RPMI-1640 medium (Hyclone) supplemented with 10%, 10%, and 20% fetal bovine serum, respectively. SKOV3 were grown in McCoy’s 5A medium (Biological Industries) supplemented with 10% fetal bovine serum. Cells were incubated with 5% CO2 atmosphere at 37°C.

Quantification of gene expression by real-time quantitative PCR

Total mRNA was purified using TRIzol reagent (TaKaRa) according to the manufacturer’s instructions. To generate cDNA, 2 μg of RNA was subjected to reverse transcription using the Hifair®II1st Strand cDNA Synthesis SuperMix (YEASEN). Real-time quantitative polymerase chain reaction (RT-qPCR) was performed with Hieff® qPCR SYBR Green Master Mix (YEASEN), using a CFX Connect Real-Time Cycler (Bio-Rad). GAPDH was set as an internal control for gene quantification. Technical and biological replicates of each gene were performed at least three times during RT-qPCR analysis. Supplementary Table 1 contains the RNA molecules evaluated on cell lines and their corresponding primers.

Survival and other statistical analysis

K-M survival analyses were tested by log rank. Receiver operating characteristic (ROC) curves were plotted by the “survivalROC” R package to evaluate the accuracy of TMBrisk and nomogram. Decision curve analysis (DCA) was conducted using “ggDCA” R package to assess the clinical outcomes of different decision strategies (49). All statistical analyses were completed by R software (version 3.6.3). R packages used in our study are mentioned above. Benjamini–Hochberg for multiple testing and false discovery rate (FDR) were used to correct the p-value. The Spearman method was used to calculate the correlation coefficient and p-value in the correlation analysis of TMBrisk and other variables. Wilcoxon rank-sum test and Kruskal–Wallis test were used for subgroup differential analyses. Two-tailed tests and p values < 0.05 for significance were used.

Results

Landscape of mutation profiles in OC samples

Mutation data of 436 OC patients were downloaded from the TCGA database and the “maftools” R package was implemented to analyze and visualize the landscape of mutation profiles. The different genetic mutations in each sample were shown in the waterfall plot, with various color annotations representing different mutation types (Figure 2A). The mutations were further classified into different categories. By comparison, missense mutations accounted for the largest fraction (Figure 2B), single-nucleotide polymorphisms (SNPs) were more frequent than insertions or deletions (Figure 2C), and C>T represented the most common type of single-nucleotide variant (SNV) (Figure 2D). Moreover, counting the number of altered bases per patient, we found that the median and maximum number of mutations were 49.5 and 724, respectively (Figure 2E). For each variant classification, the number of occurrences is shown in the box plot (Figure 2F). Furthermore, taking into consideration the total number of mutations and counting the multiple hits alone, we recalculated the top 10 mutated genes (Figure 2G), which were slightly different from the previous ones (Figure 2A). Finally, the co-occurrence and exclusive associations across mutated genes were investigated. Significant co-occurrences were observed among mutations of COL6A3 and DNAH3, APOB and MUC16, etc. (Supplementary Figure 1).

Figure 2

Calculation of TMB and its correlation with prognosis and clinical traits

A total of 436 OC patients in the TCGA database were included in our study to calculate TMB scores. The range of TMB spanned from 0.026 to 29.289 (Supplementary Table 2). For an improved accuracy of the analysis, one statistical outlier (TMB = 29.289) was removed. To assess the relationship between TMB and prognosis, patients were divided into low- and high-TMB groups with the median score as the cutoff value. According to the K-M analysis, the high-TMB group showed a trend toward better OS, PFS, and CSS than the low-TMB group (Figures 3A–C). For clinical characteristics, TMB had a significant distinct distribution in the patients presented with bilateral/unilateral tumors (Figure 3D) and patients with tumor or tumor-free patients (Figure 3E). Significant correlations also existed between TMB and histological grades (G1 vs. G2, p < 0.01; G1 vs. G3, p < 0.01) (Figure 3F), FIGO stages (SI vs. SIII, p = 0.018; SI vs. SIV, p = 0.022; SII vs. SIII, p = 0.048; SII; vs. SIV, p = 0.038) (Figure 3G), and residual disease (NO vs. 1–10 mm, p < 0.01; NO vs. >20 mm, p = 0.044) (Figure 3H), which indicated that higher TMB was associated with higher histological grades, lower clinical stages, and smaller residual tumor sizes. However, no significant association was found between TMB and age, Karnofsky performance score, or venous/lymphatic invasion (Figures 3I–L). Furthermore, TMB significantly correlated with immune scores and tumor purity inferred by ESTIMATE algorithms with positive and negative dependencies, respectively (Figures 3M, N), but there was no significant association with stromal scores (Figure 3O).

Figure 3

Evaluation of immune infiltration in the low- and high-TMB groups

In order to better understand the potential effects of TMB on immune activity and response, we compared immune infiltration between the low- and high-TMB groups. Based on the CIBERSORT algorithm, we estimated the proportion of 22 types of immune cells in each sample (Supplementary Table 3). After the filtration of samples with p < 0.05, the distribution of infiltrating immune cells in 192 samples was shown in a barplot (Supplementary Figure 2). Then, we analyzed the difference in the proportion of each type of immune cell between two TMB groups. Wilcoxon rank-sum test revealed that compared with the low-TMB group, CD8+ T cells (p = 0.043), regulatory T cells (Tregs) (p = 0.012), and activated natural killer (NK) cells (p = 0.011) exhibited lower infiltrating levels in the high-TMB group. However, samples in the high-TMB group had a significant increase in the fraction of macrophages M1 (p = 0.02) and macrophages M2 (p = 0.004) (Figure 4A).

Reportedly, immune infiltration could affect the prognosis of patients. Therefore, we performed K-M analysis on the six types of immune cells and found that high infiltration levels of CD4+ T cell and dendritic cell were associated with a favorable prognosis, while the situation was reversed in macrophage (Figure 4B).

Figure 4

Construction of co-expression modules by WGCNA

The high gene mutation rate in OC prompted us to explore the key genes that significantly influenced the overall mutational load. Through a subset-based method, differential gene expression analyses were repeated five times (Supplementary Figure 3A). The Venn diagrams showed the intersection of DEGs obtained from five analyses (Figure 5A). Between the normal and OC samples, a total of 5,120 DEGs were identified, including 2,337 upregulated genes and 2,783 downregulated genes. These DEGs were used for WGCNA. In WGCNA, no outlier samples were removed according to the sample clustering result, and the trait heatmap showed the distribution of samples according to their corresponding clinical characteristics (Supplementary Figure 3B). The optimal soft threshold value of 5 was chosen to guarantee the scale-free distribution (Supplementary Figure 3C). According to WGCNA results, 11 co-expression modules were recognized (except the gray module containing genes that were not co-expressed) (Figure 5B), and the hierarchical clustering plot revealed module eigengenes (MEs) (Figure 5C). The green module had the strongest correlation with TMB (|Cor| > 0.2, p < 0.001) (Figure 5D), which was also demonstrated by the scatter plot showing a significant correlation between module membership (MM) and gene significance (GS) (Supplementary Figure 3D).

Figure 5

GO and KEGG enrichment analyses were conducted to uncover the potential biological functions of genes in the green module. According to GO enrichment results, the genes played important roles in cell proliferation-related pathways, such as organelle fission, nuclear division, and chromosome segregation (Figure 5E, Supplementary Table 4). KEGG analysis revealed 15 significant signaling pathways, including cell cycle, oocyte meiosis, and cellular senescence (Figure 5F, Supplementary Table 5). These results implied that mutations in these genes may be potentially critical for tumorigenesis.

Furthermore, we sought to gain a deeper insight into TMB-related immune genes. Four immune genes (HDGF, NRAS, PSMD6, and NR1D1), which were significantly correlated with survival, were identified between the intersection of the immune gene set obtained from Immport and the green module (Supplementary Figure 4A). K-M curves showed that the high expression level of HDGF, NRAS, and PSMD6 was associated with an improved prognosis, while the opposite happened for NR1D1 (Supplementary Figure 4B). RT-qPCR indicated that HDGF and NRAS were overexpressed in OC, while NR1D1 and PSMD6 were underexpressed (Supplementary Figure 4C). We also used the HPA database to explore the differences in protein levels of the four immune genes between normal and tumor tissues (Supplementary Figure 4D). More importantly, the impact of copy number alteration of the four genes on immune infiltration in OC patients was evaluated. In brief, compared to patients with diploid/normal expression of HDGF, NR1D1, and NRAS, patients with high copy number amplification of HDGF or copy number deletion of NR1D1/NRAS had a lower level of immune infiltration, including CD8+ T cells, CD4+ T cells, and neutrophils. On the other hand, decreased CNV of PSMD6 upregulated immune infiltration, especially of macrophages (Supplementary Figure 5).

Establishment and validation of the TMBrisk model

Given the significant correlation between TMB and prognosis and the limitations of TMB, we considered to establish a novel risk model based on TMB-related DEGs. The 496 genes in the green module were first inputted into univariate Cox regression analysis, which identified 17 genes with p < 0.01 (Figure 6A); then, they were further screened by LASSO Cox regression (Figure 6A). The remaining 14 genes were finally enrolled in multivariate Cox analysis with stepwise regression. The model consisting of CBWD1, ST7L, RFX5-AS1, C3orf38, LRFN1, LEMD1, and HMGB3 had a minimal AIC value (AIC = 2,277.02) (Figure 6A). Therefore, the TMB-related risk model (TMBrisk) was established as follows: TMBrisk = (−0.3585) * CBWD1 + (−0.5693) * ST7L + (−0.1593) * RFX5-AS1 + (−0.2473) * C3orf38 + (0.1461) * LRFN1 + (−0.1178) * LEMD1 + (−0.2224) * HMGB3.

Figure 6

To elucidate the underlying role of TMBrisk in the development of OC, we investigated the distribution of TMBrisk in patients with different clinical characteristics. We found that tumor-bearing patients and patients with higher histological grades, advanced FIGO stages, venous invasion, or lymphatic invasion had a higher TMBrisk score (Figure 6B). Then, univariate and multivariate Cox regression analyses were conducted utilizing TCGA and three GEO datasets (GSE18520, GSE26193, and GSE63885) to examine the significance of the impact of TMBrisk on prognosis. The outcomes illustrated that in the four datasets, a higher TMBrisk score indicated poorer OS, and TMBrisk was identified as an independent prognostic factor after adjusting other clinical factors (Figure 6C). ROC curves of OS prediction in the four datasets showed that TMBrisk had a strong predictive value with the area under the curve (AUC) > 0.7 (Figure 6D, Supplementary Figures 6A–C). For 2-year OS prediction in TCGA, the AUC of TMBrisk was also significantly higher than that of other independent prognostic predictors (Figure 6E). Meanwhile, the DCA for 10-year OS prediction showed that TMBrisk had the highest net benefit (Figure 6F). Subsequently, patients in TCGA and three validation datasets were divided into low- and high-TMBrisk groups based on the median score, respectively. K-M curves showed that patients in the low-TMBrisk group exhibited improved OS, PFS, and CSS than those in the corresponding high-TMBrisk group (Figure 6G, Supplementary Figures 6D–H). The result of stratified analysis revealed significant differences in OS between low- and high-TMBrisk for subgroups with different age, FIGO stage, histological grade, residual disease, venous invasion, lymphatic invasion, and anatomic subdivision (Supplementary Figure 7). Taken together, the above results indicated that the risk model had excellent predictive ability, robust performance and potential clinical application value.

Validation of the expression of TMBrisk hub genes

To verify the reliability of the seven hub genes that made up the TMBrisk model, we validated the expression patterns of these genes through databases and RT-qPCR. According to GEPIA, mRNA expression levels of LRFN1, LEMD1, and HMGB3 were significantly higher in OC than in normal samples. The mRNA expression levels of CBWD1, ST7L, RFX5, and C3orf38 were decreased in OC, although not significantly (Figure 7A). RT-qPCR demonstrated that OC cells had higher mRNA expression of LRFN1, LEMD1, and HMGB3, and lower mRNA expression of the other five genes compared with control cells, which was consistent with the results of GEPIA (Figure 7B). In addition, according to HPA, the protein expression levels of the five genes (data for LRFN1 and LEMD1 were lacking and therefore not presented) were significantly higher in tumor samples than in normal samples (Figure 7G). Finally, we analyzed the signaling pathways and CNV of these genes in GSCA Lite. ST7L was found to be activated and inhibited in DNA damage response pathway and EMT pathway of OC, respectively. In the cell cycle pathway, RFX5, CBWD1, and HMGB3 were activated, while LEMD1 was inhibited. As for hormone-related pathways, ST7L was activated in the hormone AR/ER pathways, and HMGB3 was inhibited in the hormone ER pathway (Figure 7C). These pathways played vital roles in oncogenesis, suggesting that the hub genes may participate in OC progression. Furthermore, HMGB3 and LRFN1 had a high proportion of both heterozygous amplification and heterozygous deletion, suggesting that HMGB3 and LRFN1 may be the characteristic genes of OC. RFX5 and LEMD1 were characterized by high heterozygous amplification, while CBWD1 harbored more heterozygous deletion (Figures 7D–F).

Figure 7

Construction and validation of nomogram for survival prediction

Multivariate Cox regression analysis performed in the TCGA dataset screened the TMBrisk, age, residual disease, venous invasion, and lymphatic invasion as independent prognostic factors (Figure 6D, Table 1). To strengthen predictive power, we constructed a nomogram based on the above prognostic factors (Figure 8A). The calibration curves for the survival possibility at 1, 3, and 5 years exhibited accurate prediction ability of the nomogram (Figure 8B). Time-dependent AUC suggested that in comparison to the TMBrisk, the nomogram further enhanced the prediction ability (Figure 8C). DCA showed that the nomogram and TMBrisk provided better net benefit over different threshold probability ranges, respectively, suggesting that both might have potential applications in different contexts (Figure 8D). Taken together, the above results indicated that the nomogram had enhanced prediction efficiency, as well as potential clinical application value.

Figure 8

Association of TMBrisk with OC immune signature

Given that TMB has become an evolving biomarker in the field of immuno-oncology9, further exploration of the association of TMBrisk with OC immune signature is warranted. First, to determine the underlying mechanisms leading to the different outcomes between the low- and high-TMBrisk groups, we performed GSEA with annotations of GO (Supplementary Tables 6, 7) and KEGG gene sets (Supplementary Table 8). Figure 9A showed that cell cycle, DNA repair, purine and pyrimidine metabolism, and biosynthesis of unsaturated fatty acids were significantly enriched in the high-TMB group, while those related to immune cell migration and immune response, such as NOTCH signaling pathway, integrin binding, laminin binding, endothelial cell chemotaxis, and phosphatidylinositol signals, were enriched in the low-TMBrisk group. In addition, immune scores were significantly lower in the high-TMBrisk group, while stromal scores and tumor purity were higher (Figure 9B).

Figure 9

We also evaluated gene expression of the 48 immune checkpoints and 24 HLA family genes between the low- and high-TMBrisk groups. Our analysis showed that three immune checkpoints and eight HLA family genes were significantly modulated in the high-TMBrisk group (Figure 9C). The distribution of infiltrating immune cells between the low- and high-TMBrisk groups inferred by seven databases demonstrated that a majority of immunocytes decreased in the high-TMBrisk group, such as macrophages, NK cells, CD4+ T cells, and endothelial cells (Figure 10). However, neutrophils, and monocytes infiltrated more in the high-TMBrisk group (Figure 10).

Figure 10

TMBrisk predicts new therapeutic regimens

To find new drugs that could be incorporated into OC treatment, we estimated the IC50 value of 138 drugs in TCGA-OV patients. We found that patients in the low-TMBrisk group tended to be more sensitive to docetaxel, mitomycin C, BMS.708163, etc., while patients in the high-TMB group tended to be more sensitive to BIBW2992, camptothecin, metformin, etc. (Supplementary Figure 8A). Digging deeper, the DEGs between the low- and high-TMBrisk groups were identified, and subsequently imported into CMap. In total, 20 negatively related molecular agents (enrichment < 0 and p < 0.05) for anti-OC were identified, and 21 mechanisms of action (MOA) were shared among them (Supplementary Figure 8B). We found that a majority of drug-target MOA belong to the inhibition class. For example, indirubin and purvalanol-a shared the mechanism of CDK inhibitor; olaparib and veliparib shared the mechanism of PARP inhibitor. In summary, our study screened drugs targeting the TMBrisk-related genes. These small molecules might potentially inhibit the occurrence of OC and need to be further studied.

Discussion

Although immunotherapy is gaining prominence in the treatment of cancers, including OC, only a small percentage of patients respond to it. TMB is emerging as a potential biomarker of immune response and is linked with prognosis (18, 50). Its predictive performance has been verified in a variety of cancers, such as lung cancer (51), head and neck cancer (52), and melanoma (53). With regard to OC, the effect of TMB as a potential predictor and its relationship with immunotherapy responsiveness and immune infiltration still need to be further explored. In the current study, we evaluated the correlation between TMB and OC clinical traits as well as immune infiltration, screened TMB-related DEGs (including immune genes), developed a risk scoring system (TMBrisk) based on the TMB-related gene module and validated its performance externally, explored TMBrisk in respect of its immune signature, and finally predicted potential therapeutic regimens for OC.

By analyzing the somatic mutation profiles, we identified the two most frequently mutated genes in OC: TP53 (90%) and TTN (21%). As a tumor suppressor, TP53 is sufficient to regulate multiple cell cycle control networks, the loss of which will lead to dramatic gene expression changes (54). Consistent with our study, Costa et al. discovered that OC was characterized by great genomic instability with universal TP53 mutations (55). TTN, whose mutations are frequently identified in solid tumors, is associated with increased TMB, and patients with mutated TTN have improved outcomes in response to ICBs (56). We also found that SNP, which could explain a significant proportion of the heritable risk of cancer (57), is the main variant type in OC. The relationship between some types of SNPs and the occurrence, progression, and treatment of OC has been elucidated, such as the ALDH2∗2 polymorphism (rs671) (58), the most common SNP in Asia, but further study is still necessary.

Consistent with previous research of cutaneous melanoma (59), we discovered that OC patients in the high-TMB group had superior survival outcomes, and a higher TMB correlated with lower clinical stages. Interestingly, these findings were different from those in head and neck squamous cell carcinoma (52) and clear cell renal cell carcinoma (60), indicating that TMB may have distinct prognostic value and stage correlation in different tumor types. Thus, a pan-cancer analysis may be required to elucidate its role comprehensively. Furthermore, our research founded a positive correlation between TMB and immune score and a negative correlation between TMB and tumor purity, suggesting that patients with higher TMB presented more active immune activity, which may be because higher TMB results in more neo-antigens, increasing chances for T-cell recognition (61).

The relationship between TMB and immune infiltration was further explored. We demonstrated that Tregs and CD8+ T cells were downregulated, while macrophages M1 and M2 were upregulated in the high-TMB group. Tregs maintain immune homeostasis via suppressing excessive immune responses, but Tregs also infiltrate in the tumor microenvironment (TME) and inhibit antitumor immune activities (62). It is believed that Tregs are relevant to poor prognosis (63). As a result, some novel immunotherapies that deplete tumor-infiltrating Tregs can enhance antitumor effects (64, 65). These results partly explain why the high-TMB group, which possessed lower proportions of Tregs infiltration, had better OS and PFS. However, in OC, some studies have shown that higher CD8+ T-cell infiltration correlated with improved clinical outcomes (66, 67). Therefore, the mechanisms by which TMB and immune cell phenotypes influence the prognosis of OC patients may require more evidence and discussion. The M1/M2 macrophage paradigm plays a key role in tumor progression: M1 as pro-inflammation and antitumoral, and M2 as anti-inflammation, immunosuppression, and pro-tumoral (68–70). Thus, the M1/M2 ratio is critical for TME homeostasis and immune function. Studies have demonstrated that a high M1/M2 ratio status is associated with favorable prognosis for most solid tumors, including OC (69, 71, 72). In the context of increased M1 and M2 in the high-TMB group, it is necessary to further compare the M1/M2 ratio of the two groups in future studies, so as to better clarify the relationship between TMB and immune function of OC patients.

WGCNA harbors predominant advantages over other bioinformatics methods. Its clustering results (co-expression gene modules) have high biological significance and reliability. More importantly, WGCNA can quantify the correlation between co-expression modules and clinical features, so it can be used to search for gene modules of interest (73). Various studies (74–76) have applied WGCNA to cancer research and demonstrated its effectiveness. In this study, we identified the TMB-related gene module by WGCNA. The functional enrichment analysis revealed that genes in this module were mainly associated with cell cycle-related pathways. Cell cycle dysregulation is a hallmark of tumor progression, and the process of immune response is also dependent on the cell cycle (77). Recent studies have shown that cell cycle activity of both cancer and immune cells in TME could modulate immune function (78, 79).

From the green module, we screened four TMB-related immune genes significantly associated with the prognosis of OC patients. HDGF plays a vital role in the transformation, apoptosis, angiogenesis, and metastasis of cancer cells (80). Recent studies in OC showed that HDGF was passively released by necrotic and apoptotic cells; extracellular HDGF acted as a messenger of cellular condition and further enhanced cellular migration (81). In addition, higher HDGF expression was closely associated to poorer OS (82). NRAS, along with KRAS and HRAS, constitutes the RAS family, which is the most frequently mutated gene family in cancer (83). Similarly, OC has a high incidence of somatic mutations in MAPK pathway genes, the most common of which include NRAS, and MAPK inhibitors have shown efficacy for the treatment of NRAS-mutant OC (84). NR1D1 was reported to be a tumor suppressor in OC, which retarded the growth of cancer cells by abrogating the JAK/STAT3 signaling pathway (85). Other studies have shown that NR1D1 can also interact with PARP1 and BRCA1 to inhibit their associated DNA repair pathways, leading to the lethality of OC cells (86). The protein encoded by PSMD6 is a subunit of the 26S proteasome, which is co-located with DNA damage foci and involved in DNA damage response (DDR). The deficiency of PSMD6 delays DNA repair and leads to a decline in cell survival, which may offer new therapeutic approaches for cancer (87). Few studies on PSMD6 in OC suggest that PSMD6 may be a potential research target.

Nevertheless, TMB is an imperfect biomarker, partly because (1) it is costly to get bulk gene expression to calculate TMB scores, and results from different platforms require complex coordination and calibration (18), and (2) high TMB cannot be used as a biomarker in all solid cancer types and further tumor-type specific studies are warranted (88). Therefore, we constructed a prognostic model using only seven genes selected from the green module. Among the seven genes, except C3orf38, which was associated with a high TMBrisk, the remaining six (ST7L, LRFN1, LEMD1, HMGB3, CBWD1, and RFX5-AS1) genes were protective factors. These genes represented positive or negative regulation of TMB and immune activity. For instance, one study showed that C3orf38 was contained in a candidate tumor suppressor gene (TSG) locus of OC and may serve as a TSG (89). Similar to the ST7 tumor suppressor gene, ST7L can inhibit the proliferation, migration, and invasion of cancer cells (90). Yang et al. proved that the downregulated ST7L in OC activated the WNT/MAPK pathway, thereby promoting tumorigenic activity (91). As for LRFN1, a number of LRFN1-based prognostic models have been shown to be effective in a variety of cancers, including prostate cancer (92), kidney renal clear cell carcinoma (93), and OC (94). LEMD1 belongs to the CTA family, which is a useful target of immunotherapy (95). Guo et al. have demonstrated that LEMD1 antisense RNA 1 (LEMD1-AS1) can suppress OC progression and can be used as a biomarker to predict survival (96). HMGB proteins are closely related to DNA damage repair, among which HMGB3 is often overexpressed up to 20 times in cancer cells (97). Mukherjee et al. indicated that targeted elimination of HMGB3 reduced cisplatin resistance in OC cells, increasing tumor cell sensitivity to chemotherapy (98). However, the remaining hub genes, CBWD1 and RFX5-AS1, are poorly investigated in OC, which could be topics for further study.

The advantage of this study is that it fully evaluates the robustness of the TMBrisk model, as it has been validated on three independent datasets. ROC curves and DCA showed that TMBrisk had good accuracy in predicting OS, better than other independent predictors. Currently, FIGO staging is the most widely used system for predicting the malignancy and progression of OC. However, the weakness of this system is that it mainly considers distant metastasis, lymph node invasion, tumor size and location, ignoring the heterogeneity of other clinical features (3, 99). In this study, based on TMBrisk, by incorporating other independent prognostic factors, we developed a nomogram with higher net benefit and better predictive accuracy, which can complement the FIGO system to exert a more profound clinical value.

Regarding the relationship between TMBrisk and immune status of OC patients, our study found that high TMBrisk tended to predict an immune-suppressed status. Lower immune activities were revealed in the high-TMBrisk group, including reduced immune cell infiltration and downregulation of some HLA family genes. For example, NK cells and CD4+ T cells showed low infiltration in the high-TMB group. NK cells have unique properties, in that they can swiftly kill tumor cells and enhance antibody as well as T-cell response, which support its role as anticancer agents (100). In OC, boosting NK cell expansion and functionality has emerged as an attractive therapeutic approach and has been extensively investigated (101, 102). CD4+ T cells play a critical role in the immune system and make great contribution to the antitumor response (103). It is reported that OC patients with dense infiltration of CD4+ T cells experience favorable OS and PFS (104). In addition, three immune checkpoint genes (VSIR, TNFRSF25, and CD160) were upregulated in the high TMB group, suggesting that these immune checkpoints may be potential therapeutic targets for OC similar to PD-1/PD-L1.

Finally, we identified some untraditional antitumor compounds, such as BMS.708163, metformin, and BIBW2992, with potential advantages for patients who might not gain improvement from conventional drugs. The efficacy of these drugs for OC is expected for further investigation.

Although our model proved to be valuable in determining prognosis and guiding treatment in OC patients, it should be prospectively confirmed by large-sample clinical studies. This study is a combination of bioinformatics analysis and a small number of in vitro experiments. Consequently, more experiments are warranted to explore and verify the relationship and molecular mechanisms between TMB/TMBrisk and immune infiltration. In addition, restrained by the lack of data of OC patients undergoing immunotherapy, unfortunately, it was not possible to assess whether patients with different TMBrisk benefited differently from immunotherapy. In our further work, studies should be conducted to compare TMBrisk with current biomarkers and explore the relationship between TMBrisk and immunotherapy in OC patients.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  WebbPMJordanSJ. Epidemiology of epithelial ovarian cancer. Best Pract Res Clin Obstet Gynaecol (2017) 41:3–14. doi: 10/gbgn5g

  • Google Scholar

• 2

  FerlayJSoerjomataramIDikshitREserSMathersCRebeloMet al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer (2015) 136(5):E359–86. doi: 10.1002/ijc.29210

  • CrossRef
  • Google Scholar

• 3

  DuskaLRKohnEC. The new classifications of ovarian, fallopian tube, and primary peritoneal cancer and their clinical implications. Ann Oncol (2017) 28(suppl_8):viii8–viii12. doi: 10/gcqwsr

  • Google Scholar

• 4

  EisenhauerEA. Real-world evidence in the treatment of ovarian cancer. Ann Oncol (2017) 28(suppl_8):viii61–5. doi: 10/gcpvtp

  • Google Scholar

• 5

  FelderMKapurAGonzalez-BosquetJHoribataSHeintzJAlbrechtRet al. MUC16 (CA125): tumor biomarker to cancer therapy, a work in progress. Mol Cancer. (2014) 13:129. doi: 10/f578qp

  • Google Scholar

• 6

  ScalettaGPlottiFLuveroDCapriglioneSMonteraRMirandaAet al. The role of novel biomarker HE4 in the diagnosis, prognosis and follow-up of ovarian cancer: a systematic review. Expert Rev Anticancer Ther (2017) 17(9):827–39. doi: 10/ggdxsz

  • Google Scholar

• 7

  ArendRMartinezASzulTBirrerMJ. Biomarkers in ovarian cancer: To be or not to be. Cancer. (2019) 125 Suppl 24:4563–72. doi: 10/gpwb36

  • Google Scholar

• 8

  DochezVCaillonHVaucelEDimetJWinerNDucarmeG. Biomarkers and algorithms for diagnosis of ovarian cancer: CA125, HE4, RMI and ROMA, a review. J Ovarian Res (2019) 12(1):28. doi: 10/ggdxs6

  • Google Scholar

• 9

  Hu-LieskovanSBhaumikSDhodapkarKGrivelJ-CJBGuptaSHanksBAet al. SITC cancer immunotherapy resource document: a compass in the land of biomarker discovery. J Immunother Cancer. (2020) 8(2):e000705. doi: 10/gpwb4d

  • Google Scholar

• 10

  YarchoanMHopkinsAJaffeeEM. Tumor mutational burden and response rate to PD-1 inhibition. N Engl J Med (2017) 377(25):2500–1. doi: 10/gf7vz6

  • Google Scholar

• 11

  SamsteinRMLeeC-HShoushtariANHellmannMDShenRJanjigianYYet al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet (2019) 51(2):202–6. doi: 10/gftd4z

  • Google Scholar

• 12

  HellmannMDNathansonTRizviHCreelanBCSanchez-VegaFAhujaAet al. Genomic features of response to combination immunotherapy in patients with advanced non-Small-Cell lung cancer. Cancer Cell (2018) 33(5):843–852.e4. doi: 10/gdk9r4

  • Google Scholar

• 13

  BaciDBosiAGallazziMRizziMNoonanDMPoggiAet al. The ovarian cancer tumor immune microenvironment (TIME) as target for therapy: A focus on innate immunity cells as therapeutic effectors. Int J Mol Sci (2020) 21(9):E3125. doi: 10/gpgnfj

  • Google Scholar

• 14

  OdunsiK. Immunotherapy in ovarian cancer. Ann Oncol (2017) 28(suppl_8):viii1–7. doi: 10/gf6gqh

  • Google Scholar

• 15

  BrahmerJRTykodiSSChowLQMHwuW-JTopalianSLHwuPet al. Safety and activity of anti-PD-L1 antibody in patients with advanced cancer. N Engl J Med (2012) 366(26):2455–65. doi: 10/gfv28k

  • Google Scholar

• 16

  TopalianSLHodiFSBrahmerJRGettingerSNSmithDCMcDermottDFet al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med (2012) 366(26):2443–54. doi: 10/gcpf9b

  • Google Scholar

• 17

  NeroCCiccaroneFPietragallaADurantiSDanieleGSalutariVet al. Ovarian cancer treatments strategy: Focus on PARP inhibitors and immune check point inhibitors. Cancers (Basel). (2021) 13(6):1298. doi: 10/gpwb5s

  • Google Scholar

• 18

  ChanTAYarchoanMJaffeeESwantonCQuezadaSAStenzingerAet al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol (2019) 30(1):44–56. doi: 10/ghk6xm

  • Google Scholar

• 19

  BiFChenYYangQ. Significance of tumor mutation burden combined with immune infiltrates in the progression and prognosis of ovarian cancer. Cancer Cell Int (2020) 20:373. doi: 10.1186/s12935-020-01472-9

  • CrossRef
  • Google Scholar

• 20

  MiJLXuMLiuCWangRS. Interactions between tumor mutation burden and immune infiltration in ovarian cancer. Int J Clin Exp Pathol (2020) 13(10):2513–23.

  • Google Scholar

• 21

  KoboldtDCZhangQLarsonDEShenDMcLellanMDLinLet al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res (2012) 22(3):568–76. doi: 10/gd5963

  • Google Scholar

• 22

  MayakondaALinD-CAssenovYPlassCKoefflerHP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10/gfmnwf

  • Google Scholar

• 23

  GoodswenSJGondroCWatson-HaighNSKadarmideenHN. FunctSNP: an r package to link SNPs to functional knowledge and dbAutoMaker: a suite of Perl scripts to build SNP databases. BMC Bioinf (2010) 11:311. doi: 10/fgqwk9

  • Google Scholar

• 24

  YoshiharaKShahmoradgoliMMartínezEVegesnaRKimHTorres-GarciaWet al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10/gbd4jn

  • Google Scholar

• 25

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10/gcp3f5

  • Google Scholar

• 26

  RitchieMEPhipsonBWuDHuYLawCWShiWet al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10/f7c4n5

  • Google Scholar

• 27

  LiTFanJWangBTraughNChenQLiuJSet al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res (2017) 77(21):e108–10. doi: 10/gcjvrg

  • Google Scholar

• 28

  MermelCHSchumacherSEHillBMeyersonMLBeroukhimRGetzG. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol (2011) 12(4):R41. doi: 10/dzhjqh

  • Google Scholar

• 29

  GuoCTangYZhangYLiG. Mining TCGA data for key biomarkers related to immune microenvironment in endometrial cancer by immune score and weighted correlation network analysis. Front Mol Biosci (2021) 8:645388. doi: 10/gpwb7n

  • Google Scholar

• 30

  LangfelderPHorvathS. WGCNA: an r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10/bwdbzc

  • Google Scholar

• 31

  BotíaJAVandrovcovaJForaboscoPGuelfiSD’SaKUnited Kingdom Brain Expression Consortiumet al. An additional k-means clustering step improves the biological features of WGCNA gene co-expression networks. BMC Syst Biol (2017) 11(1):47. doi: 10/ggn2cq

  • Google Scholar

• 32

  YuGWangL-GHanYHeQ-Y. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS. (2012) 16(5):284–7. doi: 10/gdf33f

  • Google Scholar

• 33

  BhattacharyaSDunnPThomasCGSmithBSchaeferHChenJet al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. (2018) 5:180015. doi: 10/gc4bth

  • Google Scholar

• 34

  BhattacharyaSAndorfSGomesLDunnPSchaeferHPontiusJet al. ImmPort: disseminating data to the public for the future of immunology. Immunol Res (2014) 58(2–3):234–9. doi: 10/f53z4j

  • Google Scholar

• 35

  TibshiraniR. The lasso method for variable selection in the Cox model. Stat Med (1997) 16(4):385–95. doi: 10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3

  • CrossRef
  • Google Scholar

• 36

  LiCTangZZhangWYeZLiuF. GEPIA2021: integrating multiple deconvolution-based analysis into GEPIA. Nucleic Acids Res (2021) 49(W1):W242–6. doi: 10.1093/nar/gkab418

  • CrossRef
  • Google Scholar

• 37

  LiuC-JHuF-FXiaM-XHanLZhangQGuoA-Y. GSCALite: a web server for gene set cancer analysis. Bioinformatics. (2018) 34(21):3771–2. doi: 10.1093/bioinformatics/bty411

  • CrossRef
  • Google Scholar

• 38

  SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAet al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. (2005) 102(43):15545–50. doi: 10/d4qbh8

  • Google Scholar

• 39

  LiberzonABirgerCThorvaldsdóttirHGhandiMMesirovJPTamayoP. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004

  • CrossRef
  • Google Scholar

• 40

  De SimoneMArrigoniARossettiGGruarinPRanzaniVPolitanoCet al. Transcriptional landscape of human tissue lymphocytes unveils uniqueness of tumor-infiltrating T regulatory cells. Immunity. (2016) 45(5):1135–47. doi: 10/f9f4c2

  • Google Scholar

• 41

  JohnstonRJSuLJPinckneyJCrittonDBoyerEKrishnakumarAet al. VISTA is an acidic pH-selective ligand for PSGL-1. Nature. (2019) 574(7779):565–70. doi: 10/gg2d2z

  • Google Scholar

• 42

  AranDHuZButteAJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18(1):220. doi: 10/gckmjs

  • Google Scholar

• 43

  XuQXuHChenSHuangW. Immunological value of prognostic signature based on cancer stem cell characteristics in hepatocellular carcinoma. Front Cell Dev Biol (2021) 9:710207. doi: 10/gpwb7w

  • Google Scholar

• 44

  FinotelloFMayerCPlattnerCLaschoberGRiederDHacklHet al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med (2019) 11(1):34. doi: 10/gf5jwc

  • Google Scholar

• 45

  BechtEGiraldoNALacroixLButtardBElarouciNPetitprezFet al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17(1):218. doi: 10/f87sgf

  • Google Scholar

• 46

  RacleJGfellerD. EPIC: A tool to estimate the proportions of different cell types from bulk gene expression data. Methods Mol Biol (2020) 2120:233–48. doi: 10/gppzsz

  • Google Scholar

• 47

  YangWSoaresJGreningerPEdelmanEJLightfootHForbesSet al. Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res (2013) 41(Database issue):D955–961. doi: 10.1093/nar/gks1111

  • CrossRef
  • Google Scholar

• 48

  LambJCrawfordEDPeckDModellJWBlatICWrobelMJet al. The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease. Science. (2006) 313(5795):1929–35. doi: 10.1126/science.1132939

  • CrossRef
  • Google Scholar

• 49

  VickersAJElkinEB. Decision curve analysis: a novel method for evaluating prediction models. Med Decis Making. (2006) 26(6):565–74. doi: 10/cgv59z

  • Google Scholar

• 50

  RomeroD. TMB is linked with prognosis. Nat Rev Clin Oncol (2019) 16(6):336. doi: 10/gnbcsz

  • Google Scholar

• 51

  HellmannMDCiuleanuT-EPluzanskiALeeJSOttersonGAAudigier-ValetteCet al. Nivolumab plus ipilimumab in lung cancer with a high tumor mutational burden. N Engl J Med (2018) 378(22):2093–104. doi: 10/gdnqt7

  • Google Scholar

• 52

  ZhangLLiBPengYWuFLiQLinZet al. The prognostic value of TMB and the relationship between TMB and immune infiltration in head and neck squamous cell carcinoma: A gene expression-based study. Oral Oncol (2020) 110:104943. doi: 10/gmn5fw

  • Google Scholar

• 53

  ChanTAWolchokJDSnyderA. Genetic basis for clinical response to CTLA-4 blockade in melanoma. N Engl J Med (2015) 373(20):1984. doi: 10/gng7vq

  • Google Scholar

• 54

  BronderDTigheAWangsaDZongDMeyerTJWardenaarRet al. TP53 loss initiates chromosomal instability in fallopian tube epithelial cells. Dis Model Mech (2021) 14(11):dmm049001. doi: 10/gpwrxm

  • Google Scholar

• 55

  da CostaAABABaiocchiG. Genomic profiling of platinum-resistant ovarian cancer: The road into druggable targets. Semin Cancer Biol (2021) 77:29–41. doi: 10.1016/j.semcancer.2020.10.016

  • CrossRef
  • Google Scholar

• 56

  JiaQWangJHeNHeJZhuB. Titin mutation associated with responsiveness to checkpoint blockades in solid tumors. JCI Insight (2019) 4(10):127901. doi: 10/gm4pvc

  • Google Scholar

• 57

  FanfaniVCitiLHarrisALPezzellaFStracquadanioG. The landscape of the heritable cancer genome. Cancer Res (2021) 81(10):2588–99. doi: 10/gpwrx5

  • Google Scholar

• 58

  ZhangHFuL. The role of ALDH2 in tumorigenesis and tumor progression: Targeting ALDH2 as a potential cancer treatment. Acta Pharm Sin B (2021) 11(6):1400–11. doi: 10/gpwrx7

  • Google Scholar

• 59

  KangKXieFMaoJBaiYWangX. Significance of tumor mutation burden in immune infiltration and prognosis in cutaneous melanoma. Front Oncol (2020) 10:573141. doi: 10.3389/fonc.2020.573141

  • CrossRef
  • Google Scholar

• 60

  ZhangCLiZQiFHuXLuoJ. Exploration of the relationships between tumor mutation burden with immune infiltrates in clear cell renal cell carcinoma. Ann Transl Med (2019) 7(22):648. doi: 10.21037/atm.2019.10.84

  • CrossRef
  • Google Scholar

• 61

  JardimDLGoodmanAde Melo GagliatoDKurzrockR. The challenges of tumor mutational burden as an immunotherapy biomarker. Cancer Cell (2021) 39(2):154–73. doi: 10/ghz4tc

  • Google Scholar

• 62

  LeeGR. Phenotypic and functional properties of tumor-infiltrating regulatory T cells. Mediators Inflamm (2017) 2017:5458178. doi: 10/gcvpxf

  • Google Scholar

• 63

  CurielTJ. Tregs and rethinking cancer immunotherapy. J Clin Invest. (2007) 117(5):1167–74. doi: 10/d2j6vp

  • Google Scholar

• 64

  PaiC-CSSimonsDMLuXEvansMWeiJWangY-Het al. Tumor-conditional anti-CTLA4 uncouples antitumor efficacy from immunotherapy-related toxicity. J Clin Invest. (2019) 129(1):349–63. doi: 10/ggcjnk

  • Google Scholar

• 65

  SharmaNVacherJAllisonJP. TLR1/2 ligand enhances antitumor efficacy of CTLA-4 blockade by increasing intratumoral treg depletion. Proc Natl Acad Sci U S A. (2019) 116(21):10453–62. doi: 10.1073/pnas.1819004116

  • CrossRef
  • Google Scholar

• 66

  HwangW-TAdamsSFTahirovicEHagemannISCoukosG. Prognostic significance of tumor-infiltrating T cells in ovarian cancer: a meta-analysis. Gynecol Oncol (2012) 124(2):192–8. doi: 10/bfjs8z

  • Google Scholar

• 67

  SatoEOlsonSHAhnJBundyBNishikawaHQianFet al. Intraepithelial CD8+ tumor-infiltrating lymphocytes and a high CD8+/regulatory T cell ratio are associated with favorable prognosis in ovarian cancer. Proc Natl Acad Sci U S A. (2005) 102(51):18538–43. doi: 10/crfzst

  • Google Scholar

• 68

  BoutilierAJElsawaSF. Macrophage polarization states in the tumor microenvironment. Int J Mol Sci (2021) 22(13):6995. doi: 10.3390/ijms22136995

  • CrossRef
  • Google Scholar

• 69

  WangHYungMMHNganHYSChanKKLChanDW. The impact of the tumor microenvironment on macrophage polarization in cancer metastatic progression. Int J Mol Sci (2021) 22(12):6560. doi: 10.3390/ijms22126560

  • CrossRef
  • Google Scholar

• 70

  NowakMKlinkM. The role of tumor-associated macrophages in the progression and chemoresistance of ovarian cancer. Cells. (2020) 9(5):E1299. doi: 10.3390/cells9051299

  • CrossRef
  • Google Scholar

• 71

  YuanXZhangJLiDMaoYMoFDuWet al. Prognostic significance of tumor-associated macrophages in ovarian cancer: A meta-analysis. Gynecol Oncol (2017) 147(1):181–7. doi: 10.1016/j.ygyno.2017.07.007

  • CrossRef
  • Google Scholar

• 72

  MacciòAGramignanoGCherchiMCTancaLMelisLMadedduC. Role of M1-polarized tumor-associated macrophages in the prognosis of advanced ovarian cancer patients. Sci Rep (2020) 10(1):6096. doi: 10.1038/s41598-020-63276-1

  • CrossRef
  • Google Scholar

• 73

  WangYZhangQGaoZXinSZhaoYZhangKet al. A novel 4-gene signature for overall survival prediction in lung adenocarcinoma patients with lymph node metastasis. Cancer Cell Int (2019) 19:100. doi: 10/gpwr24

  • Google Scholar

• 74

  YinXWangPYangTLiGTengXHuangWet al. Identification of key modules and genes associated with breast cancer prognosis using WGCNA and ceRNA network analysis. Aging (Albany NY). (2020) 13(2):2519–38. doi: 10/gn2j65

  • Google Scholar

• 75

  LuoZWangWLiFSongyangZFengXXinCet al. Pan-cancer analysis identifies telomerase-associated signatures and cancer subtypes. Mol Cancer. (2019) 18(1):106. doi: 10.1186/s12943-019-1035-x

  • CrossRef
  • Google Scholar

• 76

  VernocchiPGiliTConteFDel ChiericoFContaGMiccheliAet al. Network analysis of gut microbiome and metabolome to discover microbiota-linked biomarkers in patients affected by non-small cell lung cancer. Int J Mol Sci (2020) 21(22):E8730. doi: 10.3390/ijms21228730

  • CrossRef
  • Google Scholar

• 77

  MuWChengXZhangXLiuYLvQLiuGet al. Hinokiflavone induces apoptosis via activating mitochondrial ROS/JNK/caspase pathway and inhibiting NF-κB activity in hepatocellular carcinoma. J Cell Mol Med (2020) 24(14):8151–65. doi: 10.1111/jcmm.15474

  • CrossRef
  • Google Scholar

• 78

  ZhangJBuXWangHZhuYGengYNihiraNTet al. Author correction: Cyclin d-CDK4 kinase destabilizes PD-L1 via cullin 3-SPOP to control cancer immune surveillance. Nature. (2019) 571(7766):E10. doi: 10/gpwsfk

  • Google Scholar

• 79

  LiJStangerBZ. Cell cycle regulation meets tumor immunosuppression. Trends Immunol (2020) 41(10):859–63. doi: 10.1016/j.it.2020.07.010

  • CrossRef
  • Google Scholar

• 80

  BaoCWangJMaWWangXChengY. HDGF: a novel jack-of-all-trades in cancer. Future Oncol (2014) 10(16):2675–85. doi: 10/f6vgrn

  • Google Scholar

• 81

  GiriKPabelickCMMukherjeePPrakashYS. Hepatoma derived growth factor (HDGF) dynamics in ovarian cancer cells. Apoptosis. (2016) 21(3):329–39. doi: 10/f8cmgf

  • Google Scholar

• 82

  LiuX-JLiuW-LYangF-MYangX-QLuX-F. Hepatoma-derived growth factor predicts unfavorable prognosis of epithelial ovarian cancer. Onco Targets Ther (2015) 8:2101–9. doi: 10/gpwsw2

  • Google Scholar

• 83

  MooreARRosenbergSCMcCormickFMalekS. RAS-targeted therapies: is the undruggable drugged? Nat Rev Drug Discovery (2020) 19(8):533–52. doi: 10/ghp85r

  • Google Scholar

• 84

  MoujaberTBalleineRLGaoBMadsenIHarnettPRDeFazioA. New therapeutic opportunities for women with low-grade serous ovarian cancer. Endocr Relat Cancer. (2021) 29(1):R1–R16. doi: 10/gpwsw6

  • Google Scholar

• 85

  WangHFuY. NR1D1 suppressed the growth of ovarian cancer by abrogating the JAK/STAT3 signaling pathway. BMC Cancer. (2021) 21(1):871. doi: 10/gpwsw9

  • Google Scholar

• 86

  BoudhraaZZaouiKFleuryHCahuzacMGilbertSTchakarskaGet al. NR1D1 regulation by ran GTPase via miR4472 identifies an essential vulnerability linked to aneuploidy in ovarian cancer. Oncogene. (2022) 41(3):309–20. doi: 10/gpwsxh

  • Google Scholar

• 87

  NarayanaswamyPBHodjatMHallerHDumlerIKiyanY. Loss of urokinase receptor sensitizes cells to DNA damage and delays DNA repair. PloS One (2014) 9(7):e101529. doi: 10/f6gppk

  • Google Scholar

• 88

  McGrailDJPiliéPGRashidNUVoorwerkLSlagterMKokMet al. High tumor mutation burden fails to predict immune checkpoint blockade response across all cancer types. Ann Oncol (2021) 32(5):661–72. doi: 10/gk43m7

  • Google Scholar

• 89

  CodyNALShenZRipeauJ-SProvencherDMMes-MassonA-MChevretteMet al. Characterization of the 3p12.3-pcen region associated with tumor suppression in a novel ovarian cancer cell line model genetically modified by chromosome 3 fragment transfer. Mol Carcinog (2009) 48(12):1077–92. doi: 10/fk2qz2

  • Google Scholar

• 90

  KatohM. Molecular cloning and characterization of ST7R (ST7-like, ST7L) on human chromosome 1p13, a novel gene homologous to tumor suppressor gene ST7 on human chromosome 7q31. Int J Oncol (2002) 20(6):1247–53.

  • Google Scholar

• 91

  YangZWangX-LBaiRLiuW-YLiXLiuMet al. miR-23a promotes IKKα expression but suppresses ST7L expression to contribute to the malignancy of epithelial ovarian cancer cells. Br J Cancer. (2016) 115(6):731–40. doi: 10/f86tph

  • Google Scholar

• 92

  NikasJBMitanisNTNikasEG. Whole exome and transcriptome RNA-sequencing model for the diagnosis of prostate cancer. ACS Omega. (2020) 5(1):481–6. doi: 10/ghv3wx

  • Google Scholar

• 93

  RuanBFengXChenXDongZWangQXuKet al. Identification of a set of genes improving survival prediction in kidney renal clear cell carcinoma through integrative reanalysis of transcriptomic data. Dis Markers. (2020) 2020:8824717. doi: 10/gpwszx

  • Google Scholar

• 94

  JiaoJJiangLLuoY. N6-Methyladenosine-Related RNA signature predicting the prognosis of ovarian cancer. Recent Pat Anticancer Drug Discovery (2021) 16(3):407–16. doi: 10/gmwzj3

  • Google Scholar

• 95

  GjerstorffMFAndersenMHDitzelHJ. Oncogenic cancer/testis antigens: prime candidates for immunotherapy. Oncotarget. (2015) 6(18):15772–87. doi: 10/f72423

  • Google Scholar

• 96

  GuoRQinY. LEMD1-AS1 suppresses ovarian cancer progression through regulating miR-183-5p/TP53 axis. Onco Targets Ther (2020) 13:7387–98. doi: 10/gpws3p

  • Google Scholar

• 97

  LiZZhangYSuiSHuaYZhaoATianXet al. Targeting HMGB3/hTERT axis for radioresistance in cervical cancer. J Exp Clin Cancer Res (2020) 39(1):243. doi: 10/gpws36

  • Google Scholar

• 98

  MukherjeeAHuynhVGainesKRehWAVasquezKM. Targeting the high-mobility group box 3 protein sensitizes chemoresistant ovarian cancer cells to cisplatin. Cancer Res (2019) 79(13):3185–91. doi: 10/gpws6g

  • Google Scholar

• 99

  Laughlin-TommasoSKHesleyGKHopkinsMRBrandtKRZhuYStewartEA. Clinical limitations of the international federation of gynecology and obstetrics (FIGO) classification of uterine fibroids. Int J Gynaecol Obstet. (2017) 139(2):143–8. doi: 10/gpws6m

  • Google Scholar

• 100

  ShimasakiNJainACampanaD. NK cells for cancer immunotherapy. Nat Rev Drug Discovery (2020) 19(3):200–18. doi: 10.1038/s41573-019-0052-1

  • CrossRef
  • Google Scholar

• 101

  MaasRJHoogstad-van EvertJSvan der MeerJMMekersVRezaeifardSKormanAJet al. TIGIT blockade enhances functionality of peritoneal NK cells with altered expression of DNAM-1/TIGIT/CD96 checkpoint molecules in ovarian cancer. Oncoimmunology. (2020) 9(1):1843247. doi: 10.1080/2162402X.2020.1843247

  • CrossRef
  • Google Scholar

• 102

  Van der MeerJMRMaasRJAGuldevallKKlarenaarKde JongePKJDEvertJSHet al. IL-15 superagonist n-803 improves IFNγ production and killing of leukemia and ovarian cancer cells by CD34+ progenitor-derived NK cells. Cancer Immunol Immunother. (2021) 70(5):1305–21. doi: 10/gjbffn

  • Google Scholar

• 103

  KimH-JCantorH. CD4 T-cell subsets and tumor immunity: the helpful and the not-so-helpful. Cancer Immunol Res (2014) 2(2):91–8. doi: 10.1158/2326-6066.CIR-13-0216

  • CrossRef
  • Google Scholar

• 104

  NelsonBH. The impact of T-cell immunity on ovarian cancer outcomes. Immunol Rev (2008) 222:101–16. doi: 10.1111/j.1600-065X.2008.00614.x

  • CrossRef
  • Google Scholar