A Risk Prediction Model of DNA Methylation Improves Prognosis Evaluation and Indicates Gene Targets in Prostate Cancer

Figures & data

Table 1. The clinicopathologic variables of patients with prostate cancer included in the study.

Characteristics	Value^â€
Patients (n)	201
Age (year), median (IQR)	62.0 (56.0–66.0)
PSA (ng/ml), median (IQR)	6.8 (4.9–9.4)
Tumor weight (g), median (IQR)	200.0 (110.0–263.0)
Risk score from the model, median (IQR)	-0.05 (-0.24–0.13)
Pathological Gleason score, n (%):
  – ≤6	31 (15.4%)
– 7 (3 + 4)	64 (31.8%)
– 7 (4 + 3)	42 (21.0%)
– ≥8	64 (31.8%)
Prior malignancy, n (%):
  – Yes	12 (6.0%)
  – No	189 (94.0%)
Race, n (%):
  – White	160 (79.6%)
  – Black or African American	35 (17.4%)
  – Asian	6 (3.0%)
Residual tumor, n (%):
  – R0	148 (73.6%)
  – R1/R2/Rx	47 (23.4%)
  – NA	6 (3.0%)
TNM (T), n (%):
  – T2	98 (48.8%)
  – T3	100 (49.8%)
  – T4	1 (0.5%)
  – NA	2 (0.9%)
TNM (N), n (%):
– N0	156 (77.6%)
– N1	14 (7.0%)
– NA	31 (15.4%)
TNM (M), n (%):
– M0	190 (94.5%)
– NA	11 (5.5%)
Operation method, n (%):
– Robotics	99 (49.3%)
– Open	102 (50.7%)
Risk based on the model, n (%):
– High risk	84 (41.8%)
  – Low risk	117 (58.2%)
Outcome, n (%):
Cancer-specific death or biochemical recurrence:	27 (13.4%)
– Disease free	174 (86.6%)

† The column ‘Value’ represents the details of clinicopathologic variables of patients in this study.

IQR: Interquantile range; NA: Not analyzed; PCa: Prostate cancer; PSA: Prostate-specific antigen; T: Tumor; TNM: Tumor, node, metastasis.

Table 2. The details of the six CpG sites in the prediction model.

CpG site	Chromosome	Start on the chromosome	End on the chromosome	Associated gene	HR	95% CI	p-value	Coefficient in the model
cg09279544	chr19	8609509	8609510	ADAMTS10	19	(0.86–430)	0.062	0.303840332
cg18390610	chr17	50201566	50201567	COL1A1	18	(1.1–290)	0.042	0.229938118
cg25741646	chr1	51729516	51729517	OSBPL9	11	(1.2–110)	0.032	0.648553628
cg21495612	chr14	95875913	95875914	TUNAR	7.4	(1.1–48)	0.036	0.147254504
cg23514619	chr6	43269592	43269593	TTBK1	6.4	(1.1–37)	0.039	0.196243553
cg23258881	chr1	150150847	150150848	PLEKHO1	0.043	(0.0031–0.6)	0.019	-1.205666942

HR: Hazard ratio; LASSO coefficient: The risk coefficients of each CpG site in the least absolute shrinkage and selection operator regression.

Figure 1. Construction of the model by least absolute shrinkage and selection operator.

(A) The minimum misclassification error was obtained by cross-validation. (B) Based on the minimum of misclassification error, the inclusion sites and their coefficients of the model were obtained. (C–E). The training group. (F–H) The validation group. (C&F) The distribution of population risk scores and risk classification based on the cut-off value. The population was divided into high-risk and low-risk groups. (D & G) The relationship between the survival outcome of patients and the cut-off value. They indicated that the probability of poor prognosis would increase when the individual’s risk score in the model was greater than 0.01. (E & H) Heat maps for six CpG sites in the model. The expression patterns of six CpG sites were different between the low-risk group and the high-risk group. 0 represents disease free, and 1 represents cancer-specific death or biochemical recurrence.

Figure 2. Verification of the effectiveness of the model.

(A) The receiver operating characteristic curve for the training group (3, 5 and 7 years). (B–D) Kaplan–Meier curve for survival analysis. (B) The curve for the training group. (C) The curve for the validation group. (D) The curve for the total 201 patients (the 201 patients were randomly divided into the training and validation groups in this study).

AUC: Area under curve; TP: True positive.

Table 3. The association between methylation status (risk) and clinicopathologic variables of prostate cancer.

Clinicopathologic variables	Risk based on the model		p-value
	Low risk (n = 117)	High risk (n = 84)
Age (year), median (IQR)	62.0 (57.0–66.0)	61.0 (55.8–66.3)	0.850^â€
PSA (ng/ml), median (IQR)	7.0 (4.9–9.8)	6.3 (4.9–8.5)	0.357^â€
Tumor weight (g), median (IQR)	200.0 (110.0–267.0)	200.0 (109.5–251.3)	0.687^â€
Pathological Gleason score, n (%):
  – ≤6	21 (17.9%)	10 (11.9%)	0.009^‡
