Identification of 4-methylation driven genes based prognostic signature in thyroid cancer: an integrative analysis based on the methylmix algorithm

Identification of MDGs in TC

The differentially methylated genes (DMGs) and differentially expressed genes (DEGs) were identified by using R package Limma and edge, respectively. The p-value <0.05, |logFC≥1| were used as the screening condition. Totally 486 genes were hypomethylated, and 397 genes were hypermethylated, while 1679 genes were downregulated, and 1751 genes were upregulated. Based on the R package methylmix, the DMGs and DEGs were reorganized as a normal methylated set, methylated cancer set, and gene expression cancer set. A p-value<0.05 and |cor|>0.3 were used as a standard criterion for the screening of the MDGs. Through a mixture model and Wilcoxon rank test, 62 MDGs were found. A heatmap was generated according to the methylation level of the 62 MDGs, with the R package pheatmap (Figure 1A). Also, the plots of correlation and the distribution map of the methylation degree were shown in Figure 1B, 1C. Representative pictures were selected in Figure 1, and the full results for MDGs were available in Table 1.

Figure 1. Overview of methylation driven genes in TC. (A) Heatmap of 62 methylation driven genes in TC. (B) Representative correlation plots of the MDGs, reflecting the correlation between expression and methylation levels of genes. (C) Representative distribution plots of MDGs, reflecting the distribution of methylation values.

Functional enrichment analysis of the MDGs

We conducted functional enrichment and pathway analysis using Metascape online on the base of 62 MDGs previously obtained. As shown in Figure 2A, the methylation-driven genes were mainly enriched in substrate-dependent cell migration, cellular response to mechanical stimulus, cell-substrate adhesion, cellular response to tumor necrosis factor, dendritic cell differentiation. Besides, the methylation gene set was mainly enriched in the regulation of JAK-STAT, activation of MAPK activity pathways (Figure 2B). Usually, hypermethylation inhibits transcription, while hypomethylation increased the expression of oncogenes. The results of the functional analysis revealed the potential mechanisms of MDGs.

Figure 2. Functional analysis of 62 MDGs based on the metascape. (A) Bar graph of enriched terms across input gene lists, colored by p-values. (B) The network of enriched terms colored by cluster-ID, where nodes that the same cluster-ID are typically close to each other. (C) Protein-protein interaction network identified in the MDGs.

Furthermore, to investigate the hub genes which play significant roles in TC, we constructed a PPI network and MCODE to identify the critical units. As shown in Figure 2C, module 1 included 11 edges and 12 nodes and involved ALDOC, PTPRC, FAS, etc. Meanwhile, module 2 included 1 edge and 2 nodes, and involved CLDN1, TACSTD2.

Construction and validation of the MDGs based on predictive model

Normalized methylation data with complete clinical information were obtained from the TCGA database. To enhance the reliability and accuracy of the predictive model based on MDGs, we treated 82 patients from our center as a validation set. Univariate cox regression analysis was first performed to identify the genes that were significantly associated with the prognosis from the 62 MDGs previously obtained. A total of 4 MDGs were then found to be prognosis-related according to their methylation β value (Figure 3A). Subsequently, Multivariate Cox regression analysis was conducted, and 4 MDGs (ALDOC, C14orf62, DVL1, PTPRC) were eventually selected to construct a predictive model. A risk score of each patient was generated as following: (3.23)× β value (ALDOC) + (2.98)× β value (C14orf62)+ (-8.96)× β value (DVL1) + (18.23)× β value (PTPRC). Then, the patients were divided into a high-risk group (n= 253) and a low-risk group (n= 254) according to the risk scores. The median survival time of patients with high risk was significantly higher than that of the patients with low risk (Figure 3B). Subsequently, the ROC curve was performed to assess the efficiency prognosis prediction of this model, and we found the AUC of the ROC curve of the predictive model exceeded that of the individual genes (Figure 3C).

Figure 3. Construction of MDGs based prognostic signature. (A) Kaplan-Meier plots of ALDOC, C14orf62, DVL1, PTPRC. (B) The Kaplan-Meier plot of the prognostic signature based on the median of the risk score. (C) The ROC curve for assessing the prediction value of the signature.

Furthermore, to evaluate the reliability and accuracy of the predictive model, we validated the power of the model in our own center. First, qMSP was performed to verify the expression of differentially methylated genes. As shown in Figure 4A, 4B, ALDOC, C14orf62, DVL1, PTPRC were hypermethylated, and their corresponding protein levels were decreased, which was consistent with the results from the TCGA cohort. Although there was no statistically significant difference between the methylation of C14orf62 and overall patient survival, the 4 MDGs joint model could successfully divide the TC patients into long-term OS group and short-term OS group, which indicated the effectiveness of the predictive model constructed based on 4 MDGs (Figure 4C, 4D). The AUC of the ROC curves of this model was 0.81 (5-year OS) (Figure 4E). Besides, as shown in Figure 5, we visualized the distribution of the 4 MDGs both on the training set and the validation set according to their risk scores by using R package ggrisk, suggesting that the higher the riskscore, the worse the prognosis. Also, the validation cohort was consistent with the TCGA cohort.

