Age-related mutational signature negatively associated with immune activity and survival outcome in triple-negative breast cancer.

Triple-negative breast cancer (TNBC) is characterized by broad genomic and transcriptional heterogeneity and results in a worse prognosis than other breast cancer types. Here, we integrated genomic and transcriptomic data combined with clinicopathologic information from 538 patients with TNBC and identified four novel significantly mutated genes (SMGs), namely, KDM6A, CD86, RBM47, and IFNGR1. A mutational signature (known as age-related signature 1) featured by enrichment of C > T mutations at NpCpG trinucleotides was associated with worse survival in TNBC (HR, 1.76 [95% CI, 1.07-2.90]; P = .025). We also analyzed gene transcriptomic profiles of TNBC samples and identified immune regulation-related gene pathways (e.g., antigen processing presentation and interferon signaling), and the cell cycle was significantly altered in samples with different signature 1 activity groups. The analysis further revealed that signature 1 was associated with decreased tumor immunogenicity and immunocyte infiltration in TNBC. This negative correlation was also observed in lung adenocarcinoma and prostate cancer samples. Furthermore, we found that patients with mutational signature 1 were markedly associated with decreased tumor mutation burden, even after controlling for age, grade, histological type, lymph node status, mutations in BRCA1/2 and ATR, and APOBEC and homologous recombination repair deficiency signatures (OR, 0.19 [95% CI, 0.11-0.32]; P < .001). Overall, this study identified previously unreported SMGs and re-annotated that age-related signature 1 was associated with a weakened immune microenvironment and predictive of poor survival in TNBC, offering opportunities to stratify patients into optimal treatment plans based on genomic subtyping.

Introduction

Triple-negative breast cancer (TNBC) accounts for 10–20% of breast carcinomas and is characterized by a lack of expression of the estrogen receptor and progesterone receptor and absence of ERBB2 (also known as HER2) amplification.[1,2] TNBCs result in worse prognosis than other types of breast cancers (BCs), with an increased likelihood of early distant recurrence and death,[1] owing mainly to the extensive molecular heterogeneity of the disease.[3,4]

Recent next-generation sequencing has become a powerful tool to elucidate potential driver genetic aberrations underlying cancer development. TNBCs are significantly associated with BRCA1 germline mutations and high levels of genomic instability; TP53 and PIK3CA are the two most significantly mutated genes (SMGs).[5] Mutational signatures are the fingerprints of endogenous and exogenous factors that have acted over the course of tumorigenesis and heterogeneity. Age-related signature 1 was dominated by C > T mutation, which is associated with spontaneous deamination of 5-methylcytosine.[6] A recent study reported that younger Asian BCs harbor a more immune-active microenvironment than older Western BCs.[7] Signature 3, characterized by T > C mutations, was associated with a deficiency in genomic integrity maintenance (homologous recombination repair deficiency [HRD]) and strongly associated with BRCA1/2 mutations in breast cancer.[8]

Clinical studies also found that patients with breast cancer with high immune infiltration, particularly TNBC and HER2-positive breast cancer, had better prognosis and treatment response.[9,10] Recently, the correlation between the specific mutational signature process (e.g., HRD and APOBEC signatures) and antitumor immune activity has been observed across cancer types;[11-13] however, the characteristics of immune infiltration and its association with mutational signatures in TNBC are still unclear. There is also evidence that immune-rich TNBCs may be under immune surveillance that continuously eliminates many immunogenic clones, resulting in lower clonal heterogeneity.[10]

Since somatic mutations can be applied to clinical practice, the assessment of larger patient series is crucial in identifying tumor genomic aberrations that are potentially predictive responses and investigating whether patients have distinct clinical outcomes. Hence, the present study aimed to identify new SMGs and characterize mutational signatures that driving the molecular subtypes with genetic prognosticators and immune activity in patients with TNBC from published multi-omic studies. We believe that these findings may be applicable for prognostic prediction and therapeutic guidance for TNBC.

Results

SMGs in TNBC

