Analysis of single-cell RNA-sequencing data identifies a hypoxic tumor subpopulation associated with poor prognosis in triple-negative breast cancer

1.   Introduction

Breast cancer (BRCA) is the most frequently diagnosed cancer and the leading cause of cancer mortality among females ^[1]. Triple-negative breast cancer (TNBC) is a subtype of BRCA with the absence of estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2), accounting for approximately 15–20% of all breast carcinomas ^[2]. Compared to other subtypes, TNBC has a more aggressive clinical process and worse prognosis due to the extremely heterogeneous tumor microenvironment (TME) ^[3]. TME is a complicated ecological system comprising vasculature, extracellular matrix (ECM), cytokines, growth factors, and lots of various cell populations ^[4], which has a strong association with tumorigenesis, tumor cells proliferation, invasion, and metastasis of BRCA ^[5]. Investigating the mechanisms by which TME influences tumor progression in TNBC could facilitate improving the prognosis of patients.

Single-cell RNA-sequencing (scRNA-seq) is an innovative technology that enables the investigation of the transcriptome of individual cells with the purpose of identification of cell subpopulations with analogical transcriptional patterns ^[6]. In contrast with traditional bulk RNA-sequencing (RNA-seq), scRNA-seq allows revealing the significant heterogeneity within TME ^[7,8]. Previous study has discovered three distinct subpopulations of breast epithelial cells and reconstructed a continuous differentiation trajectory ^[9], which provided comprehensive insights into the human mammary epithelium. Similarly, cancer-associated fibroblasts (CAFs) were revealed to exist in three spatially and functionally diverse subpopulations in the mouse breast cancer model ^[10]. Furthermore, the prognostic capacity of each CAFs subpopulation was validated in clinical cohorts, confirming the potential prognostic implications of cell subpopulations within TME.

Because of the imbalance between the rapid proliferation of tumor cells and the inadequate oxygen supply of vascular tissue, hypoxic TME is an inherent characteristic of solid tumors ^[11]. In the condition of hypoxia, tumor cells could secrete a variety of vascular growth factors that promote abnormal angiogenesis, consequently increasing the oxygen supply to meet tumor progression. Meanwhile, the invasive and metastatic capacities of tumor cells might also be greatly enhanced for the purpose of seeking a more appropriate growing environment ^[12,13]. In recent studies, prognostic prediction models have been established and validated based on hypoxia-related genes (HRGs) ^[14] and hypoxia-immune-related genes ^[15], expounding their views that HRGs within hypoxic TME were significant factors contributing to the poor prognosis of TNBC patients. However, the heterogeneity within hypoxic TME was ignored in these studies, which may play a vital role during tumor progression. To overcome this deficiency, we explored the heterogeneity within hypoxic TME utilizing scRNA-seq data and corresponding analytical methods, such as Gene Set Enrichment Analysis (GSEA) and cell-cell communication analysis in this study. Furthermore, we identified a hypoxia-related subpopulation and investigated its involved signaling pathways and prognostic implications for TNBC patients.

2.   Materials and methods

2.1.   Data acquisition

