Clinical and prognostic pan-cancer analysis of m6A RNA methylation regulators in four types of endocrine system tumors

The expression of m6A regulatory genes is related to endocrine system tumors

In view of the important role of m6A regulatory genes plays in the development of tumors, we systematically studied the relationship between 23 m6A regulatory genes and each endocrine system tumor. The expression level of each m6A regulatory gene in endocrine system tumors was presented in the form of a heat map. The results show that, except for HNRNPC gene, the expression level of other m6A regulatory genes in ACC (Figure 1A) was distinct from that in normal tissues. Except for the IGF2BP3 gene, the expression of m6A regulatory genes in PCPG (Figure 1B) generally differed from normal tissue. The expression of 23 m6A regulatory genes in THCA (Figure 1A) was generally different from that in normal tissue. Except for IGF2BP2 and FTO gene, the expression of other m6A regulatory genes in THYM (Figure 1B) was generally distinct from normal tissue.

Figure 1.

The distribution of m6A regulatory genes in endocrine system tumors. (A–C) Heatmaps defined by the expression level of 23 m6A regulatory genes (red for up-regulated, green for down-regulated and the tree on the left shows the clustering results of different genes in different samples) and vioplots visualized the differential m6A regulatory genes in endocrine system tumors (assuming blue is normal tissue and red is tumors); (D) Correlation analysis of expression of 23 m6A regulatory genes in endocrine system tumors.

Then, due to different clinical variables, we divided endocrine system tumors into “Total 1” (ACC and THCA datasets) and “Total 2” (PCPG and THYM datasets). The expression level of m6A regulatory genes of each “Total” was analyzed. In the dataset of “Total 1” (Figure 1C), 23 genes were found to be different from normal tissue. In the dataset of “Total 2” (Figure 1C), except for ZNF217, the expression of m6A regulatory genes is generally different from that of normal tissue. These results indicate that most of the expression of m6A regulatory genes are associated with endocrine tumors.

We then analyzed the correlation of 23 gene expressions in tumors. Our study found that most of expression level of the genes are negatively correlated with METTL3 expression, but there are some exceptions. With the expression level of WTAP is elevated, METTL3 expression is in high level. This can be observed in each dataset (Figure 1D). These datasets illustrate the interaction of gene expression level between m6A regulatory genes.

Consensus clustering method for m6A regulators identified two clusters of endocrine system tumors

In the view of the result that most of the expression of m6A regulatory genes are associated with endocrine tumors, we further investigated the relationship between m6A regulatory genes and the clinicopathological characteristics of the tumors. Clustering classification was performed based on the expression correlation of the m6A RNA methylation regulators, and the clustering stability was increased from k = 2 to 10. The clinicopathological features of the two subgroups clustered were cluster 1 and cluster 2 according to k=2 in ACC datasets (Figure 2A). We further used principal component analysis (PCA) to compare the expression curves between cluster 1 and cluster 2 subgroups, and found that there were significant differences between subgroups (Figure 2A). It is showed that the cluster 1 subgroup (Figure 2D) was significantly correlated with better survival status (P < 0.01), higher who grade (P < 0.01), and higher T status (P < 0.01). The cluster 2 subgroup was associated with worse survival status, lower who grade (P < 0.01), and lower T degree (P < 0.01). In addition, we observed that the overall survival rate (OS) of the cluster 1 subgroup was significantly shorter than that of the cluster 2 subgroup (Figure 2D). We then clustered the other tumors in the same way, observed some similar results and other differences. According to the THCA dataset, we observed no significant difference in OS (Figure 2D) between the cluster 1 subgroup and the cluster 2 subgroup, and it may be necessary to observe the survival in patients with the disease for a longer time. Besides, if only from the aspects of living condition and gender analysis PCPG and THYM datasets (Figure 2B, E), we found no significant difference of clinical features.

Figure 2.

