Tumor Microenvironment-Mediated Immune Profiles Characterized by Distinct Survival Outcome and Immunotherapeutic Efficacy in Breast Cancer

Background: Numerous reports have highlighted that the tumor microenvironment (TME) is closely linked to survival outcome and therapeutic efficacy. However, a comprehensive investigation of the TME feature in breast cancer (BC) has not been performed.

Methods: Here, we performed consensus clustering analysis based on TME cell expression profiles to construct TME pattern clusters and TME-related gene signature in BC. GSVA combined with CIBERSORT and ssGSEA algorithms were applied to evaluate the differences in biological pathway and immune cell infiltration level, respectively. The PCA method was employed to construct TME-score to quantify the TME-mediated pattern level in individual BC patients.

Results: We determined two distinct TME gene clusters among 3,738 BC samples, which exhibited distinct survival outcome and enriched biological processes. The TME features demonstrated that these two clusters corresponded to the established immune profiles: hot and cold tumor phenotypes, respectively. Based on TME-related signature genes, we constructed the TME-score and stratified BC patients into low and high TME-score groups. Patients with high TME-score exhibited favorable outcome and increased infiltration of immune cells. Further investigation revealed that high TME-score was also related with high expression of immunosuppressive molecules, decreased tumor mutation burden (TMB), and high rate of mutation in significantly mutated genes (SMGs) (e.g., PIK3CA and CDH1).

Conclusion: Assessing the TME-mediated pattern level of individual BC patients will assist us in better understanding the responses of BC patients to immunotherapies and directing more effective immunotherapeutic approaches.

Background

Breast cancer (BC) remains the most frequently diagnosed malignancy and the main cause of cancer-related mortality for global females, and its incidence and fatality rates are 24.2% and 15.0%, respectively (Siegel et al., 2018; Kalimutho et al., 2019). Numerous conventional therapies, including surgery, radio(chemo)therapy, endocrine therapy and targeted therapy, have been applied in clinical practice, resulting in the significant decline of BC death rate (Wang et al., 2019). Recently, immunotherapy based on immune checkpoint inhibitors (ICIs) has been identified as a promising therapeutic regimen for BC patients and could serve as a complement to traditional schemes for treatment of BC (Emens, 2018; Basu et al., 2019). However, a significant portion of BC patients fail to benefit from checkpoint-blocking antibodies and exhibit an increased risk of immune-related adverse events (Wang et al., 2017). Thus, identification of an accurate signature to predict therapeutic responsiveness of immunotherapies is required.

As a complex and continuously evolving entity, the tumor microenvironment (TME) is composed of various cell populations (immune cells, fibroblasts, endothelial cells, etc.) and extracellular constituents (extracellular matrix, growth factors, cytokines, etc.), which surround malignant cells and are supported by the vascular network (Hui and Chen, 2015; Wu and Dai, 2017; Hinshaw and Shevde, 2019). A series of studies have demonstrated that TME not only has profound effects during cancer growth and metastatic progression but also plays a pivotal role in predicting therapeutic efficacy (Deepak et al., 2020). As for BC, tumor-infiltrating lymphocyte, a TME component, has been identified as a biomarker to reflect patients’ clinical outcome and predict the potential benefits from treatment (Stanton and Disis, 2016; Denkert et al., 2018). There have been plenty of studies on evaluation of the clinical significance of TME components by computational methods, such as single-sample gene set enrichment analysis (ssGSEA), which could assess the relative abundance of each cell infiltration in the TME (Ye et al., 2019). However, few studies attempted to discover the role of comprehensive TME feature in prognosis and therapeutic responsiveness in BC.

Here, we incorporated genomic and transcriptomic data of 3,738 BC samples from TCGA, GSE, and METABRIC datasets to comprehensively evaluate the relationship between TME cell infiltration characteristics and survival outcome and therapeutic response. Two TME gene clusters were identified using consensus clustering analysis, which were characterized by cold and hot tumor phenotypes, respectively. Besides, to quantify the pattern level mediated by TME in individual BC patients, we constructed the TME-score using principal component analysis (PCA), which could predict immunotherapies based on ICIs and chemotherapeutics, suggesting that TME performed an essential role in guiding therapies for BC.

Materials and Methods

Acquisition and Preprocessing of BC Datasets

Gene-expression profiles and clinical information of BC samples were obtained from publicly available databases: TCGA cBioPortal, and GEO. A total of 3,738 BC samples were incorporated into this research, including those from TCGA-BRCA (N = 1,091), METABRIC (N = 1904), GSE20685 (N = 327) (Kao et al., 2011), GSE20711 (N = 88) (Dedeurwaerder et al., 2011), GSE42568 (N = 104) (Clarke et al., 2013), GSE58812 (N = 107) (Jézéquel et al., 2015), and GSE88770 (Metzger-Filho et al., 2013) (N = 117) (Table 1). As for TCGA dataset, we downloaded FPKM-normalized values and then transformed them into TPM format, which were more similar to microarray data and more comparable between samples (Wagner et al., 2012). The Homo_sapiens.GRCh38.104.chr.gtf from the ENSEMBLE website was used as an annotation file to map ensemble ID to gene symbol. The scale method provided by R “limma” package was applied to normalize gene expression. We discarded genes with low abundance, whose expression value of 0 accounts for >80% of total samples and calculated the average value for duplicate genes. The selection criteria for BC datasets in the GEO database were as follows: 1) datasets with a sample size larger than 80; 2) datasets created based on the GPL570 platform; and 3) datasets with necessary clinical information, especially overall survival (OS) interval and status. The GPL570 annotation file was used to map the probes. We obtained the METABRIC dataset from cBioPortal and then imputed the missing data using R “impute” package. While merging the expression matrix of these 7 BC datasets into one meta-cohort, we applied the “ComBat” algorithm from R “sva” package to achieve the batch-effect removal (Leek et al., 2012). The genomic mutation data obtained from the UCSC Xena database was used for somatic mutation and copy number variation (CNV) analysis.

