Classification of ovarian cancer associated with BRCA1 mutations, immune checkpoints, and tumor microenvironment based on immunogenomic profiling

Data

Gene expression profiles were mined from The Cancer Genome Atlas (TCGA) dataset (https://tcga-data.nci.nih.gov/tcga/) and consisted of normalized gene expression patterns for 379 ovarian cancer samples that were mapped using fragments per kilobase of transcript per million fragments. Clinical data for age, survival, stage, and tumor grade were also mined from TCGA. The somatic mutation data were obtained from single nucleotide polymorphism (SNP) data in the TCGA repository using Mutect. The expression data of the validation dataset was obtained from the Gene Expression Omnibus (GEO) repository (GSE51088), and contained 172 ovarian cancer samples. All computational and statistical analyses were performed using R software (version 3.6.1, http://www.R-project.org).

ssGSEA and clustering

We obtained 29 immune-associated gene sets, totaling 707 genes, representing different immune cell types, functions, and pathways (File S1) (He et al., 2018; Yue, Ma & Zhou, 2019). We used single sample gene set enrichment analysis (ssGSEA), using the R package GSVA (version 1.34.0), to calculate the enrichment scores of the 29 immune biosignatures for each sample in the tumor microenvironment (Barbie et al., 2009; Hänzelmann, Castelo & Guinney, 2013). ssGSEA calculated the gene signature overexpression scores by comparing the gene levels in the signature with those in all the other genes in the transcriptome. An unsupervised machine learning method was used to perform hierarchical clustering of ovarian cancer into three clusters. The clusters were further divided into three subtypes according to the immune scores: subtype 1, subtype 2, and subtype 3.

ESTIMATE and CIBERSORT

Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) (Yoshihara et al., 2013) is an approach that uses gene expression biosignatures to deduce the proportion of stromal and immune cells in tumor samples, which form the major non-tumor constituents of tumor samples. By performing ssGSEA, stromal and immune scores are calculated to estimate the levels of invading stromal and immune cells. This forms the rationale for the ESTIMATE score to deduce tumor purity in the tumor tissue. Tumor purity was calculated based on the gene expression data from the immune and stromal scores using ESTIMATE in R package. Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) (Newman et al., 2015) is a biological tool that uses the deconvolution strategy to compute the fractions of the 22 human immune cell types. We selected 1000 permutations and data with P < 0.05 as the maxim for the successful deconvolution of a sample. The Kruskal-Wallis test was used to compare the proportions of immune cell types among ovarian cancer subtypes.

Calculation of TMB scores