Somatic mutational profiles of 538 patients with TNBC from previous whole-exome sequencing (WES) genomic studies were analyzed. A median of 118 mutations per sample (range, 1–2205) in a total of 93897 coding somatic mutations were collected from previously published studies. We combined MutSigCV and OncodriveFML algorithms to identify SMGs that met the criteria of positive accumulation, clustering at a hotspot, and functional importance. We identified a total of 16 SMGs in the pooled TNBC datasets (Figure 1, Supplementary Table S1). Apart from well-known driver oncogenes (e.g., TP53, PIK3CA, FOXA1, KMT2 C/D, PTEN, RB1, BRCA1/2), we identified four new SMGs, namely, KDM6A (15 of 538, 2.8%), CD86 (8 of 538, 1.5%), RBM47 (7 of 538, 1.3%), and IFNGR1 (5 of 538, 0.9%), which were not previously implicated in TNBC. The mutation plots of these four novel SMGs are shown in Supplementary Figure S1.

Figure 1.

Mutational landscape of significantly mutated genes in triple-negative breast cancer (TNBC) patients.

We also examined the mRNA expression levels of these four novel SMGs in tumor tissues versus matched adjacent normal control tissue in a separate microarray-based TNBC gene expression dataset.[14] The analyses showed that CD86 was significantly upregulated in tumor tissues, whereas KDM6A was significantly downregulated in tumor tissues (paired t-test, P < .05, Supplementary Figure S2). We also examined the association between SMG mutation status and survival prognosis. TNBC samples with RBM47 mutations were found to be associated with short survival outcomes in the aggregated TNBC cohort (log-rank test, P < .001; Supplementary Figure S3). A previous study also identified RBM47 as a suppressor of breast cancer progression and metastasis, and patients with low expression tended to have a poor clinical outcome.[15,16]

Mutational signatures operative in the aggregated TNBC cohort

To gain further insights into the mutational processes operative in patients with TNBC, we delineated the mutation signatures from the genomic mutational profile. We first counted the number of single-nucleotide variants in the matrix of 96 possible mutations occurring in a trinucleotide context in each TNBC sample and found that the predominant mutations were the C > T and C > G transitions at TpCpW trinucleotide sites (Figure 2a). Subsequently, we extracted three mutational signatures from the TNBC genomic data with varying mutational activities and annotated them against the Catalog of Somatic Mutations in Cancer (COSMIC) signature nomenclature (Figure 2b). Cosine similarity analysis demonstrated that the three signatures showed high similarity to known COSMIC signatures 1, 3, and 13 (cosine coefficients, 0.95, 0.91, and 0.86, respectively; Supplementary Figure S4). The proportion and activities of the three extracted signatures in every tumor sample (i.e., number of somatic mutations attribute to each signature) are presented in Supplementary Figure S5. We observed that signature 1 (mutational activity, 16270/82770, 19.4%) was dominated by C > T transitions at NpCpG trinucleotides (often termed CpG dinucleotides) that was associated with the age-related accumulation of spontaneous deamination of 5-methylcytosine in most cancer types. Signature 3, which was the most prevalent signature, accounting for 56899 in 82770 mutations (68.8%), was reported to have failure of DNA double-strand break repair by homologous recombination in breast, pancreatic, and ovarian cancers. Signature 13 (9601/82770, 11.8%), characterized by C > T mutations at TpCpW (where W = A or T) trinucleotide sequences, was attributed to the activity of the AID/APOBEC family of cytidine deaminases (APOBEC) in breast, bladder, cervical, and non-small cell lung cancers (NSCLCs)[17] (Figure 2c).

Figure 2.

Mutational signatures extracted from the aggregated TNBC dataset.

Mutational signature predictive of TNBC survival

