Development a m⁶A regulators characterized by the immune cell infiltration in stomach adenocarcinoma for predicting the prognosis and immunotherapy response

The collection and pretreatment of datasets and samples

In this study, the genomics data and clinical information of STAD patients were obtained from the public TCGA database (https://cancergenome.nih.gov/), containing 350 tumor samples and 32 normal samples. The selection criteria were applied as 1) Complete clinical and OS information, including gender, age, stage, and radiation therapy, were collected for the investigation. 2) Histologically confirmed STAD. The UCSC Xena (https://gdc.xenahubs.net/) obtains mutation data. Twenty-three m⁶A regulators were gained as per previous research (Supplementary Table 1). The number of nonsynonymous and synonymous mutations was stated as the tumor mutation burden. The GSE84437 (N = 433) from Gene Expression Omnibus (GEO https://www.ncbi.nlm.nih.gov/geo/) datasets were used as the validation cohort (Table 1).

The consensus clustering of 23 m⁶A regulators by consensus cluster plus

We employed the Consensus Cluster Plus R package to elucidate the biological role of m⁶A regulators in STAD categorize patients into different m⁶A isoforms [17]. “PCA” package was used to analyze gene expression among distinct m⁶A subtypes.

Gene set variation analysis (GSVA)

The biological processes among disparate m⁶A subtypes were studied using the GSVA package [18]. We used the MSigDB and IMvigor210 CoreBiologies packages to enrich well-defined biological pathways and roles [19, 20]. The Gene Ontology (GO) annotation of m⁶A-related genes was carried out using the “ClusterProfiler” package [21].

Immune-related function and immune cell infiltration estimation via ssGSEA

We estimated the proportions of the 23 tumor-infiltrating immune cell types from each sample using the “CIBERSORT” package based on the STAD gene expression matrix. Furthermore, the association between immune-related pathways and different m⁶A subtypes was discovered by ssGSEA in TCGA-STAD profiles.

Immune response analysis

Based on the immunophenoscore (IPS) of 415 STAD patients, we created the Cancer Immunome Database (TCIA) (https://www.tcia.at/home), which could provide a comprehensive analysis of immunogenomic data. Furthermore, this study employed the Tumor Immune Dysfunction and Exclusion (TIDE) method to predict immune checkpoint blockade (ICB) response and tumor immune evasion mechanisms. The ESTIMATE package was used to assess tumor purity and tumor cellularity based on expression matrixes composed of the TIME. In addition, the summation of the immune and stromal scores from each sample was defined as the tumor purity. According to the tumor sample, which contains lower tumor purity and higher immune scores, an abundance of immune cells may have infiltrated the tumor.

DEGs associated among the m⁶A phenotypes

According to the consensus clustering algorithm, patients were divided into three m⁶A clusters and identified as differentially expressed genes (DEGs). The “limma” package was used to verify DEGs among three m⁶A clusters [22]. The screening criteria were well-defined as p-value <0.01.

Evaluation of the m⁶A gene signature

The intersected DEGs were derived from the univariate Cox regression. The number of gene clusters was then determined by adopting the consensus clustering algorithm. The prognosis-related genes were discovered to construct Principal Component Analysis (PCA) and extracted the m⁶A score [23, 24]. The advantage of this algorithm is that it focuses mainly on negatively (or positively) correlated genes. The m⁶A score formula is as follows:

$${\text{m}}^6\text{A}\text{score}=\Sigma (\text{PC}{1}_{\text{i}})+\Sigma (\text{PC}{2}_{\text{i}})$$

Lastly, the correlation between m⁶A score and TMB was performed via Spearman’s approach according to survival curve and synonymous and nonsynonymous mutation counts.

Biological pathways with m⁶A score

There is evidence that a panel of signatures associated with distinct biological pathways, including: 1) epithelial mesenchymal transition (EMT); 2) immune-checkpoint; 3) pan-fibroblast TGFb response signature; 4) CD8 T-effector signature; 5) homologous recombination; 6) Fanconi anemia pathway; 7) WNT target; 8) base excision repair; 9) mismatch repair; 10) DNA damage repair; 11) DNA replication; 12) nucleotide excision repair; 13) cell cycle regulation; 14) antigen processing; 15) cell cycle and 16) FGFR3-related genes [20, 25].

