Systematic identification of the cancer pathways and molecules related with breast cancer immunogenicity

To identify molecular features related to immunogenic activity in breast cancer (BC) and provide new targets and directions for BC immunotherapy, we firstly used ESTIMATE to evaluate the degree of immune cell infiltration of the BC patients in TCGA and METABRIC, and explore the relationship between the degree of immune cell infiltration and prognosis of BC patients. Then, we identified the cancer pathways, proteins and miRNAs related to BC immunogenicity, predicted the target genes of these miRNAs, and identified the pathways related to these target genes with KEGG pathway enrichment analysis. We also explored the correlation between PD-L1 expression level and cancer pathways and found that PD-L1 expression showed a positive association with cancer-related pathways. In this article we have successfully identified several cancer-related pathways, proteins, miRNAs, and their target genes, which provided promising new targets for BC immunotherapy. And PD-L1 blockade therapy may be more effective in BC patients with the activation of some cancer-related pathways.


Introduction
Breast cancer (BC) is the most common cancer among women, as well as the leading cause of cancer-related death in women second only to lung cancer [1]. The current treatment of BC is still limited to surgery, chemotherapy, endocrine therapy, and HER-2-targeted therapy. In recent years, immunotherapy clinical trials have shown promising results for a variety of cancers, including melanoma, nonsmall cell lung cancer (NSCLC), renal cell cancer, Hodgkin's lymphoma, bladder cancer, and head and neck cancer, especially the blockade of immune checkpoints [2]. However, current immunotherapies are only propitious to a subset of cancer patients, and fewer studies have correlated immunotherapeutic strategies and BC treatment. Although BC is not highly responsive to immunotherapy as compared to other malignant diseases, like melanoma and NSCLC, there are growing evidences suggesting the existence of variable immunogenic activity in certain BC subtypes [3].
Recently, atezolizumab in combination with paclitaxel protein bound has been accelerated approved by FDA for the treatment of unresectable adult patients with locally advanced or metastatic triple negative BC (TNBC), which was highly invasive and lacking of treatment targeting hormone or HER2 receptor. Continued approval might be contingent upon verification and description of clinical benefit in confirmatory trials [4]. This indicates a promising prospect for BC immunotherapy in future.
Concerning to immunotherapy, which has benefited only a subset of patients, some patients failed to respond to it and for the rest, there are limited responses followed by tumor progression [5]. Therefore, the systematic identification of predictive molecular biomarkers for BC immunotherapy is crucial. There are evidences indicating that PD-L1, dMMR, MSI-H, tumor mutation burden (TMB), and other biomarkers are related to cancer immunotherapy [6][7][8][9]. Immunohistochemistry (IHC) positivity of PD-L1 in tumor cells or immune cells is predictive of better response of anti-PD-(L)1 therapy [7]. A clinical phase 2 study conducted by Le et al. [8] evaluated the clinical activity of pembrolizumab, which indicated tumors with dMMR are more responsive to immunotherapy. Snyder et al. [9] have found that high 79 Page 2 of 12 mutational load in melanoma is associated with a sustained clinical benefit from ipilimumab-or tremelimumab-based CTLA-4 inhibition. Notably, the data supporting these biomarkers are mainly related to lung cancer, malignant melanoma, and colorectal cancer, and those for BC cases are still quite limited. In addition, even as the most widely accepted biomarker, there are also up to 20% patients with PD-L1negative tumors benefit from immune therapies [8,10,11]. This suggests that further screening of biomarkers for accurate prediction of clinical response is of great importance in BC immunotherapy.
In this study, we successfully identified the cancer-related signaling pathways, proteins, and miRNAs that are closely related to the immunogenicity of BC, predicted the target genes regulated by these miRNAs, and explored the relationship between PD-L1 expression level and these pathways, which provided a foundation for subsequent study on BC immunotherapy.

Activity evaluation of immune cell types and functions
We used the ssGSEA score to quantify the activity or enrichment level of different immune cells, immune functions, and cancer-related pathways in tumor samples. A higher ssGSEA score represents a higher activity. In addition, we used ESTIMATE algorithm to estimate the level of immune infiltration within tumor tissue [12]. Based on gene expression profile, ESTIMATE outputted an immune score for each sample to quantify the degree of immune cell infiltration in tumor tissue.

