2.1 Isolation of peripheral blood mononuclear cells (PBMCs)
The samples of human were collected from patients with newly diagnosed AML and healthy donors at Fujian Medical University Union Hospital, Fuzhou, China. Details on AML patient samples are shown in Table 1. The study was approved by the Committee for the Ethical Review of Research, Fujian Medical University Union Hospital and informed consent was obtained from all the patients. PBMCs were isolated through Ficoll gradient.
2.2 Processing of RNA-seq-based data
The RNA-seq-based data and clinical information of 151 AML samples were downloaded from TCGA database (https://portal.gdc.cancer.gov/) using the R package TCGAbiolinks [17], which was specifically developed for integrative analysis of GDC data. Subsequently, RNA-seq data (FPKM values) were then converted to transcripts per kilobase million (TPM) values using the limma package [18]. Normalised array data [GSE9476 (38 normal and 26 AML samples), GPL96 (HG-U133A) Affymetrix Human Genome U133A Array] and GSE23312 [GPL10107 (SMD Print_1094 Homo sapiens)] were obtained from the Gene Expression Omnibus (GEO) database. Gene expression profile matrix files of GSE9476 and GSE2331 were obtained from raw datasets using Perl [19]. Somatic mutation data were downloaded from TCGA database, and the copy number of all AML genes was obtained from the UCSC dataset.
2.3 Somatic mutation and copy number analysis of m6A regulators
To analyse the mutations and copy number variations (CNVs) of m6A regulators in AML samples, somatic mutation data were obtained from TCGA database and analysed using Perl and the R package maptools [20]. The copy number of m6A regulators was evaluated using the R packages (version 4.1.1) and Perl. The differential expression of m6A regulators in normal and AML samples were analysed using the R software. The Pearson correlation algorithm was used to analyse the correlation of m6A regulators with the occurrence and development of AML.
2.4 Unsupervised clustering of 22 m6A regulators
Only a few m6A regulators were detected using the Illumina HumanRef-8 WG-DASL v3.0 platform, 22 regulators were acquired from TCGA and GSE9476 datasets to identify different m6A modification patterns mediated by m6A regulators, including 8 writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15 and RBM15B), 1 eraser (FTO) and 13 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3 and RBMX). The established unsupervised clustering was used to determine different m6A modification patterns according to the expression of the 22 m6A regulators and classify patients for further analysis. The ConsensusClusterPlus and limma packages were used to perform the abovementioned analyses, which were repeated 1000 times to ensure the stability of the classification [21]. Subsequently, the differential expression of m6A regulators in different phenotypes were analysed using the limma [18], reshape2 and Nagpur packages [22].
2.5 Gene set variation analysis (GSVA) of the m6A regulators
Gene set variation analysis (GSVA), which is a nonparametric and unsupervised method usually used to estimate changes in pathways and biological processes [23].We used R package ‘GSVA’ to analyse differences enrichments in biological processes between m6A modification patterns. The gene set ‘C2. Cp.kegg. V6.2. Symbols was downloaded from MSigDB for GSVA. A P-value < 0.05 was considered statistically significant.
2.6 Survival analysis of m6A regulators in patients with AML
To predict the prognostic value of m6A regulators, survival analysis was performed using clinical data obtained from TCGA. And we verified some key genes in the GEO database. Kaplan–Meier (KM) survival curves of m6A regulators were plotted using the survival and survminer packages [24], and differences in survival rate were assessed using a log-rank test threshold of P-value < 0.05.
2.7 Identification of differentially expressed genes (DEGs) among distinct m6A phenotypes
The limma package was used to screen for m6A phenotype-related DEGs among the three m6A modification patterns identified from the abovementioned analyses [18], with a threshold of adjusted P-value < 0.001
2.8 Pathway and Gene Ontology enrichment
After obtaining the Entrez-ID of each DEG using the ‘org.Hs.eg.db’ R package, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of these DEGs were performed using the ‘clusterProfiler’ packages [25]. GO analysis can be divided into three categories as follows: biological processes (BPs), cellular components (CCs) and molecular functions (MFs). A P-value < 0.05 was considered significant for the Fisher’s exact test.
2.9 Tumor mutational burden (TMB) estimation and prognostic analysis
A scoring system was constructed to evaluate the m6A modification patterns of patients with AML-m6A regulator features, which was called m6Ascore. First, we screened for DEGs related to prognosis using a univariate Cox regression, and genes significantly related to prognosis were extracted for further analysis. Subsequently, principal component analysis (PCA) was performed to construct m6A regulators features. Finally, we defined the m6Ascore using a method similar to GGI [26], which is as follows:
m6Ascore =∑(PC1i + PC2i)
In this equation, i represents the expression of m6A phenotype-related genes.
Patients were divided into the high- and low-m6A-score groups based on the median score, and survival analysis was performed based on m6Ascores using the survival package. We extracted the somatic mutation data using Perl and estimated TMB values by dividing the number of somatic mutations by the total length of exons. Additionally, we performed KM survival analysis to compare OS between TMB subgroups and examined the relationship between TMB and m6Ascores via correlation analysis.
2.10 Evaluation of TME cell infiltration
Single-sample gene set enrichment analysis (ssGSEA) was used to quantify the relative abundance of infiltrating immune cells in the TME of AML using the ssGSEA package. The gene set used to label each infiltrating immune cell type of TME was obtained from a study by Charoentong [27, 28]. The gene set contains data regarding various human immune cell subtypes, including activated B cells, activated dendritic cells, T follicular helper cells, eosinophils and regulatory T cells. The enrichment fraction calculated using ssGSEA was used to represent the relative abundance of infiltrating immune cells in each sample. The Pearson correlation algorithm was used to analyse the correlation between m6A regulators and m6Ascore.
2.11 Quantitative reverse transcription polymerase chain reaction (qRT-PCR)
The TRIzol reagent (Invitrogen) was used to extract total RNA from PBMCs. cDNA was synthesised using 5×All-in-One RT MasterMix with AccuRT (abm, Canada), and Eva Green 2×qPCR MasterMix-Low ROX (abm, Canada) was used to evaluate mRNA levels. The primer sequences are listed in Table S1.
2.12 Western blot
Total proteins from PBMCs were extracted by RIPA buffer (Beyotime) containing protease and phosphatase inhibitor cocktail (Thermo). And concentration was determined using a BCA Protein Assay Kit (Thermo). Proteins were separated using SDS‑PAGE (10% gels) and transferred onto polyvinylidene difluoride (PVDF) membranes. Subsequently, the membranes were blocked with 5% BSA at room temperature for 1 h and incubated with primary antibodies overnight at 4°C. The following day, the membranes were incubated with HRP-conjugated secondary antibodies (Beyotime, #A0208, #A0216) at room temperature for 1 h. Chemiluminescence signals were visualised with a BeyoECL Star Kit (Beyotime).
2.13 Flow cytometry
All FCM studies were performed using single cell suspensions, and cells were stained in accordance with the standard protocols. Antibodies recognising CD45-eFluor506 (#69045942, Thermo), CD4-Percp/Cyanine5.5 (#30050, Biolegend), CD8-Bv786 (#563823, BD Pharmingen,), CD3-FITC, CD16/CD56-PE (#A07735, Beckman) and CD19-PerCP5 (#302227, Biolegend) were used. Samples were analyzed using a flow cytometer (BD Biosciences) and subsequent analysis was performed using FlowJo 10.1 software.
2.14 Statistical analysis
All experiments were performed in triplicate. The GraphPad Prism (version 9.2, GraphPad Software, La Jolla, CA, USA) was used for statistical analysis of all experimental data. Data were presented as mean ± SDs. The student’s t-test was performed for between-group comparisons. A P-value < 0.05 was considered statistically significant.