Genetic and Transcriptional Alterations of Histone Acetylation Regulators and m6A Regulators in LUAD
Reviewing published articles, A total of 57 genes regulated by m6A methylation and histone acetylation that were included in the analysis. Metascape analysis was performed on these genes. The biological processes that significantly enriched were mainly associated with histone modification, epigenetics and others, as summarized in Figs. 1A.Based on 57 regulators, LUAD samples were well separated from normal samples. (Figs. 1B). We evaluated the somatic mutation rates of the two groups of regulators separately to identify genetic changes of the regulators in cancer. In the LUAD cohort of TCGA, 121 of 527 samples (22.96%) experienced alterations in m6A methylation regulatory genes, while 195 (37%) experienced alterations in histone acetylation regulatory genes. The highest mutation rates were found in ZC3H13 (4%) and HDAC9 (6%) respectively (Figs. 1C, D). A survey of CNV variation frequencies showed that CNV variation was prevalent in 57 regulatory genes and mostly concentrated in copy number amplification, while CNV deletion frequencies were prevalent in genes such as ZX3H13, SIRT5 and BRDT (Figs. 1E). The location of CNV alterations in regulatory factors on chromosomes is shown in Figs. 1F. We also investigated the expression of regulatory factors in LUAD and normal lung tissues and found differences in the expression of the vast majority of genes. The above analysis showed the genetic and expression changes of 57 regulators were highly heterogeneous between LUAD and normal tissues (Figs. 2A, B), suggesting that the unbalanced expression of these regulators profoundly affects the development of LUAD. Furthermore, genes do not function in isolation; therefore, we explored the correlation between m6A regulators and histone acetylation regulators. We found that the mRNAs of both are highly correlated in expression (Figs. 2C). These imply a tight crosstalk between and within histone acetylation regulators and m6A methylation regulators. m6A methylation and histone acetylation regulators have highly heterogeneous genetic and expression patterns between LUAD and normal samples, also suggesting that regulators are critical for the development of LUAD.
Identification of Two m6A Methylation Modification Patterns and Three Histone Acetylation Modification Patterns
14 m6A methylation regulators and 37 histone acetylation regulators were used for further analysis of the corresponding expression patterns based on mRNAs from 2057 lung adenocarcinoma patient samples from TCGA-LUAD, GSE14814, GSE19188, GSE31210, GSE37745, GSE50081, GSE68465 and GSE72094. To identify the expression patterns of regulatory factors, mRNA expression data from LUAD samples were classified using ConsensusClusterPlus. Using unsupervised clustering, we identified two m6A methylation patterns with different properties (Figs. 3A), including 1322 cases of cluster A, 735 cases of cluster B. We named these patterns as MAclusterA, MAclusterB. Three histone acetylation patterns with different properties, including 170 cases of cluster A, 1071 cases of cluster B, 816 cases of cluster C (Figs. 3B), were termed as ACclusterA-C. Combined with the survival information of the samples, the prognostic analysis showed that the m6A methylation subgroup had a better survival in the MAclusterB group than MAclusterA group (Figs. 3C). Among the three subgroups with histone acetylation modification patterns, AClusterA had a better probability of survival in the first 7 years, and AClusterB had the worst probability of survival (Figs. 3D). Notably, tSNE analysis showed significant differences in transcriptional profiles among m6A methylation modification patterns and the three histone acetylation modification patterns, respectively(Figs. 3E, F), which also indicated the success of our unsupervised clustering results. Finally, we grouped 2057 samples according to the grouping of m6A methylation with histone acetylation and performed survival analysis. the common sample group of MAclusterB and AClusterB had better survival, while interestingly: Samples in the AClusterC group, which had a higher survival rate and also belonged to the MAclusterA group, which had a low survival rate, instead had a lower survival rate than those in the AClusterB group, which also belonged to the MAclusterA group and had a lower survival rate. Subsequently, we obtained similar results in TCGA-LUAD and GEO meta cohort(Supplement Figs. 2).
Molecular Background of Tumors and Infiltration Characteristics of TME Cells in Different Patterns
To identify the differences in biological behavior between the different modification patterns of m6A methylation and histone acetylation, GSVA enrichment analysis was performed (Table S4, Figs. 4A, B Supplement Fig. 3A-C). MAclusterA was enriched in cancer pathways MYC/NOTCH/PI3K/TGF-Beta, and cell cycle, while MAclusterB was enriched in Antigen processing machinery/WNT/EMT/Angiogenesis/CSCs activity pathways (Figs. 4C). In the histone acetylation pattern, the highest abundance in the cell cycle/HIPPO/MYC/NOTCH/PI3K/TGF-Beta pathway was in the AClusterB group and the lowest in the AClusterA group, while the highest abundance in the WNT/EMT/Angiogenesis pathway was in the AClusterC. interestingly, in the lowest abundance in the EMT and CSCs activity pathway was in the AClusterA group rather than in the AClusterB group, which had the worst survival rate (Figs. 4D). These results reflect that the regulatory pattern is closely related to the biological behavior of lung adenocarcinoma. Highly active regulators may be a key factor in improving the malignancy of lung adenocarcinoma.
In previous studies, immune cell infiltration of TME was significantly correlated with alterations in histone acetylation and m6A methylation. Therefore, we analyzed the role played by 57 regulators in TME. The relative abundance of infiltrating TME immune cells was quantified by ssGSEA algorithm (Supplement Fig. 3D, E) and "CIBERSOFT" package (Supplement Fig. 3F, G). And the association between regulators and TME-infiltrating immune cells was analyzed (Figs. 4E, F). Evaluation of tumor purity of samples with 'ESTIMATE' package (Supplement Fig. 3H-N). T cells gamma delta, Mast cells resting, Eosinophils and Immature. dendritic.cells were positively correlated with the majority of regulatory factors. We also analyzed the TME cell infiltration in various patterns. MAclusterB differed significantly from MAclusterA. The abundance of plasmacytoid dendritic cells was higher in MAclusterB than in MAclusterA, suggesting that MAclusterB has a higher antigen-presenting function. Natural killer cells, immature dendritic cells, etc., were in higher abundance in MAclusterB. However, activated CD8 T cells, the most potent effector in the anticancer immune system(23), as well as other important tumor-killer cells and gamma delta T cells (24), were lower in MAclusterB than in MAclusterA(Figs. 4G).In the histone acetylation pattern(Figs. 4H), plasmacytoid dendritic cells were much more abundant in the ACclusterC than in the ACclusterB group. Natural killer cells were equally abundant in AC group C than in AC group B. Similar to m6A methylation, the abundance of Activated CD8 T cell, Gamma delta T cell and other cells were higher in the ACclusterB group. What is interesting is the ACclusterA group. It has comparable levels of Regulatory T cell, Macrophage, Natural killer cell, Natural killer T cell with the ACclusterC group, and also has comparable levels of Activated CD8 T cell with the ACclusterB group. B cell abundance in the ACclusterA group was much higher than that in ACclusterB and ACclusterC, and also had the lowest abundance of Activated CD4 T cell. These results suggest that the modification patterns of both m6A methylation and histone acetylation are closely related to the tumor immune microenvironment.
Based on the analytical results of our study and the experimental findings of previous authors. We selected FTO, KAT2B, METTL3 for further analysis (Figs. 4E, F). We divided the samples into a high FTO group and a low FTO group by FTO expression. k-M plot analysis showed high FTO expression group had better survival, and we subsequently examined the expression levels of m6A methylation and histone acetylation regulators in both groups and found that the expression of regulators in the high FTO group far exceeded that in low group (Supplement Fig. 4A, B). We also obtained the same results in the corresponding analysis of two genes, KAT2B, METTL3(Supplement Fig. 4C-F). Lin et al, noted that METTL3 promotes translation in human lung adenocarcinoma. We therefore collected genome-wide mRNA expression data from GSE200649 transfected with shMETTL3 versus NC in control cells. The expression differences between m6A methylation regulators and histone acetylation regulators after transfection were analyzed. Interestingly, in this dataset, the difference in expression of stably transfected shMETTL3 at the mRNA level was not significant (Table S5). These results further suggest a tight relationship between m6A methylation regulators and histone acetylation regulators.
Identification of phenotype-related genes
Due to the differences in regulatory transcriptional profiles between m6A modification Pattern as well as histone acetylation modification Pattern, we then explored the genetic variation and mechanisms of the congregational modification patterns. Identification of differential genes between subgroups by limma R package. We identified DEGs between the MAclusterA and MAclusterB groups and between the AClusterA, AClusterB and AClusterC groups (Table S6). Subsequently, 18 and 27 genes associated with survival were selected respectively. To further deepen the understanding of m6A methylation and histone acetylation, Unsupervised clustering analysis was performed with 18 and 27 genes separately and patients were classified into different gene clusters. As with previous modification patterns, two m6A-modified gene phenotypes and three histone acetylation-modified gene phenotypes were identified. (Fig. 5A, B). We named them as MA-gene-clusterA/B, AC-gene-clusterA/B/C, respectively. These results suggest that three histone acetylation modification modes and two m6A methylation modification modes do exist in LUAD. A significant proportion of genes in the m6A gene cluster pattern are highly expressed in MA-gene-clusterB. We then performed unsupervised clustering of the TCGA cohort separately from the meta-GEO cohort (Supplement Fig. 5A-F). It was found that the samples were still clustered into two groups. The same clustering was also used to analyze the histone acetylation gene cluster patterns (Fig. 5C, D), with the AC-gene-clusterA group having the highest expression of PRMT1, RHOB and TAF10, while all other DEGs had lower expression than that both in AC-gene-clusterB and AC-gene-clusterC groups. The expression levels of DEGs in the better surviving AC-gene-clusterC were higher than those in AC-gene-clusterB group. Subsequently, we analyzed the expression of m6A and histone acetylation regulators in the gene cluster. Notably, The expression of m6A methylation regulators was also much higher in the AC-gene cluster C group than in the AC-gene-cluster-B group. The expression level of m6A methylation regulators in AC-gene-clusterC group was also significantly higher than that in AC-gene-clusterB group (Fig. 5E-G). Similarly, the expression of 57 regulators in the MA-gene-cluster was higher in MA-gene cluster B than in MA-gene cluster A(Fig. 5H-J).To reveal the role of related phenotypes in TME immunomodulation, we investigated the expression of chemokines and cytokines in gene clusters. Cytokines and chemokines were extracted from the published literature(25, 26). Interestingly: TGFb/EMT pathway-related mRNAs, immune activation transcript-related mRNAs and immune checkpoint-related mRNAs were all highly expressed in MAgA.ACgC(The group consisting of MA-gene-clusterA and AC-gene-clusterC common samples) and MAgB.ACgC (The group consisting of MA-gene-clusterB and AC-gene-clusterC common samples)groups(Supplement Fig. 5G-I).The ssGSEA analysis showed that MAgA.ACgC and MAgB.ACgC were more active in the pathways of Antigen processing machinery, EMT and CSCs activity. Consistent with our previous findings (Supplement Fig. 5J)
Construction of a Digital Model for Individual Lung Adenocarcinoma Patients
Patients with lung adenocarcinoma are heterogeneous and complex. We constructed a phenotypic gene-based scoring model, called MAscore and ACscore, respectively. The variation of patients' attributes can be clearly observed by alluvial plots (Fig. 6A, Supplement Fig. 6A, B). We analyzed the correlation between the scores and immune cells, Type 2 T helper cell, Eosinophil, Immature dendritic cell and Mast cell were significantly positively correlated with MAscore, while Gamma delta T cell and CD56dim. Natural killer cells were significantly negatively correlated with MAscore. Type 2 T helper cell, Eosinophil and Immature dendritic cell are significantly positively correlated with ACscore, while CD56dim natural killer cell and Monocyte are significantly negatively correlated with ACscore(Fig. 6B,Supplement Fig. 6C).Then, we analyzed the scores in cluster and gene-cluster. The MAclusterB and MA-gene-clusterB groups with better survival had much higher MAscore than the MAclusterA and MA-gene-clusterA groups, respectively (Fig. 6C, D). While the ACscore in the subgroup of histone acetylation, ACclusterA, which had a better survival curve among ACcluster and AC-gene-cluster, AC-gene-clusterA had the lowest ACscore score among the respective groups. The ACscore of ACclusterC and AC-gene-clusterC with better survival was higher than that of ACclusterB and AC-gene-clusterB with poorer survival, respectively (Fig. 6E, F). Cox regression analysis in the TCGA cohort showed that our constructed MAscore and ACscore can be used as independent prognostic detectors. Using the optimal cut-off point to classify patients into high and low groups, survival analysis showed that patients with high scores had a better survival prognosis than patients with low scores. Combining the high and low groups of m6A methylation and histone acetylation, the survival analysis showed that the group with both MAscore and ACscore had better survival, and encouragingly, the expression levels of m6A regulators and histone acetylation regulators in this group far exceeded those of other lung adenocarcinoma patients (Fig. 6G, H).
Models predict potential therapeutic drugs
Based on the Cancer Genome Project (CGP) database, we predicted the sensitivity of 138 drugs. We selected drugs with lower IC50 values in the low-risk group in the m6A methylation and histone acetylation subgroups, and selected drugs with higher IC50 values in both the MASH.ACSH group than in the other groups, followed by analysis of the drugs most associated with MAscore, ACscore values.GW843682X, MS.275, S.Trityl.L.cysteine. VX.680 were screened (Fig. 7A, Table S7). gw843682X is a PLK inhibitor that plays a role in regulating cell mitosis, etc. MS.275 is a histone deacetylase inhibitor that potently inhibits HDAC1 and HDAC3. s.Trityl.L.cysteine is an Eg5 inhibitor that induces mitotic arrest in HL-60 cells and VX.680 is an inhibitor of pan-Aurora, which is also a regulator of mitosis. We analyzed the expression of the target genes of the four drugs in each risk group, in which AURKB and PLK1 may be effective therapeutic targets(Fig. 7B), and then analyzed the relationship between the expression of AURKB, PLK1 and LUAD survival, survival analysis showed better survival in patients with low levels of AURKB and PLK1(Fig. 7C,D), which also reflects that VX.680, GW843682X may be a promising drug for the treatment of patients in MASL.ACSH and MASL.ACCSL groups. TIDE, which is widely used and highly recommended for assessing immune responses in cancer-related studies(27, 28). Considering the MAscore, the ACscore seems to be closely related to TME. We analyzed the differences in TIDE scores across subgroups. The results showed that the TIDE scores of MASH. ACSL were much lower than those of the MASH.ACSL and MASL.ACSL groups (Fig. 7E). However, the MASH.ACSH group did not have the lowest scores, but rather the MASL.ACSH group had the lowest TIDE scores, and it is possible that the MASL.ACSH group had a higher treatment outcome from ICI treatment. Therefore, we explored the correlation between MAscore /ACscore and immunophenoscore (IPS) in LUAD for the purpose of predicting the response of immune checkpoint inhibitors (ICIs). Correlation analysis showed that MAscore, ACscore and ips_ctla4_neg_pd1_pos, and ips_ctla4_pos_pd1_pos had the highest negative correlation (Fig. 7F). The expression of immune checkpoints were similarly higher in the MASH.ACSH group than in the others. This likewise reflects the high level of m6A methylation with high immunoreactivity of histone acetylation (Fig. 7G).
Organelle division aids m6A methylation and histone acetylation to identify LUAD patients with better prognosis
In a previous study (Fig. 3I), we found that the m6AA.AcC group, which should have survived better, instead had a worse survival curve than the m6AA.AcB group. Not coincidentally, the MAgA.ACgC group in the genome had better survival than the MAgB.ACgC group instead. We extracted these patients for further analysis. Survival analysis showed no significant difference between m6AA.AcB and m6AA.AcC (Fig. 8A), while MAgA.ACgC, MAgB.ACgC survival differed significantly (Fig. 8B). Interestingly, when we analyzed the expression of m6A methylation regulators and histone acetylation regulators in these groups, we found no significant survival differences between the m6AA.AcB and m6AA.AcC groups, but instead there were considerable differences in the expression levels of the regulators (Fig. 8C, D), while there were mostly no significant differences in the expression of the regulators between MAgA.ACgC and MAgB.ACgC group (Supplement Fig. 7A, B). This result predicts the existence of one or more pathways that effectively assist in the modification pattern of m6A methylation and histone acetylation. To explore whether such a pathway exists, we performed differential analysis of the two groups, MAgA.ACgC and MAgB.ACgC, and identified a total of 235 differential genes. We then performed GO and KEGG enrichment analysis based on these 235 genes. Most of the enrichment results were associated with cell proliferation (Table S8). We selected 6 biological features for further analysis. The activity of each sample pathway was quantified by GSVA and was significantly higher in the MAgA.ACgC group than in the MAgB.ACgC group in five biological features (Fig. 8E). Based on a total of 1152 samples from the MAgA.ACgC and MAgB.ACgC groups, we found that the lower the activity of the four pathways CELL CYCLE, CHROMOSOME SEGREGATION, MITOTIC NUCLEAR DIVISION and ORGANELLE FISSION, the higher the survival rate (Supplementary Fig. 7C-H). We then introduced ORGANELLE FISSION into the grouping of m6A methylation with histone acetylation. Survival analysis showed that ORGANELLE FISSION complemented the grouping of m6A methylation and histone acetylation (Fig. 8F, Supplementary Fig. 7I, L), and the expression levels of most m6A methylation and histone acetylation regulators were also significantly different (Fig. 8E, F).