Prediction of microRNA target gene
Firstly, TargetScan 7.2 was used to obtain the target genes predicted by miRNA. Then miRWalk2.0, miRtarBase, miRmap, and miRPathDB were used. The target genes simultaneously predicted by all the 5 databases were selected as the target genes regulated by miRNA. Finally, Cytoscape was used to construct the regulatory network between the target genes.

Survival analysis
Survival analysis of BC patients was performed based on immune score data. The Kaplan-Meier survival curve was used to represent the differences in OS and DFS between patients with high or low immune score. The median immune score was used to identify patients with high or low immune scores. The log-rank test was used to calculate the difference of survival time between the two groups of patients, and P < 0.05 represented a significant difference.

Statistics and analysis
Pearson's correlation was used to calculate the correlation between immunity and the cancer-related pathway (ssGSEA score). Spearman's correlation was used to calculate the correlation between gene expression level and cancer-related pathway enrichment level, as well as the correlation between immune score and the level of cancer-related pathway enrichment, protein expression, or miRNA expression. The error discovery rate (FDR) calculated using Benjamini and Hochberg (BH) methods was used to rectify the P value. FDR < 0.05 was considered statistically significant.

The degree of immune Cell infiltration is positively correlated with the survival prognosis of BC Patients
In the TCGA and METABRIC data set, we compared the survival of BC patients with different levels of immune cell infiltration (high or low immune score). Both of the outcomes of the two data sets showed that BC patients with high immune score had a significantly better overall survival (OS) (sequential test P < 0.05, Fig. 1), which indicates that BC patient with high immunity infiltration has better clinical outcomes.

The cancer-related pathway is positively correlated with the degree of immune Cell infiltration and immunogenicity
Using KEGG pathway enrichment analysis, we identified 23 cancer-related pathways: NOD-like receptor, Toll-like receptor, JAK-STAT, MAPK, PI3K-AKT, apoptosis, Focal adhesion, Ca2 + , Wnt, VEGF, HIF-1, TNF, Hedgehog, Notch, ECM-like receptor, NF-κB, TGF-β, glycolysis, ErbB, estrogen receptor, MMR, cell cycle, and p53 signaling pathway. We compared the relationship between the activity of cancer-related pathways and the infiltration degree of immune cells in TCGA and METABRIC data sets. The results showed that in the two data sets, almost all cancer pathways were positively correlated with the infiltration degree of immune cell, only Hedgehog and Notch were negatively correlated (Spearman's correlation, FDR < 0.05, Fig. 2a). Then we further compared the correlation between the activity of cancer-related pathways and the 8 immune characteristics. Being consistent with the degree of immune cell infiltration, almost all of the cancer-related pathways were positively correlated with the activity of the 8 immune characteristics, while Hedgehog and Notch were negatively associated (Pearson's correlation, FDR < 0.05, Fig. 2b).