TABLE 1

TABLE 1. Baseline characteristics of breast cancer patients in TCGA, METABRIC and GSE datasets.

Evaluation of Cell Infiltration in the TME

Based on the expression profile containing 547 reference genes minimally representing each cell type, CIBERSORT could characterize and accurately calculate the proportion of distinct immune cell components from bulk tumor samples using a support vector regression and deconvolution algorithm (Wagner et al., 2012). To discriminate 22 human immune cell phenotypes sensitively and specifically, the LM22 gene signature was applied to quantify T cells, NK cells, B cells, macrophages, DCs, and myeloid subsets. Gene-expression profiles processed by standard annotation files were uploaded to the CIBERSORT website (http://cibersort.stanford.edu/), with the algorithm running based on the LM22 signature and 1,000 permutations. p-value <0.05 was considered as the significance criterion.

Consensus Clustering Analysis of TME Components

Based on the expression profiles of TME-infiltrating cells, we applied the unsupervised clustering method to construct distinct TME pattern clusters and classified patients for subsequent analysis. To determine the optimal number of clusters and guarantee their stability, the consensus clustering algorithm provided by the “ConsensusClusterPlus” package was performed repeatedly for 100 times (Monti et al., 2003). Consensus clustering, as a highly useful technique in cancer research, could detect unknown groups in a dataset based on intrinsic biological features and no external information. Besides, this method could provide quantitative and visual stability evidence derived from repeated subsampling and clustering.

Identification of Differentially Expressed Genes (DEGs) Between Distinct TME Clusters

Distinct TME pattern clusters have been determined by the above consensus clustering algorithm. Subsequently, to identify DEGs between distinct TME pattern clusters, we applied the empirical Bayesian approach provided by R “limma” package (Ritchie et al., 2015) to estimate gene-expression changes and screened out DEGs using the adjusted p-value <0.05 as the significance criterion.

Gene Set Variation Analysis GSVA and ssGSEA

R “GSVA” package (Hänzelmann et al., 2013) was utilized to perform GSVA enrichment analysis, which could explore the differences in biological processes between distinct pattern clusters mediated by DEGs. The gene set of c2.cp.kegg.v7.4.symbols.gmt obtained from the MSigDB database was utilized as the well-defined signature (Subramanian et al., 2005), and adjusted p-value <0.05 was considered as the significance criterion. To investigate the immune infiltration feature between distinct pattern clusters, ssGSEA was employed to quantify the infiltration levels of 23 different types of immune cell (Ren et al., 2021). Based on a Gaussian fitting model and multidimensional scaling, we estimated the bio-similarity of immune cells, calculated the enriched score of each immune cell, and uniformly distributed the normalized score from 0 to 1.

Construction of the TME-Score

To quantify the TME-mediated pattern level of individual BC patients, we developed a scoring scheme by the following procedures. Firstly, the prognostic analysis of the overlapping DEGs identified between distinct TME pattern clusters was performed by univariate Cox regression analysis and genes with prognostic impact were extracted. Then, feature selection of these genes with prognostic value was analyzed by recursive feature elimination with random forest and the 10-fold cross-validation method provided by the “caret” package. Finally, the expression profile of the determined genes was employed to perform PCA analysis and we extracted principal component 1 and 2 as the signature score. The advantage of this method is focusing on the score on the set with the largest block of well-correlated (or anticorrelated) genes in the set, while down-weighting contributions from genes that do not track with other set members. The TME-score was defined using a formula similar to previous studies (Zhang et al., 2020; Chong et al., 2021): ∑ (PC1i + PC2i), in which i represents the expression of DEGs with prognostic value.

Quantification of Predictors Related to Immune Response: Immunophenoscore (IPS), the Estimation of Stromal and Immune Cells in Malignant Tumors Using Expression Data (ESTIMATE) and ICIs

As a superior biomarker to predict the response of anti-PD-1 and CTLA-4 therapies, IPS could calculate the determinants of tumor immunogenicity and depict the cancer antigenomes and intra-tumoral immune profiles (Charoentong et al., 2017). This scoring scheme derived from a panel of immune-related genes, which belong to four classes: suppressor cells, effector cells, immunomodulators or checkpoints, and MHC-related molecules. By averaging the samplewise Z scores of the four classes within the respective category, the sum of the weighted averaged Z score was calculated as the IPS. The ESTIMATE algorithm could take advantage of the unique properties of the transcriptional landscape to obtain the tumor cellularity and tumor purity. By using the ESTIMATE algorithm (Yoshihara et al., 2013), we calculated the immune and stromal score to predict the level of infiltrating immune and stromal cells, and these form the basis to infer tumor purity. Tumor tissues with abundant immune cell infiltration represented a higher immune score and a lower level of tumor purity. Immunotherapy-based common immunosuppressive molecules PD-1, PD-L1, and CTLA-4 have achieved tremendous success in clinical practice; novel immune checkpoint proteins such as TIGIT and LAG3 were also strongly recommended to evaluate immune response.

Statistical Analysis

All data processing was generated in R-4.1.0. For quantitative data, we performed Student’s t-tests and the Wilcoxon rank-sum test to estimate the statistical significance for normally and non-normally distributed variables, respectively. For comparisons of more than two groups, we used one-way analysis and Kruskal–Wallis test as parametric and nonparametric methods, respectively. To analyze the relationship between tumor pattern clusters and prognosis, we applied the R “Survminer” package to perform Kaplan–Meier survival analysis and the Cox proportional hazards model. The surv-cutpoint function provided by R “survival” package was used to divide patients into high and low TME-score and tumor mutation burden (TMB) groups. All statistical p-values were two-sided, with a p-value <0.05 as statistical difference and adjusted p-value calculated using the Benjamini–Hochberg correction.

Results

Landscape of Tumor-Infiltrating Immune Cells in BC

Seven independent BC datasets (TCGA-BRCA, METABRIC, GSE20685, GSE20711, GSE42568, GSE58812, and GSE88770) with completed OS information were merged into one meta-cohort. CIBERSORT was applied using the expression profiles of the meta-cohort to obtain the proportion of 22 immune cells. We classified immune cells into four clusters, and immune cells in each cluster have similar functions. Immune cells in Cluster A recognize antigens and act as messengers between the innate and adaptive immune systems, including activated dendritic cells, resting mast cells, naïve B cells, CD8 T cells, and macrophage M1. Most immune cells in Cluster B have the function of attacking and killing exogenous antigens, such as activated NK cells, T cells regulatory, T cells follicular helper, plasma cells, neutrophils, and activated mast cells. Cluster C includes macrophage M2, which decreases inflammation and encourages tissue repair, and CD4 memory resting T cells. Non-activated macrophage M0 forms Cluster D alone. Univariate Cox regression and Kaplan–Meier analysis were conducted to explore the prognostic impact of immune cells in BC samples, and p-value <0.05 was regarded as statistical difference. Figure 1A comprehensively illustrated the landscape of the interactions of the immune cells, cell lineages, and their prognostic impact in BC patients. This network demonstrated that the crosstalk of the infiltration among TME cells significantly correlated with the survival outcome of BC patients and played an indispensable role in the construction of distinct TME pattern clusters. To stratify BC samples with qualitatively different TME-mediated pattern clusters, we performed consensus clustering analysis based on the TME cell expression profiles and determined two distinct clusters, which were termed as TME cluster A, including 1,682 BC samples, and TME cluster B, containing 545 tumor ones (Figure S1). Further analysis of survival outcome revealed that TME cluster A possessed more prominent advantage in OS and progression-free survival (PFS) (Figure 1B).

FIGURE 1

FIGURE 1. Tumor microenvironment-mediated pattern clusters in breast cancer. (A) The interaction between TME-infiltrating immune cells in BC. The lines connecting immune cells represented their interaction with each other. The size of each circle represented the prognostic effect of each immune cell and scaled by p-value. Protective factors for patients’ survival were indicated by a green dot in the circle center and risk factors indicated by the purple dot in the circle center. (B) Kaplan–Meier curves of OS and PFS for BC patients in the meta-cohort and TCGA-BRCA dataset with distinct TME clusters. The numbers of patients in TME cluster A and B were 1,682 and 545, respectively. (C,D) Unsupervised clustering of immune cells in the TCGA-BRCA and METABRIC datasets. Clinicopathological information, including age, gender, ER, PR, HER2, molecular subtype, stage, grade and TME cluster were used as patient annotations. Red represented the high expression of immune cells, and blue represented the low expression of immune cells. (E, F) The fraction of tumor-infiltrating lymphocytes in distinct TME clusters using the ssGSEA and CIBERSORT algorithms. Within each group, the scattered dots represented TME cell expression values. The thick line represented the median value. The bottom and top of the boxes were the 25th and 75th percentiles (interquartile range). The statistical difference of two TME clusters was compared through the Kruskal–Wallis H test. *p < 0.05; **p < 0.01; ***p < 0.001. (G) PCA of TME-infiltrating immune cells to distinguish TME cluster A from B.

To further investigate the differences in the biological and clinicopathological characteristics underlying TME cluster A and B, we focused on the TCGA-BRCA and METABRIC datasets, because they compromised the largest BC sample size and the most comprehensive clinical features. As shown in Figures 1C, D, TME cluster A with survival benefit was significantly enriched in macrophages M1 and T cells CD8/CD4 memory activated, while T cells regulatory, macrophages M0, and M2 were mainly concentrated in cluster B. Besides, CIBERSORT and ssGSEA were applied together to estimate the differences in TME profiles between distinct TME clusters, which was highly consistent with the results obtained in the TCGA-BRCA and METABRIC datasets (Figures 1E, F). In terms of clinical features, TME cluster A exhibited a higher stromal and immune score and relatively low tumor purity, which was established by the ESTIMATE algorithm. Based on the TME cell expression profiles, PCA could separate samples with distinct TME clusters into two opposite directions (Figure 1G).

Identification of TME-Related Signature Genes

Based on TME cell expression profiles, the consensus clustering algorithm could stratify BC patients into 2 distinct TME clusters. However, the genomic alteration and transcriptomic perturbation underlying these 2 clusters were not fully elucidated. Thus, the expression changes of 13,711 genes between the two TME clusters were investigated by the “limma” package and a total of 205 DEGs were determined and regarded as the crucial TME-related gene signature (Supplementary Table S1). Further investigation of the TME-related gene signature demonstrated that these genes were associated with or specific to immune cells, such as BTLA, CD3D, CD3E, CD3G, CD8A, CD8B, CTLA4, ICOS, KLRG1, and PDCD1 for T cells; CD19, CD5, CD79A and CD79B for B cells; and CD48 for dendritic cells. One speculation is that these DEGs were determined between distinct clusters, which were identified based on infiltrating immune cells. Thus, the TME-related signature genes were mainly expressed by various types of immune cell. Then, GO enrichment analysis of these genes demonstrated that biological processes associated with immune activation, such as regulation of lymphocyte proliferation, antigen processing and presentation, and T cell activation, were prominently overrepresented (Figure 2A). KEGG analysis also revealed that these DEGs played an indispensable role in modulating the TME landscape (Figure 2B, Supplementary Table S2). Protein–protein interaction (PPI) was constructed to investigate the interactions among TME-related signature genes, and the high confidence (0.700) was set as the minimum required interaction score. We identified 202 nodes and 587 edges in the PPI network, and 16 genes with degree of interaction more than 40 were recognized as hub genes (Figures 2C, D).

FIGURE 2

FIGURE 2. Functional enrichment analysis of differentially expressed genes between distinct TME clusters. (A, B) GO and KEGG analysis of DEGs. (C) PPI of DEGs. The high confidence (0.700) was considered as the minimum required interaction score. (D) Identification of hub genes. Genes with degree of interaction more than 40 were recognized as hub genes. The X-axis label means the degree of interaction of each gene.

Landscape of Genetic Alterations of TME-Related Gene Signature in BC

Figure 3A first depicted the prevalence of somatic mutation of 205 TME-related signature genes among 986 BC samples with available variant classification and variant type information, out of which 320 (32.45%) samples experienced genetic alterations, mainly including splice site mutation, nonsense mutation, and missense mutation. Considering that COL14A1 exhibited the highest mutation frequency, we investigated the difference in the expression of TME-related signature genes between COL14A1-wild and mutant type and a total of 10 genes (ALDH1A1, CH25H, COL11A1, GJB2, HLA-DPB1, IL18R1, MMP13, P2RY10, PTGER4, and UNC5B) were found differentially expressed (Figure S2). Further exploration of 205 DEGs revealed a prevalent CNV alteration. ABCA8, CD79B, TBC1D10C, ENPP2, COL14A1, and CD7 had prevalent CNV amplification, while CD3E, IL10RA, CD3D, CD3G, ZNF683, and CD52 showed widespread CNV deletion (Figure 3B). Interestingly, we discovered that TME-related signature genes exclusively expressed in T cells, such as CD3D, CD3E, and CD3G, exhibited deleted CNV. We speculated that this phenomenon may be attributed to the lower CD3⁺ T cell infiltration. Because distinct TME clusters were determined based on the characteristics of infiltrating immune cells. Besides PCA and t-SNE based on the crucial signature gene expression profiles were performed and BC samples and normal ones could be completely distinguished (Figure 3C). To explore whether the above genetic variation affected transcriptomic expression, we calculated the differences in gene expression value between BC and normal samples and observed that HLA-DRA, HLA-DPB1, and SPP1 were significantly increased in BC tissues, while, FABP4, PLIN1, and APOD were markedly downregulated in tumor tissues (Figure 3D, Supplementary Table S3). As expected, we discovered that the CNV alteration could be the major contributor to perturbation on the expression of the TME-related gene signature. Compared with normal samples, genes with CNV amplification exhibited a remarkedly higher expression in BC samples, such as CD7, IDO1, and LYZ, while other genes with deleted CNV, such as INPP5D, were significantly decreased in tumor tissues (Figures 3B,D), indicating that there was high heterogeneity between the genomic and transcriptomic landscape of TME-related signature genes and this expression imbalance played an essential role in BC tumorigenesis and progression.

FIGURE 3

FIGURE 3. The landscape of genetic alteration of differentially expressed genes in breast cancer. (A) 320 of the 986 BC samples experienced genetic alterations, mostly including missense mutation, nonsense mutation, and splice site mutation. The number on the right indicated the mutation frequency of each DEG. Each column represented the individual BC sample. (B) The CNV mutation frequency of DEGs was prevalent. The column represented the alteration frequency. The amplification frequency, pink dot; the deletion frequency, green dot. Because of the size of the gene list (205 genes), we presented the CNV results in 3 plots: one for CNV deletion and the other two for amplified CNV. (C) PCA and t-SNE analysis based on the expression of DEGs to distinguish BC samples (red dots) from normal ones (green dots). (D) The differences in mRNA expression level of DEGs between tumor and normal samples was compared by R “limma” package. The asterisks represented the statistical p-value (*p < 0.05; **p < 0.01; ***p < 0.001). For graphical reasons, we presented the gene lists (upregulated and downregulated genes in BC samples) in 2 plots and ordered them based on the log FC.

Construction of Pattern Clusters Mediated by the DEGs

Based on the above hypotheses, we applied consensus clustering algorithm based on the TME-related gene signature to classify BC patients and obtained 2 distinct subtypes, which were defined as TME gene cluster I and II (Figures 4A, B; Supplementary Figure S3). PCA based on the signature gene expression profiles could satisfactorily separate BC samples into opposite directions (Figure 4C). The predictive ability of these 2 gene clusters in survival outcome revealed that gene cluster II exhibited a significant prognostic benefit (Figure 4D). To further explore the biological behaviors between distinct gene clusters, we conducted GSVA and found that cluster II presented enrichment pathways associated with immune activation, including allograft rejection, natural killer cell-mediated cytotoxicity, T/B cell receptor signaling pathway, Toll/NOD-like receptor signaling pathway, antigen processing and presentation, chemokine signaling pathway, and cytokine–cytokine receptor interaction (Figure 5A). Besides, we compared transcripts associated with immune activation between distinct gene clusters and observed a significant upregulation in cluster II (Figure 5B). Further analysis of the TME feature via CIBERSORT revealed that compared with gene cluster I, cluster II had elevated infiltration of antitumor immune cells, such as NK cells activated, T cells CD4 memory activated/CD8, and macrophages M1, and a decreased proportion of immunosuppressive or tumor-promoting immune cells, such as macrophages M2 and Treg cells (Figure 5C). Furthermore, we utilized ssGSEA to characterize the immune infiltration profile and observed the similar TME landscape between these 2 gene clusters (Figure 5D). We also applied the ESTIMATE algorithm to evaluate the overall infiltrating level of immune cells and stromal components between these 2 clusters. As expected, we discovered that gene cluster II exhibited an elevated immune and stromal score, indicating that cluster I had a higher tumor purity (Figure 5C). Moreover, some transcripts of immune checkpoint proteins were compared between distinct gene clusters and we observed that PD-L1, PD-1, CTLA4, TIGIT, and LAG-3 were highly expressed in gene cluster II (Figure 5E). Based on these findings, we determined that these 2 TME gene clusters were characterized by distinct immune features. Expectedly, cluster I was defined as cold tumor, characterized by the lack of tumor T cell infiltration and decreased expression of immunosuppressive molecules, and cluster II was recognized as hot tumor, characterized by abundant immune cells and elevated PD-L1 expression.

FIGURE 4

FIGURE 4. Construction of TME-related gene signature. (A, B) Unsupervised clustering of TME-related signature genes in the TCGA and METABRIC datasets. Clinicopathological information including age, gender, ER, PR, HER2, stage, grade, and molecular subtype is shown in annotations above. Red represented the high expression, while blue represented the low expression. (C) PCA based on the TME-related gene signature could satisfactorily distinguish between TME gene cluster I and II. (D) The Kaplan–Meier curves of OS and PFS for BC samples in the meta-cohort and TCGA-BRCA dataset, respectively.

FIGURE 5

FIGURE 5. Biological pathways and tumor microenvironment characteristics in distinct gene clusters. (A) Heatmap showed the GSVA score of representative hallmark pathways curated from MSigDB between distinct TME gene clusters. (B) The differences in mRNA expression of transcripts associated with immune activation between distinct gene clusters. (C, D) The fraction of tumor-infiltrating immune cells between distinct gene clusters using the CIBERSORT and ssGSEA algorithms. Within each group, the scattered dots represented TME cell expression values. The thick line represented the median value. The bottom and top of the boxes were the 25th and 75th percentiles (interquartile range), respectively. The statistical difference of the three gene clusters was compared through the Kruskal–Wallis H test. *p < 0.05; **p < 0.01; ***p < 0.001. (E) Comparison of the expression level of immunosuppressive molecules between distinct TME gene clusters.

Development of the TME-Score and Evaluation of its Clinical Significance

The above analyses have demonstrated that TME immune profiles were tightly connected with survival outcome and immune checkpoint proteins in BC. However, these findings were only based on the patient population and could not accurately predict the pattern level mediated by TME in individual tumors. Consequently, we constructed a scoring scheme based on identified TME-related signature genes and defined it as TME-score to quantify the TME-mediated pattern level in individual BC patients. Considering that the quantification of TME-score was complex, we used the alluvial diagram to illustrate the workflow of the TME-score construction (Figure 6A). The Kruskal–Wallis test demonstrated that TME cluster A and gene cluster II with survival advantage exhibited a higher TME-score (Figure 6B). Spearman analysis was performed to examine the relationship between the TME-score and immune landscape. The correlation matrix revealed that the TME-score positively correlated with tumor-infiltrating lymphocytes, including activated CD4+/CD8+ T cells, activated B cells, nature killer T cells, and macrophages M1, and negatively correlated with T cells regulatory, tumor-associated neutrophils, and macrophages M2, demonstrating the crosstalk between TME-score and tumor-infiltrating immune cells (Figure 6C). Moreover, GSEA demonstrated that biological processes associated with immune activation were significantly enriched in the high TME-score group (Figure 6D, Supplementary Table S4). The predictive ability of the TME-score for BC prognosis was estimated by stratifying patients into high- and low-score groups, according to the median value of −3.75 obtained from the METABRIC dataset (discovery cohort). As we expected, patients in the high-score group possessed a prominent survival benefit. Independent validations were performed using the TCGA-BRCA and GEO datasets (external validation cohort) to further assess the prognostic value of TME-score. The median value obtained from the discovery cohort was used to stratify BC patients in the validation cohort into high- and low-score groups. Similarly, patients with a high TME-score exhibited a promising survival outcome, compared with those with a low TME-score (Figure 6E).

FIGURE 6

FIGURE 6. Construction of the TME-score and exploration of its clinical significance. (A) Alluvial diagram of the TME cluster, TME gene cluster, and TME-score. (B) Relative distribution of TME-score between distinct TME clusters and TME gene clusters. (C) Correlations between TME-score and tumor-infiltrating immune cells using Spearman analysis. The asterisks represented the statistical p-value (*p < 0.05). (D) Representative results of KEGG pathways between high and low TME-score groups via GSEA. (E) Kaplan–Meier curves for patients with high and low TME-score groups in the meta-cohort and METABRIC, TCGA-BRCA, GSE20685, GSE58812, GSE42568, GSE88770, and GSE20711 datasets. The median TME-score obtained from the METABRIC dataset (discovery cohort) was utilized to separate BC samples into high- and low-score groups. Independent validations were then performed using TCGA-BRCA and GSE datasets (external validation cohort).

Estimation of the Role of TME-Score in Tumor Somatic Mutation Signature, Immunotherapy, and Chemotherapeutic Efficacy

Accumulating evidence has indicated that cancer-specific antigens produced by somatic mutation could influence the responsiveness to immune-checkpoint blockades. Thus, the distribution pattern of TMB was investigated between high and low TME-score groups. The results demonstrated that TME-score was negatively correlated with TMB and the high-score group exhibited lower TMB (Figures 7A,B). BC patients were categorized into high- or low-TMB group based on the median cutoff value of 0.42. Survival analysis further demonstrated that the high-TMB group had a dismal survival outcome, and as expected, the combined high-TMB and low TME-score group exhibited the worst prognosis (Figures 7C, D). Significantly mutated genes (SMGs) analysis was also conducted between different TME-score groups, and the SMGs landscape revealed that CDH1 (16.27 vs. 4.7%) and PIK3CA (37.29 vs. 26.37%) exhibited higher somatic mutation rates in the high-score group, whereas GATA3 (7.63 vs. 14.1%) had a higher somatic mutation rate in the low-score group (Figures 7E, F; Table 2).

FIGURE 7

FIGURE 7. Characteristics of the TME-score in tumor somatic mutation. (A) Correlation analysis between TME-score and TMB. (B) Relative distribution of TMB in high versus low TME-score groups. (C) Kaplan–Meier curve for high and low TMB patient groups (p = 0.001). (D) Kaplan–Meier curve for subgroup patients stratified by both TME-score and TMB (p < 0.001). (E, F) Mutational landscape of SMGs in the TCGA-BRCA dataset stratified by high (left panel) versus low TME-score (right panel) groups. Individual patients were represented in each column. The upper bar plot showed TMB, and the right bar plot showed the mutation frequency of each gene in separate TME-score groups.

TABLE 2

TABLE 2. Significantly mutated genes between high and low TME score group.

Immunotherapy based on immunosuppressive molecules has achieved a major breakthrough in antitumor response. Many well-known predictors, especially PD-L1, were extensively used to assess immune response. We discovered that the expression levels of immunosuppressive molecules, including PD-L1, CTLA4, PD-1, TIGIT, and LAG3, were pronouncedly elevated in the high TME-score group (Figure 8A), indirectly demonstrating the essential role of the TME-score in mediating immune response. In addition, we also observed that patients in the high TME-score group exhibited significant therapeutic benefits from ICIs treatment represented by CTLA4-/PD-1-, CTLA4-/PD-1+, CTLA4+/PD-1-, and CTLA4+/PD-1+, indirectly suggesting that the TME-score played an essential role in predicting response to immunotherapies (Figure 8B). In addition to immunotherapies based on ICIs, we attempted to investigate the connection between the TME-score and chemotherapeutic drugs commonly used for BC in clinical practice. We discovered that with the exception of docetaxel and lapatinib, the high TME-score group exhibited decreased IC50 of cisplatin, doxorubicin, gemcitabine, methotrexate, paclitaxel, roscovitine, vinblastine, and vinorelbine, suggesting that the high TME-score group may enjoy therapeutic benefits from these chemotherapeutic drugs (Figure 8C).

FIGURE 8

FIGURE 8. Correlation between the TME-score and immunotherapeutic benefits and chemotherapeutic efficacy. (A) The relative distribution of immunosuppressive molecules was compared between TME-score high versus low groups in the meta-cohort. (B) Relative distribution of immunotherapeutic efficacy in high TME-score versus low TME-score groups. (C) Relative distribution of IC50 for chemotherapeutic drugs commonly used in clinical practice between TME-score high versus low groups.

Discussion

Cancer immunotherapies targeting PD-1, PD-L1, and CTLA4 have achieved durable and robust responses for BC patients in clinical practice (Keenan and Tolaney, 2020; Emens, 2021). However, approximately one-third of patients fail to respond to these therapeutic agents, and numerous studies have demonstrated that microsatellite instability, PD-L1/PD-1 expression, and TMB are not efficient predictors of the efficacy of immunotherapies (Zeng et al., 2019; Tang et al., 2020). To maximize the therapeutic benefit, identification of novel predictors for the efficacy of checkpoint-blocking antibodies is essential. Mounting evidence has implied that TME not only takes on an indispensable role in survival outcome but also has profound effects on predicting immunotherapeutic responsiveness (Basu et al., 2019; Emens, 2021). In the present research, we applied several computational methodologies to quantify the infiltration level of immune cells and comprehensively explored the connections between TME-infiltrating immune cells and clinicopathological characteristics, survival outcome, and immunotherapeutic efficacy in BC.

Based on TME cell expression profiles, we employed consensus clustering analysis to determine two TME clusters, which exhibited distinct survival outcome. We speculated that this survival discrepancy may derive from the distinct immune class: TME cluster A was characterized by abundant antitumor T cells CD8 and macrophages M1, corresponding to the immune-activated subtype, while TME cluster B was characterized by immunosuppressive T cells regulatory and higher tumor purity corresponding to the immune-suppressed subtype.

Moreover, DEGs identified between immune-activated and immune-suppressed subtypes were prominently elevated in biological pathways associated with immune activation, demonstrating that these DEGs were regarded as TME-related gene signature. Based on the TME-related signature genes, two transcriptomic clusters were constructed and characterized by different survival outcome, which was in great accord with the results of TME clusters. One possible speculation for this survival difference is that TME gene clusters were characterized by distinct antitumor immunity. TME gene cluster II was characterized by immune activation and abundant antitumor immune cell infiltration, corresponding to hot tumor phenotype, while gene cluster I was characterized by immune suppression and immunosuppressive immune cell infiltration, corresponding to cold tumor phenotype (Duan et al., 2020). Several studies have demonstrated that TME components were tightly connected with aggressive properties of tumor and the likelihood of immune response (Galon and Bruni, 2019; Zeng et al., 2019). Evaluation of the density of CD3 and CD8 lymphocyte populations at the tumor center and margin was found remarkedly correlated with patients’ survival status and immunotherapeutic efficacy (Pagès et al., 2018). We also confirmed that TME gene cluster II was prominently related with antitumor lymphocyte infiltration and a high expression level of immunosuppressive molecules, indicating its potential value in predicting immunotherapeutic benefits.

Furthermore, a scoring scheme named TME-score was established to quantify the TME-mediated pattern level in individual BC patients and direct therapeutic interventions more precisely. As a result, gene cluster II characterized by a hot tumor phenotype showed a high TME-score, while gene cluster I characterized by a cold tumor phenotype exhibited a low TME-score. In addition, we observed that TME-score could serve as a prognostic biomarker in BC. Further analysis demonstrated that TME-score was markedly linked to immunosuppressive molecules and immunotherapy, implying that the TME-score could influence immunotherapeutic efficacy. Based on these findings, we believed that TME-score could be used in clinical practice to identify immune profiles and direct therapeutic strategies.

Assessment of mutated genes capable of driving tumors is one milestone toward cancer detection and therapeutic approach selection. Here, we observed that compared with the low TME-score group, PIK3CA and CDH1 exhibited elevated mutation rates in the high TME-score group, while GATA3 showed augmented mutation rates in the low TME-score group. Recent studies revealed that CDH1 and PIK3CA mutation in genetically modified mice could result in an immune-related subtype for invasive lobular carcinoma of the breast, which was characterized by enhanced immune infiltration and a strong signature for Treg cells and immunosuppressive molecule-based immune checkpoint activation (An et al., 2018). GATA3 is an essential regulator of immune cell function, and its mutation could result in the dysfunction of normal T cells (Usary et al., 2004). These TME-score mediated driver gene mutations remarkedly correlated with immune activity, highlighting the complicated connection between TME and tumor immunogenomic features.

Although, we obtained seven retrospective BC datasets (including a total of 3,738 samples) and performed a comprehensive analysis toward the correlation between TME and prognosis and immunotherapeutic efficacy in BC. Prospective datasets of BC samples were required to verify our results. Besides, it is appropriate to systematically assess the infiltrating immune cells in the tumor core and invasive margin, considering that distinct tumor regions are essential.

Conclusion

In conclusion, we comprehensively analyzed TME-mediated patterns of 3,738 BC patients based on tumor-infiltrating immune cells and systematically linked this pattern with prognosis and immunotherapeutic efficacy. This integrated analysis demonstrated that assessing the TME-mediated pattern of individual tumor will assist us in better understanding TME characteristics and directing more effective immunotherapeutic approaches.

Data Availability Statement

The datasets analyzed during the current study are available in the TCGA (https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/geo/) and cBioPortal (http://www.cbioportal.org/) databases.

Author Contributions

LX performed the statistical procedure and drafted the article. WL designed this study, and YH helped edit the manuscript. The final version of the manuscript was approved by all authors.

Funding

This work was supported by the National Scientific Foundation of China (No. 81870554) and Foundation from Renji Hospital, School of Medicine, Shanghai Jiaotong University (2019NYLYCP0102).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.840348/full#supplementary-material

Abbreviations

BC, breast cancer; ICIs, immune checkpoint inhibitors; TME, tumor microenvironment; ssGSEA, single-sample gene set enrichment analysis; PCA, principal component analysis; OS, overall survival; CNV, copy number variation; DEGs, differentially expressed genes; GSVA, gene set variation analysis; IPS, immunophenoscore; ESTIMATE, estimation of stromal and immune cells in malignant tumors using expression data; TMB, tumor mutation burden; PFS, progression-free survival; PPI, protein–protein interaction; SMGs, significantly mutated genes.

References

An, Y., Adams, J. R., Hollern, D. P., Zhao, A., Chang, S. G., Gams, M. S., et al. (2018). Cdh1 and Pik3ca Mutations Cooperate to Induce Immune-Related Invasive Lobular Carcinoma of the Breast. Cel Rep. 25 (3), 702–714. doi:10.1016/j.celrep.2018.09.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Basu, A., Ramamoorthi, G., Jia, Y., Faughn, J., Wiener, D., Awshah, S., et al. (2019). Immunotherapy in Breast Cancer: Current Status and Future Directions. Adv. Cancer Res. 143, 295–349. doi:10.1016/bs.acr.2019.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chong, W., Shang, L., Liu, J., Fang, Z., Du, F., Wu, H., et al. (2021). m6A Regulator-Based Methylation Modification Patterns Characterized by Distinct Tumor Microenvironment Immune Profiles in colon cancerA Regulator-Based Methylation Modification Patterns Characterized by Distinct Tumor Microenvironment Immune Profiles in colon Cancer [J]. Theranostics 11 (5), 2201–2217. doi:10.7150/thno.52717

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, C., Madden, S. F., Doolan, P., Aherne, S. T., Joyce, H., O’Driscoll, L., et al. (2013). Correlating Transcriptional Networks to Breast Cancer Survival: a Large-Scale Coexpression Analysis. Carcinogenesis 34 (10), 2300–2308. doi:10.1093/carcin/bgt208

PubMed Abstract | CrossRef Full Text | Google Scholar

Dedeurwaerder, S., Desmedt, C., Calonne, E., Singhal, S. K., Haibe-Kains, B., Defrance, M., et al. (2011). DNA Methylation Profiling Reveals a Predominant Immune Component in Breast Cancers. EMBO Mol. Med. 3 (12), 726–741. doi:10.1002/emmm.201100801

PubMed Abstract | CrossRef Full Text | Google Scholar

Deepak, K. G. K., Vempati, R., Nagaraju, G. P., Dasari, V. R., S, N., Rao, D. N., et al. (2020). Tumor Microenvironment: Challenges and Opportunities in Targeting Metastasis of Triple Negative Breast Cancer. Pharmacol. Res. 153, 104683. doi:10.1016/j.phrs.2020.104683

PubMed Abstract | CrossRef Full Text | Google Scholar

Denkert, C., Von Minckwitz, G., Darb-Esfahani, S., Lederer, B., Heppner, B. I., Weber, K. E., et al. (2018). Tumour-infiltrating Lymphocytes and Prognosis in Different Subtypes of Breast Cancer: a Pooled Analysis of 3771 Patients Treated with Neoadjuvant Therapy. Lancet Oncol. 19 (1), 40–50. doi:10.1016/s1470-2045(17)30904-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, Q., Zhang, H., Zheng, J., and Zhang, L. (2020). Turning Cold into Hot: Firing up the Tumor Microenvironment. Trends Cancer 6 (7), 605–618. doi:10.1016/j.trecan.2020.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Emens, L. A. (2018). Breast Cancer Immunotherapy: Facts and Hopes. Clin. Cancer Res. 24 (3), 511–520. doi:10.1158/1078-0432.ccr-16-3001

PubMed Abstract | CrossRef Full Text | Google Scholar

Emens, L. A. (2021). Immunotherapy in Triple-Negative Breast Cancer. Cancer J. 27 (1), 59–66. doi:10.1097/ppo.0000000000000497

PubMed Abstract | CrossRef Full Text | Google Scholar

Galon, J., and Bruni, D. (2019). Approaches to Treat Immune Hot, Altered and Cold Tumours with Combination Immunotherapies. Nat. Rev. Drug Discov. 18 (3), 197–218. doi:10.1038/s41573-018-0007-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data [J]. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinshaw, D. C., and Shevde, L. A. (2019). The Tumor Microenvironment Innately Modulates Cancer Progression. Cancer Res. 79 (18), 4557–4566. doi:10.1158/0008-5472.can-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

Hui, L., and Chen, Y. (2015). Tumor Microenvironment: Sanctuary of the Devil. Cancer Lett. 368 (1), 7–13. doi:10.1016/j.canlet.2015.07.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Jézéquel, P., Loussouarn, D., GUéRIN-Charbonnel, C., Campion, L., Vanier, A., Gouraud, W., et al. (2015). Gene-expression Molecular Subtyping of Triple-Negative Breast Cancer Tumours: Importance of Immune Response [J]. Breast Cancer Res. 17, 43. doi:10.1186/s13058-015-0550-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalimutho, M., Nones, K., Srihari, S., Duijf, P. H. G., Waddell, N., and Khanna, K. K. (2019). Patterns of Genomic Instability in Breast Cancer. Trends Pharmacol. Sci. 40 (3), 198–211. doi:10.1016/j.tips.2019.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kao, K.-J., Chang, K.-M., Hsu, H.-C., and Huang, A. T. (2011). Correlation of Microarray-Based Breast Cancer Molecular Subtypes and Clinical Outcomes: Implications for Treatment Optimization. BMC Cancer 11, 143. doi:10.1186/1471-2407-11-143

PubMed Abstract | CrossRef Full Text | Google Scholar

Keenan, T. E., and Tolaney, S. M. (2020). Role of Immunotherapy in Triple-Negative Breast Cancer. J. Natl. Compr. Canc Netw. 18 (4), 479–489. doi:10.6004/jnccn.2020.7554

PubMed Abstract | CrossRef Full Text | Google Scholar

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 28 (6), 882–883. doi:10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzger-Filho, O., Michiels, S., Bertucci, F., Catteau, A., Salgado, R., Galant, C., et al. (2013). Genomic Grade Adds Prognostic Value in Invasive Lobular Carcinoma. Ann. Oncol. 24 (2), 377–384. doi:10.1093/annonc/mds280

PubMed Abstract | CrossRef Full Text | Google Scholar

Monti, S., Tamayo, P., Mesirov, J., and Golub, T. (2003). Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data [J]. Machine Learn. 52 (1), 91–118. doi:10.1023/a:1023949509487

CrossRef Full Text | Google Scholar

Pagès, F., Mlecnik, B., Marliot, F., Bindea, G., Ou, F. S., Bifulco, C., et al. (2018). International Validation of the Consensus Immunoscore for the Classification of colon Cancer: a Prognostic and Accuracy Study. Lancet 391 (10135), 2128–2139. doi:10.1016/S0140-6736(18)30789-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, E.-h., Deng, Y.-j., Yuan, W.-h., Wu, Z.-l., Zhang, G.-z., and Xie, Q.-q. (2021). An Immune-Related Gene Signature for Determining Ewing Sarcoma Prognosis Based on Machine Learning. J. Cancer Res. Clin. Oncol. 147 (1), 153–165. doi:10.1007/s00432-020-03396-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2018). Cancer Statistics, 2018. CA: A Cancer J. Clinicians 68 (1), 7–30. doi:10.3322/caac.21442

CrossRef Full Text | Google Scholar

Stanton, S. E., and Disis, M. L. (2016). Clinical Significance of Tumor-Infiltrating Lymphocytes in Breast Cancer. J. Immunotherapy Cancer 4, 59. doi:10.1186/s40425-016-0165-6

CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, R., Xu, J., Zhang, B., Liu, J., Liang, C., Hua, J., et al. (2020). Ferroptosis, Necroptosis, and Pyroptosis in Anticancer Immunity. J. Hematol. Oncol. 13 (1), 110. doi:10.1186/s13045-020-00946-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Usary, J., Llaca, V., Karaca, G., Presswala, S., Karaca, M., He, X., et al. (2004). Mutation of GATA3 in Human Breast Tumors. Oncogene 23 (46), 7669–7678. doi:10.1038/sj.onc.1207966

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, G. P., Kin, K., and Lynch, V. J. (2012). Measurement of mRNA Abundance Using RNA-Seq Data: RPKM Measure Is Inconsistent Among Samples. Theor. Biosci. 131 (4), 281–285. doi:10.1007/s12064-012-0162-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P.-F., Chen, Y., Song, S.-Y., Wang, T.-J., Ji, W.-J., Li, S.-W., et al. (2017). Immune-Related Adverse Events Associated with Anti-PD-1/pd-L1 Treatment for Malignancies: A Meta-Analysis. Front. Pharmacol. 8, 730. doi:10.3389/fphar.2017.00730

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Mccullough, L. E., White, A. J., Bradshaw, P. T., Xu, X., Cho, Y. H., et al. (2019). Prediagnosis Aspirin Use, DNA Methylation, and Mortality after Breast Cancer: A Population-Based Study. Cancer 125 (21), 3836–3844. doi:10.1002/cncr.32364

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., and Dai, Y. (2017). Tumor Microenvironment and Therapeutic Response. Cancer Lett. 387, 61–68. doi:10.1016/j.canlet.2016.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, L., Zhang, T., Kang, Z., Guo, G., Sun, Y., Lin, K., et al. (2019). Tumor-Infiltrating Immune Cells Act as a Marker for Prognosis in Colorectal Cancer. Front. Immunol. 10, 2368. doi:10.3389/fimmu.2019.02368

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Li, M., Zhou, R., Zhang, J., Sun, H., Shi, M., et al. (2019). Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol. Res. 7 (5), 737–750. doi:10.1158/2326-6066.cir-18-0436

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer [J]. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar