Novel biomarkers predict prognosis and drug-induced neuroendocrine differentiation in patients with prostate cancer

Abstract

Background: A huge focus is being placed on the development of novel signatures in the form of new combinatorial regimens to distinguish the neuroendocrine (NE) characteristics from castration resistant prostate cancer (CRPC) timely and accurately, as well as predict the disease-free survival (DFS) and progression-free survival (PFS) of prostate cancer (PCa) patients.

Methods: Single cell data of 4 normal samples, 3 CRPC samples and 3 CRPC-NE samples were obtained from GEO database, and CellChatDB was used for potential intercellular communication, Secondly, using the "limma" package (v3.52.0), we obtained the differential expressed genes between CRPC and CRPC-NE both in single-cell RNA seq and bulk RNA seq samples, and discovered 12 differential genes characterized by CRPC-NE. Then, on the one hand, the diagnosis model of CRPC-NE is developed by random forest algorithm and artificial neural network (ANN) through Cbioportal database; On the other hand, using the data in Cbioportal and GEO database, the DFS and PFS prognostic model of PCa was established and verified through univariate Cox analysis, least absolute shrinkage and selection operator (Lasso) regression and multivariate Cox regression in R software. Finally, somatic mutation and immune infiltration were also discussed.

Results: Our research shows that there exists specific intercellular communication in classified clusters. Secondly, a CRPC-NE diagnostic model of six genes (HMGN2, MLLT11, SOX4, PCSK1N, RGS16 and PTMA) has been established and verified, the area under the ROC curve (AUC) is as high as 0.952 (95% CI: 0.882-0.994). The mutation landscape shows that these six genes are rarely mutated in the CRPC and NEPC samples. In addition, NE-DFS signature (STMN1 and PCSK1N) and NE-PFS signature (STMN1, UBE2S and HMGN2) are good predictors of DFS and PFS in PCa patients and better than other clinical features. Lastly, the infiltration levels of plasma cells, T cells CD4 naive, Eosinophils and Monocytes were significantly different between the CRPC and NEPC groups.

Conclusions: This study revealed the heterogeneity between CRPC and CRPC-NE from different perspectives, and developed a reliable diagnostic model of CRPC-NE and robust prognostic models for PCa.

Keywords: castration-resistant prostate cancer; cellular communication; neuroendocrine; prognosis; single-cell RNA-seq.

Figure 1

The main workflow of the study. SC-RNA, single-cell RNA; NS, normal sample; CRPC, castration-resistant prostate cancer; NEPC, Neuroendocrine prostate cancer;DEGs, differentially expressed genes; GSVA, Gene Set Variation Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ROC, receiver operating characteristic; DFS, disease free survival; PFS, progression free survival.

Figure 2

The scheme of the prognosis model.

Figure 3

Analysis of single-cell RNA seq of 3 CRPC tissues, 3 NEPC tissues and 4 NS tissues. CRPC, castration resistant prostate cancer; NEPC, neuroendocrine prostate cancer; NS, negative samples; Fibro, Fibroblast; Basal/interm, Basal/intermediate; Endoth, endothelial; (A): the number of RNA features (nFeature_RNA) and absolute UMI counts (nCount_RNA) after quality control filtering of each cell. (B): We explored 1,500 high variable genes that exhibit high cell-to-cell variation, and the names of the top 10 genes are marked. (C): Using UMAP dimensionality reduction algorithm, 12165 cells from 10 samples were displayed. (D): Cells were classified into 16 clusters via t-SNE dimensionality reduction algorithm based on the source of the cluster, each cluster was marked with the source of the cluster plus the annotated cell types. There may exist 2 same cell types in 16 clusters. (E): Heatmap depicting expressions of top 10 marker genes among 16 clusters.

Figure 4

Pseudotime and trajectory analysis revealed the tendency curve among various cell types. (A) The pseudo time is shown as the depth of the color, the darker the blue, the smaller the pseudo time, which means that the cells appear earlier. The dots above represent cells. (B, C) Pseudotime-Density diagram demonstrated cells including immune cells (such as NK-T cell, B cell, Plasma) as well as Luminal/NE cells gather around the destination. X-axis means the value of principal component 1 (the first principal direction of maximum sample change) and Y-axis means the value of principal component 2.

Figure 5

