3.1 Machine annotation combined with manual annotation used to identify cell types

A total of 7 cell clusters, epithelial cells (EPCAM), fibroblasts (DCN), vascular endothelial cells (VWF), T cells (CD3E), B cells (CD79B), mast cells (KIT), and monocytes/macrophages (CD68), were tentatively identified according to the genetic profile of the cells and the typical marker genes (NKG7 was used to validate NK cells that are subsequently identified). Subsequently, Azimuth assisted in acquiring more comprehensive annotation outcomes (Fig. 1A–C). A range of lymphocyte subtypes were identified from immune cells, including T cells, B cells, and monocytes/macrophages. Furthermore, NK cells, DC cells, and a few other blood cells were isolated from these cells (Fig. 1A–C). Ultimately, calculations were made for the cell count ratios of various cell types (Fig. 1D), and the expression of characteristic genes representative of specific cell types was depicted (Fig. 1E) (Supplementary Figure 1A). The identification of characteristic genes demonstrates the dependability in classifying cell types (Fig. 1E) (Supplementary Figure 1A). Despite the varied distribution of cell types across samples, each sample exhibits a relatively comprehensive range of cell type expressions (Fig. 1C, D).

Fig. 1

A-C t-distribution random neighbor embedding (t-SNE) plot of cells from 4 patients, colored by cell types (A), specific cell types with fine cell subpopulations (B), sample sources (C). D Relative proportions of cell types per prostate cancer patient. E t-SNE plots of marker genes for each cell subset.

3.2 Heterogeneous clustering of TAECs and their general characteristics

TAECs cells were extracted from epithelial cells, and the 16 metagene expression patterns obtained by non-negative matrix decomposition (NMF) were characterized by four gene expression patterns (Fig. 2A), of which one gene expression pattern with low confidence was eliminated. Three heterogeneous TAECs subpopulations with different gene expression patterns (TAECs1, TAECs2, and TAECs3) were finally determined (Fig. 2A, B). Hierarchical clustering reveals significant variation among the three subpopulations, with the TAECs2 and TAECs3 subpopulations showing closer affinity (Fig. 2A). Either a TSNE map or the hierarchical cluster diagram reveals clear divisions among these subpopulations, demonstrating that the outcomes of the clustering effectively depict the diversity present in TAECs (Fig. 2A, B). Furthermore, each sample contained three subsets of TAECs, indicating that these heterogeneities were not derived from individual differences (Fig. 2C).

Fig. 2

A Heat map depicting pairwise correlations of 16 tumor metagene expression patterns from 4 tumor patients. Clustering identified 4 coherent feature sets across tumors samples. B t-distribution random neighbor embedding (t-SNE) plots of tumor cells from 4 patients, colored by TAECs subpopulations. C Relative proportion of TAECs types per prostate cancer patient. D t-SNE plots of characteristic gene expression for three heterogeneous TAECs subpopulations. E GSEA enrichment analysis of TAECs1 subsets versus total normal epithelial cells. F GSEA enrichment analysis of TAECs2 subsets versus total normal epithelial cells. G GSEA enrichment analysis of TAECs3 subsets versus total normal epithelial cells.

By examining the expression of some characteristic genes, we found that the TAECs2 subpopulation overexpressed the basal cell marker genes TP63, KRT5, KRT14, and KRT19 (Fig. 2D); the TAECs1 subpopulation overexpressed hormone and secretion-related genes associated with differentiated terminal epithelial cells, such as AR, NPY, KLK3, and KLK4 (Fig. 2D), which indicates that the three subpopulations might have time-dependent differentiation development features similar to normal epithelial cells. The gene expression heatmap was used to further characterize the expression features of the three subgroups (Supplementary Figure 1B).

Enrichment analysis was performed for each subpopulation compared to total normal epithelial cells, revealing that the TAECs1 subpopulation was enriched in the androgen receptor pathway (Fig. 2E); the TAECs2 subpopulation was enriched in the epithelial-mesenchymal transition and extracellular matrix pathways (Fig. 2F); and the TAECs3 subpopulation was enriched in the interferon pathway and various antigen processing and presentation, and vesicle transport pathways (Fig. 2G). Notably, this might be the case that three subpopulations exhibit biological characteristics associated with different stages of differentiation from normal epithelial cells, rather than true tumor behavior.

3.3 Pseudo-time analysis and paired enrichment analysis

Analysis of cell differentiation and developmental trajectories using monocle 2 validates the hypothesis that TAECs follow the trajectories of TAECs2, TAECs3 to TAECs1 (Fig. 3A–C), for which TAECs3 subpopulation closely related to the TAECs2 subpopulation (Fig. 3A–C). The characteristic gene expression trends of the three subpopulations were consistent with the inferred differentiation trajectories (Fig. 3D, E).

Fig. 3