Figure 4. Validation of constructed prognostic signature. (A) The methylation level was detected by qMSP. (B) The protein level was detected by IHC. (C) The Kaplan-Meier plots of the ALDOC, C14orf62, PTPRC. (D) The Kaplan-Meier plot of the prognostic signature based on the median of the risk score in the validation set. (E) The ROC curve for assessing the prediction value of the signature in the validation set.

Figure 5. Visualization of the prognostic signature. (A) Visualization of the signature in the TCGA cohort. (B) Visualization of the signature in the validation cohort.

Construction and validation of the nomogram

To develop an effective method that could help predict an individual's prognostic risk of TC, we developed a nomogram based on the panel of these 4MDGs in the training set and validation set (Figure 6A, 6C). Calibration plots for predicting the 5-year survival rate showed that the nomograms performed well as compared with an ideal model both in the training set and validation set (Figure 6B, 6D).

Figure 6. Construction and validation of nomogram. (A) The nomogram for predicting the survival with 3- and 5-year OS. (B) The calibration plot for validation of the model. (C) The nomogram for predicting the survival with 3- and 5-year OS in the validation set. (D) The calibration plot for validation of the model in the validation set.

Knockdown of ALDOC, C14orf62, DVL1, and PTPRC significantly promoted the the proliferation, migration, and invasion of TC cells

To investigate the regulation of ALDOC, C14orf62, DVL1, and PTPRC on the proliferation, migration, and invasion of TC cells, the knockdown vectors of ALDOC, C14orf62, DVL1, and PTPRC were constructed. We found that the cell proliferation ability was remarkably increased after treatment with sh-ALDOC, sh-C14orf62, sh-DVL1, and sh-PTPRC (Figure 7A). In addition, sh-ALDOC, sh-C14orf62, sh-DVL1, and sh-PTPRC significantly promoted the cell invasion (Figure 7B, 7C) and migration ability (Figure 7D, 7E).

Figure 7. Influence of ALDOC, C14orf62, DVL1, and PTPRC on the proliferation, migration, and invasion of TC cells. (A) Influence of ALDOC, C14orf62, DVL1, and PTPRC on the proliferation of TC cells. (B) Influence of ALDOC, C14orf62, DVL1, and PTPRC on the invasion of TC cells (Scale bar=50 μm). (C) Quantitative analysis of the influence of ALDOC, C14orf62, DVL1, and PTPRC on the invasion of TC cells. (D) Influence of ALDOC, C14orf62, DVL1, and PTPRC on the migration of TC cells (Scale bar=500 μm). (E) Quantitative analysis of the influence of ALDOC, C14orf62, DVL1, and PTPRC on the migration of TC cells.

Human tissues

The samples were collected from patients that underwent surgery at the department of thyroid surgery, the affiliated Hospital of Putian University. The study was approved by the Ethics Committee of Affiliated Hospital of Putian University (Approval number: 2019-036), and all patients or their guardians signed the consent form. Totally 82 TC tissues and adjacent normal tissues were used as a validation cohort. All tissues were immediately frozen in liquid nitrogen after the resection and stored at -80° C according to the manufacturer's protocol. The inclusion criteria is that the TC patients should be diagnosed via biopsy and histological testing. The exclusion criteria is that (1) Pregnant patients; (2) Significant immunodeficiency disease patients; (3) Patients with severe underlying diseases.

DNA extraction, bisulfite modification and qMSP assay

