Research Paper Volume 13, Issue 19 pp 23210—23232
Integrated analysis of tumor-associated macrophage infiltration and prognosis in ovarian cancer
- 1 Department of Gynecology and Obstetrics, Xiangya Hospital Central South University, Changsha, Hunan, China
- 2 Department of Nephrology, Xiangya Hospital Central South University, Changsha, Hunan, China
Received: June 18, 2021 Accepted: September 28, 2021 Published: October 11, 2021https://doi.org/10.18632/aging.203613
How to Cite
Copyright: © 2021 Tan et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Ovarian cancer (OC) is a frequently lethal gynecologic malignancy, characterized by a poor prognosis and high recurrence rate. The immune microenvironment has been implicated in the progression of OC. We characterized the immune landscape in primary and malignant OC ascites using single-cell and bulk transcriptome raw OC data acquired from the Gene Expression Omnibus and The Cancer Genome Atlas databases. We then used the CIBERSORT deconvolution algorithm, weighted gene co-expression network analysis, univariate and multivariate Cox analyses, and the LASSO algorithm to develop a tumor-associated macrophage-related gene (TAMRG) prognostic signature, which enabled us to stratify and predict overall survival (OS) of OC patients. In addition, inter- and intra-patient heterogeneity of infiltrating immune cells was characterized at single-cell resolution. Tumor-infiltrating macrophages with an M2 phenotype exhibited immunosuppressive activity. M1 macrophages positively correlated with OS, whereas activated mast cells, neutrophils, M2 macrophages, and activated memory CD4+ T cells were all negatively correlated with OS. A total of 219 TAMRGs were identified, and a novel 6-gene signature (TAP1, CD163, VSIG4, IGKV4-1, CD3E, and MS4A7) with independent prognostic value was established. These results show that a TAMRG-based signature may be a promising prognostic and therapeutic target in OC.
With 13,770 estimated deaths in the United States in 2021, ovarian cancer (OC) has emerged as a highly lethal gynecologic malignancy . Almost 80% of OC cases are diagnosed in advanced stages, and its 5-year survival rate is about 48% . Although the standard platinum-based chemotherapy and cytoreductive surgery result in complete remission, cancer tends to relapse and disseminate to distant organs in most patients .
Recent studies have shown the contribution of the tumor microenvironment (TME) to OC metastasis . Compared to other solid tumors, the TME of epithelial OC is unusual because the cancer cells frequently escape the primary tumor site to create a microenvironment in the peritoneal cavity, known as malignant ascites . Several therapeutic approaches based on angiogenesis, tumor-associated macrophages, cancer-associated fibroblasts, and immune checkpoint blockade are being devised to target the TME [6, 7].
The diversified TME in OC is a manifestation of its high heterogeneity. Jiménez-Sánchez et al. reported both regressing and progressing metastases, with different immune molecular patterns, in the same patient with high-grade serous ovarian cancer (HGSOC) who had undergone multiple chemotherapies . They further reported the co-existence of inflammatory and immune cell-excluded microenvironments in patients with untreated HGSOC, suggesting widespread variation in infiltrating immune cells . Therefore, identifying the subset of highly metastatic cancer cells and OC cell type diversity at the single-cell level is necessary for developing effective clinical biomarkers and treatment strategies for OC.
A TME comprises myeloid-derived suppressor cells (MDSCs), lymphocytes, macrophages, mast cells, neutrophils, and dendritic cells (DCs) and contributes to tumor growth, invasion, and metastasis [10, 11]. Moreover, abundant tumor-associated macrophage (TAM) infiltrates are associated with poor prognoses . The ovarian cancer TME is mostly populated with TAM that exhibits functional plasticity and polarizes into M2-like cells [5, 13]. The immunosuppressive M2-like phenotype induces inflammation, tissue remodeling, and tumor angiogenesis .
Identifying prognostic markers of OC using large-scale public gene expression data may provide novel therapeutic targets for OC. We integrated single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq to develop and validate a TAM-based prognostic signature for OC. The TME of HGSOC ascites samples, identified using scRNA-seq data, was highly enriched in M2-like macrophages. TAP1, CD163, VSIG4, IGKV4-1, CD3E, and MS4A7 were identified as the major OS-predicting gene signatures using the bulk RNA-seq data. Finally, a prognostic model was validated using a GEO validation dataset.
Tumor cell heterogeneity in ovarian cancer
Figure 1 shows the schematic illustration of the study design. A total of 9,609 cells from ascites samples of patients with OC were obtained for the analysis after quality control (Figure 2A and 2B). Variance analysis showed 1,500 highly variable genes among 10,048 genes. The top 10 genes were identified as IGLL5, IGJ, WT1-AS, CCL17, GNLY, CCL5, MZB1, HBB, and NKG7 (Figure 2C). The principal component analysis (PCA) showed a mixed representation of intra-and inter-patient cells (Figure 2D). The p-value of the first 20 principal components (PCs) was less than 0.05 (Figure 2E). The highly related genes in the top four PCs are shown in Supplementary Figure 1. Afterward, t-distributed stochastic neighbor embedding (tSNE) was performed with PC 1–20 at a resolution of 0.3 to group the OC cells into 13 separate clusters (Figure 2F). A total of 6,453 marker genes were identified and the top 10 differentially expressed genes from the 13 clusters were displayed in the heatmap (Figure 2G).
Figure 2. Heterogeneity in patients with ovarian cancer (OC) was identified using single-cell RNA-seq data. (A and B) A total of 9,609 cells from six OC patients were included in the analysis. (C) Scatter plots displayed 10,048 corresponding genes in all cells from OC samples. (D) The principal component analysis (PCA) revealed unclear separations of OC cells. (E) The first 20 principal components with a p-value < 0.05 were generated by PCA. (F) OC cells were categorized into 13 clusters using the tSNE algorithm with the first 20 principal components. (G) The heatmap shows the top 10 differential marker genes of each cluster. A total of 124 unique genes were identified after removing the same marker genes among the clusters.
Identification of immune cells and GSVA
We next annotated the cell clusters as epithelial, immune, or stromal (fibroblasts) cells (Figure 3A) (Supplementary Table 1) using previously established cell surface markers. Clusters 9 and 10, containing 419 cells, were classified as epithelial cells; clusters 2 and 7, containing 1,895 cells, were classified as stromal cells; and clusters 0, 1, 3–8, 11, and 12, containing 6,706 cells, were classified as immune cells. We next extracted immune cells (n = 6,706) separately and annotated them as B cells, macrophages, dendritic cells (DCs), and T cells (Figure 3B). Clusters 0, 1, 3, and 4, containing 5,930 cells, were classified as macrophages; cluster 6, containing 384 cells, was classified as B cells; cluster 8, containing 320 cells, was annotated as T cells; and cluster 12, containing 72 cells, was classified as DCs.
Figure 3. Clustering of immune cell populations and GSVA enrichment scores. (A) Cell-types visualized using tSNE dimensionality reduction revealed the clustering of tumor-stroma immune cells. (B) Immune cells in the tumor microenvironment were annotated into four subpopulations. (C) Hierarchical clustering was used using GSVA enrichment scores for gene set for B cell signature, DC signature, macrophage signature, M1 up signature, M2 up signature, and T cell signature. (D) Violin plot indicates the genes corresponding to immune cell subpopulations.
Gene set variation analysis (GSVA) was performed to analyze the B cell functional status, including cytokine production, naïveness, anti-apoptotic, proliferation, pro-apoptotic functions, and germinal center characteristics. The functional status related to pro-apoptosis showed heterogeneity in the same patient (Supplementary Figure 2). For T cells, we estimated the cytotoxic, native, regulatory, exhausted, co-stimulatory, and G1/S- and G2/M-related gene signatures. The functional status of cytotoxicity in patient 2 was highly enriched (Supplementary Figure 3). Figure 3C shows a comprehensive analysis of the four kinds of immune cells. Macrophages constituted the largest proportion (88.4%) of immune cells. The M2 up signature was highly enriched in macrophages (Figure 3C). Macrophages were marked by CD163, CSF1R, MS4A7, and VSIG4; B cells were marked by CD79B and MZB1; T cells were marked by CD2 and CD3D; and DCs were marked by CD1E and CD83 (Figure 3D). In summary, B cells, T cells, and macrophages showed heterogeneity in the same patient and among different patients.
High M2 and low M1 phenotypes were correlated with poor survival in ovarian cancer patients
The CIBERSORT-based Nu support vector regression algorithm was used to evaluate the immune score of the transcribed dataset of bulk RNA-seq obtained from a TCGA-OC cohort. The estimated proportion of 22 immune cell types is shown in Figure 4A and 4B. Among them, four kinds of tumor-infiltrating immune cells were negatively correlated with the OS of OC patients, including M2 macrophages (p = 0.031), activated mast cells (p = 0.0033), activated memory CD4+ T cells (p = 0.04), and neutrophils (p = 0.027), whereas M1 macrophages (p = 0.00042) were positively associated with the OS of OC patients (Figure 4C).
Figure 4. Tumor-infiltrating immune cell profile of OC samples and survival analysis. (A) Boxplot and barplot (B) display the proportion of 22 infiltrating immune cell types in OC samples. (C) Kaplan–Meier analysis of overall survival with respect to M2 macrophages, M1 macrophages, activated mast cells, and neutrophils in TCGA-OC.
Identification of TAM-related genes
The expression values of 4,043 genes were used to build a gene co-expression network using the WGCNA R package. Pearson’s correlation values and average linkage values were used to cluster the OC samples. A hierarchical clustering tree was constructed using the dynamic hybrid cutting model. Each leaf on the tree represented a single gene. Genes with similar expression patterns were clustered into branches to form a gene module. Furthermore, 10 modules were constructed (Figure 5A), among which the turquoise module exhibited a strong relationship with the black and pink modules (Figure 5B).
Figure 5. WGCNA analysis of co-expression modules. (A) Dendrogram of gene modules based on the dynamic hybrid cutting model. Ten modules were constructed. (B) Heatmap and hierarchical clustering of adjacencies in module eigengenes. (C) Heatmap of the correlation between module eigengenes and the proportion of tumor-infiltrating immune cells. (D) Scatter plot of M1 macrophage and M2 macrophage module eigengenes in three modules.
We correlated the gene module and CIBERSORT fraction and found that the black module exhibited a high correlation with M1 macrophages (R2 = 0.52, p = 2e-27) and naïve B cells (R2 = 0.5, p = 5e-25). The pink module was highly correlated with M1 macrophages (R2 = 0.54, p = 2e-30) and CD8+ T cells (R2 = 0.36, p = 1e-12), whereas the turquoise module was highly correlated with M1 (R2 = 0.47, p = 1e-22) and M2 macrophages (R2 = 0.31, p = 1e-09; Figure 5C). The scatter plots illustrated a strong positive correlation between module membership and gene significance in the black module (cor = 0.9, p = 4.9e-32), turquoise module (cor = 0.73, p = 8.9e-124), and pink module (cor = 0.79, p = 3.6e-17; Figure 5D), indicating that these modules were highly correlated with TAMs. Because we were specifically interested in macrophages, we focused on the black, pink, and turquoise modules that showed a correlation with macrophages; these were identified as TAM-related modules. According to the cut-off criteria, 219 genes in three modules were defined as tumor-associated macrophage-related gene (TAMRG) signatures (Supplementary Table 2).
Molecular subtypes based on TAMRG signature
The RNA-seq data of 375 patients with OC in the TCGA database were clustered using the unsupervised k-means-based clustering method and the expression patterns of 219 TAMRG signatures. The cumulative distribution function (CDF) plot showed k = 2 (2–6) as the optimal number of clusters (Figure 6A). The consensus heatmap divided all OC patients into two clusters: 151 (40.0%) in cluster 1 and 224 (60.0%) in cluster 2 (Figure 6B). The heatmap revealed differentially expressed genes between the two molecular subtypes (Figure 6C). The Kaplan–Meier survival analysis indicated that patients in cluster 2 had worse survival than those in cluster 1 (p = 0.0071) (Figure 6D). The violin plot showed that cluster 1 had higher M1 macrophage scores than cluster 2 (p = 6.8e-13; Figure 6E). The GSEA revealed the enriched pathways between the two groups. Some of the upregulated pathways included apoptosis, antigen processing and presentation, angiogenesis, epithelial–mesenchymal transition (EMT), and M1 macrophage upregulation in cluster 1 (Figure 6F). We compared six main immune checkpoints between the two molecule subtypes of OC patients. PDCD1 (PD1; p = 6.2e-10), CTLA4 (p = 6.6e-16), CD274 (PDL1; p = 5.4e-06), CD80 (p = 7.7e-07), PDCD1LG2 (PDL2; p = 1.9e-0), and CD86 (p = 1.6e-05) were highly expressed in cluster 1 (Figure 6G).
Figure 6. Identification of the molecular subtypes based on 219 TAMRG signatures. (A) CDF plot of the consensus score (k = 2–6). (B) Consensus clustering matrix for k = 2. (C) Heatmap showing differentially expressed genes between the two groups. (D) Kaplan–Meier analysis of overall survival for clusters 1 and 2. (E) The TAM abundance of two groups is shown in violin plots. (F) Upregulated pathways in GSEA. (G) The expression of six immune checkpoints between the two molecular subtypes.
Identification of the best TAMRG-based gene signature for predicting survival in OC patients
We identified 25 prognosis-associated TAM-based signatures using univariate Cox analysis in the TCGA training cohort (Supplementary Figure 4A). We identified a 6-gene signature using least absolute shrinkage and selection operator (LASSO) algorithm followed by multivariate Cox analysis (Supplementary Figure 4B and 4C): CD163 (hazard ratio [HR] = 1.19, 95% confidence interval [CI] = 1.05–1.34, p < 0.01), transporter 1 (TAP1, HR = 0.74, 95% CI = 0.63–0.87, p < 0.001), V-set and immunoglobulin-domain containing 4 (VSIG4, HR = 1.15, 95% CI = 1.04–1.28, p < 0.01), immunoglobulin kappa chain variable 4–1 (IGKV4-1, HR = 0.64, 95% CI = 0.47–0.86, p < 0.01), CD3E (HR = 0.82, 95% CI = 0.70–0.96, p = 0.016), and membrane spanning 4-domains A7 (MS4A7, HR = 1.16, 95% CI = 1.02–1.32, p = 0.023). We explored the expression of the 6-gene signature in the scRNA-seq set (Supplementary Figure 5A; MS4A7, VSIG4, and IGKV4-1 were not provided). CD163, MS4A7, and VSIG4 were upregulated in macrophages, whereas CD3E was significantly upregulated in T cells. TAP1 was upregulated in macrophages and B cells. The expression of IGKV4-1 was not detected. Furthermore, Gene Expression Profiling Interactive Analysis (GEPIA) was performed to study the expression of the 6-gene signature in 88 normal (GTEx) samples and 426 OC (TCGA) samples. We noticed that compared with normal samples, CD163 was upregulated, whereas TAP1, CD3E, and IGKV4-1 were downregulated in OC (Supplementary Figure 5B).
Generation and validation of the 6-gene signature-based prognostic risk score model
The risk score was calculated using the following formula: risk score = 0.0284 × ExpCD163 + (–0.0297) × ExpCD3E + (–0.0007) × ExpIGKV4 + (–0.0130) × ExpTAP1 + 0.0009 × ExpVSIG4 + 0.0487 × ExpMS4A7. Patients were dichotomized into low-risk or high-risk groups according to the median risk score (Figure 7B). The Kaplan–Meier survival curve revealed that the OS of the low-risk group was higher than that of the high-risk group (log-rank, p = 0.00016; Figure 7A and 7B). The concordance index (C-index) for OS prediction was 0.614 (95% CI = 0.593–0.636). The receiver operating characteristic (ROC) curve demonstrated that areas under the curve (AUC) of 0.624, 0.68, and 0.718 revealed a predictive ability of 3-, 5- and 10-year OS, respectively (Figure 7C). The relationship between the proportion of six infiltrating immune cell types and the risk score was analyzed to validate the effect of the 6-gene signature on TAMRG. The scatter plot showed a highly negative correlation between the risk score and the proportion of infiltrating M1 macrophages (R = 0.38, p = 1.4e-14). The proportion of infiltrating M2 macrophages was positively correlated with the risk score (R = 0.29, p = 1.3e-08; Supplementary Figure 5C).
Figure 7. Prognostic analysis of the 6-gene signature in OC patients. Kaplan–Meier OS curve of low-risk and high-risk groups in the TCGA training dataset (A) and GEO validation dataset (D). Risk score distribution in TCGA (B) and GEO (E) datasets. Upper panel: The curve of risk score. Middle panel: patients’ overall survival status and time. Bottom panel: Heatmaps of gene expression profiles. ROC curve for OS prediction in the TCGA (C) and GEO (F) datasets.
Next, the GEO cohort was used to validate the prognostic predictive performance. As shown in Figure 7D and 7E, 185 OC patients were divided into low-risk and high-risk groups. Consistent with the TCGA training cohort, the Kaplan–Meier survival curve suggested that the OS of the high-risk group was considerably lower than that of the low-risk group (log-rank, p = 0.0026;). The C-index was 0.611 (95% CI = 0.584–0.638). Moreover, the ROC suggested AUC values of 0.674 and 0.704, indicating that the model could predict 3- and 5-year OS, respectively (Figure 7F).
Comparison with clinical characteristics and other signatures
We next compared the predictive accuracy of the 6-gene signature with clinical characteristics and published prognostic signatures of OC. Among the 13 survival predictors, TAMRG-based signature exhibited the best mean C-index (0.613) for age (0.593), grade (0.532), stage (0.519), residual (0.556) (Supplementary Table 3), and other signatures (0.516–0.584; Supplementary Table 4). These results revealed that the TAMRG-based signature was an independent prognostic indicator that effectively predicted the prognosis of patients with OC.
Solid tumors are characterized by a unique microenvironment formed by malignant and several non-malignant cells that can modify tumor characteristics . We analyzed the heterogeneity of tumor-infiltrating immune cells in OC using scRNA-seq and bulk RNA-seq data. At the single-cell level, the majority of non-malignant cells were immune cells, identified as four distinct clusters of B cells, T cells, DCs, and macrophages. Both inter-and intra-patient heterogeneity was observed for tumor-infiltrating immune cells in HGSOC. Macrophages and T cells exhibited immunosuppressive characteristics: macrophages with an M2 phenotype and T cells with an exhausted phenotype. M2 macrophages express the ligand receptors for PD-1 and CTLA-4, whose activation inhibits T cell proliferation and cytotoxic function, thus contributing to the formation of an immunosuppressive TME [14, 15]. These findings suggested that a distinct immune system status and dynamic immune cell interactions in the surrounding microenvironment contributed to the wide tumor heterogeneity of ovarian cancer transcriptome.
The bulk RNA-seq data were analyzed using CIBERSORT to estimate tumor-infiltrating immune cell subsets and their correlation with prognosis. Activated mast cells, neutrophils, M2 macrophages, and activated memory CD4+ T cells were negatively correlated with OS of OC patients, whereas M1 macrophage infiltration indicated better clinical outcomes. Similarly, CD4+CD25+FOXP3+ Treg cell infiltration was associated with high mortality and decreased survival in 104 individuals affected with OC . Although the increased density of M2-like TAM is known to be associated with poor OS [17–19], the relationship between M2 density and OS in OC remains controversial. For instance, Zhang et al.  reported no correlation between these two factors, whereas Lan et al.  found a negative correlation, consistent with our results. These differences in the findings could be ascribed to varying types of tumor tissues.
A total of 219 TAMRG signatures were established by WGCNA. OC patients with M1-like TAM reported a better prognosis, suggesting that TAMRG-based patient classification could effectively predict patient survival. Although GSEA revealed that M1-like TAM was enriched in immune-related, angiogenesis, EMT, and JAK-STAT3 pathways, M2 macrophages constituted the dominant population in the TAM of OC and were implicated in tumor angiogenesis, invasion, metastasis, and early recurrence [12, 22, 23]. The simple dichotomy of M1/M2 macrophages is insufficient to explain the complexity of TAM heterogeneity . Classically, M1/M2 phenotypes are extremes of a continuum of activation states [24, 25], whereas the TAM subset shares characteristics of both M1 and M2 phenotypes . For example, a recent study revealed that TAMs simultaneously express M1/M2 markers, with early-stage TAMs co-expressing T cell co-stimulatory and co-inhibitory receptors . Similarly, Müller identified phenotypic differences in TAMs from distinct lineages at single-cell resolution in human gliomas. These results indicate that TAMs frequently co-express M1/M2 markers in single cells, making it difficult to isolate M1 and M2 phenotypes . Thus, TAMs exhibit functional plasticity and intermediate states, resulting in reversible changes in their distribution and functional states under different microenvironment stimuli [24, 29–31].
Immunotherapy has emerged as a promising treatment strategy for cancer that alleviates the immunosuppression status of the tumor. For instance, an anti-PD-1 antibody is known to prolong the progression-free survival and OS of patients with melanoma  and non-small cell lung cancer . Unfortunately, clinical trials involving immune checkpoint inhibitors either as a single agent or in combination with other therapeutic modalities have been unsuccessful in OC. Our results revealed that cluster 1 with an increased macrophage M1 density expressed high levels of PD1/PDL1/PDL2 and CTLA4/CD80/CD86 molecules, indicating that patients in cluster 1 may benefit more from anti-PD1 and anti-CTLA4 therapies than those in cluster 1.
A novel 6-gene signature risk score module was successfully established and validated using a GEO dataset. Of the six genes, a high expression of CD163, VSIG4, and MS4A7 was related to poor prognosis, whereas that of CD3E, IGKV4, and TAP1 was associated with a favorable prognosis. Some of these genes have been implicated in cancers. For example, CD163(+) TAMs correlate with poor prognosis, OS, and metastasis of various malignancies. Chen et al. reported that CD163 contributes to gliomagenesis via CK2, and its high expression is associated with an unfavorable patient prognosis . Similarly, CD163(+) TAMs were associated with poor OS and increased microvessel density in gastric cancer . VSIG4 is overexpressed in OC . Agnes et al. demonstrated that underexpressed TAP1 resulted in low infiltration of macrophages and poor prognosis in patients with colorectal cancer . Univariate and multivariate Cox regression analyses revealed that the 6-gene signature could be applied as an independent prognostic factor. We believe this prognostic module is the first to incorporate a TAM-related signature to predict survival in patients with OC. Although the absolute value of the prognostic module C-index was low, it was superior to the traditional clinical characteristics. Therefore, it can be used to estimate the prognosis of patients and classify them into distinct subgroups for effective treatment.
Our study had certain limitations. The data of a few patients acquired from TCGA and GEO databases were incomplete in terms of grade and medical history (e.g., unavailable for multivariate Cox analysis). In addition, this was a retrospective study; further prospective, large-scale trials are warranted to verify its clinical application.
We constructed a novel TAMRG-based, the 6-gene prognostic signature for patients with OC. The TAMRG-based signature could serve as a potential target for the prognosis and predicting the therapeutic response in patients with OC.
Materials and Methods
Dataset collection and preprocessing
We used both the bulk RNA-seq data and scRNA-seq data of human OC samples. The scRNA-seq data (GSE146026) consisted of 9,609 cells from six ascites samples of HGSOC and were acquired from the GEO database. The single-cell library was constructed on a 10× genomics platform and read on an Illumina NextSeq 500 sequencing system. The bulk RNA-seq data of patients with OC were downloaded from the GEO and TCGA databases. We included 375 OC samples with available clinical characteristics from the TCGA database as the training dataset and 185 OC samples from the GEO database (GSE26712) as the validation dataset.
Single-cell RNA-seq analysis
The R4.0.5 software was used to perform all analyses. The Seurat 3.0 package was used for scRNA-seq data quality control, filtering, statistical analysis, and subsequent analysis . First, cells with < 200 detected genes, genes detected in < 3 cells, and mitochondrial genes ≥ 5% were used as the filtering criteria. In total, 9,609 cells were included in the study. Next, a linear regression model was used to normalize the gene expression. Data were analyzed using PCA to visualize the available dimensions (p-value < 0.05). The first 20 PCs were used for tSNE with a resolution of 0.3 for dimension reduction and clustering. The differentially expressed genes were identified using the Wilcoxon test and FindAllMarkers function of Seurat. Marker genes were identified using the following cut-off criteria: |log2[fold change (FC)]| > 0.5 and p-value < 0.05. Afterward, three cell types were annotated using the established cell surface marker genes . To further characterize the immune cells, we annotated the cells using the singleR package and corrected with the CellMarker database. The list of cell surface markers is shown in Supplementary Table 1.
Gene set variation analysis of immune cells
GSVA is an unsupervised gene set enrichment method to estimate the variations in pathway activity within a sample population . The enrichment scores of the gene sets were evaluated using the GSVA package in R. Immune cell-related gene sets were derived from supplementary documents provided by Chung et al. .
Estimate of the abundance of tumor immune infiltration and survival analysis
The CIBERSORT R script was used to calculate the abundance of 22 tumor-infiltrating immune cell types in the TCGA OC cohort based on the bulk RNA-seq dataset . The median of each immune cell enrichment score was used to separate the samples into two groups to compare the OS. The R package “survival” and “survminer” were used for computing the survival.
Construction of weighted correlation network analysis
WGCNA analysis was implemented using the R package “WGCNA” . To construct a signed, scale-free co-expression gene network, the scale-free topology fitting index R2 > 0.85 and power of β = 3 were selected as soft-threshold parameters. We used the dynamic hybrid cutting method to classify genes with similar expression patterns; the minimum size cut-off of the module was 30. Module eigengenes were used to perform the component analysis of each module. We estimated the correlation between tumor-infiltrating immune cell enrichment score and module eigengenes to determine the significance of modules using Pearson’s test. The selected macrophage subtypes and modules were used for follow-up analysis. The major genes were determined by calculating the gene significance (GS) and module membership (MM). The genes in the module with |GS| > 0.2 and MM > 0.8 were considered significant.
Identification of molecular subtypes
The R package “ConsensusClusterPlus” was used to perform k-means-based unsupervised consensus clustering based on the expression patterns of TAM-related gene signature . We obtained consensus matrices and cumulative distribution function (CDF) plots with a set of parameters, including 1,000 iterations and an 80% resampling rate in Pearson’s correlation. Next, we compared OS between different clusters using the Kaplan–Meier survival analysis. GSEA was used to explore the potential molecular mechanisms. In addition, we compared six immune checkpoints (PDCD1, CTLA4, CD274, CD80, PDCD1LG2, and CD86) between different clusters.
Establishment and estimation of the prognostic risk score model
First, the univariate Cox regression analysis was used to assess the associations between TAM-related genes and survival in the TCGA training set. Second, the LASSO algorithm and multivariate Cox regression analysis were applied to identify prognosis-related genes with p < 0.05. Subsequently, a prognostic risk model was established based on the major prognosis-related TAMRG-based genes. The risk score was calculated by a linear method to assemble the Cox coefficient and prognostic gene expression using the following formula: Risk score = β1 × Expgene1 + β2 × Expgene2 + … + βn × Expgenen, where “β” and “Exp” represent the regression coefficient and the expression of prognostic genes, respectively.
The patients in the TCGA training set were divided into high-risk and low-risk groups according to the median risk score. Kaplan–Meier survival analysis was used to compare the OS of these two groups. Time-dependent ROC curve analysis and Harrell’s C-index were applied using the “survivalROC” and “survminer” packages in R to evaluate the prediction accuracy of the prognostic risk model. Finally, a validation cohort was obtained from the GEO database to validate the prognostic value of the risk model. In addition, we performed a comparison of the C-index between the prognostic model and age, stage, grade, residual, and eight published signatures.
Student’s t-test was used to compare the mean values between the two groups. Cox proportional hazard models were applied to assess the associations between factors and OS. Survival curves were drawn using the Kaplan–Meier method and compared using log-rank tests. The C-index and the ROC curve were estimated using survival, survminer, and survivalROC packages in the R software.
Statistical analyses were performed using the R version 4.0.5 (The R Foundation). The p-values were two-tailed, and p < 0.05 was considered significant.
QT and FD initiated the comprehensive study; QT, HL, and FD wrote the main manuscript; and JX and YM analyzed the data and prepared all the figures. All authors reviewed the manuscript.
We are grateful to the TCGA and GEO databases for providing the data sources.
Conflicts of Interest
The authors declare no conflicts of interest related to this study.
This research was supported by the Nature Scientific Foundation of China (grant number: 81472434).
- 1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021; 71:7–33. https://doi.org/10.3322/caac.21654 [PubMed]
- 2. Kuroki L, Guntupalli SR. Treatment of epithelial ovarian cancer. BMJ. 2020; 371:m3773. https://doi.org/10.1136/bmj.m3773 [PubMed]
Bowtell DD, Böhm S, Ahmed AA, Aspuria PJ, Bast RC
Jr, Beral V, Berek JS, Birrer MJ, Blagden S, Bookman MA, Brenton JD, Chiappinelli KB, Martins FC, et al. Rethinking ovarian cancer II: reducing mortality from high-grade serous ovarian cancer. Nat Rev Cancer. 2015; 15:668–79. https://doi.org/10.1038/nrc4019 [PubMed]
- 4. Luo Z, Wang Q, Lau WB, Lau B, Xu L, Zhao L, Yang H, Feng M, Xuan Y, Yang Y, Lei L, Wang C, Yi T, et al. Tumor microenvironment: The culprit for ovarian cancer metastasis? Cancer Lett. 2016; 377:174–82. https://doi.org/10.1016/j.canlet.2016.04.038 [PubMed]
- 5. Nowak M, Klink M. The Role of Tumor-Associated Macrophages in the Progression and Chemoresistance of Ovarian Cancer. Cells. 2020; 9:1229. https://doi.org/10.3390/cells9051299 [PubMed]
- 6. Jiang Y, Wang C, Zhou S. Targeting tumor microenvironment in ovarian cancer: Premise and promise. Biochim Biophys Acta Rev Cancer. 2020; 1873:188361. https://doi.org/10.1016/j.bbcan.2020.188361 [PubMed]
- 7. Yang Y, Yang Y, Yang J, Zhao X, Wei X. Tumor Microenvironment in Ovarian Cancer: Function and Therapeutic Strategy. Front Cell Dev Biol. 2020; 8:758. https://doi.org/10.3389/fcell.2020.00758 [PubMed]
- 8. Jiménez-Sánchez A, Memon D, Pourpe S, Veeraraghavan H, Li Y, Vargas HA, Gill MB, Park KJ, Zivanovic O, Konner J, Ricca J, Zamarin D, Walther T, et al. Heterogeneous Tumor-Immune Microenvironments among Differentially Growing Metastases in an Ovarian Cancer Patient. Cell. 2017; 170:927–38.e20. https://doi.org/10.1016/j.cell.2017.07.025 [PubMed]
- 9. Jiménez-Sánchez A, Cybulska P, Mager KL, Koplev S, Cast O, Couturier DL, Memon D, Selenica P, Nikolovski I, Mazaheri Y, Bykov Y, Geyer FC, Macintyre G, et al. Unraveling tumor-immune heterogeneity in advanced ovarian cancer uncovers immunogenic effect of chemotherapy. Nat Genet. 2020; 52:582–93. https://doi.org/10.1038/s41588-020-0630-5 [PubMed]
- 10. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011; 144:646–74. https://doi.org/10.1016/j.cell.2011.02.013 [PubMed]
- 11. Condeelis J, Pollard JW. Macrophages: obligate partners for tumor cell migration, invasion, and metastasis. Cell. 2006; 124:263–66. https://doi.org/10.1016/j.cell.2006.01.007 [PubMed]
- 12. Pollard JW. Tumour-educated macrophages promote tumour progression and metastasis. Nat Rev Cancer. 2004; 4:71–78. https://doi.org/10.1038/nrc1256 [PubMed]
- 13. Sica A, Allavena P, Mantovani A. Cancer related inflammation: the macrophage connection. Cancer Lett. 2008; 267:204–15. https://doi.org/10.1016/j.canlet.2008.03.028 [PubMed]
- 14. Yang Q, Guo N, Zhou Y, Chen J, Wei Q, Han M. The role of tumor-associated macrophages (TAMs) in tumor progression and relevant advance in targeted therapy. Acta Pharm Sin B. 2020; 10:2156–70. https://doi.org/10.1016/j.apsb.2020.04.004 [PubMed]
- 15. Noy R, Pollard JW. Tumor-associated macrophages: from mechanisms to therapy. Immunity. 2014; 41:49–61. https://doi.org/10.1016/j.immuni.2014.06.010 [PubMed]
- 16. Curiel TJ, Coukos G, Zou L, Alvarez X, Cheng P, Mottram P, Evdemon-Hogan M, Conejo-Garcia JR, Zhang L, Burow M, Zhu Y, Wei S, Kryczek I, et al. Specific recruitment of regulatory T cells in ovarian carcinoma fosters immune privilege and predicts reduced survival. Nat Med. 2004; 10:942–49. https://doi.org/10.1038/nm1093 [PubMed]
- 17. Väyrynen JP, Haruki K, Lau MC, Väyrynen SA, Zhong R, Dias Costa A, Borowsky J, Zhao M, Fujiyoshi K, Arima K, Twombly TS, Kishikawa J, Gu S, et al. The Prognostic Role of Macrophage Polarization in the Colorectal Cancer Microenvironment. Cancer Immunol Res. 2021; 9:8–19. https://doi.org/10.1158/2326-6066.CIR-20-0527 [PubMed]
- 18. Cersosimo F, Lonardi S, Bernardini G, Telfer B, Mandelli GE, Santucci A, Vermi W, Giurisato E. Tumor-Associated Macrophages in Osteosarcoma: From Mechanisms to Therapy. Int J Mol Sci. 2020; 21:5207. https://doi.org/10.3390/ijms21155207 [PubMed]
- 19. Zheng X, Weigert A, Reu S, Guenther S, Mansouri S, Bassaly B, Gattenlöhner S, Grimminger F, Pullamsetti S, Seeger W, Winter H, Savai R. Spatial Density and Distribution of Tumor-Associated Macrophages Predict Survival in Non-Small Cell Lung Carcinoma. Cancer Res. 2020; 80:4414–25. https://doi.org/10.1158/0008-5472.CAN-20-0069 [PubMed]
- 20. Zhang M, He Y, Sun X, Li Q, Wang W, Zhao A, Di W. A high M1/M2 ratio of tumor-associated macrophages is associated with extended survival in ovarian cancer patients. J Ovarian Res. 2014; 7:19. https://doi.org/10.1186/1757-2215-7-19 [PubMed]
- 21. Lan C, Huang X, Lin S, Huang H, Cai Q, Wan T, Lu J, Liu J. Expression of M2-polarized macrophages is associated with poor prognosis for advanced epithelial ovarian cancer. Technol Cancer Res Treat. 2013; 12:259–67. https://doi.org/10.7785/tcrt.2012.500312 [PubMed]
- 22. Reinartz S, Schumann T, Finkernagel F, Wortmann A, Jansen JM, Meissner W, Krause M, Schwörer AM, Wagner U, Müller-Brüsselbach S, Müller R. Mixed-polarization phenotype of ascites-associated macrophages in human ovarian carcinoma: correlation of CD163 expression, cytokine levels and early relapse. Int J Cancer. 2014; 134:32–42. https://doi.org/10.1002/ijc.28335 [PubMed]
- 23. Yin M, Li X, Tan S, Zhou HJ, Ji W, Bellone S, Xu X, Zhang H, Santin AD, Lou G, Min W. Tumor-associated macrophages drive spheroid formation during early transcoelomic metastasis of ovarian cancer. J Clin Invest. 2016; 126:4157–73. https://doi.org/10.1172/JCI87252 [PubMed]
- 24. Ostuni R, Kratochvill F, Murray PJ, Natoli G. Macrophages and cancer: from mechanisms to therapeutic implications. Trends Immunol. 2015; 36:229–39. https://doi.org/10.1016/j.it.2015.02.004 [PubMed]
- 25. Mantovani A, Sozzani S, Locati M, Allavena P, Sica A. Macrophage polarization: tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol. 2002; 23:549–55. https://doi.org/10.1016/s1471-4906(02)02302-5 [PubMed]
- 26. Qian BZ, Pollard JW. Macrophage diversity enhances tumor progression and metastasis. Cell. 2010; 141:39–51. https://doi.org/10.1016/j.cell.2010.03.014 [PubMed]
- 27. Singhal S, Stadanlick J, Annunziata MJ, Rao AS, Bhojnagarwala PS, O'Brien S, Moon EK, Cantu E, Danet-Desnoyers G, Ra HJ, Litzky L, Akimova T, Beier UH, et al. Human tumor-associated monocytes/macrophages and their regulation of T cell responses in early-stage lung cancer. Sci Transl Med. 2019; 11:eaat1500. https://doi.org/10.1126/scitranslmed.aat1500 [PubMed]
- 28. Müller S, Kohanbash G, Liu SJ, Alvarado B, Carrera D, Bhaduri A, Watchmaker PB, Yagnik G, Di Lullo E, Malatesta M, Amankulor NM, Kriegstein AR, Lim DA, et al. Single-cell profiling of human gliomas reveals macrophage ontogeny as a basis for regional differences in macrophage activation in the tumor microenvironment. Genome Biol. 2017; 18:234. https://doi.org/10.1186/s13059-017-1362-4 [PubMed]
- 29. Stout RD, Jiang C, Matta B, Tietzel I, Watkins SK, Suttles J. Macrophages sequentially change their functional phenotype in response to changes in microenvironmental influences. J Immunol. 2005; 175:342–49. https://doi.org/10.4049/jimmunol.175.1.342 [PubMed]
- 30. Okabe Y, Medzhitov R. Tissue-specific signals control reversible program of localization and functional polarization of macrophages. Cell. 2014; 157:832–44. https://doi.org/10.1016/j.cell.2014.04.016 [PubMed]
- 31. Kim J, Bae JS. Tumor-Associated Macrophages and Neutrophils in Tumor Microenvironment. Mediators Inflamm. 2016; 2016:6058147. https://doi.org/10.1155/2016/6058147 [PubMed]
- 32. Robert C, Schachter J, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, Larkin J, Lorigan P, Neyns B, et al, and KEYNOTE-006 investigators. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N Engl J Med. 2015; 372:2521–32. https://doi.org/10.1056/NEJMoa1503093 [PubMed]
- 33. Reck M, Rodríguez-Abreu D, Robinson AG, Hui R, Csőszi T, Fülöp A, Gottfried M, Peled N, Tafreshi A, Cuffe S, O'Brien M, Rao S, Hotta K, et al, and KEYNOTE-024 Investigators. Pembrolizumab versus Chemotherapy for PD-L1-Positive Non-Small-Cell Lung Cancer. N Engl J Med. 2016; 375:1823–33. https://doi.org/10.1056/NEJMoa1606774 [PubMed]
- 34. Chen T, Chen J, Zhu Y, Li Y, Wang Y, Chen H, Wang J, Li X, Liu Y, Li B, Sun X, Ke Y. CD163, a novel therapeutic target, regulates the proliferation and stemness of glioma cells via casein kinase 2. Oncogene. 2019; 38:1183–99. https://doi.org/10.1038/s41388-018-0515-6 [PubMed]
- 35. Park JY, Sung JY, Lee J, Park YK, Kim YW, Kim GY, Won KY, Lim SJ. Polarized CD163+ tumor-associated macrophages are associated with increased angiogenesis and CXCL12 expression in gastric cancer. Clin Res Hepatol Gastroenterol. 2016; 40:357–65. https://doi.org/10.1016/j.clinre.2015.09.005 [PubMed]
- 36. Byun JM, Jeong DH, Choi IH, Lee DS, Kang MS, Jung KO, Jeon YK, Kim YN, Jung EJ, Lee KB, Sung MS, Kim KT. The Significance of VSIG4 Expression in Ovarian Cancer. Int J Gynecol Cancer. 2017; 27:872–78. https://doi.org/10.1097/IGC.0000000000000979 [PubMed]
- 37. Ling A, Löfgren-Burström A, Larsson P, Li X, Wikberg ML, Öberg Å, Stenling R, Edin S, Palmqvist R. TAP1 down-regulation elicits immune escape and poor prognosis in colorectal cancer. Oncoimmunology. 2017; 6:e1356143. https://doi.org/10.1080/2162402x.2017.1356143 [PubMed]
- 38. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018; 36:411–20. https://doi.org/10.1038/nbt.4096 [PubMed]
- 39. Maynard A, McCoach CE, Rotow JK, Harris L, Haderk F, Kerr DL, Yu EA, Schenk EL, Tan W, Zee A, Tan M, Gui P, Lea T, et al. Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing. Cell. 2020; 182:1232–51.e22. https://doi.org/10.1016/j.cell.2020.07.017 [PubMed]
- 40. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
- 41. Chung W, Eum HH, Lee HO, Lee KM, Lee HB, Kim KT, Ryu HS, Kim S, Lee JE, Park YH, Kan Z, Han W, Park WY. Single-cell RNA-seq enables comprehensive tumour and immune cell profiling in primary breast cancer. Nat Commun. 2017; 8:15081. https://doi.org/10.1038/ncomms15081 [PubMed]
- 42. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
- 43. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9:559. https://doi.org/10.1186/1471-2105-9-559 [PubMed]
- 44. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–73. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]