Identification of consensus clusters by m6A RNA methylation regulators in endocrine system tumors. (A–C) Consensus clustering matrix for k = 2, and consensus clustering cumulative distribution function (CDF) for k = 2 to 9, and relative change in area under CDF curve for k = 2–9, and principal component analysis (PCA) of total RNA expression profiles in tumors data from the cancer genome atlas (TCGA) dataset; (D, E) Heatmaps, and Kaplan–Meier overall survival (OS) curve and clinicopathologic features of two clusters defined by consistent expression of the m6A regulatory genes (clusters1/2).

Then we conducted a cluster analysis on “Total”. Due to incomplete information of age and tumor metastasis in ACC clinical variables, these two variables will not be analyzed. Both the clinicopathological characteristics of cluster 1 and cluster 2 subgroups clustered by k=2 (Figure 2C), principal component analysis (PCA) was further used to compare the expression curves between cluster 1 and cluster 2 subgroups (Figure 2C). we found that the cluster 1 subgroup was associated with lower survival status (P < 0.01), higher who grade (P < 0.01), greater T (P < 0.01), and more lymph node metastasis (P < 0.01) in “Total 1” (Figure 2D). We observed cluster 1 subgroup of OS is not different from cluster 2 subgroup, only from the aspects of living condition and gender analysis, found the cluster 1 subgroup and gender (P < 0.01) in “Total 2” (Figure 2E). In general, the results showed that the expressions of m6A regulatory genes were closely related to the clinicopathological characteristics of the tumors.

The prognostic role of m6A regulators in endocrine system tumors

In order to evaluate the prognostic effect of the m6A regulatory genes in tumors, we performed univariate Cox regression analysis on the expression level of m6A regulatory genes in the dataset. The results showed that IGF2BP1, HNRNPA2B1, METTL14 and IGF2BP3 were significantly correlated with OS of ACC patients (Figure 3A). METTL14 was the protective gene for HR = 1, and IGF2BP1, HNRNPA2B1 and IGF2BP3 were the risk genes for HR = 1. For genes with prognostic values in PCPG, RBMX, HNRNPC, IGF2BP3 and HNRNPA2B1 were significantly correlated with the OS of PCPG patients (Figure 3B), and were all risk genes with HR = 1. For genes with prognostic value in THCA, YTHDF3 and FTO were significantly correlated with OS of THCA patients (Figure 3A), and were all risk genes with HR = 1. For genes with prognostic values in THYM, CBLL1, RBM15B, KIAA1429, WTAP were significantly correlated with OS of THYM patients (Figure 3B), all of which were HR = 1 protective genes.

Figure 3.

Risk signatures with three m6A RNA methylation regulators in endocrine system tumors. (A–C) The coefficients calculated by multivariate Cox regression using LASSO are shown, and Kaplan–Meier overall survival (OS) curve of patients was divided into high-risk and low-risk groups according to the risk score.

We also analyzed the prognostic value of each m6A RNA methylation regulator in “Total 1” (Figure 3C). The results showed that among the 23 tested genes, METTL14, WTAP, RBM15B, METTL16, FTO and IGF2BP1 were significantly correlated with OS (P <0.05). Among the six genes, IGF2BP1 is the risk gene of HR = 1, METTL14, WTAP, RBM15B, METTL16 and FTO are the protective genes of HR = 1. In “Total 2” (Figure 3C), WTAP and METTL16 were significantly correlated with OS (P<0.05), and they were both protective genes of HR = 1. In general, WTAP, METTL16 and IGF2BP3 had a significant influence in the occurrence and development of endocrine system tumors. The reduction expression of WTAP and METTL16 can cause poorer survival in endocrine system tumors. IGF2BP3 is danger genes.

The prognostic risk score is closely related to the clinical characteristics of endocrine system tumors