The immune-checkpoint cohorts analysis

Three independent anti-PD-L1 cohorts that could apply to genomic and clinical information were acquired from the IMvigor210 cohort [20] and GSE78220 cohort [26] to calculate the predictive value of the m⁶A score for immunotherapy.

Evaluation of the sensitivity of chemotherapeutic drugs

Genomics of Drug Sensitivity in Cancer (GDSC) is the largest public pharmacogenomics database that could predict sensitivity between high and low m⁶A score subtypes [27]. The pRRophetic package could process the half-maximal inhibitory concentration (IC₅₀) via the ridge regression model [28].

Immunohistochemistry (IHC)

IHC experiment detected IRGs protein levels in STAD and corresponding normal tissues. First, each sample was subjected to 10% formalin fixation, paraffin embedding, and processing up to 4-μm consecutive sections. After that, the sections were treated with methanol, followed by BSA incubation and primary antibody staining. Then, we stained the samples with a secondary antibody after washing with PBS. Finally, each section was observed and photographed using a microscope. The primary antibodies were obtained from Abcam.

Statistical analyses

R software (version 3.6.3) was applied for this study. Wilcoxon test was computed using the difference comparison of the DEGs between high and low m⁶A score groups. Association between tumor-infiltrating immune cells and m⁶A score was applied for Spearman’s correlation. Kaplan-Meier methods were visualized for the difference in overall survival (OS) between the high-risk and low-groups. In addition, the screening criterion was p-value <0.01.

Data availability statement

