Genetic variation of m6A regulators in TC
Figure 1A showed a protein-protein interaction (PPI) network of 23 m6A regulators, including 8 writers, 13 Readers, and 2 Erasers, using different colors to distinguish the functions of the regulators through the STRING database. Further analysis in Figure 1B revealed CNV mutations in m6A regulatory factors. The copy numbers of YTHDC2, YTHDF2, RBM15, LRPPRC, METTL14 and FMR1 increased significantly, while IGFBP2, ZC3H13 and METTL16 decreased significantly. The CNV changes of the 23 m6A regulatory factors on chromosomes were shown in Figure 1C. Figure 1D showed the expression of 23 m6A methylation regulators in TC tumors and normal tissues. In the heat map, red indicates high expression and blue indicates low expression. It could be seen that 23 m6A regulators had different expression characteristics in tumor and normal tissues. In order to determine whether the above genetic variations affect the expression of m6A regulatory factor in TC patients, we investigated the mRNA expression levels of regulatory factor in normal and TC samples. The differential expression of m6A gene in TC was analyzed through TCGA database. In the boxplot 1E, red represents TC tissue and blue represents normal tissue. It can be seen in Figure 1E that most m6A genes are differentially expressed in cancer and normal tissue. Pearson correlation analysis shows the correlation between 23 m6A genes. In Figure 1F, red and blue represent positive correlation and negative correlation respectively. It can be seen from the figure that VIRMA has high correlation with YTHDF3, METTL14 and YTHDF2, and their correlation coefficients are 0.92, 0.89 and 0.82, respectively. YTHDC1 was also highly correlated with METTL14, with a correlation coefficient of 0.88.
Table 2
The expression levels of 23 m6A regulators in TCGA-THCA and normal tissues
Gene
|
Normal (median)
|
Tumor (median)
|
logFC
|
pValue
|
P symbol
|
METTL3
|
24.4677
|
20.21352
|
-0.27556
|
2.83E-06
|
***
|
METTL14
|
18.33431
|
13.11829
|
-0.48297
|
8.53E-15
|
***
|
METTL16
|
20.69898
|
18.64735
|
-0.15059
|
0.004616
|
**
|
WTAP
|
59.9992
|
47.9017
|
-0.32487
|
5.25E-09
|
***
|
VIRMA
|
19.19407
|
15.63908
|
-0.2955
|
5.78E-07
|
***
|
ZC3H13
|
37.88995
|
27.2297
|
-0.47663
|
5.57E-11
|
***
|
RBM15
|
9.226988
|
6.034931
|
-0.61252
|
4.87E-22
|
***
|
RBM15B
|
38.27544
|
38.53054
|
0.009583
|
0.825583
|
ns
|
YTHDC1
|
57.38977
|
41.05247
|
-0.48332
|
2.95E-15
|
***
|
YTHDC2
|
11.11029
|
8.937243
|
-0.314
|
4.15E-07
|
***
|
YTHDF1
|
78.59671
|
67.79571
|
-0.21327
|
4.28E-06
|
***
|
YTHDF2
|
58.7508
|
55.69421
|
-0.07708
|
0.088826
|
*
|
YTHDF3
|
42.34521
|
34.01845
|
-0.31588
|
4.54E-07
|
**
|
HNRNPC
|
153.0058
|
159.6669
|
0.061479
|
0.101208
|
ns
|
FMR1
|
33.6059
|
34.29579
|
0.029317
|
0.798682
|
ns
|
LRPPRC
|
40.77588
|
35.37785
|
-0.20487
|
0.000152
|
***
|
HNRNPA2B1
|
320.5051
|
261.1527
|
-0.29545
|
6.37E-10
|
***
|
IGFBP1
|
0.072435
|
0.221925
|
1.615309
|
0.000409
|
***
|
IGFBP2
|
9.695936
|
26.5383
|
1.452624
|
2.15E-18
|
***
|
IGFBP3
|
27.98718
|
86.16594
|
1.622352
|
1.20E-08
|
***
|
RBMX
|
112.6987
|
102.968
|
-0.13027
|
0.014141
|
*
|
FTO
|
18.36648
|
13.08476
|
-0.48919
|
1.49E-12
|
***
|
ALKBH5
|
134.8082
|
113.6899
|
-0.2458
|
6.41E-08
|
***
|
ns P > 0.05; * P < 0.05; ** P < 0.01; ***P < 0.001
Relationship between m6A related genes and prognosis of TC
In Figure 2A, m6A survival analysis was performed using TCGA database, and a total of 9 m6A genes related to the prognosis of TC patients were screened out. In Table 3, HR is the risk value, HR > 1 means the risk increases, HR < 1 means the risk decreases, and HR.95L and HR.95h are the risk fluctuation range. According to the relationship between m6A related genes and prognosis, the prognostic network was drawn in Figure 2B, which reflected the comprehensive situation of m6A regulator interaction, regulator connection and its prognostic significance for TC patients. We found that not only m6A regulators in the same functional category showed significant correlation in expression, but also between writers, readers and erasers. In Figure 2C, the protein expression of m6A molecule related to prognosis in TC tissue was analyzed by human protein Atlas (https://www.proteinatlas.org/). Immunohistochemical images showed that the expression levels of RBM15, VIRMA, YTHDC2, YTHDF2 and METTL14 proteins in tumor cells were significantly higher than those in adjacent normal gland cells in the same tissue sample, while IGFBP1 was lower in tumor cells.
Table 3
Survival analysis of 23 m6A regulators in TC
Gene
|
HR
|
HR.95L
|
HR.95H
|
pvalue
|
P symbol
|
METTL3
|
0.991493
|
0.920541
|
1.067914
|
0.179112
|
ns
|
METTL14
|
1.063449
|
0.938832
|
1.204607
|
0.035013
|
*
|
METTL16
|
1.017094
|
0.933827
|
1.107786
|
0.103822
|
ns
|
WTAP
|
1.003511
|
0.965689
|
1.042815
|
0.248033
|
ns
|
VIRMA
|
1.094501
|
0.98445
|
1.216854
|
0.017295
|
*
|
ZC3H13
|
0.996324
|
0.953178
|
1.041423
|
0.102423
|
ns
|
RBM15
|
1.288094
|
0.99871
|
1.661329
|
0.000185
|
***
|
RBM15B
|
1.010316
|
0.96265
|
1.060343
|
0.143754
|
ns
|
YTHDC1
|
1.006299
|
0.967866
|
1.046259
|
0.182387
|
ns
|
YTHDC2
|
1.155641
|
0.968981
|
1.378258
|
0.019884
|
*
|
YTHDF1
|
1.01795
|
0.997272
|
1.039057
|
0.000819
|
***
|
YTHDF2
|
1.018298
|
0.979898
|
1.058202
|
0.01272
|
*
|
YTHDF3
|
1.039552
|
0.998568
|
1.082217
|
0.00662
|
**
|
HNRNPC
|
0.997087
|
0.985051
|
1.00927
|
0.092627
|
ns
|
FMR1
|
1.003451
|
0.965708
|
1.04267
|
0.309718
|
ns
|
LRPPRC
|
0.994555
|
0.955562
|
1.035138
|
0.103564
|
ns
|
HNRNPA2B1
|
0.999429
|
0.992254
|
1.006656
|
0.195462
|
ns
|
IGFBP1
|
0.853222
|
0.223739
|
3.25374
|
0.003011
|
**
|
IGFBP2
|
1.003457
|
0.989907
|
1.017193
|
0.023581
|
*
|
IGFBP3
|
0.99757
|
0.990476
|
1.004716
|
0.198068
|
ns
|
RBMX
|
1.002753
|
0.985517
|
1.020289
|
0.24036
|
ns
|
FTO
|
1.007692
|
0.904929
|
1.122124
|
0.140034
|
ns
|
ALKBH5
|
0.995477
|
0.979609
|
1.011602
|
0.086091
|
ns
|
ns P > 0.05; * P < 0.05; ** P < 0.01; ***P < 0.001
Determination of m6A modification mode
The category discovery tool "consensesclusterplus" was used to uniformly cluster the data of TC patients based on m6A methylation regulator. Figure 3A showed that between k=2 and 9, the most stable clustering results can be obtained when k=2. PCA principal component analysis in figure 3B showed that this classification method can effectively distinguish samples. In Figure 3C, heatmap analysis was used to show the distribution of different clinical features in two m6A grouped samples. In Figure 3D, histogram and bubble diagram were drawn through GO enrichment analysis to observe the functional types mainly involved in differential genes. It can be seen in the figure that the functions mainly focus on cell cycle regulation and cell mitosis. In boxplot 3E, the abscissa represents the m6A related gene and the ordinate represents the gene expression level. It can be seen that m6A molecules were differentially expressed in genotype, and most m6A regulatory factors were highly expressed in genotype-A.
Construction of m6A gene subgroup
Although consistent clustering algorithm based on expression of m6A regulatory factor divided TC patients into two m6A-modified phenotypic patterns, the underlying genetic changes and expression disturbances in these phenotypes remain unclear. Therefore, we applied empirical Bayesian method to screen the overlapping differentially expressed genes (DEGs) between two m6A clusters, and performed unsupervised consistent cluster analysis, so as to obtain two different m6A gene characteristic subgroups, which were defined as gene Clusters-A and B respectively (Figure 4A). Figure 4B showed PCA principal component analysis, indicating that this classification method could effectively distinguish samples. Heatmap 4C was drawn according to genotype, in which the abscissa represents the sample, the ordinate represents the gene, the blue and yellow represent the gene group, and the clinical traits are above the group. It could be seen in the figure that the two component types had different expressions. Figure 4D was the survival curve of TC patients, and the results showed that there were significant differences between the two subtypes of m6A gene cluster A and B. In boxplot 4E, the abscissa represents m6A related genes and the ordinate represents gene expression level. It could be seen in the figure that m6A molecules were differentially expressed in genotype, and most m6A regulatory factors were highly expressed in genotype-B.
Construction of m6A score system
PCA was used to obtain the m6A score of each sample based on prognostic genes, and the optimal cut-off value was selected to divide the patients into high score and low score, and the survival curve was drawn. As can be seen in Figure 5A, there was a significant difference between the high score and low score of m6A, indicating that the patients in the low m6A score group had a relatively good prognosis. In Figure 5B, it can be seen that there was a significant difference in m6A score between the two groups of m6A subtypes, which has the highest score in m6A subtype A. There was no significant difference in m6A score between genotypes, but it seemed that the trend of type A was slightly higher than that of type B. Figure 5C was drawn based on m6A cluster, genotype, m6A score and survival status, and the distribution of different genotypes in other genotypes could be seen in the figure. Two waterfall plots were drawn according to the m6A score, with red indicating high m6A score and blue indicating low m6A score. By comparing the two plots, it could be seen that mutation frequency of each gene in high and low rating groups was generally different.
M6A score predicts immunotherapy benefit
The abscissa of histogram 6A is m6A high and low score groups, and the ordinate is the percentage of survival status within 5 years. It could be seen from the figure that the survival rate of patients in the low-score group was relatively higher. According to the 5-year survival, TC patients were divided into two groups: survival and death. The m6A scores of the two groups were compared, and boxplot 6B showed that the m6A scores of dead patients were higher than those of surviving patients. In figure 6C, further analysis of the relationship between m6A score and clinical features showed that the survival rate of the low-grade m6A group was better than that of the high-grade T-stage group regardless of the early stage (T1-2) or the late stage (T3-4). The correlation between m6A score and immune cells could be observed through immune correlation analysis. Figure 6D showed that m6A score was negatively correlated with most immune cells. We analyzed the relationship between m6A score and Tumor Mutation Burden (TMB). In boxplot 6E, m6A score was divided into two groups with high score and low score, and the results showed that TMB was significantly overexpressed in the low-score m6A group. We conducted survival analysis on Tumor Mutation Burden, and the results in figure 6F showed that the survival rate of patients with low TMB was significantly better than that of patients with high TMB. In combination with TMB and m6A score, the results in figure 6G showed that the survival rate of patients with low TMB and low m6A score was significantly higher than that of patients with high TMB and high m6A score. Next, through the correlation analysis of immunotherapy, we could observe which group of patients with high and low rating was more suitable for immunotherapy. In figure 6H, the abscissa is the m6A score and the ordinate is the immunotherapy score. Different expression levels of CTLA4 antibody and PD-1 antibody in the two groups were detected when they were positive (POS) and negative (NEG) respectively. It could be seen that the immunotherapy score of m6A low evaluation group was generally higher, suggesting that low m6A score group was better than high m6A score group for immunotherapy benefit (P<0.05).