Based on the above data analysis results that m6A regulated gene expression was closely related to the clinical characteristics and prognosis of tumors, we further studied the relationship between prognostic risk score and endocrine tumors. To better determine the clinical outcome of endocrine system tumors using m6A regulatory genes, we applied the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to 23 prognostic genes in the TCGA dataset, obtained genes significantly associated with prognosis to establish risk characteristics based on the minimum criteria, and obtained the coefficients from the LASSO algorithm to calculate the risk scores of the validation dataset. To investigate the effect of genetic risk on prognosis, patients with endocrine system tumors in the dataset were separated into the low-risk group and the high-risk group according to the median risk scores. The significant differences in OS between the two groups were observed (P<0.05).

In our analysis of the tumors, we found that there were significant differences in survival status, T status and who grade between the high-risk group of ACC (Figure 4A) and the low-risk group, and significant differences in survival status between the high-risk group and the low-risk group of PCPG (Figure 4A), THCA (Figure 4A) and THYM (Figure 4A). Besides, we observed that the OS of high-risk patients in each disease dataset was much shorter than that of low-score patients.

Figure 4.

Relationship between risk score, m6A regulatory genes and clinicopathologic characteristics. (A, B) Heatmaps showed the expression levels of different m6A regulatory genes in low-risk and high-risk patients. The distribution of clinicopathological features in low-risk and high-risk groups was compared. (C) ROC curve shows the signature of risk prediction efficiency in endocrine system tumors.

We also analyzed the expression of six selected m6A RNA methylation regulators in the “Total 1” dataset in high-risk and low-risk patients. We observed significant differences in “Total 1” (Figure 4B) between the high-risk group and the low-risk group, and the high-risk group associated with lower survival status (P < 0.01), higher who grade (P < 0.01), greater T (P < 0.01), and more lymph node metastasis (P < 0.01). We then analyzed the “Total 2” dataset (Figure 4B) and failed to find the difference between the high-risk group and the low-risk group in the clinical information of gender and survival status. But OS of the high-risk patients in both “Total” datasets was much shorter than the OS of the low-score patients.

ROC curves showed that the risk score was a perfect predictor of 5-year survival in ACC patients (AUC = 80.1%) (Figure 4C), PCPG (AUC = 88.8%) (Figure 4C), THCA (AUC =72.8%) (Figure 4C) and THYM (AUC = 81.4%) (Figure 4C), “Total 1” (AUC =68.3%) (Figure 4C), and “Total 1” (AUC =70.9%) (Figure 4C).

We then performed univariate and multivariate Cox regression analyses for the dataset to determine whether the risk characteristics were independent predictors. Univariate regression analysis of ACC data showed that risk scores, who rating and T status were all correlated with OS (P <0.01, Table 1). When these factors were included in multivariate Cox regression, T status (P <0.01) and risk scores (P =0.01) remained significantly correlated with OS (Table 1). Univariate regression analysis of PCPG data showed that the risk score was correlated with OS (P =0.004, Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P =0.002) remained significantly correlated with OS (Table 1). Univariate regression analysis of THCA data showed that risk scores (P =0.014), grade (P =0.003), T status (P =0.028), and age (P <0.01) were all associated with OS (Table 1). When these factors were included in multivariate Cox regression, age (P <0.01) remained significantly correlated with OS (Table 1). A univariate regression analysis of THYM data shows that the risk score (P <0.01) is OS related. In “Total 1”, risk scores (P <0.001), who rating (P <0.001), and T status (P <0.001) were all associated with OS by univariate analysis (Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P <0.001), who grade (P =0.002), and T status (P =0.012) remained significantly correlated with OS (Table 1). In “Total 2”, the risk scores (P =0.012) was correlated with OS by univariate analysis (Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P =0.013) remained significantly correlated with OS (Table 1). In general, risk scores could be used to greatly predict the prognosis of patients with endocrine system tumors. Gender was not associated with the prognosis of endocrine tumors.

Data sources