To identify the extracted mutational signature activities associated with TNBC clinical features, we applied unsupervised hierarchical clustering for activities of signatures[18] and identified two distinctive clusters, C1 and C2 (Supplementary Figure S6a). Survival analysis showed that their association with prognosis was marginally statistically significant (log-rank test, P = .075; Cox regression, HR, 0.57 [95% CI, 0.28–1.14]; Supplementary Figure S6b–c). Therefore, we followed the previously recommended method and criteria and stratified the extracted mutational signatures into binary variable (i.e., no and yes) model.[8,19] We observed that patients with mutational signature 1 were markedly associated with worse survival outcomes (log-rank test, P = .025; HR, 1.76 [95% CI, 1.07–2.90]; Figure 3a). To rule out the possibility that associations between signature 1 and prognosis were affected by confounding factors, we incorporated age, grade, histology type, and lymph node status into the multivariate Cox regression model. Associations between signature 1 and the survival of patients with TNBC remained statistically significant after controlling for these factors (HR, 1.85 [95% CI, 1.08–3.18]; P = .024; Figure 3b). We also investigated the correlation between signature 1 activities and patient age and found a positive correlation using the Spearman regression model (Spearman r = 0.197, P < .001; Supplementary Figure S7), further confirming previous findings.[6] We also examined the difference between mutational activity and signature 1 with respect to the mutational status of these 16 SMGs. Increased mutational activities of signature 1 were observed in samples with mutations in CBFB, FOXA1, NF1, PIK3CA, and TP53 (Wilcoxon rank-sum test, all P < .05; Supplementary Figure S8).

Figure 3.

Association of mutational signature 1 with prognosis in TNBC cohort.

Immune-related pathways and immunocyte infiltration associated with mutational signature 1

In this analysis, we investigated the potential mechanism behind mutational signature 1, which may regulate TNBC prognosis. Gene set enrichment analysis (GSEA) on REACTOME sets revealed that the enrichment of genes involved in the immune system, cell cycle, and class I MHC antigen processing presentation were significantly altered in samples with different signature 1 grouping (q < 0.001, Figure 4a, Supplementary Figure S9). We further analyzed the determinants of tumor immunogenicity with respect to mutational signature 1 using The Cancer Immunome Atlas (TCIA)-recommended immunophenoscore (IPS) scheme. A significant negative correlation between signature 1 activity and immunogenicity level was identified across the TNBC samples (Spearman r = −0.13, P = .001, Figure 4b). We also conducted mutational signature analysis on other age-related tumor types from the TCGA dataset, including lung adenocarcinoma (LUAD) and prostate adenocarcinoma (PRAD) (Supplementary Figure S10a-b), and found that signature 1 was inversely correlated with the IPS (Spearman correlation test, LUAD, r = −0.23; PRAD, r = −0.14; P < .001 Supplementary Figure S10c-d).

Figure 4.

Significantly enriched pathways and immune infiltration alteration with signature 1 activities.

The T cell-inflamed gene expression profile (GEP) and immune cytolytic activity (CYT) were also analyzed, and a lower score was found in the signature 1 group (GEP scores, 0.25 vs −0.30, P = .004, Figure 4c; CYT, 15.7 vs 11.2, P = .042, Supplementary Figure S11; Wilcoxon rank-sum test). We further compared clonal heterogeneity measured using the mutant allele tumor heterogeneity (MATH) score, which quantified the dispersion of variant allele frequencies in TNBC samples. Accordingly, patients with signature 1 had significantly higher MATH scores, indicating stronger clonal genomic heterogeneity (MATH scores, 43.6 vs 52.2, Wilcoxon rank-sum test, P = .007; Figure 4d). Moreover, we evaluated (using the CIBERSORT algorithm) the abundance of tumor-infiltrating immunocytes in the TNBC microenvironment using the gene expression profile data. We found that CD8+ T cells, monocytes, macrophages, and antitumor leukocytes were less clustered in the mutational signature 1 group (Figure 4e).

Tumor mutation burden (TMB) associated with mutational signature 1