A–C TAECs cell developmental trajectories inferred using monocle2 colored by TAECs subpopulations (A), inferred differentiation state (B), pseudo time (C). D Curve plots showing expression changes of characteristic genes related to differentiation state along pseudo time. Point colors correspond to cell state colors in B. E The distribution of TAECs is shown along pseudo time, Heatmap showing dynamic expression changes of selected genes and related pathways along pseudo time, clustering was carried out according to the dynamic expression changes of genes with pseudo-time. F Differential enrichment analysis (GSEA) of the TAECs1 subpopulation with its matched N 1 normal epithelium. G Differential enrichment analysis (GSEA) of the TAECs2 subpopulation with its matched N 2 normal epithelium. H Differential enrichment analysis (GSEA) of the TAECs3 subpopulation with its matched N 3 normal epithelium.

Given possible problems with previous enrichment analysis, following earlier findings on dimensionality reduction clustering from anchor point integration, clusters of normal epithelial cells, identical in differentiation to TAECs, were pinpointed. (N1, N2, N3) (Supplementary Figure 2A-C). The paired enrichment study was conducted on groups of three TAECs subtypes along with their matched normal epithelial cells, subpopulations in each pair were at an identical stage of differentiation (Fig. 3F–H). Findings indicated that each of the three subtypes possessed more potent anti-apoptosis properties than their corresponding normal epithelial cells, as demonstrated by the inverse enhancement in the P53, apoptosis, and NFKB-driven tumor necrosis factor signaling pathways (Fig. 3F–H). Additionally, the TAECs1 subpopulation shows an down-regulated enhancement in epithelial-mesenchymal transition, interferon signaling, and specific immune response routes (Fig. 2F), indicating its function as a final differentiation phase and immune avoidance. The TAECs2 subpopulation, resembling basal cells, show enhanced epithelial-mesenchymal transition activity compared to matched normal cells (Fig. 3G), suggesting possible stromal infiltration, though this is not histologically verified. TAECs3 subpopulations exhibit enhanced secretion capabilities (Fig. 3H).

3.4 Cell communication analysis

TAECs cells are the major senders of MIF signals (Fig. 4A–C). Compared with matched normal epithelial cells at the same differentiation stage, they show higher levels of expression of EGFR, ERBB2, and interferon receptors (Fig. 4D–F). In addition, TAECs express VEGF signaling at an earlier stage of differentiation than normal epithelial cells and might undergo a partial expression shift from VEGFA to VEGFB signaling (Fig. 4G).

Fig. 4

A Cell-to-cell communication mediated by the MIF signaling pathway. B Expression of all signaling genes associated with the MIF signaling pathway. C Heatmap visualization of the calculated centrality score to identify the main signaling role of the cell population in the MIF pathway. D Cell-to-cell communication mediated by EREG_EGFR_ERBB2 signaling pathways. E Expression of all signaling genes associated with the EGFR signaling pathway. F Expression of all signaling genes associated with the interferon signaling pathway. G Expression of all signaling genes associated with the VEGF signaling pathway. H Transcription factor activity heat map of tumor subsets and their paired normal epithelial cells at the same TSNE neighbor plot location based on anchor integration.

Tumor cells abnormally expressed cell adhesion related signals (SEMA4A, OCLN, NECTIN, MPZL1, JAG, ALCAM, CD46, CD99) (Supplementary Figure 3A), which were closely related to tumor survival and invasion, angiogenesis, and tumor immunity. Compared to matched normal epithelial cells, the expression of MHC molecules related to antigen presentation is higher in TAECs2 and TAECs3 subpopulation (Supplementary Figure 3A).

Compared to matched normal epithelial cells, the TAECs2 and TAECs3 subpopulations expressed higher levels of SEMA3C, CXCL16 (Supplementary Figure 4A-C), and WNT5B signal (Supplementary Figure 4D); the TAECs1 and TAECs3 subpopulations showed higher levels of BAG signal (Supplementary Figure 4E, F).

The TAECs2 subpopulation showed higher expression of LAMB3 and LAMA5 molecules compared to matched normal epithelial cells (Supplementary Figure 3B). Also, it, as well as TAECs3 subpopulation, exhibited higher expression of matrix-related receptors (SDC1, SDC4, ITG) (Supplementary Figure 3C), indicating a closer association with the extracellular matrix.

Compared to matched normal epithelial cells, the TAECs1 subpopulation had higher expression of FZD5 and LRP6 (Supplementary Figure 5A, B), which could serve as receptors for WNT10A signal in plasma cells; the TAECs2 subpopulation specifically expressed ADM signal, acting on endothelial cells (Supplementary Figure 5C); the TAECs3 subpopulation showed higher expression of TNFRSF15, acting on MAIT cells and ILC cells (Supplementary Figure 5D).

3.5 Transcription factor regulatory network analysis

Using scenic to infer the transcription factor regulatory network in epithelial cells (Supplementary Table 1), we calculated the transcription factor activity score (Supplementary Table 2) for three TAECs subpopulations with reference to matched normal epithelial cells at the same stage of differentiation. We extracted the top 20 transcription factor regulatory units with high and low transcriptional activity to draw the transcription factor activity heatmap (Supplementary Table 3) (Fig. 4H). We found that in the three TAECs subpopulations, the transcriptional activity of the transcription factors YY1, NKX3-1, and EHF was higher than in the matched normal epithelium at the same differentiation stage (YY1 regulatory unit ranked in the top 10) (Supplementary Table3).