The counts matrix of TNBC scRNA-seq was collected from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) with accession number GSM4476488 ^[16]. The BRCA RNA-seq expression profiles of The Cancer Genome Atlas (TCGA) were downloaded through the UCSC Xena platform (https://xena.ucsc.edu/) ^[17], consisting of tumor samples of 1050 patients with corresponding clinical features. Subtyping by R package genefu ^[18] with the PAM50 model obtained 188 TNBC patients. Furthermore, we retrieved the data of 107 TNBC patients from GEO with accession number GSE58812 ^[19].

2.2.   Analysis of scRNA-seq data

Seurat package ^[20] was used to analyze the raw counts matrix of scRNA-seq data containing 33,694 genes and 532 cells. Imported into R by CreateSeuratObject with genes detected in at least one cell, cells with the proportion of mitochondrial genes higher than 15% and doublet cells detected by DoubletFinder package were further filtered. After data normalization through NormalizeData, the top 2000 highly variable genes were identified using FindVariableFeatures for downstream analysis. Subsequently, we performed ScaleData to scale the data and remove the effects of certain variables such as the proportion of mitochondrial genes. Furthermore, we conducted principal component analysis (PCA) to reduce the dimensionality. The first 15 significant principal components were applied by t-Distributed Stochastic Neighbor Embedding (t-SNE) for nonlinear dimensionality reduction and FindNeighbors for cells clustering. Finally, the resolution of FindClusters was set to 0.8 thus identifying seven clusters. Cell types were annotated by typical marker genes such as EPCAM, VWF, PTPRC and PDGFRA along with the SingleR package ^[21]. The genomic copy number profiles of cells were inferred to distinguish tumor cells from non-malignant (i.e., normal) cell types through the CopyKAT package with default parameters ^[16]. For the epithelial cells extracted, the resolution was adjusted to 1.3 after re-clustering. FindMarkers was utilized to identify marker genes with the log₂FC (fold change) > 0.2 and p-value < 0.05. Several plotting functions built-in the Seurat package were applied for visualization such as DimPlot, VlnPlot, and FeatureScatter.

2.3.   Single-sample Gene Set Enrichment Analysis

Single-sample Gene Set Enrichment Analysis (ssGSEA) is an extension of GSEA that calculates separate enrichment scores for each paired sample and gene set. Each ssGSEA enrichment score represents the extent to which genes in a particular gene set are coordinately up- or down-regulated within single sample. Acquiring the collection of hallmark gene sets from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) ^[22] as an input, which includes 50 gene sets representing specific well-defined biological states or processes, we performed ssGSEA and visualized the enrichment scores using R package escape.

2.4.   Cell-cell communication analysis

Analysis of cell-cell communication was conducted based on R package CellChat ^[23] built-in human ligand-receptor interactions database, consisting of 1939 interactions among 427 ligands and 357 receptors. Through the integration of the expression levels with prior known ligand-receptor interactions, we compute the communication probability and infer intercellular communication networks. Afterwards, we summarized the communication probabilities of all interactions associated with each signaling pathway, and identified outputting communication patterns primarily regulating different signaling pathways.

2.5.   Tumor immune estimation resource (TIMER) analysis tool