TNBC is genomically heterogeneous with varying TMBs, and its level is associated with the survival and immune infiltration of multiple cancer types.[20-22] Therefore, we compared signature 1 with TMB and found a significantly lower TMB in patients with signature 1 (log2 TMB, median, 7.1 vs. 6.3, P < .001; Wilcoxon rank-sum test, Figure 5a). Furthermore, MHC Class I-associated tumor neoantigen counts (burden) in the TCGA dataset also showed a similar tendency (log2 TNB, median, 6.7 vs 5.4, P < .001; Figure 5b). TMB was largely attributed to genomic instability. To rule out the influence of confounding factors, we incorporated genomic integrity maintenance-associated genes (mutations in BRCA1/2, ATR, and SMGs in TNBC), extracted mutational signatures (1, 3, 13), and clinical features (age, grade, histology type, and lymph node status) into the Bayesian logistic regression model. The association between signature 1 and TMB remained statistically significant after adjusting for such factors (OR, 0.19 [95%CI, 0.11–0.32]; P < .001; Figure 5c).

Figure 5.

Association of signature 1 with tumor mutation burden in pooled TNBC cohort.

Discussion

In this study, a genomic meta-analysis of 538 TNBC samples from five published studies was performed and identified several novel less frequently mutated SMGs that were not recognized previously. We revealed that patients with age-related mutational signature 1 were associated with worse prognosis in TNBC and confirmed this finding in a multivariable Cox regression model. We further identified that this age-related signature was associated with lower immune activities and decreased TMB, suggesting that cancers with signature 1 probably escaped from immune surveillance.

Immune infiltration and surveillance in determining the prognosis of various types of cancers are increasingly recognized.[10,23,24] More than 70% of breast cancers contain at least some tumor-infiltrating lymphocytes and eliminate neoplastic cells by antitumor immunity. In this study, we identified that tumors with signature 1 were negatively associated with cytotoxic cell infiltration and antitumor immunogenicity, further supported by less enrichment among genes related to antigen processing presentation and interferon signaling. These findings suggest that the tumor immune escape mechanism of mutational signature 1 may predict the relapse-free survival of patients with TNBC. Moreover, signature 1 activities are negatively correlated with the anti-CTLA-4/PD-1 therapies immune response predictor IPS (determinant of tumor immunogenicity)[25] and neoantigen counts, suggesting that signature 1 may regulate treatment response to immunotherapy.

Comprehensive knowledge of the oncogenes underlying human cancers is a critical foundation for cancer diagnostics, therapeutics, and selection of rational combination therapies.[26] Here, we combined MutSigCV and OncodriveFML algorithms, followed by filter criteria, to re-annotate mutations and identified four unreported novel SMGs in TNBC. KDM6A (lysine demethylase 6A, mutated in 3.0% of patients) is located in the X chromosome and catalyzed the demethylation of tri/dimethylated histone H3.[27] Mutations in KDM6A could promote carcinogenesis and confer drug resistance in tumors.[28,29] T-lymphocyte activation antigen CD86 is a receptor involved in the costimulatory signal essential for T-lymphocyte proliferation and interleukin-2 production[30] and plays a vital role in anti-CTLA-4 blockade therapy immune response.[31] The RNA-binding protein RBM47 has been shown to strongly inhibit tumor progression and metastasis and was found to be substantially downregulated during the epithelial-to-mesenchymal transition in breast cancer. Another study demonstrated that RBM47 functions as a regulator of the p53/p21-signaling axis and plays an important role in antitumor effects. In our study, patients with mutations in RBM47 were associated with poor clinical outcome, probably due to the non-synonymous functional mutation causing gene expression aberration. The ligand-binding chain (alpha) of the gamma interferon receptor encoding gene IFNGR1, acted in IFN-γ pathways, and its mutation was reported to regulate the immune response to PD-1 treatment.[13,32]

Mutations patterns of signature 1 were commonly found in most cancer types (BRCA, LUAD, PRAD, etc.), known as clock-like mutational process properties, which showed a correlation between the number of mutations and age at diagnosis.[6] In this meta-analysis, the age at diagnosis was also correlated with signature 1 mutation activity in TNBCs. Indeed, older patients with breast cancer harbored a less active immune environment than younger patients with breast cancer.[7] Recent studies also noted a cancer type-specific (BRCA, PRAD) negative correlation of immune cell infiltration and mutational signature 1,[33] and further investigation is required to understand how signature 1-related mutated clones escaped from immune recognition.