TMB is the overall enumeration of mutations per million bases in tumor tissue. It details the mutation density of tumor genes, i.e., the enumeration of mutations in the tumor genome, including the total number of genetic coding errors, base substitutions, and gene insertions or deletions. We computed the mutation frequency with the number of variants/the length of exons (38 million) for each sample using Perl (v5.30.1, https://www.perl.org/).

Survival analyses

We obtained the follow-up information of patients from their clinical data and calculated the significance of survival time using the log-rank test and their differences using a threshold of P < 0.05. The median was used as the cut-off value to divide the samples into high- or low scores to obtain the relationship between the related immune gene sets score and prognosis. We plotted the Kaplan–Meier curves to indicate the differences in the survival periods.

Gene-set enrichment analysis

Gene-set enrichment analysis of TCGA datasets was conducted in the GSEA (R package) (Subramanian et al., 2005). The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were used to assess the functional role of the differentially expressed genes between subtype 1 and subtype 3. Differential gene set enrichment was inspected by the limma R package. P < 0.05 was used as the cut-off value.

Immunogenomic profiling identifies three ovarian cancer subtypes

We obtained enrichment scores for each sample of the 29 immune-associated gene sets in the tumor microenvironment using ssGSEA. Ovarian cancer cases were then hierarchically clustered into three classes according to their immune scores. These classes were clearly defined and classified as subtype 1, subtype 2, and subtype 3, which represented low immunity, median immunity, and high immunity, respectively (Fig. 1A). We determined the immune, stromal, ESTIMATE scores, and tumor purity using the ESTIMATE algorithm. Our results showed that the immune and stromal scores were significantly different (P < 0.001) and were highest in subtype 3 and lowest in subtype 1; tumor purity was the highest in subtype 1 and lowest in subtype 3 (Figs. 1B–1D). These results indicated that immune and stromal cells are most prevalent in subtype 3, while tumor cells are most prevalent in subtype 1.

We found that most of the human leukocyte antigen (HLA) gene expression levels were lowest in subtype 1 and highest in subtype 3 (P < 0.001) (Fig. S1A). The expression levels of many immune cell subgroup biomarker genes, such as FOXP3 [regulatory T cell (Treg)], CD45RO (memory T cell), CD8A (cytotoxic T cell), CD20 (B cell), CD1A [immature dendritic cell (iDC)], CXCR5 (Tfh cell), and IL3RA [plasmacytoid dendritic cell (pDC)] were remarkably higher in subtype 3 and markedly lower in subtype 1 (Fig. S1B).

Three subtypes show differential expression of immune checkpoint genes

We analyzed the expression levels of the checkpoint receptors in the three ovarian cancer subtypes responsible for decreasing T cell bioactivity, including PDCD1 (PD1), CTLA4, LAG-3, and TIM-3. We then analyzed the PDCD1 ligand CD274 (PD-L1), PDCD1LG2 (PD-L2), CTLA4 ligand CD86, and CD80. We found that the expression levels of these 8 immune checkpoint genes were all remarkably lower in subtype 1 and remarkably higher in subtype 3 (P < 0.001) (Figs. 2A–2H). This indicated that the immunophenotype of our hierarchical clusters was clearly distinguished, and that ovarian cancer subtype 3 may respond more effectively to checkpoint inhibitor therapy.

Analysis of the TMB and BRCA mutations among the three subtypes of ovarian cancer

TMB is a predictor of tumor behavior and immunological response in a diverse range of cancers (Goodman et al., 2017). In general, tumors with a high TMB have elevated levels of neoantigens, which play an important role in immunotherapy activities (Goodman et al., 2017; Schumacher & Schreiber, 2015). We mined the somatic mutation profiles of 436 ovarian cancer patients from the SNP data in TCGA using Mutect, and calculated the TMB using the enumeration of mutation events per million bases. We then analyzed the TMB between the three subtypes of ovarian cancer and found that the three subtypes were not significantly correlated with TMB (p = 0.726) (Fig. 3A).

An increasing number of studies have documented that targeted therapies can stimulate the immune response of the host. The relationship between BRCA mutations and immunity is of great interest at present. We analyzed the connection linking BRCA1 and BRCA2 mutations in the three subtypes. The data from BRCA1 and BRCA2 mutations were mined from the SNP data in TCGA via Mutect. There were 23 patients with the BRCA1 mutation and 20 patients with the BRCA2 mutation out the 436 ovarian cancer patients sampled. We found 13 patients with the BRCA1 mutation and 13 patients with the BRCA2 mutation out of 274 ovarian cancer patients using the intersection between the mutation data and the immunity cluster data samples. Surprisingly, all BRCA1 mutation patients were in subtype 3, and the difference was significant (χ2 test, P = 0.0016) (Fig. 3B). The BRCA2 mutation ratio was greater in subtype 2 and subtype 3 compared with subtype 1, but this difference was not statistically significant (χ2 test, P = 0.577) (Fig. 3C).

Different immune cells among the 3 subtypes of ovarian cancer

CIBERSORT can deduce 22 types of human immune cells, such as B cells, myeloid subset cells, T cells, NK cells, macrophages, and DCs, according to the gene expression data, using the gene-based deconvolution algorithm method (Newman et al., 2015). We set 1000 permutations and P < 0.05 as the maxim for the successful deconvolution of a sample and found that CD8 T cells, CD4 memory activated T cells, Tregs, macrophages M1, and resting dendritic cells were all at the highest levels in subtype 3 and the lowest levels in subtype 1 (P < 0.01). However, activated dendritic cells showed an opposite trend (Fig. 4).

Prognostic analysis of the ovarian cancer subtypes and immune-associated gene sets

Survival analyses indicated significantly different prognoses among the three ovarian cancer subtypes. Subtype 2 had the worst survival prognosis among the three subtypes and there was no significant difference in survival between subtypes 1 and 3 (Fig. 5A). We analyzed the prognostic value of the different immune gene set expression scores for predicting patient survival. We found that the high expression level of the check-point, major histocompatibility complex (MHC) class I, APC co-inhibition, T cell co-inhibition, Th1 cells, Th2 cells, Tfh, inflammation-promoting, and Tregs was associated with a significantly better prognosis (P < 0.05) compared with low expression levels (Figs. 5B–5J).

Identification of the ovarian cancer subtype-specific pathways and GO

GSEA identified 628 GO terms and 56 KEGG terms in subtypes 1 and 3. GO analysis indicated that the immunoglobulin complex, circulating immunoglobulin complex, MHC class II protein complex, immunoglobulin receptor binding, and the MHC protein complex were the top five significantly enriched biological processes in subtype 3. In addition, glucuronidation, metabolic processes, and methyl-CpG binding were the most enriched terms in subtype 1 (Figs. 6A and 6B). GSEA showed that the immune-associated pathways were most active in subtype 3, and consisted of Th17 cell differentiation, NF-κB signaling axis, the B cell receptor signaling cascade, the T cell receptor signaling cascade, PD-L1 expression and the PD-1 checkpoint axis in cancer, the IL-17 signaling cascade, and the tumor necrosis factor (TNF) signaling axis. Our results verified that immune activity was elevated in subtype 3. However, subtype 1 was enriched in pathways, including maturity onset diabetes of the young, ascorbate and aldarate metabolism, pentose and glucuronate interconversions, fat digestion and absorption, and porphyrin and chlorophyll metabolism (Figs. 6C and 6D). These cascades may be inversely linked to ovarian cancer immunity.

Validation of external datasets

We used the same method to hierarchically cluster ovarian cancer in the GSE51088 dataset, which included 172 ovarian cancer samples. Interestingly, there was a similar clustering result, with three distinct clusters (Fig. 7A). The immune and stromal scores were remarkably higher in subtype 3 and were markedly lower in subtype 1, while tumor purity showed an opposite result (Figs. 7B–7D). Most HLA genes and CD8A, CD1A, CD45R, and IL3RA expression levels were significantly lower in subtype 1 and significantly higher in subtype 3, which was consistent with the TCGA datasets (Figs. S2A and S2B). The expressions of the immune checkpoint genes, including PDCD1, CD274, CTLA4, CD80, CD86, TIM-3, and LAG-3 were all remarkably lower in subtype 1 and were significantly higher in subtype 3 (Figs. 7E–7L). These results were similar to results from previous research. Our results suggest that there are varied immune status subtypes in ovarian cancer that may affect the treatment of immune checkpoints differently.