Integrated transcriptome and single-cell sequencing reveals M2-like tumor-associated macrophages-related prognostic signatures in breast cancer

doi:10.21203/rs.3.rs-5236065/v1

Download PDF

Research Article

Integrated transcriptome and single-cell sequencing reveals M2-like tumor-associated macrophages-related prognostic signatures in breast cancer

https://doi.org/10.21203/rs.3.rs-5236065/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: M2-like tumor-associated macrophages (M2-like TAMs) have crucial functions in the tumor microenvironment (TME) and cancer development. This study explored M2-like TAMs-related prognostic signatures in breast cancer (BRCA) by combining transcriptome and scRNA-seq.

Methods: All datasets (TCGA-BRCA, GSE20685, GSE176078) were downloaded from UCSC xena and GEO database. AUCell score of immune-related genes (IRGs) was calculated using R package. M2-like TAMs-related genes were screened by WGCNA. Signature genes predictive of BRCA prognosis were identified by univariate Cox and LASSOregression analyses. Then, RiskScore model was constructed and validated in external dataset. Nomogram was established by integrating RiskScore and clinical characteristics. Correlation analysis between RiskScore and TME or chemotherapeutic drugs was conducted.

Results: Macrophage exhibited the highest IRGs AUCell score in single-cell transcriptomic atlas of BRCA. Macrophage was split into C1 and C2 subtypes, and Macrophage C1 highly expressed the marker genes of M2 macrophages. Then, 903 M2-like TAMs-related genes were acquired. ARHGAP26, RILP, KLRB1, CSTA, KLHDC7B, PSMB8, KYNU, RNASE1, LONRF3, and TRPM2were screened as prognostic signatures to establish RiskScore model, exhibiting good robustness. Nomogram was constructed by combining stage, Age, and RiskScore, with good predictive performance. Besides, most immune cells showed a negative correlation with RiskScore. Eight drugs were notably correlated with RiskScore, and the high-risk group had higher half maximal inhibitory concentration (IC50) of Ribociclib_1632, conversely indicating that the BRCA patients with low RiskScore were more sensitive to Ribociclib_1632.

Conclusion: We identified 10 M2-like TAMs-related prognostic signatures, providing potential therapeutic targetsfor BRCA.

breast cancer

M2 macrophages

tumor microenvironment

immune-related genes

prognosis

therapeutic targets

Breast cancer (BRCA) is a prevalent malignant tumor in women (Wang et al., 2024a). Increasing incidence rate and mortality of BRCA in recent years has severely threatened female health (Cao et al., 2023). BRCA exhibits highly heterogeneous, with multiple genomes, gene expression profiles, as well as molecular markers (Zhou et al., 2024). BRCA is divided into Luminal A, Luminal B, HER2+, and triple-negative breast cancer (TNBC) by the expressions of estrogen-receptor (ER), progesterone-receptor (PR), and human epidermal growth factor receptor 2 (HER2), respectively (Hu et al., 2023). Different molecular subtypes of BRCA appear disparate clinical features and prognosis (Yang et al., 2023a). TNBC has high metastasis, strong invasion, together with poor prognosis, and lacks efficacious therapeutic targets (Fu et al., 2023; Liu et al., 2024b). In clinical practice, the common treatment methods for BRCA chiefly contain breast-conserving therapy, mastectomy, drug therapy, radiotherapy, chemotherapy, etc (Tan et al., 2023). Although great advancement has been made in early diagnosis and treatment, which extends the survival period of BRCA patients, 5-year overall survival (OS) of some patients remains very low due to recurrence, metastasis or drug resistance (Chen et al., 2024; Wang et al., 2024c). Hence, developing more accurate and reliable prognosis signatures and models that can help to grade the risk level and provide potential therapeutic targets for BRCA patients.

The TME of BRCA is a complicated micro-ecosystem, which has been proved to significantly impact the development and prognosis of BRCA by the interactions between tumor cells and stromal cells or immune cells of TME (Guo et al., 2023). TAMs are the immune cells with the highest abundance in breast TME that can appear different phenotypic polarization to perform specific functions (Li et al., 2021). Generally, M1 macrophages exert beneficial effects in TME, such as pro-inflammatory, cytotoxic, and anti-tumor (Nie et al., 2023). Noteworthily, M2 macrophages are deemed to be "tumor promoters", can stimulate angiogenesis by secreting VEGF, and supply nutrition for tumor cells, facilitating the progression of BRCA (Huang et al., 2022; Lu et al., 2023). Moreover, M2-like TAMs are also indispensable in promoting tumor metastasis and enhancing tumor drug resistance (Wang et al., 2024b). It has been proved that the high infiltration level of TAMs is connected with bad prognosis and low survival probability of many cancers (Kim et al., 2022; Wang et al., 2022a). Recent research of Cheng suggested that glycyrrhetinic acid could inhibit the M2-like polarization of macrophages by activating JNK1/2 signaling pathway, decrease the level of VEGF, and thus efficiently restrain the development and metastasis of BRCA (Cheng et al., 2023). Therefore, targeting M2-like TAMs may become a new direction for controlling tumor metastasis and overcoming drug resistance of BRCA. However, to our knowledge, dependable and specific gene signatures related to M2-like TAMs have not been found to predict the prognosis of BRCA patients.

The scRNA-seq approach can be applied to understand the molecular characteristics of various cell clusters in TME (Kodous et al., 2023; Huang et al., 2024). We downloaded the Cancer Genome Atlas (TCGA)-BRCA, GSE20685, GSE176078 datasets from the public database. Then, the single-cell transcriptomic atlas of BRCA and heterogeneity of Macrophage were analyzed in GSE176078 dataset. M2-like TAMs-related genes in TCGA-BRCA dataset were filtered by WGCNA. Subsequently, prognostic signatures were identified to construct RiskScore model for BRCA, the robustness of which was validated in GSE20685 dataset. Nomogram was established to investigate whether RiskScore was independent of other clinicopathological features. Furthermore, the correlation between RiskScore and TME, as well as chemotherapeutic drugs was analyzed with the expectation of contributing to personalized medicine for BRCA.

2.1 Data collection and preprocessing