Recent advances have reported that genomic mutational signatures are associated with clinical prognosis and treatment response. Xing et al. revealed that mutational signature 18 was significantly enriched in tumors with cadherin 1 (CDH1) mutations and associated with poor prognoses in gastric cancer.[34] There is evidence that the APOBEC mutational signature is a potential predictive marker for PD-1 immunotherapy response in NSCLC.[13,35] In our genomic meta-analysis, we identified that signature 1 was associated with shortened survival time in patients with TNBC and suggested a significant association with tumor escape from immunological surveillance. A recent study indicated that short-term chemotherapy (e.g., doxorubicin and cisplatin) may induce a more favorable tumor microenvironment and increase the likelihood of response to PD-1 blockade in TNBC.[36] Therefore, we speculated that patients with signature 1 after induction treatments of doxorubicin and cisplatin may improve the clinical benefit of immune checkpoint inhibitor therapy.

The main limitation of this study was the use of the public datasets from different cohorts, which introduced a heterogeneity in data processing and could have resulted in a batch effect in the final mutation lists. Although we utilized multiple genomic and transcriptomic datasets for analysis, the dataset with both WES and RNA sequencing was only available in TCGA and FUSCC cohorts. As a result, the association between mutational patterns and gene expression, including analysis of immune cell infiltration and oncogenic pathways, needs further validation. In addition, considering that SMG non-synonymous mutations causing amino acid changes in protein function are huge and complicated, we will focus on these areas in future studies.

To our knowledge, this is the largest number of samples of a genomic study that combined multiple WES and RNA-Seq datasets, with the ultimate goal of providing novel genomic-driven therapeutic strategies for TNBC. We identified a putative clock-like mutation signature and molecular biomarkers of response to TNBC treatment, demonstrating the complex interplay of host and tumor in immune surveillance. However, the mechanisms underlying the association between signature 1 and prognostic outcome remain unclear, and further studies are warranted.

Materials and methods

Genomic and clinical data

