Thyroid cancer risk prediction model using m6A RNA methylation regulators: integrated bioinformatics analysis and histological validation

Abstract

Background: Epigenetic reprogramming has been reported to play a critical role in the progression of thyroid cancer. RNA methylation accounts for more than 60% of all RNA modifications, and N6-methyladenosine (m6A) is the most common modification of RNAs in higher organisms. The purpose of this study was to explore the related modification mode of m6A regulators construction and its evaluation on the clinical prognosis and therapeutic effect of thyroid cancer.

Methods: The levels of 23 m6A regulators in The Cancer Genome Atlas (TCGA) were analyzed. Differentially expressed genes (DEGs) and survival analysis were performed based on TCGA-THCA clinicopathological and follow-up information, and the mRNA levels of representative genes were verified using clinical thyroid cancer data. In order to detect the effects of m6A regulators and their DEGs, consensus cluster analysis was carried out, and the expression of different m6A scores in Tumor Mutation Burden (TMB) and immune double antibodies (PD-1 antibody and CTLA4 antibody) were evaluated to predict the correlation between m6A score and thyroid cancer tumor immunotherapy response.

Results: Different expression patterns of m6A regulatory factors were detected in thyroid cancer tumors and normal tissues, and several prognoses related m6A genes were obtained. Two different m6A modification patterns were determined by consensus cluster analysis. Two different subgroups were established by screening overlapping DEGs between two m6A clusters, with cluster A having the best prognosis. According to the m6A score extracted from DEGs, thyroid cancer patients can be divided into high and low score subgroups. Patients with lower m6A score have longer survival time and better clinical features. The relationship between m6A score and Tumor Mutation Burden (TMB) and its correlation with the expression of PD-1 antibody and CTLA4 antibody proved that m6A score could be used as a potential predictor of the efficacy of immunotherapy in thyroid cancer patients.

Conclusions: We screened DEGs from cluster m6A and constructed a highly predictive model with prognostic value by dividing TCGA-THCA into two different clusters and performing m6A score analysis. This study will help clarify the overall impact of m6A modification patterns on thyroid cancer progression and formulate more effective immunotherapy strategies.

Keywords: N6-methyladenosine (m6A); The Cancer Genome Atlas (TCGA); geneCluster; m6Acluster; m6Ascore; thyroid cancer.

Figure 1

Genetic variation of m6A regulators in thyroid cancer. (A) The protein-protein interaction (PPI) network of 23 m6A regulators. (B) The mutation frequency. (C) The location of the change of m6A regulator CNV on chromosome. (D) m6A waterfall plot. The right vertical coordinate represents m6A regulators, and the left vertical coordinate represents the mutation rate of m6A regulators in thyroid cancer. (E) Pearson correlation analysis shows the correlation of 23 m6A methylation modification regulators in thyroid cancer. (F) m6A methylation regulators expression in thyroid cancer. The figure shows the expression of 23 m6A regulators in thyroid cancer tumors and normal specimens. (G) The difference of mRNA expression levels of 23 m6A regulators between normal and thyroid cancer samples. The asterisk indicates statistical p value (^*P < 0.05, ^** P < 0.01, ^***P < 0.001).

Figure 2

Relationship between m6A related genes and prognosis of thyroid cancer. (A) In the survival curve, the abscissa is the survival time (years) and the ordinate is the survival rate. (B) The m6A prognosis network shows the expression and interaction of 23 m6A regulators in thyroid cancer. (C) Human protein Atlas (https://www.proteinatlas.org/) is used to analyze the protein expression of some m6A molecules related to prognosis in thyroid cancer tissue.

Figure 3

Determination of m6A modification mode. (A) According to the expression similarity of m6A RNA methylation regulator, 506 thyroid cancer patients in TCGA cohort were divided into m6A Cluster A and B. (B) PCA analysis shows that m6A related genes can distinguish the two groups of m6A genotyped samples. (C) The Heatmap shows an unsupervised cluster of 23 m6A regulators in TCGA-THCA. (D) GO enrichment analysis was performed on the difference genes screened by comparison between the two groups of m6A cluster to observe the functions of these genes. The ordinates of the histogram and bubble diagram represent the name of GO, which can be divided into three categories: BP (biological process), CC (Cell Component), and MF (Molecular function). (E) The expression of 23 m6A regulators in the two groups of m6A cluster. The asterisk indicates statistical P value (^*P < 0.05; ^**P < 0.01; ^***P < 0.001). (F) Survival analysis for RFS among two m6Aclusters. Kaplan–Meier curves and log-rank P values are shown in the graph, and the numbers at risk are shown at the bottom.

Figure 4

Construction of m6A gene subgroup. (A) By screening the overlapping DEGs between two m6A clusters and conducting unsupervised consensus cluster analysis, the samples are classified into two types according to the internal correlation. The types 1 and 2 correspond to gene cluster-A and gene cluster-B respectively. (B) PCA analysis shows that m6A related DEGs can distinguish two groups of m6A cluster samples. (C) Heatmap is drawn for m6A cluster of the two groups according to different types. The abscissa in the figure represents samples and the ordinate represents m6A related genes. (D) Kaplan-Meier curve is used to evaluate the survival of phenotypic m6A related gene characteristics, and the results show that the prognosis of genotype A is significantly better than that of genotype B (P = 0.023). (E) Expression of 23 m6A regulators in three gene clusters. The top and bottom of the box represent the quartile range of values, the lines in the box represent the median, and the colored dots represent outliers. The asterisk indicates the statistical p value (^*P < 0.05; ^**P < 0.01; ^***P < 0.001).

Figure 5

Construction of m6A score system. (A) Survival curve shows that the prognosis of thyroid cancer patients in m6A low rating group is significantly better than that in high rating group (P < 0.01). (B) Histogram shows the proportion of patients who survived or died within 5 years in the low or high m6A group. Comparison of survival and death: 98% and 2% in the low m6A score group, and 86% and 15% in the high m6A score group, respectively. (C) The abscissa in the boxplot represents the survival and death groups, and the ordinate is the m6A score. It can be seen that the m6A score in the death group is significantly higher than that in the survival group (P = 0.046). (D) There is a significant difference in m6A score between m6A cluster A and B (P = 0.026), while m6A score shows no significant difference between genotypes (P = 0.39). (E) Alluvial diagram is drawn based on m6A cluster, genotype, m6A score and patient survival status, which shows the distribution of different genotypes.

Figure 6

m6A score predicts the benefits of immunotherapy. (A) The correlation between m6A score and immune cells can be observed by immune correlation analysis. (B) In the waterfall plot, the abscissa is the sample, the ordinate is the mutation related gene, different colors represent different mutation types, and different base changes are shown below the graph. (C) Correlation analysis of m6Ascore and TMB value in thyroid cancer was performed through Spearman correlation analysis. (D) The survival curve shows that patients with low TMB had significantly better survival than those with high TMB (P < 0.001). (E) TMB and m6A score were compared in the survival curve, and the results shows that the survival rate of patients with low TMB and low m6A score is significantly higher than that of patients with high TMB and high m6A score (P < 0.001). (F) Box plot of PD-L1 expression in the low and high m6Ascore groups. The P value is shown in box plot. (G) Box plot of CTLA4 expression in the low and high m6Ascore groups. The P value is shown in box plot. (H) The expression levels of CTLA4 antibody and PD-1 antibody in high m6A score group and low m6A score group were compared.

Figure 7

Clinical evaluation of m6A score. Survival analysis of different m6Ascore groups among thyroid cancer patients with Age (A), Gender (B), tumor Stage (C), T stage (D), N stage (E), and M stage (F).