Differential Expression of Necroptosis-Related Genes with regulating lncRNA in Triple Negative Breast Cancer to Predict Prognosis and Immunotherapy Response

Differential Expression of Necroptosis-Related


Introduction
According to the global cancer statistic GLOBOCAN published in 2020, female breast cancer (BC) was the most common malignancy, with approximately 2.3 million newly diagnosed cases (11.7%).It took the fourth place (6.9%) regarding cause of cancer-specific death of female patients globally [1].Triple-negative BC (TNBC) is one of four main subtypes of BC.Different from the other three types, TNBC lacks of expression of estrogen and progesterone receptors (ER, PR) as well as the expression of human epidermal growth factor receptor 2 (HER-2).This characteristic makes it the most aggressive and refractory BC subtype, especially in younger patients [2,3].Traditional endocrine treatment (tamoxifen and aromatase inhibitors) and targeting anti-HER-2 therapy trastuzumab are ineffective in patients with TNBC [4].Although neoadjuvant systemic therapy was used in earlystage TNBC to improve survival, a series of age-related adverse events occurred, including cognitive damage and premature ovarian failure, severely affecting life quality [5][6][7].What's worse, heterogeneity of TNBC with unknown mechanism of carcinogenesis prompts more individualized anticancer treatments [4].
Necroptosis is termed as a kind of caspase-independent, regulated, and programmed cell death [8].Different from apoptosis, necroptosis is featured as plasma membrane rupture, followed by releasing intracellular components which stimulating the immune system [9].Necroptosis is associated with pathology of various types of diseases, including non-tumor and tumor diseases.The former part includes cardiovascular, cerebrovascular, neurodegenerative diseases, a series of infectious diseases, and some special types of ischemic morbidities [8,[10][11][12].The latter part consists of hematological and solid malignancies [13][14][15][16][17]. Receptor-interacting serine/threonine-protein kinase 1/3 (RIPK1/3) and mixed lineage kinase domain-like proteins (MLKL) are pivotal molecules in the mechanism of necroptosis in cancers [18][19].Via a cascade of signaling pathways involving interferon receptors (IFNR), T cell receptor (TCR), tumor necrosis factor superfamily (TNFSF), Toll-like receptors (TLR), and various types of physicochemical agents, RIPK1 was activated, which is the process of necroptosis induction [18,[20][21][22].In tumor biology, necroptosis serves as a double-edged sword.On one hand, it is capable of inhibiting proliferation, aggression, and metastasis of malignancy; on the other hand, it stimulates production of reactive oxygen species (ROS) and maintains tumor microenvironment (TME) [23][24][25].Bench studies have found that the function and survival of luminal or non-luminal BC cells could be impacted by some specific chemical elements which were capable of regulating necroptosis [26][27][28].Hence, regulation of necroptosis could be a brand-new approach for anti-TNBC therapy.
Long non-coding RNAs (lncRNAs) are a large family of endogenous RNAs specifically enriched in exosomes and each lncRNA has more than 200 nucleotides [29].LncRNAs do not code for protein production [30].LncRNAs in human tissues have been increasingly explored via bioinformatics and secondgeneration sequencing techniques [31].In TME, lncRNAs regulate gene expression, affecting variety of biological process of cancer cells including proliferation, differentiation, migration, cell death, invasion and metastasis, and drug resistance [30,32,33].Previous studies found that lncRNA could also indirectly affect necroptosis by different pathways.In a study from Germany, researchers found that an intergenic lncRNA, Linc00176, exerted a protective function for hepatocellular (HCC) cells by maintaining the capability of proliferation and survival of cancer cells [34].In a non-tumor model, scientists have shown that a lncRNA called necrosis-related factor (NRF) was capable of targeting RIPK1/ RIPK3 and miR-873 to control necroptosis.miR-873 inhibited necrotic death of cardiomyocytes via inhibiting protein production of RIPK1/RIPK3 [35].Nevertheless, whether lncRNAs contribute to TNBC invasion and proliferation is still unclear.Further intense investigation of necroptosis-related lncRNAs is necessary for more precise and individualized management for TNBC.In this study, we aim to explore remarkable necroptosis-related genes (NRGs) of TNBC as well as related lncRNAs.Subsequent risk model was built based on these genes to predict prognosis of BC patients.

Data downloading and pre-processing
Transcriptome RNA-seq datasets [HTseq-Counts and Fragments Per Kilobase of exon model per Million mapped fragments (FPKM)] of female BC were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/).Corresponding clinical data were acquired from cbioportal (https:// www.cbioportal.org/).Gene expression matrices from HTseq-Count FPKM files were created via the strawberry perl algorithm (perl 5, version 30).The Count matrix was used to filter the differentially expressed coding genes between TNBC and NTNBC samples, while the FPKM matrix was used for identification of DEG-correlated lncRNAs as well as other analysis.After excluding 12 male samples, 112 normal samples, and 183 samples with unknown state of BC immunohistochemistry (ER, PR, HER-2), 901 samples were used for DEG analysis.For prognosis and immune response studies, 779 female BC samples without unclear clinicopathological information (age, race, staging, subtype, surgery, survival time, etc) were finally included (Supplementary Figure 1).

Prognosis-associated lncRNAs and establishment of a prognostic prediction model
We merged the clinical data and data of expression of necroptosis-related lncRNAs.The merged file had 864 BC samples.Then we divided these samples into two subgroups with 1:1 ratio.One group was called as "training set" and the other one was called as "testing set".In the training set, univariate Cox regression analysis was preliminarily performed to select prognosticclassified lncRNAs with p value <0.1.Least absolute shrinkage and selection operator (LASSO) regression was subsequently performed to filter prognosis-related lncRNAs with 10-fold cross-validation and 1,000-cycle running as well as a p value of 0.05.Then the prediction model was created by multivariate Cox regression analysis in the training set, and then the model was used to predict the risk of each sample in the training and testing set, respectively.The risk score was calculated by the following formula: Risk score= (1), in which "expr" stands for the expression of lncRNA in each sample while "coef" stands for the coefficient of lncRNAs correlated with survival.Median risk scores in both sets were calculated to evaluate the risk of each sample as either "high" or "low" risk so that the samples would be stratified into low-risk or high-risk groups.Clinicopathological variables including age, race, subtype, and type of surgery were evaluated by multi-Cox regression analysis to select potential predictors for nomogram preparation.Receiver operating characteristic (ROC) curves were created to measure the quality of prediction of model.The above analysis was completed via R packages "survival", "rms", "foreign", "caret", "glmnet", "survminer", and "timeROC".Survival analysis between high-and low-risk groups was conducted with R packages "survival" and "survminer".

Construction of nomogram and calibrated plots
Nomogram was produced according to the results of multivariate Cox regression analysis.Clinicopathological features with p value <0.05 were included in the nomogram.Based on each point derived from every variable, a total point of each sample could be acquired so that the nomogram could predict 1-year, 3-year, 5-year, and 10-year survival rates.To evaluate the predictive accuracy of the nomogram, 1-year, 3-year, 5-year, and 10-year calibrated plots were illustrated.R packages "rms", "survival", and "foreign" were used in this analysis.

GSEA analysis
According to high-risk and low-risk group in all 864 available BC samples, we classified differentially expressed KEGG pathways into different groups through Gene Set Enrichment Analysis software (GSEA, version 4.1.0,Broad Institute, Inc., Massachusetts Institute of Technology, and Regents of the University of California, USA) [36][37][38][39][40]. Amongst, 91 biological functions and pathways were active in high-risk group, while 87 biological functions and pathways were active in low-risk group.We selected each five significantly enriched biological process and pathways in high-and low-risk group, respectively, based on p value< 0.05 and FDR< 0.25.

Comparison of immune landscape between high-risk and lowrisk groups
For correlation between immune cell infiltration and risk score, we firstly downloaded the immune cell infiltration file from TIMER (http://timer.cistrome.org/)[41][42][43].This file includes data of tumor-infiltrating immune cells derived from several types of algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, and Cibersort).Using Spearman correlation test, we calculated correlations of different immune cells with different risk scores, which is illustrated in a bubble chart.Immune-cell infiltration was furtherly assessed by single sample gene set enrichment analysis (ssGSEA) and major immune functions were also evaluated between high-risk and low-risk group.ssGSEA and its visualization was performed with the help of R packages "GSVA", "limma", "GSEABase", "ggpubr", and "reshape2".

TME and immunotherapy response between high-risk and low-risk groups
Based on Estimate algorithm, we evaluated the TME status by calculating immune score, stromal score, and ESTIMATE score (immune score + stromal score).Then we compared these scores between high-and low-risk groups [44].We collected immune checkpoint molecules from ImmPort website (https://www.immport.org/shared/genelists)and extracted their expression levels from FPKM file [45,46].Then we compared their expressions between patients with different risk levels.Finally, potential efficacy of immunotherapy [anti-programmed death-1 (PD-1) and anti-cytotoxic T lymphocyte-associated antigen-4 (CTLA-4) monoclonal antibodies] in the two groups were evaluated by using data of immunophenoscore (IPS) from the Cancer Immunome Atlas (TCIA, https://tcia.at/home)[47,48].IPS is an effective predictor of treatment response to immune checkpoint inhibitor (ICI), anti-CTLA-4/anti-PD-1 antibodies.Its calculation is based on four categories of immunogenicity: expression of effector cells (CD4 + T cells, activated CD8 + T cells, and memory CD4/ CD8 + T cells), immunosuppressive cells (regulatory T cells and myeloid-derived suppressive cells), major histocompatibility complex (MHC) molecules, and immunomodulators.Previous results showed that IPS was associated with several kinds of solid malignancies including BC [49].
The complete process of statistical analysis was performed by R software 4.0.4(R Foundation for Statistical Computing, Vienna, Austria).

Construction of a prognostic prediction model for risk assessment
The complete process of model construction was conducted in the training set (N=432).Twelve NRG-related lncRNAs were significantly associated with prognosis according to univariate Cox regression analysis (Figure 2A).Seven of them were indicative of poor prognosis while the other five suggested a protective element of prognosis.Ten lncRNAs with significantly association with prognosis were acquired by LASSO analysis (Figure 2B-2C).Ultimate six prognosis-related lncRNAs were resulted from multivariate Cox analysis, in which expression of three lncRNAs [Z68871.1,SPRY4-AS1, AC011247.1,hazard ratio (HR)>1] were high-risk elements for prognosis and the other three played a protective role in prognosis (COL4A2-AS1, CCDC144NL-AS1, ADAMTS9-AS2, HR<1) (Figure 2D).Based on the results of multi-Cox regression, we obtained a risk score of each sample in the training set using the formula below: Then we used this formula to calculate the risk score of samples in the testing set.In both groups, samples were stratified into highrisk and low-risk levels according to the median risk score (Supplementary Figure 2A-C).Survival analysis suggested that patients in high-risk group had a significant poorer overall survival than those of low-risk group, no matter in the training or testing set (Figure 3A-3C, Supplementary Figure 2D-2F).Heatmap reaffirmed the results that lncRNAs Z68871.1,SPRY4-AS1, and AC011247.1 were highly expressed in high-risk groups while the other three lncRNAs COL4A2-AS1, CCDC144NL-AS1, and ADAMTS9-AS2 were highly expressed in low-risk groups, no matter in training, testing or the complete datasets, suggesting their prognostic meaning (Supplementary Figure 2G-2I).The association of clinicopathological factors other than risk scores with overall survival was analyzed and the results showed that risk score together with age, tumor stage, subtype, and type of surgery were independent risk factors of overall survival.Multi-Cox regression analysis demonstrated that older age, higher staging of BC, TNBC, and more invasive operation of BC were indicative of a poorer prognosis (Figure 3D-3E).ROC curves illustrating the area under curve (AUC) in terms of 1-, 3-, 5-, and 10-year were 0.718, 0.623, 0.670, and 0.720 (Figure 4A).Regarding 10-year ROC curves of the model, the AUC of the variable "risk score" was larger than any other clinicopathological variables (AUC=0.720),illustrating a relatively powerful predictive capability (Figure 4B).The HR of risk score in uni-and multivariate Cox model was 1.208 and 1.236, respectively (both p value<0.0001).

Nomogram construction and calibration plot assessment
Based on the results of the multivariate Cox model for survival analysis, we brought variables "risk score", "age", "stage", "subtype", and "type of surgery" into the nomogram, a visualization tool for prediction of 1-year, 3-year, 5-year, and 10-year survival rate (Figure 5A).Calibration plot demonstrates that the nomogram exhibited an excellent concordance with prediction (Figure 5B).

Biological pathway analysis
Through GSEA analysis, we summarized different KEGG biological process and pathways in high-risk (N=417) and low-risk (N=447) groups.Each five of significantly enriched biological functions and/or pathways were collected from high-and low-risk groups, respectively.Amongst, high-risk groups included DNA mismatch repair, DNA replication, nucleotide excision repair, RNA degradation, and cell cycle regulation; low-risk groups included cytokine receptor interaction, calcium signaling pathway, TGF-β signaling pathway, regulation of actin cytoskeleton, and cancer pathway (Figure 6).

Comparison of immune signature, TME between high-and low-risk groups
Through ssGSEA, significant discrepancy of immune cell abundance could be found between high-and low-risk groups, especially cytotoxic T cells, B cells, natural killer cells (NK cells), dendritic cells (DCs), and neutrophils.Low-risk group had a higher immune cell infiltration than its counterpart (Figure 7A), which is in line with the result of the activity of major immune functions (cytolytic activity, T-cell costimulation, and interferon response, etc.) between the two groups (Figure 7B).Correlation analysis between immune cell infiltration and risk score was also performed and data were derived from several algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, and CIBERSORT, Figure 8A).Take TIMER algorithm for instance, scatter plots demonstrate that T cells (CD4 + , CD8 + ), B cells, DCs, and neutrophils were all negatively correlated with the risk score which corroborated the results of ssGSEA above (Figure 8B-8F).

Comparison of immune checkpoint molecules and potential immunotherapy response between different risk levels
47 immune checkpoint molecules were collected and the expression of these genes were extracted from the original FPKM file of female BC samples.We compared the expression level of these genes between high-and low-risk group, which is demonstrated in Figure 10.Results showed that expression of genes in low-risk group were higher than that of high-risk group except CD80.IPS of patients in high/low-risk groups were obtained.Score 2 (IPS-PD-1/PD-L1/PD-L2 blocker), score 3 (IPS-CTLA-4 blocker), and score 4 (IPS-PD-1/PD-L1/PD-L2 and CTLA-4 blocker) in low-risk group were significantly higher than those in high-risk group (Figure 11A-11D), suggesting a better expected response of ICIs (anti-PD-1/anti-CTLA-4 monoclonal antibodies) in the former group.The results were consistent with the counterpart of immune landscape and TME mentioned in the previous section.

Discussion
In our study, seven DEGs were found related to necroptosis between NTNBC group and TNBC group.Amongst, PLA2G4E and HDAC9 were upregulated in TNBC samples while H2AC 4/12/14/20 were downregulated in TNBC samples.PLA2G4E belongs to the cytosolic phospholipase A2 group IV family.This group of genes regulate membrane tubule-mediated transport.Proteins produced by this gene get involved in recreating tubules to move specific cargo molecules back to cell surface (see details in https://www.ncbi.nlm.nih.gov/gene/123745).Previous studies have found that PLA2G4E was related to necroptosis.Inhibition of PLA2G4E by microRNA Mir504-5p reduced necroptosis of skin and prolonged survival of skin cells [50].Previous studies have found that inactivated function of RIPK1, RIPK3, and MLKL, especially by the necroptosis inhibitor Nec-1, protected neurons from being damaged, degenerated, and keeping the viability of neural cells in these neurodegenerative diseases [51][52][53][54][55].In an animal model of Alzheimer´s disease, PLA2G4E served as a protective role for the survival of cerebral tissue.Overexpression of PLA2G4E could resume memory and cognition, indicating a potential role in treatment of neural degenerative disease [56].It is worthwhile to further explore relationship between this gene and necroptosis in neurodegenerative diseases.HDAC9 is a member of the class IIa family of histone deacetylases.This gene family get involved Volume 6; Issue 01 in gene silencing by modifying chromatin into a transcriptionally repressed conformation.HDAC9 are expressed in many human tissues and organs, especially in brain and skeletal muscles [57,58].HDAC9 negatively regulates muscle differentiation by interacting with myocyte enhancer factor (MEF2) and inhibiting its transcriptional activity [59].Similar to PLA2G4E, HDAC9 also plays a protective role in the nervous system.One study has found that reduced expression of HDAC9 caused death of neural cells [60].HDAC9 possesses several functions in carcinogenesis as well.It has been found to be upregulated in various types of tumors including BC, hepatocellular carcinoma, pancreatic cancer, osteosarcoma, and hematological malignancies [4,[61][62][63][64][65].In BC, researchers found that overexpression of HDAC9 was associated with metastasis and poor prognosis [62].HDAC9 upregulated NF-κB expression and its downstream inflammatory pathways, stabilizing TME and leading to recurrence of BC [61].Another study led by Salgado et.al found that HDAC9 was overexpressed in TNBC, which was consistent with our results [66].HDAC contributed to increased invasion and angiogenesis in TNBC tissues via VEGF and MAPK3 signaling pathway, which is negatively affected by miR-206.A bench study from China using genome-wide CRISPR screening found that HDAC9 was a pivotal molecule in regulating paclitaxel-resistance of TNBC.Inhibition of HDAC9 caused tubulin instability and cell cycle arrest, leading to cancer cell death and increasing the sensitivity to chemotherapy [67].Other related studies also found that HDAC9 suppressed the expression of Erα and endowed antiestrogen-resistance to BC cells [68].These results indicated a potential HDAC-targeting therapy in the treatment of TNBC to suppress its aggression and improve the sensitivity of TNBC to systemic therapy.Different from the upregulated genes, the significantly downregulated NRGs in TNBC samples were all H2A clustered genes.Proteins encoded by these genes were all belong to replication-dependent histone H2A family (see https://www.ncbi.nlm.nih.gov/gene/8338).A study led by Su found that H2AC was upregulated in ER-positive BC.Suppression of this replication-dependent histone H2A isotype caused abnormal estrogen signaling and cell cycle arrest of BC [69].In another study, the authors also furtherly found that H2AC was related to telomeres in human cells and interacted with telomeres protein 1 (POT1) as well as telomere repeat factor 2 (TRF2).Depletion of H2AC would dysfunction the telomere, destabilize chromosome structure, and cause DNA damage [69].
According to the expression and coefficient of lncRNAs, we constructed the model and measured the risk scores for each sample in the training set and testing set.Samples in both sets were stratified into two subgroups: high-risk group and low-risk group.Survival analysis in both sets were similar that patients in high-risk group had a worse prognosis compared with low-risk counterparts.In addition to risk score, age, tumor stage, subtype, and type of surgery were all independent variables for prognosis.Nomogram including these variables was created to predict 1-, 3-, 5-, and 10-year survival rates for BC patients.Predictive efficacy was testified by ROC curves and calibration plots.Results showed that risk score was the best independent predictor of overall survival.
Based on the results of GSEA, we found that significantly enriched biological pathways were focused on cell cycle, cellular metabolism, and cellular proliferation as well as other special pathways including TGF-β and Ca2 + signaling pathways.These are important pathways that necroptosis gets involved in.For instance, tumor-induced endothelial necroptosis triggered tumor cell extravasation and metastasis.TGF-β-activated kinase 1 (TAK1) could inhibit necroptosis of endothelial cells [82,83].Studies also found that dysfunction or deficiency of TAK1 enhanced the necroptosis in which RIPK3 was upregulated [84].In another study, scientists found that thymoquinone, an apoptosis-inducing agent with anticancer activity, could also induce necroptosis through elevating cytosolic calcium concentration by endoplasmic reticulum calcium depletion and activation of store-operated calcium entry (SOCE) in lymphoma cell [85].
We compared the immune signatures of BC samples between high-risk and low-risk group.Results showed that low-risk group had higher level of immune cell infiltration, higher immune, and stromal scores, and higher expression level of immune checkpoint molecules.However, expression of CD80 was lower in low-risk group than that in high-risk group.Also named as B7-1, CD80 could not only bind to immune-stimulating signaling molecule CD28 but also combine with immunosuppressive molecule CTLA-4 [86].As reported, the affinity of CD80 to CTLA-4 was higher than that to CD28, providing its immune-inhibiting feature [87].A study led by Arutha Kulasinghe found that patients with refractory TNBC had a higher level of CD80 [88], corroborating with our results.More investigation should be performed to validate this conclusion.Correlation analysis illustrates that the abundance of immune cells (DCs, innate, and adaptive immune cells) was all inversely proportional to the risk score.The above results were consistent with the immunotherapy response prediction that lowrisk group had a potentially better response to PD-1 and CTLA-4 monoclonal antibodies.Necroptosis not only directly eliminates tumor cells by promoting cell death, but also indirectly suppress tumor growth and proliferation by inducing immune response in TME [89].During necroptosis, rupture of the plasma membrane led to outpouring damage associated molecular patterns (DAMPs) into intercellular space, then necroptotic cells recruited immune and inflammatory cells to clear the cell debris and helped tissue repair [89,90].Animal studies showed that CD8 + T cells and NKT cells played major roles in necroptosis-inducing immune response [91,92].Necroptotic cells offered tumor antigens and inflammatory stimuli to antigen-presenting cells (APC), which then activated CD8 + T cells by antigen cross-priming [93].CD8 + T cells exerted antitumor cytotoxic immunity by producing and secreting multiple cytokines.The above biological process was highly dependent on RIP1 and NF-κB signaling [91].Different from CD8 + T cells, NKT-induced immune responses were RIPK3dependent, of which downstream signaling mediated by pGAM5 and Drp1 [92].Counterintuitively, some researchers suggested that inflammation derived from necroptosis was a pro-cancer factor in TME [94].In the process of carcinogenesis, inflammation could promote tumor cell growth, proliferation, aggressiveness, invasion, and metastasis by providing bioactive molecules [90,95,96].ROS and other mutagenic chemicals also facilitated evolution of tumor cells by genetic alterations [90].Given that the above conclusions were controversial, the contribution of immune and inflammatory responses induced by necroptosis in tumor prevention and development needs to be further elucidated.
Despite a comprehensive analysis of necroptosis in TNBC and its relationship with prognosis and immune landscape, some limitations still existed in our study.Firstly, the results were not further validated in basic studies (in vitro and in vivo).The biological functions of prognosis-related lncRNAs were needed to be clarified.Secondly, different extents of biases during data processing and case inclusion/exclusion were inevitable because of a retrospective feature of this study.Thirdly, the expression of PD-1 (CD279) and PD-L1 (CD274) were not included in the original transcriptome data file of TCGA-BRCA.Hence, we failed to directly compare the expression level of these two biomarkers between high-and low-risk groups though data of anti-PD-1 and PD-L1 response could be retrieved from TCIA database.
To conclude, seven differentially expressed NRGs were selected between NTNBC and TNBC groups.Six NRG-related lncRNAs were significantly associated with prognosis and were brought into risk model construction.Risk score is an excellent predictor for overall survival of BC and the model was of high predictive quality.Compared with low-risk group, high-risk group has a gloomier prognosis, less abundance immune cell infiltration, and less responsive ICI treatment.The results pave way for further exploration of precision therapy for TNBC.

Figure 1 :
Figure 1: Network between differentially expressed NRG and NRG-related lncRNA; lncRNA long non-coding RNA, NRG necroptosisrelated gene.Lines between two nodes indicate significant correlations between lncRNA and NRG.

Figure 3 :
Figure 3: Prognostic analysis of risk model in different sample sets.(a-c) Kaplan-Meier survival curves of overall survival of breast cancer patients between high-and low-risk group in training, testing, and entire samples sets, respectively; (d-e) Forest plots visualizing independent clinicopathological variables of overall survival in the entire sample set.

Figure 4 :
Figure 4: ROC curves of the risk model in different comparisons.(a) 1-, 3-, 5-, 10-year ROC curves of the model in complete sample set; (b) 10-year ROC curves for different clinicopathological variables; AUC area under curve, ROC receiver operating characteristic.

Figure 6 :
Figure 6: GSEA analysis of the top five KEGG biological pathways in high-and low-risk groups, respectively.GSEA Gene Set Enrichment Analysis, KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 7 :
Figure 7: ssGSEA analysis of immune landscapes in high-and low-risk group in the entire sample set.(a) Enrichment level of major immune cells in samples from high-and low-risk groups; (b) Activity level of typical immune functions or pathways between high-and low-risk groups; aDC activated dendritic cell, APC antigen presenting cell, CCR chemokine receptor, DC dendritic cell, HLA human lymphocyte antigen, iDC immature dendritic cell, IFN interferon, MHC major histocompatibility complex, NK cell natural killer cell, pDC plasmacytoid dendritic cell, ssGSEA single sample Gene Set Enrichment Analysis, Tfh follicular helper T cell, Th helper T cell, TIL tumor-infiltrating lymphocyte, Treg regulatory T cell; ***, p<0.001; **, p<0.01; *, p<0.05.

Figure 8 :
Figure 8: Correlation of major innate and adaptive immune cell infiltration with different risk scores in the entire sample set.(a) Bubble plot demonstrating the correlation between immune cell infiltration with risk score via different various algorithms; (b-f) correlation of infiltration scores of CD4 + T cells (b), CD8 + T cells (c), B cells (d), myeloid Dendritic cells (e), and neutrophils (f) with risk scores.Furthermore, we calculated immune scores and stromal scores for each BC patients and the boxplots show that patients with higher risk had lower immune/stromal/ESTIMATE scores than patients with lower risk (Figure9A-9C).

Figure 9 :
Figure 9: Comparison of tumor microenvironment of samples in entire sample sets between groups with different risk levels.(a) Comparison of immune cell infiltration between high-and low-risk group; (b) Comparison of stromal cell infiltration between high-and low-risk group; (c) Comparison of ESTIMATE score between high-and low-risk group.ESTIMATE score is the sum of the immune and stromal score.