The data of TCGA-BRCA dataset in FPKM format were acquired from the University of California, Santa Cruz (UCSC) xena database (https://xena.ucsc.edu/). Excluding samples with no clinical follow-up data or survival time less than 10 days, a total of 1053 BRCA samples and 113 normal samples were obtained. Then, the Ensembl was converted to Gene symbol, and the maximum value was taken to indicate the expressions of multiple gene symbols. The expression matrix was transformed to Trusted Platform Module (TPM) format and the log2 conversion was performed. The TCGA-BRCA dataset was the training set.

The scRNA-seq data and clinical data of GSE20685 dataset, containing 327 BRCA samples, were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20685) database. The probe was mapped to the gene based on the annotation information. When multiple probes matched one gene, the maximum value was to indicate the expression level of the gene. The GSE20685 dataset was the validation set.

The scRNA-seq data of GSE176078 dataset was collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE176078), including 26 BRCA samples (10 TNBC, 11 ER+, and 5 HER2+). The gene expressed in at least 3 cells and the cell expressing at least 200 genes were selected, and the cells with mitochondrial gene ratios <15% and gene numbers between 500-6000 were retained. The SCTransform function was applied to normalize data (Zhang et al., 2020), and the RunPCA function was employed for principal component analysis (PCA) (Sun et al., 2024a). Next, the batch effects between different samples were eliminated using the Harmony R package (Zulibiya et al., 2023).

2.2 Cell clustering analysis and marker genes identification

The top 30 principal components were subjected to the Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) for dimensionality reduction in GSE176078 dataset (Ba et al., 2024). The FindNeighbors and FindClusters functions were utilized for cell clustering analysis, with a resolution of 0.1 for all BRCA cells and 0.3 for Macrophages. Then, each cell type was annotated using the marker genes offered by the CellMarker2.0 database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/). The marker genes that highly expressed in each cell subpopulation were identified using the FindAllMarkers function, with the filtering conditions of logfc.threshold=0.5 and min.pct=0.25.

2.3 Calculation of AUCell enrichment scores