To determine the correlation between gene expression and immune infiltration in TNBC, we applied TIMER 2.0 (http://timer.cistrome.org/) web platform, comprising of more 10,000 TCGA samples across 32 cancer types. We investigated the infiltration levels of six types of immune cells, including CD8+ T cells, CD4+ T cells, B cells, neutrophils, macrophages, and dendritic cells. Furthermore, we utilized the Gene_DE and Gene_Corr modules of TIMER to explore the expression of ARTN and L1CAM in TNBC, along with their correlation.

2.6.   Construction of the risk score model for TNBC prognosis

Utilizing the TNBC patients extracted from TCGA-BRCA cohorts as the train set, univariate Cox proportional hazards regression analysis was performed to identify prognostic genes (p-value < 0.05) associated with overall survival (OS). Whereafter, multivariate Cox proportional hazards regression analysis was implemented to construct a prognostic model for the calculation of risk score as Eq (1), where the β_i and Exp_i represented the multivariate Cox regression coefficient and expression level of gene i respectively. The number of prognostic genes was denoted by n. To ascertain how well the risk score model predicts the survival of patients at the time point of one, three and five years, receiver operating characteristic (ROC) curves were drawn by survivalROC package to obtain the area under curve (AUC) values. Afterward, TNBC patients of the GSE58812 dataset were divided into high- and low-risk groups according to the median risk score as the cutoff value. The performance of the model was further validated through Kaplan-Meier curves with the log-rank test.

3.   Results and discussion

3.1.   Investigation of the heterogeneity within hypoxic TME

A total of 483 cells and 22,237 genes were remained after quality control, consisting of 196 epithelial cells, 250 fibroblasts, 13 endothelial cells, and 24 monocytes. On the basis of the genomic copy number profiles distinguishing cellular malignant status, we divided the epithelial cells into the tumor and normal ones. Cell numbers of the former (157) were nearly four times that of the latter (39), suggesting that a large proportion of the epithelial cells were malignant. The final annotation results of cell types were shown in Figure 1A, while the expression levels of typical marker genes were illustrated in Figure A1. To further explore the heterogeneity within epithelial cells, we extracted them and identified two normal subpopulations along with four tumor subpopulations (Figure 1B). For each subpopulation, marker genes were arranged in descending order according to log₂FC, of which the top five genes were illustrated in Figure 1C.

Figure 1.  Dimensional reduction plot by t-SNE for (A) all cell types after quality control; (B) subpopulations of epithelial cells. (C) Heat map of the top five genes ordered by log₂FC among all subpopulations of epithelial cells.

DownLoad: Full-Size Img PowerPoint

Notably, the co-expression of homologous genes was detected in some subpopulations of epithelial cells. TFF1 and TFF3, both belonging to the same trefoil factor family, were overexpressed in the normal subpopulation 1 of epithelial cells (epi_nor1). It has been demonstrated that TFF1 was strongly expressed under estrogen transcriptional control in a type of human breast cancer cells containing estrogen receptor, but at low levels in normal mammary tissue ^[24,25]. In line with previous studies, the expression level of estrogen receptor 1 (ESR1) was up-regulated in epi_nor1 compared to other subpopulations (Figure A1). Consistently, ssGSEA results further confirmed that epi_nor1 was intensely involved in the response to estrogen (Figure 2A). Furthermore, the Pearson correlation coefficient among TFF1, TFF3, and ESR1 in epithelial cells reached 0.97, 0.85, and 0.86 respectively (Figure 2B), which indicates that the expression of TFF3 might be under estrogen regulation similar to that of TFF1. Moreover, S100 calcium binding protein A8 (S100A8) and S100 calcium binding protein A9 (S100A9) were up-regulated in epi_tum2, which can form a heterodimer with biological activity ^[26]. The expression of S100A8 and S100A9 was regulated by tamoxifen, an estrogen receptor inhibitor, in breast cancer tissue ^[27]. Meanwhile, mammaglobin A (SCGB2A2) and lipophilin B (SCGB1D2), members of the secretoglobin superfamily, were identified as the marker genes of epi_tum1. These two genes were reported to be co-expressed in breast cancer ^[28] and several malignant tumors of the female genital tract ^[29], where their proteins form a covalent compound ^[30,31]. It has been also proved that SCGB1D2 was expressed in estrogen receptor-positive tumors with more frequency ^[28]. These results suggest that the co-expression of above mentioned gene pairs could be regulated by ESR1-positive epi_nor1.

Figure 2.  Scatter plot showing (A) the ssGSEA enrichment score of early and late response to estrogen in each subpopulation of epithelial cells (Each hexagon represents a certain number of cells, indicated by color shades); (B) the co-expression among TFF1, TFF3 and ESR1 in epi_nor1, with the Pearson correlation coefficient reached 0.97, 0.85 and 0.86 respectively.

DownLoad: Full-Size Img PowerPoint

As the ssGSEA results characterized diverse biological states for subpopulations of epithelial cells shown in Figure 3, there were no hallmark gene sets enriched in epi_nor2, indicating its comparatively weak biological behaviors contrary to other subpopulations. Moreover, gene sets associated with cell division (G2/M checkpoint, DNA repair, and mitotic spindle assembly) as well as target genes of MYC and E2F transcription factors family were significantly enriched in epi_tum3 and epi_tum4. It has been proved that the transcription factors coded by MYC ^[32,33,34] and E2F family ^[35] target genes related to DNA synthesis and cell cycle, thus regulating cell division. Although the enrichment scores of epi_tum2 on these gene sets are slightly less than epi_tum3 and epi_tum4, there appear to be similarities among them in some gene sets such as DNA repair and target genes of MYC. These results suggest that these three subpopulations might be in the process of mitosis and proliferation, of which epi_tum3 and epi_tum4 share similar malignant states. It was further confirmed that genes encoding proteins involved in oxidative phosphorylation were enriched significantly in epi_tum2, epi_tum3, and epi_tum4 due to excess energy consumption by cell proliferation. In contrast to them, up-regulation of genes associated with hypoxia in epi_tum1 indicates it was a vital subpopulation within hypoxic TME. To clarify how the hypoxic tumor subpopulation (epi_tum1) influences the behaviors of other subpopulations through intercellular interactions, we further performed cell-cell communication analysis.

Figure 3.  Heat map showing the ssGSEA enrichment scores of certain gene sets in all subpopulations of epithelial cells, where each column represents a cell.

DownLoad: Full-Size Img PowerPoint

3.2.   Effects of the hypoxic tumor subpopulation on tumor progression

The cell-cell communication network was illustrated in Figure 4A, where the width of the edges represents the strength of the communication. As shown more intuitively in Figure 4B, the outputting/inputting communication strength of all cell subpopulations demonstrated that fibroblasts and endothelial cells served as the dominant signal sender/receiver, though monocytes have the highest expression levels of ligands/receptors (Figure A1). Furthermore, we identified statistically significant signaling pathways with p < 0.05, of which 39 signaling pathways involved by epi_tum1 as a signal sender were illustrated in Figure 5. There were six outputting communication patterns coordinating with each other to participate in these signaling pathways. Particularly, epi_nor2 was filtered out since it has relatively weak communications with other subpopulations at the level of signaling pathways. Distinct from the other three tumor subpopulations, epi_tum1 has a unique outputting communication pattern corresponding to eight signaling pathways, where ligands and receptors involved were listed in Table 1. In MPZ and OCLN signaling pathways, the communication way was cell-cell contact and the ligands and receptors were identical, which primarily mediated cell adhesion and tight junction between cells respectively ^[36,37].

Figure 4.  (A) Network diagram of cell-cell communication with the width of edges representing the strength of the communication, where the arrow points from the ligand cells to the receptor cells. (B) Scatter plot showing the outputting and inputting communication strength of all cell subpopulations.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Alluvial diagram presenting the correspondence between six outputting communication patterns with epi_tum1 involved 39 signaling pathways. Each color of line represents a distinct outputting communication pattern, and whether the dashed or solid line represents the communication strength of the signaling pathway involved.

DownLoad: Full-Size Img PowerPoint

Table 1.  Signaling pathways corresponding to the epi_tum1 outputting communication pattern.

Signaling pathway	Ligand	Log2FC	p-value	Receptor	Communication way
ANGPTL	ANGPTL4	1.017	0.002	ITGA5_ITGB1	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	CDH5	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	CDH11	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	SDC1	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	SDC2	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	SDC3	Secreted signaling
ANGPTL	ANGPTL4	1.017	0.002	SDC4	Secreted signaling
CALCR	ADM	1.358	< 0.001	CALCRL	Secreted signaling
EDN	EDN1	0.435	< 0.001	EDNRB	Secreted signaling
GDNF	ARTN	0.354	< 0.001	GFRA1_RET	Secreted signaling
L1CAM	L1CAM	0.240	< 0.001	ITGA4_ITGB7	Cell-cell contact
L1CAM	L1CAM	0.240	< 0.001	L1CAM	Cell-cell contact
MPZ	MPZL1	0.355	< 0.001	MPZL1	Cell-cell contact
OCLN	OCLN	0.371	< 0.001	OCLN	Cell-cell contact
VEGF	PGF	-0.781	0.042	FLT1	Secreted signaling
VEGF	VEGFA	1.963	< 0.001	FLT1	Secreted signaling
VEGF	VEGFA	1.963	< 0.001	KDR	Secreted signaling
VEGF	VEGFA	1.963	< 0.001	FLT1_KDR	Secreted signaling
VEGF	VEGFB	0.103	0.027	FLT1	Secreted signaling

---

 | Show Table

DownLoad: CSV

Angiogenesis-related signaling pathways were revealed to be participated predominantly by epi_tum1 such as CALCR, VEGF, EDN, and ANGPTL (Figure 6A–D). Calcitonin receptor-like receptor (CALCRL) was one receptor of adrenomedullin (ADM), a 52 amino acid peptide with 24% homology to calcitonin gene-related peptide (CGRP) initially discovered in 1993 ^[38]. ADM could promote proliferation and migration of endothelial cells thus facilitating vasodilation and angiogenesis ^[38,39], which was functionally similar to vascular endothelial growth factor A (VEGFA) ^[40,41]. Except for being synthesized and secreted by endothelial cells ^[42], ADM has also been reported as an autocrine growth factor within TME ^[43]. As shown in Figure 6A, B, endothelial cells are the only signal receivers of CALCR and VEGF signaling pathways. Compared to other cell subpopulations, epi_tum1 is the dominant signal sender of these two signaling pathways promoting the proliferation and migration of endothelial cells, proving that the autocrine growth factor ADM was secreted primarily by the hypoxic tumor subpopulation within TME. Moreover, endothelin 1 (EDN1) and angiopoietin like 4 (ANGPTL4) are ligands secreted by epi_tum1 involved in the EDN and ANGPTL signaling pathways, whose roles in the regulation of angiogenesis have also been demonstrated in previous studies ^{[44,45,46,47]}. As illustrated in Figure 6C–D, fibroblasts are also participating in both signaling pathways. In summary, the hypoxic tumor subpopulation within TME could accelerate abnormal angiogenesis and tumor progression through the above signaling pathways due to its up-regulation of ADM, VEGFA, EDN1, and ANGPTL4, which were shown in Figure A1.

Figure 6.  Chord diagram visualizing the cell-cell communication in (A) CALCR; (B) VEGF; (C) EDN; (D) ANGPTL; (E) GDNF and (F) L1CAM signaling pathways. The color of the arrow represents the ligand cells and the direction points to the receptor cells.

DownLoad: Full-Size Img PowerPoint

Besides, epi_tum1 serves as the primary signal sender of GDNF and L1CAM signaling pathways (Figure 6E–F). Artemin (ARTN) was a member of the glial cell derived neurotrophic factor (GDNF) family ligands ^[48], the expression of which was regulated by estrogen ^[49]. It has been demonstrated that ARTN was highly expressed in breast cancer ^[50] and could enhance the metastasis and invasion of estrogen receptor-positive cells ^[51,52]. Consistent with previous studies, ARTN secreted by epi_tum1 mainly targets the ESR1-positive subpopulation (epi_nor1) within TME as shown in Figure 6E. L1 cell adhesion molecule (L1CAM) and ANGPTL4 have been demonstrated to accelerate the vascular metastasis of tumor cells from the breast to the lungs ^[53], whose expression was mediated by hypoxia inducible factor 1 subunit alpha (HIF1A), especially under hypoxia ^[53]. The expression levels of ARTN and L1CAM were illustrated in Figure A1. These results suggest that the hypoxic tumor subpopulation contributes to the metastasis and invasion of tumors.

3.3.   Prognostic implications of the hypoxic tumor subpopulation

Totally 10 ligands secreted by epi_tum1 involved in above eight signaling pathways were listed in Table 1, of which PGF and VEGFB were filtered out due to not meeting the criteria for marker genes of epi_tum1 (log₂FC > 0.2 and p-value < 0.05). For the rest eight genes, we investigated their potential correlations with immune cells infiltration in TNBC using TIMER. As shown in Figure 7, the infiltration levels of dendritic cells had significant positive correlations with the expression levels of L1CAM (R = 0.196, p = 9.58e-03), EDN1 (R = 0.272, p = 2.89e-04), VEGFA (R = 0.162, p = 3.24e-02), and ADM (R = 0.202, p = 7.65e-03). The infiltration levels of macrophages had significant positive correlations with the expression levels of ANGPTL4 (R = 0.15, p = 4.79e-02), OCLN (R = 0.182, p = 1.60e-02), and MPZL1 (R = 0.249, p = 9.42e-04). The infiltration levels of CD4+ T cells had significant positive correlations with the expression levels of L1CAM (R = 0.296, p = 7.50e-05) and EDN1 (R = 0.253, p = 7.57e-04). The infiltration levels of neutrophils had significant positive correlations with the expression levels of MPZL1 (R = 0.177, p = 1.98e-02). On the other hand, the infiltration levels of B cells had significant negative correlations with the expression levels of ARTN (R = -0.175, p = 2.11e-02), ANGPTL4 (R = -0.273, p = 2.70e-04), VEGFA (R = -0.236, p = 1.69e-03), and ADM (R = -0.284, p = 1.5e-04). These results suggested that these genes modulate infiltration of immune cells into tumor tissues in TNBC.

Figure 7.  Correlation of expression levels of eight marker genes of epi_tum1 with immune infiltration levels in TNBC.

DownLoad: Full-Size Img PowerPoint

Then, we performed univariate Cox regression analysis to identify whether these eight genes were associated with OS of TNBC patients. As shown in Figure 8A, ARTN and L1CAM were considered as risk factors associated with prognosis due to hazard ratio (HR) > 1 and p-value < 0.05, then used for the construction of the prognostic model. The risk score was calculated as follows: Risk score = 0.3893 × Exp_ARTN + 0.1824 × Exp_L1CAM. The time-dependent ROC curves were shown in Figure 8B, where the AUC values reached 0.802, 0.728, and 0.637 at the time point of one, three, and five years respectively, indicating the reliability of the model in terms of predicting OS of TNBC patients. Moreover, as the Kaplan-Meier curves are shown in Figure 8C, the high-risk group was associated with worse OS, suggesting high expression levels of ARTN and L1CAM imply poor prognosis.

Figure 8.  (A) Forest plot showing the results of univariate Cox proportional hazards regression analysis. (B) Time-dependent ROC curve for predicting the OS of TCGA patients. (C) Kaplan-Meier curve showing the OS of GEO patients stratified by median risk score.

DownLoad: Full-Size Img PowerPoint

Considered as metastasis-related marker genes of the hypoxic tumor subpopulation in this study, the roles of ARTN and L1CAM in angiogenesis have also been studied. ARTN has been reported to stimulate de novo tumor angiogenesis in ER-positive BRCA partially mediated by up-regulation of VEGFA ^[54]. Furthermore, the sixth Ig-like domain of L1CAM (L1Ig6) has been demonstrated as a pro-angiogenic factor in vivo ^[55,56,57]. A recent study has developed and verified a risk score model based on TNBC-specific differential HRGs (dHRGs) for prognostic prediction in TNBC patients ^[14]. Through differential expression analysis and survival analysis, they identified 48 dHRGs associated with prognosis. The eventual 3-gene dHRGs signature consisted of ALDOA, PFKL, and PGK1, which performed well in predicting prognosis and distinguishing OS between high- and low-risk groups. All three genes (ALDOA, PFKL, and PGK1) considered as risk factors in TNBC patients by them were members of the hypoxia gene set containing 200 HRGs. In contrast to the previous study ^[14], neither ARTN nor L1CAM of the prognostic model constructed in this study belongs to the hypoxia gene set. However, both of them were transcriptionally activated by HIF1A according to previous studies ^[53,58]. Specifically, ARTN was identified as a responsive factor to hypoxia that promotes tumor progression in hepatocellular carcinoma ^[58]. Therefore, the relationship between ARTN and hypoxia in TNBC deserves to be further investigated. Figure 9A showed the expression levels of ARTN and L1CAM in tumor tissues of different BRCA subtypes and normal tissues, which were investigated using TIMER. It indicated that both ARTN and L1CAM were up-regulated in TNBC (i.e., BRCA-Basal) compared to normal tissues. We further analyzed their expression correlation in four subtypes of BRCA. As illustrated in Figure 9B, ARTN and L1CAM had a significant positive correlation in TNBC, suggesting that their crucial roles in tumor progression. In summary, the hypoxic tumor subpopulation within TME tends to overexpress angiogenesis- and metastasis-related genes in response to hypoxia, thus presenting a self-perpetuating state. In addition, the marker genes of the hypoxic tumor subpopulation ARTN and L1CAM were indicators of poor prognosis for TNBC patients.

Figure 9.  (A) The expression levels of ARTN and L1CAM in tumor tissues of different BRCA subtypes and normal tissues; (B) The correlation of expression levels between ARTN and L1CAM in TNBC.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In conclusion, we revealed the considerable heterogeneity within TME based on scRNA-seq data of TNBC, focusing on the roles of different cell subpopulations, especially epi_nor1 and epi_tum1 with high expression levels of ESR1 and HRGs respectively. Subsequently, the intercellular communication network was inferred to further explore how distinct communication patterns were coordinately involved in diverse signaling pathways. Particularly, the hypoxic tumor subpopulation was identified as the predominant signal sender of angiogenesis- and metastasis-related signaling pathways such as CALCR, VEGF, GDNF, and L1CAM. Finally, ARTN and L1CAM were applied for the construction of the risk score model, both of which were risk factors during tumor progression, indicating the hypoxic tumor subpopulation with self-perpetuating capacity could cause a dismal prognosis.

Acknowledgments

This work was supported by the Beijing Natural Science Foundation (No. 2202002) and the Chinese Natural Science Foundation project (No. 21173014).

Conflict of interest

The authors declare that they have no conflicts of interest.

Appendix

Figure A1.  Violin plot presenting the relative expression levels of certain genes among all cell subpopulations.

DownLoad: Full-Size Img PowerPoint