As previously noted, the MIF signaling pathway is overexpressed in various TAECs subpopulations (Fig. 4A–C). The scenic-inferred transcription factor regulatory network suggests that YY1 acts as a dependable transcriptional upstream regulator for MIF (Supplementary Table 1), with its target genes containing diverse genetic activity elements, underscoring its vital function in tumor growth and progression. Since it was previously found that the functional enrichment of the three TAECs subpopulations were weakened in the apoptotic related pathway (Fig. 3F–H), we further intersected the differential genes enriched by the apoptotic and P53 pathway of the three subpopulations (Supplementary Tables 4), and found that ATF3 was the upstream regulator of most of these genes (Supplementary Tables 4), and the transcriptional activity of ATF3 in normal epithelial cells was higher than that of its paired TAECs (Supplementary Table3), suggesting that TAECs might inhibit ATF3 by some mechanism to escape apoptosis.

3.6 Correlations of prostate cancer immune infiltration with YY1, EHF, NKX3-1 and ATF3 expressions

Correlation analysis showed that the expression of YY1 was significantly positively correlated with the immune infiltration of CD4 resting memory T cells, resting dendritic cells and M1 macrophages (P<0.001, Cor= 0.30000817, 0.22245352, 0.18262690) (Supplementary Figure 6 A).

The expression of NKX3-1 was significantly positively correlated with the immune infiltration of resting mast cells and M1 macrophages (P<0.001, Cor= 0.163012839, 0.154186876). The expression of NKX3-1 was significantly negatively correlated with the immune infiltration of suppressor T cells (Tregs) and large chemical cells (P<0.001, Cor= 0.246266240, 0.218359890) (Supplementary Figure 6B).

The expression of EHF was significantly positively correlated with the immune infiltration of CD4 resting memory T cells and M1 macrophages (P<0.001, Cor= 0.298198031, 0.154686150). It was negatively correlated with the immune infiltration of inhibitory regulatory T cells (Tregs), CD8 T cells and activated natural killer cells (P<0.001, Cor= 303334670, 0.210507210, 0.155223607) (Supplementary Figure 6C).

The expression of ATF3 was significantly positively correlated with the immune infiltration of large chemical cells, neutrophils, naïve B cells and follicular helper T cells (P<0.001, Cor= 0.491918316, 0.241133868, 0.218357002, 0.177083614). It was significantly negatively correlated with the immune infiltration of resting mast cells and memory B cells (P<0.001, Cor= 0.391263842, 0.154324034) (Supplementary Figure 6D).

3.7 Pan-cancer analysis revealing the expression and prognosis of YY1, NKX3-1, EHF and ATF3

In order to further explore the prognosis of several transcription factors with aberrant transcriptional activity in prostate cancer, we subsequently performed a pan-cancer analysis (Supplementary Figure 7), and the results showed that YY1 and NKX3-1 were significantly overexpressed in tumor tissues of prostate cancer and ATF3 was significantly underexpressed in tumor tissues relative to adjacent tissues, which was consistent with the results obtained by single-cell analysis. Among them, the hazard ratio of YY1 in overall survival of prostate cancer was 11.9 (95% confidence interval (1.76−81.2, P = 0.0112), the hazard ratio of ATF3 in prostate cancer was 0.791 (95% confidence interval 0.631−0.991, P = 0.0411), and the risk ratio of ATF3 in prostate cancer was 0.88 (95% confidence interval 0.778−0.995, P = 0.0413). In addition, these four transcription factors were associated with an increased or decreased risk of poor prognosis in a variety of cancers, among which YY1, EHF, and NKX3-1 all showed a significant increase in the risk of consistent adverse prognosis in pancreatic cancer (YY1: DFI HR 8.11 95% confidence interval 2.18−30.1, P = 0.00177, OS HR 3.17 95% confidence interval 1.68−5.97, P = 0.000357, PFI). HR 2.63 (95% CI 1.41−4.89), P = 0.00232; EHF: DFI HR 1.57 95% CI 1.07−2.32, P = 0.0224), OS HR 1.26 95% CI 1.05−1.5, P = 0.0106, PFI HR 1.24 95% CI 1.06−1.46, P = 0.00875; NKX3-1: DFI HR 2.46 95% CI 1.41−4.28, P = 0.00145, OS HR 1.3 95% The confidence interval was 1.05−1.61, P = 0.0152, PFI HR 1.61, 95% confidence interval was 1.3−1.99, P = 0.106*10-6). <0.001), high expression of YY1 showed a consistent increase in the risk of poor prognosis of renal papillary cell carcinoma (DFI HR 3.71 95% confidence interval 1.27−10.9, P = 0.0166, OS HR 9.1 95% confidence interval 3.78−21.9, P = 0.832*10-8, PFI HR 4.19 95% confidence interval 2−8.78, P = 0.000149) (Supplementary Figure 8).