Somatic mutations and clinical data (a total of 93897 coding somatic mutations in 538 samples) of the aggregated TNBC cohort were curatedfrom recent publications, including 65 samples from BCCRC cohort,[37] 13 from Broad cohort,[38] 14 from Sanger cohort,[39] 279 from FUSCC cohort,[40] and 167 from TCGA (https://portal.gdc.cancer.gov). All patients were female, and the extracted DNA and RNA for sequencing were obtained from primary untreated tumor tissues. We re-annotated all previously called somatic mutations using the Oncotator[41] against the hg19 reference genomics database. Transcriptomic data were obtained from TCGA and FUSCC cohorts. Both peptides resulting from wild-type and mutated sequences for predicted binding affinities scores below 500 nM to patient HLA are defined as neoantigens.[42] Somatic copy number alteration was obtained from Davoli et al.[43] The gene microarray data for 33 paired tumor-normal tissues were downloaded from GEO (ID: GSE76250) published by Jiang et al.[14] Detailed clinical information, including age, histological type, lymph node status, relapse-free survival, and status, was also collected from these studies and is presented in Supplementary Table S2. The overall clinical characteristics with specific signature grouping are summarized in Supplementary Table S3.

Identification of SMGs

We identified SMGs by combining MutSigCV[26] and OncodriveFML[44] algorithms, followed by further filtering procedures. MutSigCV measures the significant enrichment of non-silent somatic mutations in a gene by addressing mutational context-specific background mutation rates. OncodriveFML was designed to identify genes with highly clustered and nonrandomly distributed functional mutations. Candidate novel SMGs were required to meet these criteria: statistically significance in two algorithms (MutSigCV and OncodriveFML, both q < 0.1) and expression in human Cancer Cell Line Encyclopedia (https://portals.broadinstitute.org/ccle)[45] and TCGA pan-cancer dataset.[46]

Deciphering the mutational signature operative in the genome

We used the framework suggested by Kim et al.[47] to extract mutational signatures from aggregated TNBC samples’ (n = 538) genomic data. The method was based on Bayesian variant non-negative matrix factorization and can automatically determine the optimal number of mutational signatures. The mutation portrait matrix was factorized into two nonnegative matrices and , where represents mutational processes and represents the corresponding mutational activities. Mutational activity in every sample refers to the number of somatic mutations attributed to the extracted relevant signature. The number of columns in matrix is the number of mutational signatures. The rows of matrix are the 96 mutational contexts, and its columns are the 538 TNBC samples. The 96 mutational contexts are derived from combinations of 6 mutational types (i.e., C > A, C > G, C > T, T > A, T > C, and T > G) and their 5′- and 3′-adjacent bases. The extracted mutational portrait of TNBC was compared and annotated by cosine similarity analysis against the COSMIC.[8,46]

Extracted signature with TMB

The extracted mutational signatures were stratified as binary variables (i.e., no and yes) in the multivariate model. The classification method is based on the previous method, in which a signature was considered significant if it contributed to more than 100 substitutions or more than 25% of the total mutation activities.[8] The mutational signature activities and designation of each tumor sample were also indicated in Supplementary Table S2. As mutations in BRCA1/2 and ATR and some mutational signatures increase mutation rates in the cancer genome, we recognized them as confounding factors and incorporated them into the generalized linear models and fit proportional hazards regression to analyze associations between mutation signatures and TMB.

GSEA and network analysis

The R package limma[48] was used to evaluate the differential expression of each gene in TNBC samples with different signature 1 activities (no vs yes). Specifically, gene expression data were normalized by voom and then fed to lmFit and eBayes functions in the R limma package. Differential expression statistics were used as input to perform GSEA[49] based on the REACTOME gene set (downloaded from the MSigDB database version 6.3). The fast GSEA algorithm implemented in the Bioconductor R package fgsea was used.

Tumor infiltrating lymphocyte cell analysis

The deconvolution approach CIBERSORT (http://cibersort.stanford.edu/) was used to estimate the abundance of 22 distinct leukocyte subsets[50] with the gene expression profile in TNBC.

Tumor immunogenicity analysis

IPS was used to determine the tumor immunogenicity, which was based on TCIA-recommended method.[25] The scoring scheme was developed from a panel of immune-related genes belonging to four classes: MHC-related molecules, checkpoints or immunomodulators, effector cells (activated CD8 + T cells and CD4 + T cells, Tem CD8+ and Tem CD4+ cells), and suppressor cells (Tregs and MDSCs). For each class, a sample-wise Z score from the gene expression data was calculated. Then, the weighted averaged Z score was calculated by averaging the Z scores within the respective category leading to four values and the sum of the weighted averaged Z score of the four categories was recognized as the IPS.

Quantifying the T cell-inflamed gene expression profiling (GEP) and CYT

T cell-inflamed gene expression profiling (GEP) proposed by Ayers et al.[51] were applied to assess T cell infiltration and quantify the GEP scores. GEP is composed of 18 inflammatory genes associated with chemokine expression, antigen presentation, and adaptive immune resistance. T cell-inflamed GEP scores were calculated as a weighted mean of normalized expression values for these genes. The CYT was calculated as the geometric mean of GZMA and PRF1 transcript levels.

MATH score

The MATH score was calculated as the median absolute deviation of each somatic mutation’s allelic fraction from the median allelic fraction of all mutations in the tumor, divided by the median variant allelic fraction (VAF).[10] We calculated MATH values for TCGA samples from mutant variant allele frequencies (VAF) of all genes.

Statistical analysis

The statistical analysis in this study was generated using R 3.6.1. Quantitative data were presented as median. Continuous variables between groups were compared using the Wilcoxon rank-sum test or paired t-test depending on data distribution. The Spearman correlation coefficient was used to analyze the correlation between two quantitative variables. Kaplan–Meier survival analysis and Cox proportional hazards model were used to analyze the association between mutational signatures and prognosis using the R survival package (survminer 0.4.6). We used stan_lm from the R package rstanarm (version 2.19.1) to perform multivariate Bayesian logistic regression analyses. All comparisons were two-sided with an alpha level of 0.05, and the Benjamini–Hochberg method was applied to control the false discovery rate for multiple hypothesis testing.[52]

Click here for additional data file.

Click here for additional data file.