Key proteins associated with immune score and immunogenicity
Eight proteins were identified being positively correlated with the immunity score: LCK, cleaved_capase-7, Syk, Axl, PI3K, Annexin-1, PKC-α, and ATM. Five proteins were negatively correlated with immune score: ER-α, E-Cadherin, GATA3, HER3, and Claudin-7 (Spearman's correlation, |R|> 0.25, FDR < 0.05, Fig. 3a). In further analysis, the correlation between these proteins and the activities of the 8 immune characteristics were almost the same (Pearson's correlation, FDR < 0.05, Fig. 3b). These results suggest that the expression of these up-or down-regulated key proteins is associated with the increase or decrease of immunogenicity in BC.
The correlation between microRNA (miRNA) and immune score miRNA is a conserved noncoding RNA segment of 22-25 nucleotide sequences. By targeting specific mRNA, miRNA leads to mRNA degradation or inhibits the mRNA translation. As a result, it regulates gene expression at the transcriptional level. This kind of regulation is critical to cellular function. More than one-third of mRNAs might be the target of miRNAs [11]. Recent studies have shown that miRNA also plays a key role in the regulation of immune function, including innate and adaptive immune responses, development and differentiation of immune cells, as well as prevention of autoimmunity [13].
We further used TargetScan 7.2 to predict target genes of these 4 miRNAs which had a correlation coefficient of immune score above 0.5: hsa-mir-155, hsa-mir-150, hsa-mir-146a, and hsa-mir-223. In order to reduce falsepositive rate and ensure the reliability, miRWalk 2.0, miRPathDB, miRTarBase, and miRNAMap were simultaneously used for prediction. The intersection of above databases and software was taken and the results are shown in Table 1.
We further compared the correlation between these target genes and the infiltration degree of immune cells (Spearman's correlation, FDR < 0.05), and analyzed the regulation network of the selected target genes (Fig. 5a).
Finally, we performed KEGG pathway enrichment analysis on these genes. The results showed these genes were involved in multiple cancer-related pathways in BC, including NOD-like receptor, Toll-like receptor, JAK-STAT, MAPK, apoptosis, Focal adhesion, Ca2 + , Wnt, VEGF, Hedgehog, Notch, TGF-β, glycolytic, ErbB, cell cycle, and p53 signaling pathway. These genes were also involved in several immune-relating pathways, including B cell receptor signaling pathway, T cell receptor signaling pathway, chemokine signaling pathway, NK cell-mediated cytotoxicity signaling pathway, and cytokine-cytokine receptor interaction pathway (Fig. 5b). This indicates that miRNA plays an important regulatory role in BC as well as the immunotherapy of BC, which may be achieved by regulating the target genes.

Correlation between PD-L1 and the activity of cancer-related pathway
The expression of PD-L1 on tumor cells plays a key role in tumor immune escape [14]. We analyzed the correlation between PD-L1 expression and the activity of cancerrelated pathway. It was found that in the two data sets the expression level of PD-L1 was significantly and positively correlated with the activity of 11 cancer-related pathways (NOD-like receptor, Toll-like receptor, NF-κB, JAK-STAT, TNF, cell apoptosis, HIF-1, VEGF, ErbB, and p53 pathway), while negatively correlated with the activity of Notch pathway (Spearman's correlation, P < 0.05) (Fig. 6). This reveals that PD-1 or PD-L1 blocking therapy may be more effective for BC subtypes in which these cancer pathways were activated.

Discussion
In order to identify the molecular characteristics associated with immunogenicity of BC, we firstly compared the survival of BC patients with different degrees of immune cell infiltration. The results showed that BC patients with high immune infiltration had better clinical prognosis.
Then we identified 23 common cancer-related pathways in BC from KEGG, and analyzed the correlation between cancer-related pathway activity and the degree of immune cell infiltration or 8 immune signatures. We found that almost all of these cancer-related pathways were positively correlated with the degree of immune cell infiltration, while Hedgehog signaling pathway and Notch signaling pathway were negatively correlated.
Otsuka A et al. has shown that inhibition of the Hedgehog signaling pathway could change tumor microenvironment, transform the chemokine and cytokine network, increase the amount of CD8 + cytotoxic T cells invading tumor cells, and promote adaptive immune response of basal cell carcinoma [15]. Notch signal negatively regulates the inflammatory response that activated by Tolllike receptor (TLR), an important molecule in nonspecific immunity, which plays a connecting role in nonspecific immunity and specific immunity [16]. Therefore, inhibition of Hedgehog or Notch pathway may enhance the efficacy of immunotherapy in BC, indicating that the combination of Hedgehog or Notch pathway inhibitor and immunotherapy might be a potential solution for the treatment of BC.
At the protein level, we identified 13 proteins that were most relevant to immune score and immune signatures. Among them, 8 were identified being positively correlated with the immunity score--LCK, cleaved_capase-7, Syk, Axl, PI3K, Annexin-1, PKC-α, and ATM--while the other 5 proteins were negatively correlated--ER-α, E-Cadherin, GATA3, HER3, and Claudin-7. Among the first 8 proteins, LCK is one of the molecules from protein tyrosine kinase (SRC) family that plays a key role in the T cell receptor (TCR) signaling pathway. There had been study shown that the inhibitory effect of PD-1 on TCR signal molecule activation is significantly enhanced by PD-1-mediated LCK inhibition [17]. Capase-7 is a key protein in the cell apoptotic pathway, which regulates immune response by mediating inflammation [18]. As an adaptive immune receptor signal, the role of Syk is crucial in innate immune recognition [19]. Axl overexpresses in a variety of tumor tissues and cells. Wei et al. [20] combined Axl with chimeric antigen receptor T cells, namely, AXL-CAR-T, and found that it exhibits obvious anti-tumor effect in TNBC xenotransplantation model, implying Axl could be a potential therapeutic target for TNBC [20]. Annexin-1, a 1 3 79 Page 6 of 12 B A significant anti-inflammatory protein controlling the interleukin-mediated immune response [21], has both innate and adaptive immune function and mainly participates in inflammatory cell regulation. As adherence junction protein, E-cadherin adjusts a number of intracellular signal transduction pathways. Recent studies have revealed that E-cadherin affects immunity via regulating mononuclear macrophages [22]. Down-regulation of Claudin-7 induces metastasis and invasion of colorectal cancer by promoting epithelial-mesenchymal transition [23]. The zinc finger transcription factor GATA-3, a major regulator of T helper 2 cell (Th2) differentiation, modulates the expression of IL-4, IL-5, and IL-13 [24]. GATA-3 facilitates type II immunity by adjusting the development and function of type 2 innate lymphocyte cells (ILC-2) [25]. HER3 is crucial heterodimer ligand of EGFR and HER2 and its high expression accelerates the invasion of BC. Therefore, HER3 inhibition could be widely used in the treatment of EGFR-and HER2-driven tumors. Osada T et al. [26] reported that HER3-targeted vaccine activates HER3-specific T cells and induces anti-HER3-specific antibodies, which changes both the T cell infiltration inside tumor and the response to immune checkpoint inhibition. The expression of these up-or down-regulated protein molecules we identified is closely related to the increase or decrease of immunogenicity in BC. On this basis, we can provide a direction for immunotherapy of BC from the protein level.
miRNA is a group of small ultra-conserved noncoding RNAs that regulates gene expression at transcriptional level. They target specific mRNA and consequently induce their degradation or inhibit their translation. By this way, the miRNA plays a quite important role in immunization, cell differentiation, tumor genesis, and cell death [27]. We identified miRNAs that are significantly associated with immunity based on Spearman's analysis. The result indicated that the expression levels of these up-or down-regulated miRNAs are also related to the increase or decrease of immunogenicity in BC. It has been reported that the serum expression of miR-155 has increased more than 2 times in BC patients, which is positively correlated with the clinical stage and Ki-67 expression, while negatively correlated with p53 status [28]. It also regulates the production of mature B cells [11]. Another important miRNA, miR-150, targets at NOTCH3 and impacts both the physiology and development of T cell [29]. miR-146a inhibits the proliferation and metastasis of TNBC cells by regulating SOX5 [30], and also adjusts immune response in the microenvironment of melanoma [31]. miR-223 regulates immune response via adjusting granulocyte production and inhibiting of macrophage differentiation [32]. Subsequently, we predicted the target genes that regulated by the 4 miRNAs with a correlation coefficient greater than 0.5 through prediction software and database, and then analyzed the regulatory network and enrichment of KEGG pathway of these target genes. Regulatory network analysis showed these target genes greatly correlated with immunogenicity. KEGG pathway enrichment analysis showed that these genes involved in BC-relating signal pathways that identified by us, and also participated in several immune-relating pathways. These results indicate that the miRNAs we identified play a crucial regulatory role in BC and its immunotherapy, the basis of which might be regulating the target genes we have identified. Therefore, these miRNAs have large potential to be a therapeutic target for BC, suggesting a promising direction for immunotherapy of BC at the miRNA level.
The expression of PD-L1 is of great importance in tumor immune escape [14]. We also analyzed the correlation between PD-L1 expression and cancer pathway activity. The results showed that the expression level of PD-L1 was positively correlated with the activities of 11 cancer pathways (NOD-like receptor, Toll-like receptor, NF-κB, JAK-STAT, TNF, apoptosis, HIF-1, VEGF, ErbB, and p53 pathway). We have already found these cancer pathways were positively correlated with the degree of immune cell infiltration or their immunogenicity, indicating the activation of these cancer pathways will strengthen the immunogenicity in BC. However, the expression level of PD-L1 was negatively correlated with the activity of NOTCH pathway. Considering our previous finding that inhibition of Notch pathway activity enhances the immunogenicity of BC, we assume that the increased activity of NODlike receptor, Toll-like receptor, NF-κB, JAK-STAT, TNF, cell apoptosis, HIF-1, VEGF, ErbB, p53 pathway, and the inhibited activity of Notch pathway may make BC patients more sensitive to PD-1 / PD-L1 blocking treatment.
(b) KEGG pathway (cancer-related and immune-related) enrichment analysis of the genes targeted by the miRNAs whose expression is significantly associated with immune signatures