Bioinformatics reveal macrophages marker genes signature in breast cancer to predict prognosis

Data source and acquisition

ScRNA-seqdata in the form of RSEM normalized counts from six primary triple-negative breast cancers was downloaded fromGSE118389 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE118389) [18], and used to screen macrophages marker genes. Normalized RNA sequencing data in the form of log₂ (FPKM + 1) was downloaded from UCSC Xena (https://xenabrowser.net/datapages/) for further survival-related genes screening and model construction. Individual patient files and mRNA expression raw data were obtained from TCGA data portal (https://portal.gdc.cancer.gov/). To validate the prognostic ability of the constructed model, normalized RNA sequencing data in the form of log₂ (FPKM + 0.1) was downloaded from GSE96058 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96058) [24,25]. DNA methylation raw data were obtained from TCGA data portal.

Identification of macrophages marker genes and functional analysis

ScRNA-seq data analysis was performed by using “Seurat”, “SingleR” packages [26]. Cells with more than 5% of mitochondrial gene were removed [27]. Cells with number of gene mapped less than 200 and clusters with cell counts less than 5 were moved. We performed principal component analysis (PCA) using the most 1500 variable genes in the dataset in order to visualize transcriptional variability over the complete scRNA-seq dataset. T-distributed Stochastic Neighbour Embedding (t-SNE) was used for further dimensional reduction of the significant principal components [27]. Genes that exhibited a |log₂ (fold change)|>0.8 and adjusted p value < .01 were considered as the marker genes. Kyoko Encyclopaedia of Genes and Genomes (KEGG) pathway enrichments and gene ontology (GO) analysis were conducted by using “ClusterProfiler” [28], “org.Hs.eg.db”, “GOplot” [29], “enrichplot” packages.

Construction and validation of macrophages marker genes prognosis risk model

Cases with follow-up time more than 30 days in the TCGA database were included to build the risk score model. Macrophages marker genes expression data downloaded from UCSC Xena were merged with overall survival (OS) time and status for each case. To screen out the most significant genes account for OS, p value was set as less than .01 in the univariate Cox analysis. Statistically significant genes in the univariate Cox analysis were included to build multivariate Cox proportional hazards model. Further, we detected whether the genes in the signature model differentially expressed between the tumour tissue and tumour adjacent normal tissue in the TCGA dataset. X-tile software (version 3.6.1) was used to define the optimum cut-off value for risk scores based on the association with OS [30].”Survival”, “survminer” packages were used to construct Kaplan–Meier survival curves to evaluate survival differences between high-risk and low-risk groups both in the training and validation set. Then the capacity of the risk score to predict OS was tested by adjusting for classical clinical variables in the multivariate Cox model in both training and validation groups. Furthermore, time-dependent receiver operating characteristic (ROC) curve, area under the curve (AUC) and the concordance index (C-index) were performed to evaluate the discrimination power of this model [31]. Then, the heatmaps of these model genes in both training and validation groups were performed.

Correlation between gene methylation and expression

Perl 64 was used to merge gene methylation and gene expression data. Then we analysed the correlation between methylation and gene expression by using correlation test in R [32].

Statistical analysis

GSE96058 dataset normalized RNA sequencing data were transformed into the same format as TCGA data before model validation. All figure construction in this study was conducted by using R package software (version 4.0.3). Univariate and multivariate Cox regression model were performed by using “survival”, “survminer” packages. Wilcox test was used to determine statistical differences of categorical variables. Genes included in the multivariate analysis were selected through stepwise Cox regression analysis with both directions. Bootstrap method was used to perform internal validation in the TCGA dataset. The survival curves were measured by the Kaplan–Meier method and the significance of disparity was assessed by log-rank test. “timeROC”, “boot” packages were used to evaluate the prediction capacity of the prognostic model. Pearson correlation coefficient |r|>0.3 and p < .05 were defined as closely correlated.

Identification of macrophages marker genes expression profiles and function annotation

According to the screening criteria, a total of 1205 cells from 6 samples were analysed to identify and characterize cell populations. We first reduced the dimensionality of the data by PCA by using the 1500 variable genes (Figure 1(a)). We then assigned cell-type identities by cross-referencing differentially expressed genes in each cluster with previously reported cell-type-specific marker genes (with |log₂ (fold change)|>0.8 as well as adjusted p value <.01) [18,33], and identified 14 cell clusters using Seurat. Cells in cluster 7 expressing macrophage-specific markers were classified as macrophages (Figure 1(b)). We also found that this cluster had distinct gene expression profiles with a subset of genes differentially expressed between the 14 clusters (Figure 1(c)). As a result, we found out breast cancer-related 314 macrophage marker genes.

A total of 314 macrophage marker genes were screened from GSE118389 dataset according to the screening criteria. Then GO analysis (Figure 2(a)) and KEGG pathway enrichment (Figure 2(b)) were performed to explore the biological function of these marker genes. We found that macrophages marker genes were mostly involved in the process of neutrophil activation, phagocytosis, macrophage activation, etc. (Figure 2(a)). The KEGG pathway enrichment demonstrated that macrophages marker genes primarily participated in the pathway of phagosome, antigen processing and presentation, complement and coagulation cascades, etc.

Construction of prognostic macrophage marker genes signature (MMGS)

In order to identify a macrophages gene prognostic signature, we used TCGA cohort as the training set. According to the screen criteria previously described, there were a total of 1034 patients included in the training cohort with follow-up time ranging from 31 to 8605 days, with a median follow-up of 865 days. The clinicopathological characteristics of included cases were shown in Supplemental Table S1.

The univariate Cox proportional hazards regression analysis of TCGA database revealed that10macrophage marker genes with significantly associated with OS were selected as candidate genes of MMGS (Figure 2(c)). These 10 marker genes were then incorporated into a multivariate Cox proportional hazards regression model to determine the genes and their coefficients. Ultimately, 7 macrophage marker genes were included in the MMGS model (Figure 2(d)). The risk score which was used to predict prognosis was given as follows: MMGS risk score = SERPINA1expression × (−0.155) + CD74expression × 0.222 + STX11expression × (−0.572) + ADAM9expression × 0.190 + CD24expression × 0.107 + NFKBIAexpression × (−0.371) + PGK1expression × 0.360. The relative expression of the 7 marker genes in various clusters were shown in Figure 2(e), and the marker genes except CD24 were relatively higher in the macrophages cluster (cluster 7) compared to the other clusters as a whole.

Since macrophages in TME were activated to pro-tumoral phenotype, we assumed the MMGS risk score would discriminate macrophage characters between tumour and normal tissues. We detected the relative marker genes expression between tumour and tumour adjacent normal tissue. As shown in Supplemental Figure 1, ADAM9 and SERPINA1 expression between tumour and tumour adjacent normal tissue were not significantly different (p = .634, .491, respectively). CD24, CD74, PGK1 showed significantly higher expression in tumour tissue compared with tumour adjacent normal tissue, while STX11 and NFKBIA showed the opposite trend (p < .001 for all). Risk score showed significantly higher expression in tumour tissue compared with tumour adjacent normal tissue (p < .001).

MMGS risk score was calculated for each individual patient from the TCGA cohort (Figure 3(a)). The heatmaps of model genes in the training set and the validation set are shown in Supplemental Figure S2(a, b), respectively. The cut-off of risk score (2.8425) generated by X-tile software was set to divide patients into high- and low-risk groups (Figure 3(b)). Patients in the high-risk group had a significantly shorter OS than those in the low-risk group (hazard ratio: 3.077 [95%confidence interval (CI): 1.912–4.954], p < .001) (Figure 3(c)). A time-dependent ROC analysis demonstrated that the AUC for 3-year and 5-year OS of this classifier were 0.662 and 0.701, respectively, indicating MMGS possesses good sensitivity and specificity in predicting the prognosis in training set (Figure 4(a)). The C-index of MMGS for OS prediction in training set was 0.666 (95%CI: 0.609–0.723).

Validation of the prognostic value of MMGS

We first performed internal validation by bootstrap method (B = 1000), and the C-index of bootstrap validation of the prediction model was 0.667 (95% CI: 0.589–0.745), indicating good distinguishing capacity of this model. In an external independent cohort from the GSE96058 validation dataset, MMGS risk score was calculated for each individual patient (Figure 3(d)), and patients were classified into high- and low-risk subgroups based on the same risk score cut-off value previously mentioned (Figure 3(e)). The OS in the high-risk group was significantly shorter than that of the low-risk one (hazard ratio: 1.808 [95%CI: 1.458–2.243], p < .001) in the validation set (Figure 3(f)). AUC for 3-year and 5-year time-dependent ROC were 0.646 and 0.599, respectively (Figure 4(b)). The C-index of MMGS for predicting OS invalidation set was 0.618 (95%CI: 0.586–0.651).

The correlation of MMGS risk score with classical clinical variables

To investigate the correlation of MMGS risk score with clinicopathological factors, including age, tumour size, lymph node status, ER status, PR status, HER2 status, we calculated MMGS risk scores distribution in patients from TCGA database stratified by each clinical risk factors. We found the risk score was lower in the ER and PR positive groups while higher in HER2 positive group in training set (p < .001) (Figure 5). The same results were observed in the validation set (Figure 6). While the correlations of MMGS risk score with age, tumour size, lymph node status were inconsistent in training set and validation set (Figures 5 and 6). Then we explored the prognostic value of MMGS in different subgroups. In the training set, a significant poor survival probability was associated with high MMGS risk score in patients with different clinicopathological factors except age (≤40) and HER2 positive subgroups (Figure 7(a)). In the validation set, MMGS risk score effectively predict patients’ prognosis in hormone receptor positive and HER2 negative subgroups as well as all age and lymph node status subgroups (Figure 7(b)).

Next, we determine whether MMGS can be employed to independently predict the survival of breast cancer patients. MMGS and clinical characteristics, including age, tumour size, lymph node status, ER status, PR status, HER2 status, were included in univariate and multivariate analysis. MMGS remained an independent prognostic factor after adjusting for clinical characteristics in the multivariate analysis with the hazard ratio of 3.185 (95%CI: 1.983–5.115, p < .001) in training set (Table 1) and 1.362 (95%CI: 1.055–1.757, p = .018) in the validation cohort (Table 2).

The association between gene methylation and expression in MMGS

It is said that gene transcription could be regulated by DNA methylation level [22,23]. To lay the foundation for further research of these model genes, we obtained the gene methylation data from TCGA database and investigated the association between gene methylation and expression levels in the MMGS model. As shown in Figure 7(c), CD74 and SERPINA1 promoter hypermethylation is significantly correlated with lower gene expression (p < .001, cor= −0.532 for CD74 and p < .001, cor= −0.429 for SERPINA1), while STX11, ADAM9, NFKBIA, PGK1 promoter hypermethylation is not significantly correlated with gene expression (|r|<0.3, Supplemental Figure 3). CD24 promoter hypermethylation data were not obtained in the TCGA dataset, so the relationship remains unclear.