The IRGs were collected from the ImmPort database (https://www.immport.org/shared/), obtaining 2483 IRGs. Then, the AUCell enrichment scores of IRGs for each cell type were computed using the AUCell R package (Weiliang et al., 2024).

2.4 Pathway enrichment analysis

The highly expressed marker genes in Macrophages were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the clusterProfiler R package (Song et al., 2023). The pathways with p<0.05 were recognized as significant enrichment.

2.5 Single sample gene set enrichment analysis (ssGSEA)

The M2-like TAMs scores between normal and tumor samples in TCGA-BRCA dataset were compared by ssGSEA (Fan et al., 2024), and the specific high expression genes of M2-like TAMs was utilized as the background gene set. According to the optimal critical value of M2-like TAMs scores, the patients were split into high- and low- risk groups, and the OS probability of which was assessed.

2.6 Weighted correlation network analysis (WGCNA)

The key genes relevant to M2-like TAMs of BRCA patients in TCGA dataset were screened using the WGCNA R package (Zhang et al., 2023). The optimal soft threshold (power) was determined by the pickSoftThreshold function, and the minimum gene number in each module was set as 80 to obtain the co-expression modules. Next, the module-trait relationships between modules and M2-like TAMs scores were analyzed through Spearman’s correlation analysis, and the key module with the highest correlation coefficient was screened. Then, the correlation between gene significance (GS) for M2-like TAMs scores and module membership (MM) in key module was conducted to identify the M2-like TAMs related genes, with the criteria of |GS|>0.3 and |MM|>0.3.

2.7 Construction and validation of RiskScore model

The prognosis genes of BRCA were identified from the M2-like TAMs related genes by univariate Cox regression analysis (p<0.05). The LASSO Cox regression analysis was further performed using the glmnet R package to reduce gene number (Wang et al., 2023a). Then, the prognostic signatures associated with M2-like TAMs were obtained by stepwise regression analysis, and RiskScore model was established using 3-fold cross validation. The RiskScore for each sample in TCGA-BRCA dataset was calculated by the following formula (Sun et al., 2024b):

βi represents the coefficient of gene in Cox regression model, and ExPi indicates the expression level of gene.

The BRCA patients were separated into high- and low-risk groups based on the optimal critical value of RiskScore, and the OS rate was compared between different risk groups. Meanwhile, the receiver operating characteristic (ROC) curve was plotted using the timeROC R package (Tang et al., 2022). Moreover, the robustness of RiskScore model was validated in the GSE20685 validation set.

2.8 Establishment and evaluation of nomogram

To further investigate whether RiskScore was independent of other clinicopathological features (Age, Gender, stage, pathologic_M, pathologic_N, pathologic_T), the univariate and multivariate Cox regression analysis was conducted in the TCGA-BRCA dataset. Then, a nomogram was established by integrating stage, Age, and RiskScore using the rms R package (Zhang et al., 2022b). The prediction accuracy and reliability of nomogram were evaluated by the calibration curve and decision curve analysis (DCA), respectively.

2.9 Correlation analysis between RiskScore and TME

The relationship between RiskScore and TME was assessed by ESTIMATE and MCP-counter algorithms. The ImmuneScore, StromalScore, and ESTIMATEScore of high/low-risk groups in TCGA-BRCA dataset were calculated using the ‘ESTIMATE’ R package by performing ssGSEA (Mo et al., 2023). Additionally, the infiltration of 10 immune cells was evaluated using the MCP-counter algorithm (Yang et al., 2021), and the correlation between RiskScore and 10 immune cells infiltration was assessed by Spearman method.

2.10 Correlation analysis between RiskScore and chemotherapeutic drugs

The oncoPredict R package was utilized to calculate the IC50 value of each drug in the TCGA-BRCA samples (Maeser et al., 2021). Then, the correlation between RiskScore, key model genes and 8 chemotherapeutic drugs was evaluated by Pearson correlation analysis. P<0.05 and |cor|>0.3 was deemed to be significant correlation. Besides, the IC50 values of Ribociclib_1632 between the two risk groups were compared.

2.11 Statistical analysis

All the statistical data was performed in R language (version 4.2.0). The OS rate was assessed by the Kaplan-Meier method, and the survival difference was determined by the log-rank test. The difference between two groups of continuous variables was compared by Student’s t-test or Wilcoxon rank sum test. A P<0.05 was considered to be statistically significant.

3.1 Macrophage exhibited the highest IRGs AUCell score in single-cell transcriptomic atlas of BRCA

The scRNA-seq data of 26 BRCA samples in GSE176078 was downloaded, after cell filtration, standardization, dimension reduction, and cell clustering analysis, 13 cell clusters were obtained (Fig. 1A). According to the expressions of classical maker genes in each cell cluster, 6 cell types were classified (Fig. 1B), including T cell (marked with CD2, CD3D, CD3E, and CD69), Epithelial cell (marked with EPCAM, KRT18, and KRT8), Macrophage (marked with CD68, CD163, AIF1, FCER1G, LYZ, and TYROBP), Endothelial cell (marked with PECAM1, PLVAP, and VWF), Fibroblast (marked with COL1A1, COL1A2, DCN, and LUM), and B cell (marked with CD79A, DERL3, JCHAIN, MZB1, and MS4A1) (Fig. 1C). Furthermore, the abundance of each cell type in 26 BRCA samples was compared, showing that T cell accounted for a high proportion in the samples (Fig. 1D). In addition, Macrophage exhibited the highest IRGs AUCell score among 6 cell types (Fig. 1E).

Fig. 1 The single-cell transcriptomic atlas of breast cancer (BRCA) in GSE176078 dataset. (A) UMAP plot of different cell clusters. (B) UMAP plot of different cell types. (C) Bubble plot of classical marker genes in each cell type. (D) The proportion of cell types in 26 BRCA samples. (E) AUCell enrichment scores of immune-related genes (IRGs) in different cell types.

3.2 Macrophage C1 highly expressed the marker genes of M2 macrophages

The Macrophage was re-clustered by UMAP clustering and separated into two subpopulations: C1 and C2 (Fig. 2A), which was marked as Macrophage C1 (CD163, FGL2, and MSR1) and Macrophage C2 (CIB1 and KYNU) (Fig. 2B). The IRGs AUCell score of Macrophage C1 was higher than that of Macrophage C2 (Fig. 2C). Moreover, Macrophage C1 highly expressed the marker genes of M2 macrophages such as CD163, FGL2, and MSR1 (Fig. 2D). Thus, Macrophage C1 was selected as M2 macrophages for further study. KEGG enrichment analysis indicated that the specific high expression genes of M2 macrophages were notably enriched in Phagosome, IL-17 signaling pathway, Chemokine signaling pathway, NK-kappa B signaling pathway, Antigen processing and presentation, TNF signaling pathway, Apoptosis, and so on (Fig. 2E).

Fig. 2 The heterogeneity of Macrophage. (A) UMAP plot of Macrophage. (B) The marker genes with high expression in Macrophage C1 and Macrophage C2. (C) AUCell scores of IRGs in Macrophage C1 and Macrophage C2. (D) Violin plot of expression levels of marker genes in M2 macrophages. (E) KEGG pathway enrichment analysis of marker genes in M2 macrophages.

3.3 903 genes highly relevant to M2-like TAMs were identified by WGCNA

The ssGSEA showed that the M2-like TAMs score of tumor group was noticeably lower than that of normal group in TCGA-BRCA dataset (Fig. 3A), manifesting that M2-like TAMs were dysregulated in the tumor samples. According to the Kaplan-Meier survival curve, the OS rate of BRCA patients in high M2-like TAMs score group was higher than that in low M2-like TAMs score group within 10 years, yet the OS rate after 10 years showed the opposite outcome (Fig. 3B). Further, WGCNA was applied to identify the key module correlated with M2-like TAMs of BRCA patients in TCGA-BRCA dataset. The optimal soft threshold (power) was determined to be 6 by pickSoftThreshold function (Fig. S1) and 9 modules were acquired (Fig. S2). The number of genes contained in each module was displayed by lollipop plot (Fig. 3C). The heatmap of module-trait relationships showed that the black module exhibited a remarkable positive correlation with M2-like TAMs score (cor:0.88, p:0.00) (Fig. 3D). Then, 903 genes highly relevant to M2-like TAMs were screened based on the |GS|>0.3 and |MM|>0.3 (Fig. 3E).

Fig. 3 Identification of M2-like tumor-associated macrophages (TAMs) related genes by WGCNA in TCGA-BRCA dataset. (A) M2-like TAMs scores in tumor and normal samples by ssGSEA. ** indicates p < 0.01. (B) OS survival curves between high and low M2-like TAMs score groups. (C) The number of genes included in each module identified by WGCNA. (D) The heatmap of module-trait relationships. (E) Scatter plot of module membership (MM) in black module versus gene significance (GS) for M2-like TAMs score.

3.4 RiskScore model based on 10 M2-like TAMs-related prognostic signatures showed robust performance

First, the prognosis genes were screened from the 903 M2-like TAMs related genes by univariate and LASSO Cox regression analysis, and 17 genes were selected for subsequent study with lambda=0.0172 (Fig. 4A, B). Next, 10 prognostic signatures associated with M2-like TAMs were obtained, including 6 “Protective” genes (ARHGAP26, RILP, KLRB1, CSTA, KLHDC7B, and PSMB8) and 4 “Risk” genes (KYNU, RNASE1, LONRF3, and TRPM2) (Fig. 4C). Then, the RiskScore model was constructed, namely “RiskScore=0.42*TRPM2+0.295*RNASE1-0.266*RILP-0.146*PSMB8+0.353*LONRF3+0.197*KYNU-0.252*KLRB1-0.153*KLHDC7B-0.166*CSTA-0.287*ARHGAP26”. Furthermore, high-risk group had lower OS rate and shorter time of survival than the low-risk group in TCGA-BRCA dataset (Fig. 4D). The area under ROC curve (AUC) analysis suggested that the Riskscore model was reliable with 1-year AUC of 0.73, 3-years AUC of 0.76, 5-years AUC of 0.72, and 7-years AUC of 0.73 (Fig. 4D). The robustness of RiskScore model was validated in the validation set GSE20685, and the results were similar to the TCGA-BRCA dataset (Fig. 4E).

Fig. 4 Establishment and validation of RiskScore model. (A) The coefficients of each independent variable changing with lambda. (B) The confidence interval under lambda. (C) The prognostic signatures of RiskScore model. (D) The ROC and OS survival curves of RiskScore model in TCGA-BRCA dataset. (E) The ROC and OS survival curves of RiskScore model in validation set GSE20685.

3.5 Nomogram exhibited good performance for predicting the prognosis of BRCA

The results of univariate and multivariate Cox regression analyses showed that stage, Age, and RiskScore were the independent prognostic factors (Fig. 5A, B). Then, a nomogram was developed by integrating stage, Age, and RiskScore (Fig. 5C). Further, the predictive calibration points for 1, 3, 5, and 7 years of the calibration curve were close to the standard curve (Fig. 5D), manifesting that the nomogram had strong predictive performance. The net benefit of nomogram was higher than the extreme curve by DCA curve (Fig. 5E), which demonstrated that the nomogram was dependable.

Fig. 5 Construction and validation of nomogram. (A) Univariate Cox regression analysis on RiskScore and clinicopathological features in TCGA-BRCA dataset. (B) Multivariate Cox regression analysis on RiskScore and clinicopathological features in TCGA-BRCA dataset. (C) Construction of nomogram by combining stage, Age, and RiskScore. (D) The calibration curve of nomogram. (E) The DCA curve of nomogram.

3.6 The BRCA patients with low RiskScore were more sensitive to chemotherapeutic drugs

The StromalScore, ImmuneScore, and ESTIMATEScore of high-risk group were considerably lower than the low-risk group in TCGA-BRCA dataset (Fig. 6A), suggesting that the decrease of immune scores may be related to the tumor progression and poor prognosis of patients. Moreover, the correlation between 10 immune cells and RiskScore was assessed by MCP-counter algorithm, founding that most immune cells exhibited a negative correlation with RiskScore (Fig. 6B). In addition, there was a remarkable correlation between RiskScore and IC50 of 8 chemotherapeutic drugs (FDR<0.05 and |cor|>0.3) (Fig. 6C), such as Ribociclib_1632 (R=0.34) (Fig. 6D). The IC50 of Ribociclib_1632 in high-risk group was higher compared to low-risk group (Fig. 6E), which conversely indicated that the BRCA patients with low RiskScore were more sensitive to Ribociclib_1632.

Fig. 6 Correlation analysis between RiskScore and TME or chemotherapeutic drugs in TCGA-BRCA dataset. (A) StromalScore, ImmuneScore, and ESTIMATEScore of high/low-risk groups assessed by ESTIMATE algorithm. (B) Correlation analysis between RiskScore and 10 immune cells infiltration calculated by MCP-counter algorithm. (C) Correlation analysis between RiskScore, key model genes and 8 chemotherapeutic drugs. (D) Correlation analysis between RiskScore and Ribociclib_1632. (E) IC50 values of Ribociclib_1632 between high- and low-risk groups. **** indicates p < 0.0001; *** indicates p < 0.001; ** indicates p < 0.01; * indicates p < 0.05.

M2-like TAMs play a critical part in the development and metastasis of BRCA via regulating immune response and TME, which is correlated with a worse prognosis of BRCA patients (Yang et al., 2018). This study constructed the single-cell transcriptomic atlas of BRCA, and found that Macrophage exhibited the highest IRGs AUCell score. Macrophage was further separated into C1 and C2 subpopulations, and Macrophage C1 showed the characteristics of M2 macrophages. By WGCNA, 903 M2-like TAMs-related genes were obtained. Then, 10 prognostic signatures were identified, including 6 “Protective” genes (ARHGAP26, RILP, KLRB1, CSTA, KLHDC7B, and PSMB8) and 4 “Risk” genes (KYNU, RNASE1, LONRF3, and TRPM2). RiskScore and nomogram model was established, which exhibited good performance for predicting the prognosis of BRCA. Moreover, RiskScore model showed strong correlation with immune cells infiltration and drugs sensitivity.

Numerous studies SHOWED that TAMs-related genes can be served as prognostic signatures of BRCA (Singh et al., 2018). Previous research by Long developed 45 gene signatures from TAMs-regulated BRCA cell MCF7 transcriptome, which were verified to possess favorable performance in predicting the OS and prognostic outcomes for patients suffering from BRCA, and were almost relevant to all clinical features of BRCA (Long et al., 2022). Additionally, Cassetta screened 37 gene signatures associated with human breast TAMs that were significantly enriched in the aggressive BRCA subtypes and related to shorter disease-specific survival (Cassetta et al., 2019). M2-like TAMs were dominant in the breast TME, markedly influencing the progression and metastasis of BRCA (Yang et al., 2023b). Chang screened 722 M2-like TAMs genes using WGCNA, identified two modular subtypes of M2-like TAMs in BRCA, and constructed a prognostic signature for evaluating the prognosis and immunotherapy efficacy of BRCA (Chang et al., 2024). In our present study, we found that Macrophage exhibited the highest IRGs AUCell score, which was further separated into C1 and C2 subpopulations, and Macrophage C1 showed the characteristics of M2 macrophages. Then, 903 M2-like TAMs-related genes were acquired via WGCNA, based on which, 10 prognostic signatures were screened to establish RiskScore model. According to RiskScore, BRCA patients were split into high- and low-risk groups, with those showing a high risk having lower OS rate and poor prognosis. Besides, RiskScore showed strong correlation with immune cells infiltration and drugs sensitivity.

Furthermore, the 10 M2-like TAMs-related prognostic signatures for BRCA contained 6 “Protective” genes (ARHGAP26, RILP, KLRB1, CSTA, KLHDC7B, and PSMB8) and 4 “Risk” genes (KYNU, RNASE1, LONRF3, and TRPM2). ARHGAP26 belongs to a Rho GTPase-activating protein that can transform GTPases into an inactive form, and thereby suppressing the signaling transduction (Zhang et al., 2022a). ARHGAP26 is involved in the progression of numerous cancers, which is considered to be a putative tumor suppressor, such as ovarian cancer (Chen et al., 2019). RILP, a Rab-interacting lysosomal protein, fulfills a crucial fuctions in regulating the transport of lysosomal (Wang et al., 2022b). RILP has been evidenced to inhibit the proliferation and invasion of BRCA (Wang et al., 2015), prostate cancer, and lung cancer cell lines (Lin et al., 2021). KLRB1 is a vital member of killer cell lectin-like receptor family (Xia and Zhou, 2024). It has been proved that KLRB1 is an independent factor for BRCA prognosis, and exhibits low expression in BRCA patients, predicting poor survival outcomes (Huang et al., 2023a; Liu et al., 2024a). CSTA is a cysteine protease inhibitor and an important estrogen-regulated gene in BRCA cells (Lee et al., 2012). The deletion of CSTA in breast tumors can facilitate extracellular matrix remodeling, inducing the invasion and metastasis of BRCA (John Mary et al., 2020). KLHDC7B is a Kelch domain-containing protein and exerts crucial roles in protein degradation, gene expression, and extracellular communication (Ding et al., 2023). KLHDC7B was recognized as an epigenetic marker, highly methylated at the promoter but up-regulated expression in BRCA (Jeong et al., 2018). PSMB8 encodes a critical subunit of the special immunoproteasome complex (Chen et al., 2021). The higher mRNA expression level of PSMB8 was reported to be relevant to the better prognosis of basal-like and TNBC patients (Geoffroy et al., 2023). In addition, KYNU serves as an enzyme that participates in tryptophan catabolism (Liu et al., 2019). The expression of KYNU was associated with the malignancy grade of BRCA and usually indicated poor outcomes (Li et al., 2023). RNASE1 pertains to the pancreatic RNase A superfamily and involves in hemostasis, anti-inflammatory, and innate immunity (Wang et al., 2023b). The overexpression of RNASE1 in BRCA predicts poor survival (Chen et al., 2023). LONRF3 encodes a ring-finger domain-containing protein. LONRF3 showed unfavorable prognostic significance in pancreatic cancer (Feng et al., 2021), yet the role of LONRF3 in BRCA has not been reported. TRPM2 is a key channel of plasmalemma and lysosome in immune and tumor cells (Sumoza-Toledo et al., 2016). TRPM2 was highly expressed in many cancers, such as ovarian cancer, neuroblastoma, and BRCA, suggesting that TRPM2 may facilitate the development of cancers (Huang et al., 2023b). Therefore, the above results manifested that M2-like TAMs-related prognostic gene signatures can be regarded as potential targets for BRCA therapies.

Nevertheless, our current study also has some limitations. Firstly, the data used in the present work is only based on the public database, with a limited sample size, thus further validation is needed by adding clinical samples and detailed experiments data. Second, the prognostic genes involved in the progression of BRCA, such as ARHGAP26 and LONRF3, are rarely reported in BRCA. Thus, the biological functions and molecular mechanisms of the 10 prognostic genes in BRCA should be investigated by more basic experiments and large-scale clinical trials in the future.

In summary, by integrating transcriptome and single-cell sequencing analysis, 10 M2-like TAMs-related prognostic signatures for BRCA were identified, including ARHGAP26, RILP, KLRB1, CSTA, KLHDC7B, PSMB8, KYNU, RNASE1, LONRF3, and TRPM2. RiskScore model was established, with good robustness, which could be an effective tool for the prediction of BRCA prognosis. Moreover, RiskScore showed strong association with immune cell infiltration and drug sensitivity. The present work discovered potential targets and drugs selection for BRCA therapies in clinical practice.

Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Competing Interests The authors have no relevant financial or non-financial interests to disclose.

Author Contributions Z.J.H. and Z.Q.contributed to the study conception and design. Material preparation, data collection and analysis were performed by X.Y.H.C. and L.P.Y.. The first draft of the manuscript was written by L.P.Y. and then Z.Y. modified it. All authors discussed the results and revised the manuscript. All authors reviewed and approved the submission of final manuscripts.

Availability of data and materials The data analysed during the current study are available in UCSC xena database (https://xena.ucsc.edu/) and GEO databases (https://www.ncbi.nlm.nih.gov/geo/).

Ethics approval and consent to participation This manuscript is not a clinical trial, hence the ethics approval and consent to participation are not applicable.

Consent for publication Not applicable.

Ba, Q., Wang, X., Hu, H., and Lu, Y. (2024) Single-Cell RNA Sequencing Analysis Reveals Metabolic Changes in Epithelial Glycosphingolipids and Establishes a Prognostic Risk Model for Pancreatic Cancer. Diagnostics (Basel) 14(11). http://doi.org/10.3390/diagnostics14111094.
Cao, Q., Mushajiang, M., Tang, C.Q., and Ai, X.Q. (2023) Role of hypoxia-inducible factor-1α and survivin in breast cancer recurrence and prognosis. Heliyon 9(3), e14132. http://doi.org/10.1016/j.heliyon.2023.e14132.
Cassetta, L., Fragkogianni, S., Sims, A.H., Swierczak, A., Forrester, L.M., Zhang, H., et al. (2019) Human Tumor-Associated Macrophage and Monocyte Transcriptional Landscapes Reveal Cancer-Specific Reprogramming, Biomarkers, and Therapeutic Targets. Cancer Cell 35(4), 588-602.e510. http://doi.org/10.1016/j.ccell.2019.02.009.
Chang, K., Yue, Q., Jin, L., Fan, P., Liu, Y., Cao, F., et al. (2024) Comprehensive Molecular Analyses of an M2-Like Tumor-Associated Macrophage for Predicting the Prognosis and Immunotherapy in Breast Cancer. J Immunother 47(6), 205-215. http://doi.org/10.1097/cji.0000000000000517.
Chen, D., Jin, C., Dong, X., Wen, J., Xia, E., Wang, Q., et al. (2021) Pan-cancer analysis of the prognostic and immunological role of PSMB8. Sci Rep 11(1), 20492. http://doi.org/10.1038/s41598-021-99724-9.
Chen, X., Chen, S., Li, Y., Gao, Y., Huang, S., Li, H., et al. (2019) SMURF1-mediated ubiquitination of ARHGAP26 promotes ovarian cancer cell invasion and migration. Exp Mol Med 51(4), 1-12. http://doi.org/10.1038/s12276-019-0236-0.
Chen, X., Tang, P., Kong, Y., Chen, D., and Tang, K. (2024) Identification and validation of Golgi apparatus-related signature for predicting prognosis and immunotherapy response in breast cancer. J Cancer Res Clin Oncol 150(2), 61. http://doi.org/10.1007/s00432-024-05612-w.
Chen, X., Yang, C., Wang, W., He, X., Sun, H., Lyu, W., et al. (2023) Exploration of prognostic genes and risk signature in breast cancer patients based on RNA binding proteins associated with ferroptosis. Front Genet 14, 1025163. http://doi.org/10.3389/fgene.2023.1025163.
Cheng, Y., Zhong, X., Nie, X., Gu, H., Wu, X., Li, R., et al. (2023) Glycyrrhetinic acid suppresses breast cancer metastasis by inhibiting M2-like macrophage polarization via activating JNK1/2 signaling. Phytomedicine 114, 154757. http://doi.org/ 10.1016/j.phymed.2023.154757.
Ding, J., Ji, X., Liu, L., Chen, D.Z., Luo, N., Yu, X.T., et al. (2023) A prognostic and immunological analysis of 7B-containing Kelch structural domain (KLHDC7B) in pan-cancer: a potential target for immunotherapy and survival. J Cancer Res Clin Oncol 149(10), 7857-7876. http://doi.org/10.1007/s00432-023-04738-7.
Fan, M., Lu, L., Shang, H., Lu, Y., Yang, Y., Wang, X., et al. (2024) Establishment and verification of a prognostic model based on coagulation and fibrinolysis-related genes in hepatocellular carcinoma. Aging (Albany NY) 16(9), 7578-7595. http://doi.org/10.18632/aging.205699.
Feng, L., Wang, J., Zhang, J., Diao, J., He, L., Fu, C., et al. (2021) Comprehensive Analysis of E3 Ubiquitin Ligases Reveals Ring Finger Protein 223 as a Novel Oncogene Activated by KLF4 in Pancreatic Cancer. Front Cell Dev Biol 9, 738709. http://doi.org/10.3389/fcell.2021.738709.
Fu, Y., Yang, Q., Yang, H., and Zhang, X. (2023) New progress in the role of microRNAs in the diagnosis and prognosis of triple negative breast cancer. Front Mol Biosci 10, 1162463. http://doi.org/10.3389/fmolb.2023.1162463.
Geoffroy, K., Araripe Saraiva, B., Viens, M., Béland, D., and Bourgeois-Daigneault, M.C. (2023) Increased expression of the immunoproteasome subunits PSMB8 and PSMB9 by cancer cells correlate with better outcomes for triple-negative breast cancers. Sci Rep 13(1), 2129. http://doi.org/10.1038/s41598-023-28940-2.
Guo, S., Ma, Y., Li, X., Li, W., He, X., Yuan, Z., et al. (2023) Identification of stromal cell proportion-related genes in the breast cancer tumor microenvironment using CorDelSFS feature selection: implications for tumor progression and prognosis. Front Genet 14, 1165648. http://doi.org/10.3389/fgene.2023.1165648.
Hu, J., Lv, H., Zhao, S., Lin, C.J., Su, G.H., and Shao, Z.M. (2023) Prediction of clinicopathological features, multi-omics events and prognosis based on digital pathology and deep learning in HR(+)/HER2(-) breast cancer. J Thorac Dis 15(5), 2528-2543. http://doi.org/10.21037/jtd-23-445.
Huang, G., Xiao, S., Jiang, Z., Zhou, X., Chen, L., Long, L., et al. (2023a) Machine learning immune-related gene based on KLRB1 model for predicting the prognosis and immune cell infiltration of breast cancer. Front Endocrinol (Lausanne) 14, 1185799. http://doi.org/10.3389/fendo.2023.1185799.
Huang, W., Wu, Y., Luo, N., Shuai, X., Guo, J., Wang, C., et al. (2023b) Identification of TRPM2 as a prognostic factor correlated with immune infiltration in ovarian cancer. J Ovarian Res 16(1), 169. http://doi.org/10.1186/s13048-023-01225-y.
Huang, X., Cao, J., and Zu, X. (2022) Tumor-associated macrophages: An important player in breast cancer progression. Thorac Cancer 13(3), 269-276. http://doi.org/10.1111/1759-7714.14268.
Huang, Y., Zhang, Y., Duan, X., Hou, R., Wang, Q., and Shi, J. (2024) Exploring the immune landscape and drug prediction of an M2 tumor-associated macrophage-related gene signature in EGFR-negative lung adenocarcinoma. Thorac Cancer. http://doi.org/10.1111/1759-7714.15375.
Jeong, G., Bae, H., Jeong, D., Ham, J., Park, S., Kim, H.W., et al. (2018) A Kelch domain-containing KLHDC7B and a long non-coding RNA ST8SIA6-AS1 act oppositely on breast cancer cell proliferation via the interferon signaling pathway. Sci Rep 8(1), 12922. http://doi.org/10.1038/s41598-018-31306-8.
John Mary, D.J.S., Sikarwar, G., Kumar, A., and Limaye, A.M. (2020) Interplay of ERα binding and DNA methylation in the intron-2 determines the expression and estrogen regulation of cystatin A in breast cancer cells. Mol Cell Endocrinol 504, 110701. http://doi.org/10.1016/j.mce.2020.110701.
Kim, S., Choi, I., Han, I.H., and Bae, H. (2022) Enhanced Therapeutic Effect of Optimized Melittin-dKLA, a Peptide Agent Targeting M2-like Tumor-Associated Macrophages in Triple-Negative Breast Cancer. Int J Mol Sci 23(24). http://doi.org/10.3390/ijms232415751.
Kodous, A.S., Balaiah, M., and Ramanathan, P. (2023) Single cell RNA sequencing – a valuable tool for cancer immunotherapy: a mini review. 25(6), 635-639. http://doi.org/10.1515/oncologie-2023-0244.
Lee, S., Stewart, S., Nagtegaal, I., Luo, J., Wu, Y., Colditz, G., et al. (2012) Differentially expressed genes regulating the progression of ductal carcinoma in situ to invasive breast cancer. Cancer Res 72(17), 4574-4586. http://doi.org/10.1158/0008-5472.Can-12-0636.
Li, Y., Wang, M., Zhao, L., Liang, C., and Li, W. (2023) KYNU-related transcriptome profile and clinical outcome from 2994 breast tumors. Heliyon 9(6), e17216. http://doi.org/10.1016/j.heliyon.2023.e17216.
Li, Y., Zhao, X., Liu, Q., and Liu, Y. (2021) Bioinformatics reveal macrophages marker genes signature in breast cancer to predict prognosis. Ann Med 53(1), 1019-1031. http://doi.org/10.1080/07853890.2021.1914343.
Lin, J., Zhuo, Y., Yin, Y., Qiu, L., Li, X., and Lai, F. (2021) Methylation of RILP in lung cancer promotes tumor cell proliferation and invasion. Mol Cell Biochem 476(2), 853-861. http://doi.org/10.1007/s11010-020-03950-0.
Liu, X., Cui, Q., and Qin, N. (2024a) Low expression of KLRB1 predicts poor survival outcomes and is associated with immune infiltration in breast cancer. Transl Cancer Res 13(3), 1225-1240. http://doi.org/10.21037/tcr-23-1231.
Liu, X., Sun, K., Yang, H., Zou, D., Xia, L., Lu, K., et al. (2024b) Molecular subtype identification and prognosis stratification based on lysosome-related genes in breast cancer. Heliyon 10(4), e25643. http://doi.org/10.1016/j.heliyon.2024.e25643.
Liu, Y., Feng, X., Lai, J., Yi, W., Yang, J., Du, T., et al. (2019) A novel role of kynureninase in the growth control of breast cancer cells and its relationships with breast cancer. J Cell Mol Med 23(10), 6700-6707. http://doi.org/10.1111/jcmm.14547.
Long, M., Wang, J., and Yang, M. (2022) Transcriptomic Profiling of Breast Cancer Cells Induced by Tumor-Associated Macrophages Generates a Robust Prognostic Gene Signature. Cancers (Basel) 14(21). http://doi.org/10.3390/cancers14215364.
Lu, X., Yuan, Q., Zhang, C., Wang, S., and Wei, W. (2023) Predicting the immune microenvironment and prognosis with a anoikis - related signature in breast cancer. Front Oncol 13, 1149193. http://doi.org/10.3389/fonc.2023.1149193.
Maeser, D., Gruener, R.F., and Huang, R.S. (2021) oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 22(6). http://doi.org/10.1093/bib/bbab260.
Mo, B., Zhao, X., Wang, Y., Jiang, X., Liu, D., and Cai, H. (2023) Pan-cancer analysis, providing a reliable basis for IDO2 as a prognostic biomarker and target for immunotherapy. 25(1), 17-35. http://doi.org/10.1515/oncologie-2022-1026.
Nie, X., Fu, L., Cheng, Y., Wu, X., Lv, K., Li, R., et al. (2023) Garcinone E suppresses breast cancer growth and metastasis by modulating tumor-associated macrophages polarization via STAT6 signaling. Phytother Res 37(10), 4442-4456. http://doi.org/10.1002/ptr.7909.
Singh, R., Dagar, P., Pal, S., Basu, B., and Shankar, B.S. (2018) Significant alterations of the novel 15 gene signature identified from macrophage-tumor interactions in breast cancer. Biochim Biophys Acta Gen Subj 1862(3), 669-683. http://doi.org/10.1016/j.bbagen.2017.12.004.
Song, Z., Yu, J., Wang, M., Shen, W., Wang, C., Lu, T., et al. (2023) CHDTEPDB: Transcriptome Expression Profile Database and Interactive Analysis Platform for Congenital Heart Disease. 18(6), 693--701.
Sumoza-Toledo, A., Espinoza-Gabriel, M.I., and Montiel-Condado, D. (2016) Evaluation of the TRPM2 channel as a biomarker in breast cancer using public databases analysis. Bol Med Hosp Infant Mex 73(6), 397-404. http://doi.org/10.1016/j.bmhimx.2016.10.001.
Sun, H., Chang, Z., Li, H., Tang, Y., Liu, Y., Qiao, L., et al. (2024a) Multi-omics analysis-based macrophage differentiation-associated papillary thyroid cancer patient classifier. Transl Oncol 43, 101889. http://doi.org/10.1016/j.tranon.2024.101889.
Sun, J., Wang, Y., Zhang, K., Shi, S., Gao, X., Jia, X., et al. (2024b) Molecular subtype construction and prognosis model for stomach adenocarcinoma characterized by metabolism-related genes. Heliyon 10(7), e28413. http://doi.org/10.1016/j.heliyon.2024.e28413.
Tan, X.J., Cheor, W.L., Cheng, E.M., Rahman, K.S.A., Muhamad, W.Z.A.W., and Leow, W.Z. (2023) Breast cancer status, grading system, etiology, and challenges in Asia: an updated review. 25(2), 99-110. http://doi.org/10.1515/oncologie-2022-1011.
Tang, X., Qi, C., Zhou, H., and Liu, Y. (2022) A novel metabolic-immune related signature predicts prognosis and immunotherapy response in lung adenocarcinoma. Heliyon 8(8), e10164. http://doi.org/10.1016/j.heliyon.2022.e10164.
Wang, H., Wang, X., and Xu, L. (2024a) Transforming growth factor-induced gene TGFBI is correlated with the prognosis and immune infiltrations of breast cancer. World J Surg Oncol 22(1), 22. http://doi.org/10.1186/s12957-024-03301-z.
Wang, R., Zhang, X., He, C., and Guo, W. (2023a) An effective prognostic model for assessing prognosis of non-small cell lung cancer with brain metastases. Front Genet 14, 1156322. http://doi.org/10.3389/fgene.2023.1156322.
Wang, S., Jiang, H., Zheng, C., Gu, M., and Zheng, X. (2022a) Secretion of BMP-2 by tumor-associated macrophages (TAM) promotes microcalcifications in breast cancer. BMC Cancer 22(1), 34. http://doi.org/10.1186/s12885-021-09150-3.
Wang, S., Wang, J., Chen, Z., Luo, J., Guo, W., Sun, L., et al. (2024b) Targeting M2-like tumor-associated macrophages is a potential therapeutic approach to overcome antitumor drug resistance. NPJ Precis Oncol 8(1), 31. http://doi.org/10.1038/s41698-024-00522-z.
Wang, Y.N., Lee, H.H., Jiang, Z., Chan, L.C., Hortobagyi, G.N., Yu, D., et al. (2023b) Ribonuclease 1 Enhances Antitumor Immunity against Breast Cancer by Boosting T cell Activation. Int J Biol Sci 19(10), 2957-2973. http://doi.org/10.7150/ijbs.84592.
Wang, Y.Y., Cheng, K.H., Hung, A.C., Lo, S., Chen, P.Y., Wu, Y.C., et al. (2024c) Differential impact of cytoplasmic vs. nuclear RAD51 expression on breast cancer progression and patient prognosis. Int J Oncol 64(2). http://doi.org/10.3892/ijo.2023.5600.
Wang, Z., Zhou, Y., Hu, X., Chen, W., Lin, X., Sun, L., et al. (2015) RILP suppresses invasion of breast cancer cells by modulating the activity of RalA through interaction with RalGDS. Cell Death Dis 6(10), e1923. http://doi.org/10.1038/cddis.2015.266.
Wang, Z., Zhou, Y., Nie, D., Tan, Y., Zhao, S., Wang, G., et al. (2022b) RILP inhibits proliferation, migration, and invasion of PC3 prostate cancer cells. Acta Histochem 124(7), 151938. http://doi.org/10.1016/j.acthis.2022.151938.
Weiliang, W., Yinghao, R., Weiliang, H., Xiaobin, Z., Chenglong, Y., Weimiao, A., et al. (2024) Identification of hub genes significantly linked to tuberous sclerosis related-epilepsy and lipid metabolism via bioinformatics analysis. Front Neurol 15, 1354062. http://doi.org/10.3389/fneur.2024.1354062.
Xia, J., and Zhou, X. (2024) Necroptosis-related KLRB1 was a potent tumor suppressor and immunotherapy determinant in breast cancer. Heliyon 10(6), e27294. http://doi.org/ 10.1016/j.heliyon.2024.e27294.
Yang, H., Zhao, L., Zhang, Y., and Li, F.F. (2021) A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma. Cancer Med 10(16), 5696-5711. http://doi.org/10.1002/cam4.4117.
Yang, M., Li, Z., Ren, M., Li, S., Zhang, L., Zhang, X., et al. (2018) Stromal Infiltration of Tumor-Associated Macrophages Conferring Poor Prognosis of Patients with Basal-Like Breast Carcinoma. J Cancer 9(13), 2308-2316. http://doi.org/10.7150/jca.25155.
Yang, P., Yang, Z., Dong, Y., Yang, L., Peng, S., Yuan, L., et al. (2023a) Clusterin is a biomarker of breast cancer prognosis and correlated with immune microenvironment. Transl Cancer Res 12(1), 31-45. http://doi.org/10.21037/tcr-22-1882.
Yang, Y., Lu, T., Jia, X., and Gao, Y. (2023b) FSTL1 Suppresses Triple-Negative Breast Cancer Lung Metastasis by Inhibiting M2-like Tumor-Associated Macrophage Recruitment toward the Lungs. Diagnostics (Basel) 13(10). http://doi.org/10.3390/diagnostics13101724.
Zhang, L., Zhou, A., Zhu, S., Min, L., Liu, S., Li, P., et al. (2022a) The role of GTPase-activating protein ARHGAP26 in human cancers. Mol Cell Biochem 477(1), 319-326. http://doi.org/10.1007/s11010-021-04274-3.
Zhang, W., Zhang, P., Jiang, J., Peng, K., Shen, Z., and Kang, M. (2022b) Development and validation of a prognostic model related to pyroptosis-related genes for esophageal squamous cell carcinoma using bioinformatics analysis. J Thorac Dis 14(8), 2953-2969. http://doi.org/10.21037/jtd-22-948.
Zhang, X., Chen, T., Qian, X., and He, X. (2023) Bioinformatics Analysis of Immune Cell Infiltration and Diagnostic Biomarkers between Ankylosing Spondylitis and Inflammatory Bowel Disease. Comput Math Methods Med 2023, 9065561. http://doi.org/10.1155/2023/9065561.
Zhang, Y., Ma, Y., Huang, Y., Zhang, Y., Jiang, Q., Zhou, M., et al. (2020) Benchmarking algorithms for pathway activity transformation of single-cell RNA-seq data. Comput Struct Biotechnol J 18, 2953-2961. http://doi.org/10.1016/j.csbj.2020.10.007.
Zhou, Y., Xue, W., Meng, X., Bhandari, A., Zeng, H., Kc, R., et al. (2024) GNPNAT1 is a Biomarker That Predicts a Poor Prognosis of Breast Cancer. Breast Cancer (Dove Med Press) 16, 71-89. http://doi.org/10.2147/bctt.S451054.
Zulibiya, A., Wen, J., Yu, H., Chen, X., Xu, L., Ma, X., et al. (2023) Single-Cell RNA Sequencing Reveals Potential for Endothelial-to-Mesenchymal Transition in Tetralogy of Fallot. 18(6), 611--625.

No competing interests reported.

SupplementaryInformation.pdf

Download PDF

Version 1

posted

You are reading this latest preprint version

Integrated transcriptome and single-cell sequencing reveals M2-like tumor-associated macrophages-related prognostic signatures in breast cancer

Status:

Version 1

Abstract

Figures

1 Introduction

2 Materials and methods

2.1 Data collection and preprocessing

2.2 Cell clustering analysis and marker genes identification

2.3 Calculation of AUCell enrichment scores

2.4 Pathway enrichment analysis

2.5 Single sample gene set enrichment analysis (ssGSEA)

2.6 Weighted correlation network analysis (WGCNA)

2.7 Construction and validation of RiskScore model

2.8 Establishment and evaluation of nomogram

2.9 Correlation analysis between RiskScore and TME

2.10 Correlation analysis between RiskScore and chemotherapeutic drugs

2.11 Statistical analysis

3 Results

3.1 Macrophage exhibited the highest IRGs AUCell score in single-cell transcriptomic atlas of BRCA

3.2 Macrophage C1 highly expressed the marker genes of M2 macrophages

3.3 903 genes highly relevant to M2-like TAMs were identified by WGCNA

3.4 RiskScore model based on 10 M2-like TAMs-related prognostic signatures showed robust performance

3.5 Nomogram exhibited good performance for predicting the prognosis of BRCA

3.6 The BRCA patients with low RiskScore were more sensitive to chemotherapeutic drugs

4 Discussion

5 Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1