– 7 (3 + 4)	46 (39.3%)	18 (21.4%)
– 7 (4 + 3)	19 (16.2%)	23 (27.4%)
– ≥8	31 (26.6%)	33 (39.3%)
Prior malignancy, n (%):
– Yes	9 (7.7%)	3 (3.6%)	0.224^â€
– No	108 (92.3%)	81 (96.4%)
Race, n (%):
– White	90 (76.9%)	70 (83.3%)	0.585^§
– Black or African American	23 (19.7%)	12 (14.3%)
– Asian	4 (3.4%)	2 (2.4%)
Residual tumor, n (%):
– R0	86 (73.5%)	62 (73.8%)	0.736^‡
– R1/R2/Rx	26 (22.2%)	21 (25%)
– NA	5 (4.3%)	1 (1.2%)
TNM (T), n (%):
– T2	66 (56.4%)	32 (38.1%)	0.009^§
– T3	48 (41.0%)	51 (60.7%)
– T4	1 (0.9%)	1 (1.2%)
– NA	2 (1.7%)	0 (0.0%)
TNM (N), n (%):
– N0	94 (80.3%)	62 (73.8%)	0.074^‡
– N1	5 (4.3%)	9 (10.7%)
– NA	18 (15.4%)	13 (16.5%)
Operation method, n (%):
– Robotics	61 (52.1%)	38 (45.2%)	0.335^‡
– Open	56 (47.9%)	46 (54.8%)
Outcome, n (%):
– Cancer-specific death or biochemical recurrence	5 (4.3%)	22 (26.2%)	<0.001^‡
– Disease free	112 (95.7%)	62 (73.8%)

† p-values were calculated by the Mann–Whitney U test.

‡ χ2 test.

§ Fisher’s exact test.

IQR: Interquantile range; NA: Not analyzed; PSA: Prostate-specific antigen; T: Tumor; TNM: Tumor, node, metastasis.

Table 4. Univariate and multivariate COX regression for clinical variables and the risk prediction model.

Prognostic index	Univarite Cox regression			Multivariate Cox regression
	HR	95% CI	p-value	HR	95% CI	p-value
Age (year)	1.00	0.98–1.10	0.252
PSA (ng/ml)	1.00	1.00–1.00	0.100
Tumor weight (g)	1.00	1.00–1.00	0.140
Risk score from the model	19.00	6.00–62.00	<0.001	10.30	2.43–43.30	0.002
Pathological Gleason score	1.90	1.30–2.90	0.002	1.20	0.68–2.00	0.562
Prior malignancy	NA	NA	NA
Race	0.80	0.35–1.80	0.600
Residual tumor	1.60	0.73–3.60	0.231
TNM (T)	4.10	1.90–8.60	<0.001	3.00	0.93–9.90	0.066
TNM (N)	6.20	2.50–16.00	<0.001	3.70	1.43–9.50	0.007
Operation method	0.85	0.39–1.90	0.686

HR: Hazard ratio; NA: Not analyzed; PSA: Prostate-specific antigen; T: Tumor.

Figure 3. Multivariate time-dependent receiver operating characteristic curves for the 201 patients.

(A) The multivariate time-dependent receiver operating characteristic curve (multi-tdROC) at 3 years. (B) The multi-tdROC curve at 5 years. (C) The multi-tdROC curve at 7 years.

AUC: Area under curves; N: TNM (N); ROC: Receiver operating characteristic.

Figure 4. Results of weighted correlation network analysis.

(A) Gene modules divided by weighted correlation network analysis. (B) The heatmap for the topological overlap measure adjacent matrix. The darker the color, the higher the genetic similarity. (C) Identification of the key module.

Figure 5. Enrichment analysis in the Metascape database.

(A) The top 20 enrichment terms. (B) The network of the enrichment terms.

Figure 6. Hub genes selection.

(A) Correlation between MM and GS for the risk scores and the yellow module. (B) The enrichment clusters related to PCa. (C) The candidate hub genes. (D) The expression of hub genes between normal and tumor samples in the Gene Expression Profiling Interactive Analysis database.

GS: Gene significance; MM: Module membership; N: Normal sample; PCa and PRAD: Prostate cancer; T: Tumor sample.

Figure 7. In vitro, the knockdown of FOXD1 inhibits cell proliferation, migration and invasion of prostate cancer.

(A) The relative expression of FOXD1 was investigated by quantitative real-time PCR in PCa cell lines LNCaP, DU145, PC-3, 22RV1 and a normal human prostate cell line WPMY-1. GAPDH was used as an internal control. (B) The expression of FOXD1 at protein level was evaluated by Western blot. (C) The relative expression of FOXD1 was investigated by qRT-PCR in PCa cell lines PC-3 and DU145 after lentivirus transfection. (D) The expression of FOXD1 at protein level was evaluated by Western blot after lentivirus transfection. (E) The effects of FOXD1 on the proliferation of PC-3 and DU145 cells were evaluated by Cell Counting Kit-8 tests. (F) The effects of FOXD1 on the proliferation of PC-3 and DU145 cells were evaluated by colony formation tests. (G) Statistical analysis of the results from colony formation assays. (H) The effects of FOXD1 on the migration of PC-3 and DU145 cells were evaluated by wound healing assays. (I) Statistical analysis of the results from wound healing assays. (J) Representative images of transwell assays from control, sh-control and sh-FOXD1 groups of PC-3 and DU145 cells. (K) Statistical analysis of the results from transwell assays. *p < 0.05, **p < 0.01, ***p < 0.001 vs the sh-control group.

Data sharing statement

The authors used DNA methylation data and RNA-seq data from patients with prostate cancer in the TCGA database to train the model and explore therapeutic targets. This study is the first in prostate cancer to construct a DNA methylation risk model and explore therapeutic targets based on the risk prediction model. The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.