Genomic DNA was extracted from the 82 TC tissues and paired normal tissues by using the TIANamp Genomic DNA Kit (#DP304, TIANGEN, Beijing, China). EZ Methylation-Direct Kit (#D502, Zymo Research, Irvine, CA, USA) was used for bisulfite conversion of 1000 ng from each sample according to the requirement of the instruction. The purified DNA was stored at −20° C. The purity/quality of DNA was measured using OD 260 nm/280 nm method, and the was OD 260 nm/280 nm value was 1.96 detected using NanoDrop 8000 (Thermo, Waltham, Massachusetts, USA).

qMSP was used to quantify methylation levels of the MDGs. The qMSP reaction consisted of a 1×ABI master mix containing Taq polymerase, dNTPs, 10 ng bisulfite-converted DNA, SYBR green dye and ROX as a passive dye (Thermo, Waltham, Massachusetts, USA) and 200 nM of specific primers. Specific primers for the promoter region of target genes were designed (Supplementary Table 1). Moreover, the ACTB gene, which did not contain any CpG dinucleotide, was used as a reference. The relative methylation level of each MDG was calculated using 2^{(ΔCttarget- ΔCtreference)}.

Data acquisition and pretreatment

The normalized mRNA expression and methylation data of TC were downloaded from the TCGA database. All data were transferred from the Genome Data Commons (GDC, https://gdc.cancer.gov/). The data includes mRNA expression from 568 specimens, including 510 tumor specimens and 58 normal specimens, as well as methylation level from 571 specimens, including 515 tumor specimens and 56 normal specimens. The R package edge was used to process the downloaded data to obtain levels of normalized expression and methylation. In addition, 507 samples with TC, which were accompanied by both clinical and expression information, were obtained from TCGA. The study was performed under the publication guidelines of TCGA and no restrictions were imposed in this research.

Methylation-driven genes in TC

The methylmix is an efficient and accurate algorithm used for automatically analyzing the aberrantly methylated genes and the correlation between gene expression and methylation level. Before the application of the methylmix, preparation of documents for three datasets: cancer methylation data (METCancer), normal methylation data (METnormal), and matched gene expression data (GECancer), were normalized and obtained from the previous step. The 3 mentioned files were served as input files and then submitted to the methylmix according to the requirements of the algorithm. Briefly, three main steps were carried out: The first step was to identify the methylation status of genes by using the β-mixed model, which was applied for avoiding overfitting according to the Bayesian information criterion. Second, the Wilcoxon rank-sum test was performed to compare the methylation state between tumor and normal tissues to identify the significant difference on the basis of a Q-value of 0.05, which was performed using P-value multiple testing correction with false discovery rate (FDR). Finally, linear regression was used to model the expression of genes in terms of their DNA methylation. Considering the numbers of samples, |cor|> 0.3 was used as a standard criterion.

Functional and pathway analysis

Furthermore, to explore the underlying mechanisms of these MDGs, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed via the Metascape database (http://metascape.org). Additionally, the protein-protein interaction (PPI) network was used to identify the densely connected regions further. Here, P-value < 0.05 was used as a cutoff.

Prognostic risk model construction

To further search for MDGs with prognostic value, Cox regression analysis was applied to assess univariate and multivariate associations of risk factors with the development of TC. Univariate Cox regression analysis was first performed to identify survival-associated MDGs. MDGs with a p-value < 0.05 were considered to be target genes. Multivariate Cox regression analysis was then applied to eliminate non-independent prognostic predictor genes. The prognostic model was established based on these target genes weighted by their estimated regression coefficients. Subsequently, Areceiver operating characteristic (ROC) was constructed to evaluate the predictive value of the prognostic model. The area under curve (AUC) was used to determine the efficiency of the model. The R package survival was applied in this section based on the Rstudio (Version 1.3.1073). Here, the data obtained from TCGA was used as a training set, and the 82 samples from our center were employed as the validation set.

Construction and validation of the nomogram

The risk target genes obtained mentioned above were used for the nomogram model building to generate the predictive probability of 3-year and 5-year overall survival (OS). The C-index and calibration curves were applied to evaluate the internal validation of the nomogram. Besides, external validation of the predictive model was evaluated in a validation cohort. The R package rms was used to plot nomograms and calibration curvesbased on the Rstudio (Version 1.3.1073).

Immunohistochemistry validation

Immunohistochemistry was performed according to the manufacturer's instruction. Sections were incubated with primary antibodies against aldolase C (AldoC), C14orf62, dishevelled 1 (DVL1), and protein tyrosine phosphatase receptor type C (PTPRC) overnight at 4° C. The image was captured at an appropriate magnification in the microscope (Nikon Microsystems, Shanghai, China).

CCK8 assay

Cells were plated into the 96-well plate. After treatment with sh-ALDOC, sh-C14orf62, sh-DVL1, and sh-PTPRC for 48 h, the CCK-8 kit (#C0037, Beyotime, Beijing, China) was applied to detect cell proliferation. 20 μL regent was added to each well, OD at 450 nm was measured after 1 h incubation.

Wound healing method

Cells were seeded into 6-well plates firstly. After reaching 70% confluence, 200 μL pipette tip was used to drawn a line in the middle plate. The medium was replaced with new medium. Then, the cells were incubated on the condition of 5% CO₂ and 37° C. The cells were captured at 0 h and 48 h, and the relative migrated distance was analyzed.

Transwell assay

The transfected cells were seeded into the upper chamber of a Transwell plate. The lower chambers were supplemented with 15% FBS (Gibco, Langley, OK, USA). After incubation (48 h), the cells were fixed using 70% ethanol (30 min). After staining with 0.2% crystal violet (10 min), the cells in the lower chamber were counted using a microscope.

Ethics approval and consent to participate

The study was approved by the Ethics Committee of Affiliated Hospital of Putian University, and all patients or their guardians signed the consent form.

Availability of data and materials

All data generated or analysed during this study are included in this published article.