Publicly available datasets were analyzed in this study. The datasets generated during the current study are available in the Cancer Genome Atlas (TCGA) public dataset (https://portal.gdc.cancer.gov/) and the Gene-Expression Omnibus (GEO) public dataset (https://www.ncbi.nlm.nih.gov/geo/).

Genetic expression of 23 m⁶A regulators in STAD

Figure 1 shows the flow chart of this study. The TCGA dataset identified 23 m⁶A regulators (eleven “readers,” eight “writers,” and two “erasers”), which were summarized the biological functions and processes via the Metascape database (Figure 2A). Figure 2B shows the prevalent deletions of m⁶A regulators in copy number, while VIRMA, YTHDF1, and FMR1 featured an overall frequency of CNV amplification. Subsequently, as per Figure 2C, the location of CNV of all m⁶A regulators was listed in the chromosomes. Figure 2D shows that ZC3H13 showed the highest mutation frequency (8%), followed by RBM15B and YTHDC1, whereas VIRMA, METTL16, ALKBH5, FTO, and METTL3 did not display any mutations. Based on our analysis of the most mutated ZC3H13, nine of the other 22 m⁶a regulators are of interest (Supplementary Figure 1). Further investigation demonstrated that the expression of most m⁶a genes was significantly up-regulated in tumor samples except IGFBP2 (Figure 2E). m⁶A regulators (like VIRMA, YTHDF1, FMR1, and ZC3H13) with amplificated CNV showed higher expression than normal tissues. (Figure 2B, 2E). Following the obtained findings, we could demonstrate that m⁶A regulators exhibited significant transcriptome-altering landscapes and heterogeneity in genomics between STAD and normal samples.

Figure 1. Research flow chart.

Figure 2. Genetic topography of the 23 m⁶A regulators in STAD. (A) Metascape network of 23 functionally enriched m⁶A regulators. Annotations denoted by different circles varied. (B) CNV map for 23 m⁶A regulators, where the column stood for the alteration frequency, and the green and pink dots, respectively, indicated the deletion and amplification of CNV. (C) CNV alteration sites for the cellular m⁶A regulators. (D) Among 433 patients, varying genetic alterations were noted in 94 patients, such as missense, nonsense, and splice-site mutations. (E) The different expression levels of 23 m⁶A regulators between normal and STAD.

Identification of m⁶A subgroups mediated via 23 m⁶A regulators

The TCGA-STAD dataset with available clinical and survival information was enrolled into the training cohort. By applying the R package “Consensus ClusterPlus”, 350 STAD patients were stratified into three distinct following the expression of 23 m⁶A genes (Figure 3A, 3B, and Supplementary Figure 2). Variable selection was performed by LASSO Cox regression and the parameter λ indicated that the most suitable model to predict survival included RBM15, IGFBP1, RBMX, FTO, and ALKBH5 with coefficients of -0.306, 0.125, 0.042, 0.543, and -0.293, respectively (Supplementary Figure 3). The regulator network comprehensively represented the prognostic significance of 23 m⁶A regulators and their whole interactions (Figure 3C, and Supplementary Figure 4). FTO was found to be a risk factor and a favorable factor for ALKBHS in the eraser gene, and all other writers, except for IGFBP2 and IGFBP3 genes, were also favorable. In addition, we found that the 23 m⁶A regulators are positively correlated except for IGFBP2 and IGFBP3. Based on these results, the crosstalk among the 23 m⁶A regulators might play a key role in the pathogenesis and progression of individual tumors and in the formation of distinct m⁶A modifications. Thus, unsupervised clustering was adopted to classify samples into distinct m⁶A clusters. Furthermore, we could distinguish one m⁶A cluster from others completely in line with PCA (Figure 3D). Accordingly, three distinct m⁶A clusters were ultimately detected in the TCGA dataset, containing 161 in m⁶A cluster A, 133 in m⁶A cluster B, and 143 in m⁶A cluster C (Figure 3E). Among the above clusters, m⁶A cluster A, m⁶A cluster B, and m⁶A cluster C, m⁶A cluster A showed poor prognosis in the TCGA cohort, whereas patients in the m⁶A cluster C had an advantage in overall survival (p = 0.003). Similar results could also be found in the validation cohort (GEO cohort (Figure 3F).

Figure 3. Patterns of m⁶A methylation modification. (A) CDF (cumulative distribution function; k = 2–9) in the right panel. (B) Depending on the consensus clustering matrix (k = 3), the patients with STAD were classified into three clusters. (C) Interactions among 22 m⁶A regulators in STAD. Circles in varying colors were used to represent differing RNA modifications, where red indicated Erasers, orange indicated readers, and gray indicated writers. Besides, green and purple circles, respectively, referred to the favorable and risky factors. (D) PCA (principal component analysis)-based map depicting prominent differences among the three m⁶A clusters. (E) Kaplan–Meier OS (overall survival) plots of 3 m⁶A clusters for the TCGA cohort (p = 0.003). (F) Kaplan–Meier OS plots of three m⁶A clusters for the GEO cohort (p < 0.001). OS was worse among those in the m⁶A cluster C as compared to the other 2 clusters.

The distinct immune landscapes of TIME in m⁶A clusters

We performed the heatmap and visualized the expression of 23 m⁶A regulators (Figure 4A) to investigate the correlation of these m⁶A clusters with clinical features. Most of the m⁶A regulators were up-expressed in cluster B, while IGFBP1 and IGFBP2 were down-regulated in cluster A. To our surprise, most people over 65 years significantly increased in cluster B. TIME and CIBERSORT package were applied to investigate the subsets of immune cells and feature the immune cell infiltration according to the expression file. Anti-tumor lymphocyte cells, like NK cells and activated CD8+ T cells, were primarily engaged in the m⁶A cluster A (Figure 4B). Cluster C was enriched with 2/17 T helper cells and gamma delta T cells. We demonstrated that m⁶A cluster A represented the highest immune scores, followed by m⁶A cluster C and m⁶A cluster B (Figure 4C). Accordingly, m⁶A cluster A displayed a higher tumor purity than m⁶A cluster C and m⁶A cluster B, suggesting that more stromal cells and immune cells surround tumors in m⁶A cluster C and B (Figure 4D). Previous studies had suggested a tumor stroma with multiple immune cells could retain these cells rather than the parenchyma of an immune-excluded phenotype.

Figure 4. TIME properties in the three m⁶A clusters. (A) Thermogram depicting the consensus clustering outcome for the TCGA cohort, where age, gender, T&N stage, radiation, survival, and stage constituted the clinical traits. (B) CIBERSORT-based immunocyte infiltration in the three m⁶A clusters. *p < 0.05;**p < 0.01; ***p < 0.001. ns stood for no significance. (C–E) Assessment results of the (C) immune scores, (D) estimate scores, as well as (E) stromal scores across three m⁶A clusters.

The m⁶A-related DEGs in STAD

1018 DEGs were overlapped via the Venn diagram to identify the biological behaviors of these m⁶A clusters (e.g., expression perturbations and genetic alterations), and then, we examined transcriptional expression changes associated with three m⁶A clusters in STAD (Figure 5A). According to GO functional annotation, BPs associated with ribonucleoprotein complex biogenesis and RNA localization were significant bio-functions (Figure 5B). As per Figure 5C, three m⁶A gene clusters displayed distinct clinicopathological characteristics. Further, patients of late stages were mostly clustered into gene cluster C. Furthermore, all those three m⁶A gene clusters showed significantly different prognoses among the STAD samples based on survival analysis. Generally, samples of m⁶A cluster B had dismal prognostic outcomes, while those of m⁶A cluster C had superior prognostic outcomes (Figure 5D).

Figure 5. The construction of m⁶A gene clusters and relevant functional annotations. (A) Venn plot depicting 1,028 DEGs (differentially expressed genes) among three m⁶A clusters. (B) GO enrichment findings of 1,028 intersecting genes. (C) Intersect gene-based consensus clustering result of patients into three separate gene clusters. (D) Kaplan–Meier plots for the three m⁶A gene clusters (p < 0.001).

Prognostic signature establishment and clinical feature exploration

According to those findings above, m⁶A regulators have essential effects on modulating TIME and prognostic outcome. However, the above analysis is based on the general population alone, whereas the heterogeneous and complicated m⁶A regulators are not explained separately. Based on those discovered m⁶A genes, this work established one scoring system referred to as m⁶A score to quantify scores of different individuals using those discovered m⁶A genes.

The alluvial diagram represents quantitative alterations in STAD samples (Figure 6A). These findings indicated that m⁶A clusters A/C showed increased m⁶A scores, while cluster B showed decreased m⁶A scores. In addition, we performed Spearman’s correlation analysis to illustrate m⁶A regulator patterns. Based on the heatmap, the m⁶A score showed immature B cell and activated CD4 T cell (Figure 6B). Next, this study evaluated the m⁶A score’s significance as a prognosis predictor for patients. We divided patients as low-or high m⁶A score subgroups as per the threshold. As expected, high m⁶A score patients were related to significantly poorer outcomes (Figure 6C). In this study, we also analyzed the relationship between m⁶A score and PD-L1 expression; as a result, high m⁶A score patients had increased PD-L1 expression compared to low m⁶A score patients (Figure 6D). Increasing evidence has illustrated that TMB is related to immunotherapeutic response. Thus, this work analyzed TMB distribution between high and low m⁶A score patients. According to our results, low m⁶A score patients were related to a decreased TMB frequency (Figure 6E). As shown in Figure 6G, the m⁶A score showed a marked positive correlation with TMB (R = 0.35, p = 6.8e-12). Further, we also examined somatic mutation gene distributions between both subgroups. According to Figure 6F, patients with high m⁶A scores had more somatic mutations than those with low scores.

Figure 6. The m⁶A score construction and relevant genetic trait assessment. (A) Alluvial chart of m⁶A clusters regarding gene cluster and score of m⁶A, as well as patient survival. (B) Spearman correlations of m⁶A score with immunocytes. (C) Survival findings for the TCGA cohort patients marking high and low m⁶A scores. (D) PD-L1 level comparison between patients marking high and low m⁶A scores. (E) TMB (tumor mutation burden) distribution comparison between patients marking high and low m⁶A scores. (F) Mutational waterfall plot for the TCGA cohort patients marking high (left panel) and low (right panel) m⁶A scores. All patients were represented by individual columns. (G) Diagram illustrating significant positive association of m⁶A score with TMB (R = 0.35, p = 6.8–12). (H) Kaplan–Meier plots for patients exhibiting high (H) and low (L) TMBs. (P< 0.001). (I) Kaplan–Meier plots for the TCGA cohort patients stratified by both m⁶A score and TMB. H, high; L, Low; TMB, tumor mutation load (P = 0.003).

Additionally, cases with increased TMB frequency had survival benefits (p<0.001; Figure 6H), whereas those with high m⁶A scores did not exhibit any difference in survival advantage between high and low TMB frequency (Figure 6I). Therefore, cases with an increased TMB frequency possibly gained more survival benefits from radiotherapy than those with a decreased frequency.

This work also examined the m⁶A score’s prognostic significance under diverse clinical factors. As a result, low m⁶A score cases were associated with superior prognosis to low m⁶A score patients from the T-stage subgroup (Supplementary Figure 5). Besides, low m⁶A score cases possibly gained more advantages from signatures than high m⁶A score cases.

m⁶A score’s effect on estimating immunotherapeutic responses

This study examined differential expression of immune checkpoints (ICPs) to predict immunotherapeutic responses among STAD patients. IPS for the 140 TCIA- derived PAAD cases were identified as the favorable factor that predicted anti-cytotoxic T lymphocyte antigen-4 (CTLA-4) and anti-programmed cell death protein 1 (anti-PD-1) antibody responses, and it was significantly different between the two groups that they were both positive (Figure 7A–7D). Additionally, the other two immunotherapy cohorts, IMvigor210 cohort and GSE78220, were used to analyze the role of m⁶Acore in predicting anti-PD-L1 therapeutic response in patients. As for IMvigor210 and GSE78220 cohorts, patient survival did not consistently respond between low and high m⁶A score groups (Supplementary Figure 6). In conclusion, the results of this study suggest a close relationship between m⁶A regulators and TIME, which is responsible for immunotherapeutic responses.

Figure 7. The m⁶A score predicts immunotherapeutic benefits. (A–D) TCIA database-based correlations of m⁶A score with IPS among the PAAD patients: (A) CTLA4(−) PD1(−) (B) CTLA4(−) PD1(+) (C) CTLA4(+) PD1(−) (D) CTLA4(+) PD1(+). (E) Distribution of m⁶A score in different MSI groups. (F) The ratio of MSI (microsatellite instability) in STAD samples with MSS and MSI-L and MSI-H.

External validation of five key prognostic m⁶A RNA modulators

IHC staining of five paired STAD and adjacent non-tumorous tissues was used to validate the protein expression of the five prognosis genes (Figure 8). The compared with adjacent non-tumorous tissues, we found a higher level of RBM15, IGFBP1, RBMX, FTO, and ALKBH5 in all five STAD tissues, RBM15, FTO and ALKBH5 were overexpressed in STAD tumor tissues, while the expression of IGFBP1, RBMX was not changed in STAD tumor tissues.

Figure 8. IHC analysis of five m⁶A modulators in stomach adenocarcinoma tissues and adjacent normal tissues.