RNA-seq transcriptome data and the corresponding clinicopathological and prognostic information were obtained for 79 ACC patients,183 PCPG patients, 510 THCA patients, 119 THYM patients from the TCGA database (https://cancergenome.nih.gov/). And TCGA normal group data were obtained for 0 ACC patients, 3 PCPG patients, 58 THCA patients, 2 THYM patients. Since the number of TCGA controls was insufficient, we combined the Genotype-Tissue Expression (GTEx) (http://www.bio-info-trainee.com/3705.html) database with a small number of TCGA normal group data to form the control group. We added GTEx_adrenal gland to the ACC control group (n=128), added GTEx_adrenal gland to the PCPG control group (n=128), added GTEx_Thyroid to the THCA control group (n=279), and added GTEx_Heart (Atrial Appendage, Left Ventricle) to the THYM control group (n=377) (Table 2).

In order to study the common characteristics of the expression value of m6A regulatory genes between endocrine system tumors, we divided the four groups of tumors for different clinical variables into two groups. “Total 1” is the dataset of ACC and THCA, mainly analyzing the clinical information of their common survival time, survival status, gender, grade, tumors range and lymph node metastasis. “Total 2” is the dataset of PCPG and THYM, mainly analyzing the only clinical information of survival time, survival status and gender (Figure 5).

Figure 5.

Clinical data processing and subsequent studies in the pan-oncogenic analysis of m6A regulatory genes in endocrine system tumors. Total 1 is the dataset of ACC and THCA, and Total 2 is the dataset of PCPG and THYM. The figure shows from TCGA and GTEx after downloading data in the database of datasets, Exclude patients with incomplete information after the rest of the sample size and the process of analysis of m6A regulatory genes.

Selection of m6A regulatory genes

We compiled a list of 23 tumor-related m6A regulatory genes [46–49] based on published literature, and then restricted the list to the genes with available RNA expression data in TCGA data. We obtained mRNA expression data of 23 m6A-related genes [50] from the TCGA database and compared the relationship between the expression of 23 m6A-related genes and the clinicopathological variables of endocrine system tumors.

Bioinformatic analysis

We used the “R” package (R v3.6.2) for analysis (Figure 5). And the Practical Extraction and Report Language (Perl) was used to accurately handle text format that require R package analysis. Firstly, we used the “limma” package to conduct gene difference analysis, used the “pheatmap” and the “vioplot” packages to visualize the expression of 23 genes, and then used the “corrplot” package to analyze the correlation of each gene expression in the tumors. Due to the expression similarity, the “ConsensusClusterPlus” package was used to cluster endocrine system tumors into different groups, and the “ggplot2” package was used to analyze the gene expression of each cluster group. Moreover, we removed patients with incomplete clinical information and used the “survival” package to determine the survival prognosis, and used the “pheatmap” package to visualize the difference expression of 23 genes between cluster 2 and cluster 1.

Next, we used the “forestplot” package to construct univariate COX proportional regression model, and evaluated different m6A expression values and used the “glmnet” package to obtain m6A regulatory genes. Moreover, we use the “survival” package and the “survivalROC” package to analyze the survival of the cluster, and used the “pheatmap” package to visualize the relationship between the expression of relative genes and the different clinical factors. Finally, we used the “forestplot” and the “survival” packages to screen out m6A regulatory genes and clinical factors with a significant influence on prognosis.

Statistical analysis

In order to determine m6A RNA methylation regulating factors in the prognostic role of endocrine tumor patients, we use the Cox univariate analysis data from the TCGA and GTEx database. From this, we selected value genes significantly associated with survival (P < 0.05), which we chose for further analysis of the function and developed a potential risk signature with the LASSO Cox regression algorithm [51]. Risk score calculation is as follows:

$\text{Risk\hspace{0.17em}score\hspace{0.17em}}={\displaystyle \sum }_{i=1}^n{\text{Coef}}_i\times {\text{\hspace{0.17em}x}}_i$

where Coef_i is the coefficient, and x_i is the expression value of each selected gene. This formula was used to calculate the risk scores for each patient.