Identification and functional analysis of CRPC-NE featured markers. (A) Marker gene expression for epithelial cells (KRT5, KRT19, KRT8, KRT18, AR) and neuroendocrine characteristic cells (CHGB, ENO2, LMO3, EZH2, SOX2, SIAH2), in which dot size and color represent percentage of marker gene expression and the averaged scaled expression value, respectively. (B) 102 genes with higher expression in NEPC than that in CRPC, and the latter is higher than that in NS were screened out. Then we selected genes shared between the SC-RNA data (102 genes) and Bulk-RNA data (1529 genes). (C, D) GO enrichment and KEGG pathway enrichment analysis of differentially expressed 102 genes. (E, F) Heatmap illustrating the differential KEGG pathways (upper panel) between CRPC_Luminal cluster and NEPC_Luminal/NE cluster at the single cell RNA-seq level, and discrepant KEGG pathway (lower panel) from the aspect of Bulk-RNA seq. The color indicates the level of pathway expression.

Figure 6

Intercellular ligand–receptor prediction among differernt clusters. (A) The chord diagram shows the expression of ANGTP, IL16, CSF, LIFR and OSM pathways among different cell clusters. In the peripheral ring, different colors represent different cells, Cells that send the arrow express the ligand, and cells that the arrow points to express the receptor, the more ligand-receptor pairs, the thicker the line. (B) Relative contribution of each ligand-receptor pair to the signal pathway, which may affect the overall communication network of the signaling pathway. CSF, IL16, LIFR and OSM pathways are shown in turn. (C) The extensive ligand-receptor mediated cellular interaction between different cell clusters of CRPC and NEPC has been further explored and demonstrated. The color gradient indicates the probability of cellular communication.

Figure 7

Identification of the markers to establish CRPC-NE diagnostic model. (A) The figure shows the weight of 12 genes to elucidate the importance of genes to disease classification (CRPC-NE or CRPC). The larger the”MeanDecreaseGini”index, the more likely this gene is to be classified as a characteristic gene. (B) Heatmap visualizing the expression levels of the six CRPC-NE diagnostic genes in the Cbioportal training cohort. (C) Results of neural network visualization: six CRPC-NE diagnostic genes were selected as the input nodes. Positive weights are connected by black lines, negative weights are connected by gray lines, and the thickness of the lines reflects the weight value. (D) The receiver operating characteristic (ROC) curves of 6-gene CRPC-NE diagnostic model in training cohort and validation cohort.

Figure 8

The distribution landscape of immune cell, and TMB pattern between CRPC and CRPC-NE (A) The difference of 22 immune infiltration between CRPC and CRPC-NE groups, red color indicates the abundance of immune cells in the latter, blue color indicates the abundance in the former. (B) Pseudotime trajectory analysis elucidated luminal/NE cluster and immune cells like NK, T, B and Plasma cells moved towards the termini of the trajectory. (C) Waterfall plots summarize the mutation landscape of 12 CRPC-NE featured genes in CRPC and CRPC-NE samples, showing that the mutation rate of these genes is low except CAMTA1.

Figure 9

A 2-gene prognosis model for DFS (NE-DFS signature) in the TCGA PanCancer training cohort and GSE21035 validation cohort in PCa. (A) Three genes significantly correlated with DFS were identified through LASSO regression analysis and ten-fold cross-validations for screening of the optimal parameter lambda (B) Kaplan–Meier curves displayed that high-risk group exhibited worse DFS than low risk group in TCGA PanCancer training group (n=276) and GSE21035 group (n=138). (C, D) Receiver operating characteristic (ROC) curves of the NE-DFS signature had better Predictive effectiveness than age, tumor stage and lymph node status to evaluate the predictability of DFS at 1-, 2- and 3- year in the TCGA PanCancer training cohort, similar phenomena were observed in the GSE21035 validation group.

Figure 10

Construction and validation of the prognosis model for PFS in the TCGA PanCancer cohort. (A) Four genes correlated with PFS were selected for multivariate analysis by LASSO regression analysis. (B) Kaplan-Meier plots evaluate the predictive ability of the constructed prognostic model in the TCGA PanCancer training cohort and internal validation cohort, respectively. (C, D) NE-PFS signature exhibited better predictive ability than other clinical features as displayed, the 1-, 2- and 3- year AUC for PFS was 0.700 (95% CI: 0.587−0.814), 0.659 (95% CI: 0.566−0.752), and 0.707 (95% CI: 0.622−0.792) in the TCGA PanCancer training cohort.

Figure 11

Nomogram construction and calibration plot validations for DFS and PFS prediction in PCa. (A, B) The composite nomogram consists of the DFS- or PFS- signature and clinical features of the individual patient, by adding the points from variables listed together, the 1-,3- and 5-year survival (DFS or PFS) probability can be inferred by the clinician. (C, D) Calibration curves for validation the consistence between 1-, 3- and 5-year (blue, red and orange color, respectively) inferred DFS and actual data in TCGA cohort and GSE21035 cohort. The dashed line represents the best match between the nomogram-predicted probability and the actual data evaluated by Kaplan–Meier analysis.