Abstract
Background/Aim: Prediction of response to azacitidine (AZA) treatment is an important challenge in hematooncology. In addition to protein coding genes (PCGs), AZA efficiency is influenced by various noncoding RNAs (ncRNAs), including long ncRNAs (lncRNAs), circular RNAs (circRNAs), and transposable elements (TEs). Materials and Methods: RNA sequencing was performed in patients with myelodysplastic syndromes or acute myeloid leukemia before AZA treatment to assess contribution of ncRNAs to AZA mechanisms and propose novel disease prediction biomarkers. Results: Our analyses showed that lncRNAs had the strongest predictive potential. The combined set of the best predictors included 14 lncRNAs, and only four PCGs, one circRNA, and no TEs. Epigenetic regulation and recombinational repair were suggested as crucial for AZA response, and network modeling defined three deregulated lncRNAs (CTC-482H14.5, RP11-419K12.2, and RP11-736I24.4) associated with these processes. Conclusion: The expression of various ncRNAs can influence the effect of AZA and new ncRNA-based predictive biomarkers can be defined.
- Noncoding RNAs
- circular RNAs
- transposable elements
- myelodysplastic syndrome
- acute myeloid leukemia
- azacytidine
Myelodysplastic syndromes (MDS) are a heterogenous group of clonal malignancies characterized by ineffective hematopoiesis resulting in peripheral cytopenia and hypercellular bone marrow dysplasia. MDS patients have an increased tendency to transform to acute myeloid leukemia (AML). When the disease progresses to leukemia, prognosis is clearly unfavorable, with an estimated survival of less than one year (1). The prognosis of MDS patients is assessed based on the Revised International Prognostic Scoring System (IPSS-R) and depends on the number and severity of cytopenia, percentage of bone marrow (BM) blasts and cytogenetics (2).
Azacitidine (AZA) is a hypomethylating agent (HMA) that is currently used as a standard care for patients with higher-risk MDS or AML with 20-30% BM blasts when not eligible for intensive chemotherapy and allogeneic transplantation. AZA is a cytidine analogue that, at low doses, inhibits DNA methyltransferases, resulting in DNA hypomethylation. At high doses, AZA is directly cytotoxic to abnormal BM hematopoietic cells because of its incorporation into DNA and RNA (3). Although AZA is extensively used, 30-40% of patients fail to respond or relapse after treatment (4). The International Working Group (IWG) proposed standardized response criteria for evaluating clinically significant responses in MDS (5). Stratification of patients who will benefit from those who will not respond to AZA is essential. In addition to the high costs associated with the treatment, an important issue lies in prolonged exposure to AZA as current recommendations are for a minimum of 6 months of treatment before deeming it a failure. Therefore, proper prediction biomarkers that could potentially prevent prolonged unnecessary adverse effects in AZA-nonresponsive patients are needed.
There have been many efforts to characterize clinical, biological, or molecular parameters predictive of AZA treatment outcome (6-12). Because DNA hypermethylation can contribute to tumorigenesis by silencing tumor suppressor genes, it was supposed that the methylation pattern may be predictive of AZA response. Although both global and gene-specific hypermethylation have been performed in MDS (13, 14), little relation between the degree of demethylation following hypomethylating treatment and hematologic response has been evidenced (15). Furthermore, due to the number of somatic mutations found in regulators of DNA methylation, it has been hypothesized that such mutations may impact the response to AZA. Some studies have suggested the utility of DNMT3A (7, 8) and TET2 (9, 10) mutations as predictive markers for AZA response, yet these datasets have not been largely confirmed. To date, no clinical, cytogenetic or molecular markers of AZA treatment outcome to support clinical decisions have been validated (16).
Gene expression profiling has repeatedly been performed in MDS and has identified many protein coding genes (PCGs) and signaling pathways associated with disease pathophysiology (17) and treatment (11). To provide biological insights into the contribution of noncoding RNAs (ncRNAs) to the AZA response or resistance and to propose novel appropriate molecular markers able to predict the response, whole transcriptome RNA sequencing (RNA-seq) was performed in CD34+ BM cells in patients with higher-risk MDS and AML with myelodysplasia-related changes (AML-MRC) before AZA therapy. A more detailed analysis of PCG data was not the major aim of this study, as they have already been addressed in a previous study (18). This study focused on the noncoding transcriptome and paid special attention to the examination of various classes of ncRNAs as predictive markers of the AZA response and compared their power to predict the response with the predictive power of PCGs. In addition to regular analysis of RNA-seq data that standardly examines the expression of PCGs and several classes of noncoding transcripts (mainly long noncoding RNAs, lncRNAs), the NGS data with regard to circular RNAs (circRNAs) and transposable elements (TEs) was also analysed. To the authors’ knowledge, this is the first study addressing and integrating transcriptional data of these molecules in MDS.
Materials and Methods
Patients. The study included 26 patients with a diagnosis of MDS or AML-MRC and 9 healthy controls. BM samples were obtained from the patients before the first administration of AZA during routine clinical assessment at the Institute of Hematology and Blood Transfusion and the First Department of Internal Medicine of the General Faculty Hospital in Prague. The study included only patients with no known history of previous malignancy, chemotherapy or radiation therapy, and none of the patients received hematopoietic stem cell (HSC) transplantation. The patient’s diagnoses were assessed based on the standard WHO 2016 classification criteria (19), and all the patients were classified according to the IPSS-R categories (2) at the time of sample collection. AZA was administered at 75 mg/m2 subcutaneously for 7 days and repeated every 4 weeks. The response status of the patients was evaluated according to IWG criteria (5). Patients with complete remission (CR), partial remission (PR), marrow CR with hematological improvement (mCR with HI), or stable disease with HI (SD with HI) were considered responders (RESPs, N=11). Patients with SD without HI and with progressive disease (PD) accounted for nonresponders (nRESPs, N=15). The control group contained nine hematologically healthy donors (four males and five females; mean age 41 years; range=27-69 years). Written informed consent was obtained from all tested subjects and the study was approved by the Scientific Board and the Ethics Committee of the Institute of Hematology and Blood Transfusion and the General Faculty Hospital, and performed in accordance with the ethical standards of the Declaration of Helsinki. The detailed clinical and laboratory characteristics of the cohort, including the patient classification into subgroups, IPSS-R categories, AZA response, blood counts, cytogenetics, and mutational status, are summarized in Table I.
Patient characteristics.
Cell separation and nucleic acid extraction. Mononuclear cells (MNCs) were separated from BM aspirates by Ficoll-Paque density gradient centrifugation (GE Health care, Munich, Germany). CD34+ cells were isolated from MNCs using magnetic cell separation system from Miltenyi Biotec (Bergisch Gladbach, Germany). The acid-guanidine-phenol-chloroform method was used to extract total RNA and DNase I (Qiagen, Hilden, Germany) treatment was applied to prevent genomic DNA contamination. Each RNA was quantified using Qubit 3 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was assessed using the TapeStation 4200 HS RNA kit (Agilent Technologies, Santa Clara, CA, USA).
RNA sequencing. For RNA-Seq, a rRNA depletion-based approach (RiboCop rRNA Depletion Kit v1.2, Lexogen, Vienna, Austria) was used followed by library construction with NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). Each library was marked with unique indexed primers. Library quality and size were assessed using a TapeStation 4200 HS D1000 kit and quantified with a Qubit 3 Fluorimeter. Libraries were subsequently equimolarly pooled and 2 x 125 bp paired-end sequenced with a HiSeq SBS V4 kit on a HiSeq 2500 instrument (Illumina, San Diego, CA, USA).
Bioinformatical data processing. FASTQ files were subjected to initial quality control by FastQC (20), and adaptor trimming and low-quality sequence removal were performed by Trimmomatic (21). For identification of “regular” transcripts (i.e., PCGs and lncRNAs), the reads were aligned to the human genome version hg38 using STAR (22). Genomic features in mapped data were then counted using the StringTie tool (23). For circRNA identification, the reads were mapped to the human genome version hg19 using BWA (24) because genomic coordinates of circRNAs in the circBase database are in hg19. The mapped data were further processed by CIRI2 software (25). The SalmonTE tool (26) was used for the detection of TEs from raw FASTQ files and for further processing of TE data. All other subsequent statistical analyses were performed in the R statistical environment with appropriate Bioconductor packages. The pcaMethods (27) was applied for sample clustering. Differential expression analysis (DEA) was performed with edgeR (28) and DESeq2 (29). For DEA, only those transcripts that had at least 5 counts-per-million (CPM) in at least this number of libraries, which corresponded to half of the smallest sample group, were retained. In standard settings, only the transcripts with |logFC|>1 (absolute value of binary logarithm of fold change) and FDR<0.05 (false discovery rate; p-value adjusted for multiple testing) were considered significantly deregulated.
Three different supervised machine learning methods were used for the identification of a subset of transcripts applicable as the best prognostic classifiers of patient response. The methods that can cope with a large number of features and are not prone to overfitting were selected. All the methods had been applied to RNA-seq data classification. In particular, the following methods were employed: support-vector machines with recursive feature elimination (SVM-RFE) (30), lasso and elastic-net regularized generalized linear models (GLMNET) (31), and variance modeling at the observational level (voom) extensions of the nearest shrunken centroids (voomNSC) (32). Considering unbalanced classes (the number of RESP and nRESP patients differs) the area under the ROC curve (AUC) was used for the evaluation of classifier performance. To assess the generalization ability of the classifiers, the evaluation was performed with the aid of 10 times repeated 10-fold cross-validation.
To define a formula calculating the prediction score of a patient based on the weight of the individual selected transcripts, GLMNET was applied to fit generalized linear models via penalized maximum likelihood. When the response was binary (RESP from nRESP patients in our case), it fit a logistic regression model:
Prediction score=logPr(RESP|X=x)Pr(nRESP|X=x)=β0+βTx
where β0 stands for intercept, β for regression coefficients of individual transcripts (T), and x for an expression vector of a patient. The prediction score is equal to a logarithm of patient response odds; thus, a score >0 predicts future AZA response, and a score <0 suggests AZA failure.
Pathway analysis was conducted using the Gene Set Enrichment Analysis (GSEA) method (33). As a reference, the c5 (c5.all.v7.0.symbols.gmt [Gene Ontology, GO]) gene set from the Molecular Signatures Database was utilized. Computational modeling was applied for the construction of coexpressional networks as previously described (34). Briefly, ncRNA network analysis was implemented as follows. The GO terms significantly enriched in deregulated PCGs were identified. Then, ncRNAs (from all three classes of noncoding transcripts: lncRNAs, circRNAs, and TEs) were found whose expression profiles were most correlated with the PCGs associated with the given GO term. Specifically, ncRNAs were scored with preranked GSEA, and the enrichment score statistic was calculated based on absolute Spearman correlation between the ncRNA under consideration and the set of PCGs associated with the GO term. Eventually, a network was displayed that included nodes corresponding to the highest scoring ncRNAs and the core PCGs representing the given GO term. The core PCGs were those that most often contributed to the enrichment signal for the highest scoring noncoding transcripts. The edges between the network nodes represented the upper quartile correlations (bold line) and the above median correlations (regular line).
Results
Characterization of transcriptomic data. Deep RNA-seq was performed on 35 samples obtained from CD34+ BM cells isolated from 26 patients (clinical details are summarized in Table I) and 9 healthy controls and, on average, 75,8 M (range=56.8-99.9 M) of raw pair-end reads per sample were detected. The output data were aligned to hg38 and, on average, 87% reads per sample uniquely mapped to the genome. Overall, 19,873 PCGs and 38,343 ncRNAs were detected with at least one read across all the samples. Using the CIRI2 software, the data were further realigned to hg19 and mapped the sequences to the circBase database containing 58,216 circRNAs. On average, 3.7% of reads per sample mapped to circRNAs and 49,280 circRNAs were identified across all the samples. Finally, data analysis with the SalmonTE tool showed that 5.7% of raw reads mapped to TE sequences. Of the 687 elements utilized by the tool, 627 TEs were detected across all samples.
Only the transcripts that had ≥5 CPMs at least in any 5 samples were considered as expressed (12,330 PCGs, 2,059 ncRNAs, 6,015 circRNAs, and 554 TEs) and were therefore used for subsequent data analyses. The majority of the ncRNAs were composed of antisense RNAs (20%, N=418) or long intergenic noncoding RNAs (lincRNAs, 15%, N=375). Therefore, for the purpose of simplification, the “lncRNA” term was used for this category of transcripts within the paper. Regarding circRNAs, the majority of them were of exonic origin (94%, N=5,628) whereas only a minority of them were intronic (4%, N=233) or intergenic (2.5%, N=145). Within TE data analysis, the expression of 410 LTR (68%; e.g., 235 ERV1, 41 ERV2, and 133 ERV3) and 196 non-LTR (32%; e.g., 166 L1 and 70 SINE) retrotransposons was detected.
Disease-specific transcriptome. In the initial step, the patient-specific transcriptome must be characterized. Although this was not a major aim of this study, it was necessary to explore the general molecular background in disease pathogenesis to better understand AZA-related changes. Additionally, there was almost no information about circRNA and TE transcription in MDS. This data was therefore of particular interest.
Four PCA clustering analyses (i.e., separate analyses for PCGs, lncRNAs, circRNAs, and TEs) were performed to compare the contribution of different transcript groups to the disease pathogenesis. Based on PCA, the degree of variability was evaluated in gene expression among the samples. All four PCA analyses showed that healthy controls grouped into a compact homogenous sample cluster whereas expression in patient samples was more heterogeneous and future response status cannot be predicted based on any type of genome-wide transcriptomic data. Interestingly, patients were most clearly differentiated from healthy controls based on the levels of PCGs and lncRNAs, yet information about circRNA or TE expression was not sufficient for stratification of patients and controls (Figure 1).
Principal component analysis (PCA). Clustering of samples was performed separately with four different types of RNA-seq output data. PCGs: Protein coding genes, lncRNAs: long noncoding RNAs, circRNAs: circular RNAs; TEs: transposable elements; green – healthy controls; red – responders, yellow – patients with stable disease without hematological improvement, blue – patients with progressive disease.
In the next step, DEA was performed between samples obtained from patients and those from normal individuals. All patient samples were included in these analyses, irrespective of their future AZA response status. Table II summarizes the most deregulated transcripts of each RNA category in MDS/AML-MRC patients compared to healthy controls.
The top deregulated transcripts between patients with myelodysplastic syndromes/acute myeloid leukemia with myelodysplasia-related changes and healthy controls from each RNA category defined by differential expression analysis.
Regarding PCGs, DEA revealed 457 significantly deregulated transcripts (|logFC|>1, FDR <0.05). In patients, 331 PCGs were up-regulated (based on logFC, the 12 most up-regulated genes were COL4A5, DSC2, CASS4, SLITRK5, LGALSL, CD5L, TPSB2, LOXL4, HBA1, MFAP3L, HBA2, and CANCB4) and 126 PCGs were down-regulated (the 12 most down-regulated genes were ARPP21, PAX5, EBF1, AKAP12, POU2AF1, RAG1, MME, AREG, CD79A, NEIL1, CD24, and LEF1). These differences were comparable to those published in earlier gene expression profiling studies (34-37). For example, down-regulation of LEF1 significantly correlated with clinical outcome and might be promising for assessing prognosis in MDS (37). Importantly, the similarity of the PCG data presented in this study with the previous findings proved the validity of our RNA-seq outputs. GSEA further revealed that the PCG deregulation in the patients was predominantly associated with epigenetic regulation of gene expression (normalized enrichment score (NES)=2.02, p=0.002) and chromatin silencing (NES=1.97, p=0.002). Core enrichment genes for these two pathways included multiple histone proteins, especially those from the HIST1H and HIST2H families. Other significantly deregulated pathways included DNA ligation involved in DNA repair (NES=1.82; p<0.001; core enrichment genes: HMGB1/2, LIG1/3, PARP1/2, XRCC1, etc.) and regulation of alternative mRNA slicing via spliceosomes (NES=1.82; p=0.002; core enrichment genes: splicing factors such as SRSF2-4/6/9-10 and SF3A1, etc.). All these pathways were suppressed in patient samples compared to healthy controls, suggesting that transcriptional regulation and reparatory processes were disrupted at multiple levels in myelodysplasia. GSEA output is summarized in Figure 2. Because PCG profiling has repeatedly been performed in MDS and many signaling pathways associated with disease pathophysiology have been described (34-37), the expression of various classes of noncoding transcripts was a focus of this study.
Pathway analysis of differentially expressed genes in patients with myelodysplastic syndromes/acute myeloid leukemia with myelodysplasia-related changes compared to healthy controls. A plot of gene set enrichment analysis results shows p-values (logarithmically scaled) and normalized enrichment scores (NES) of the top 10 most enriched pathways in each group of samples sorted according to descended NES values. The four most significant gene sets (marked by asterisks; negative regulation of gene expression - epigenetic; chromatin silencing; DNA ligation involved in DNA repair; and regulation of alternative mRNA slicing via spliceosome) were selected for detailed visualization of the data. Enrichment plots and expression heatmaps of core enrichment genes are shown for these four gene sets.
Among the lncRNAs identified by standard RNA-seq data analysis, 80 differentially expressed transcripts (|logFC|>1, FDR<0.05; 64 up-regulated/16 down-regulated, Figure 3A) were found. MEG3 and MIRLET7BHG were the most up-regulated lncRNAs and RP13-685P2.5 and RP11-424C20.2 were the most down-regulated lncRNAs (Figure 3B). Interestingly, increased expression of MEG3 has already been associated with poor outcome in high risk MDS (38). MIRLET7BHG is a precursor gene for let7b miRNA which is also up-regulated in MDS (39).
Differentially expressed lncRNAs in patients with myelodysplastic syndromes/acute myeloid leukemia with myelodysplasia-related changes compared to healthy controls. (A) Heatmap of differentially expressed transcripts (|logFC|>1, FDR<0.05, N=80). (B) Boxplots of the four most deregulated lncRNAs.
Furthermore, circRNA data were analyzed by DEA. Of the 6,015 circularized transcripts used for the analysis, a significant change (|logFC|>1, FDR<0.05) was observed in the level of 23 circRNAs (3 up-regulated/20 down-regulated). Up-regulation was detected of circRNAs spliced from DPP10, ZEB1, and FNDC3B genes and down-regulation of circRNAs processed from ARPP21, METRNL, EBF1, SKI, ABCC1, and AKAP12 genes. More detailed information about the deregulated circRNAs, including their sequences and localizations in the genome, is shown in Table III.
Differentially expressed circRNAs between patients with myelodysplastic syndromes/acute myeloid leukemia with myelodysplasia-related changes and healthy controls (|logFC|>1, FDR<0.05).
DEA performed on TE data showed that the most significantly deregulated clades of TEs between patient and healthy control samples were ERV1 (mean logFC=0.14, p=9.34E-09), ERV3 (mean logFC=0.15, p=5.74E-04), and SINE (mean logFC=–0.12, p=1.27E-03) (Figure 4A). At the level of individual elements, only four of them were identified with significant deregulation (|logFC|>1 and FDR<0.05): up-regulation of LOR1b_LTR (ERV1 clade), MLT2F (ERV3 clade), X5A_LINE (CR1 clade), and MIR3 (SINE clade (Figure 4B). However, 80 TEs were significant at the cutoff level FDR <0.05, and the illustrative heatmap in Figure 4C shows their levels in the patient and control samples.
Deregulation of transposable elements in patients with myelodysplastic syndromes/acute myeloid leukemia with myelodysplasia-related changes and controls. (A) Boxplot showing deregulation (logFC) of different clades of TEs. (B) List of the four most significantly (|logFC|>1, FDR<0.05) deregulated transposable elements and (C) their expression boxplots. (D) Illustrative heatmap of 80 transposable elements deregulated with the cutoff FDR<0.05.
Pre-treatment transcriptional differences related to future AZA response. Following the collection of BM samples, all tested patients started to receive AZA therapy. To assess expression differences in various groups of transcripts with relation to their predictive power, patients were divided into two groups of samples, RESP (N=11) and nRESP (SD and PD, N=16) and DEA was performed. Table IV shows lists of the most deregulated transcripts from each category. Overall, DEA between RESPs and nRESPs identified significant (|logFC|>1, FDR <0.05) deregulation of 202 PCGs (111 up-regulated/91 down-regulated in RESP), 34 lncRNAs (34 down-regulated), 21 circRNAs (4 up-regulated/17 down-regulated), and down-regulation of one TE. The three PCGs most up-regulated in RESPs were SEC14L5, TUBB1, and TREML1, and the most down-regulated PCGs were ADGRB3, LRAT, and ARHGAP20. Interestingly, several of the most up-regulated PCGs were involved in megakaryocyte differentiation (e.g., TUBB1, TREML1, PPBT, ITGB3, and THBS1). Detailed pathway analysis of PCG results is described below in a separate chapter. Among the lncRNAs with significantly altered levels between RESPs and nRESPs (details included in Table V), four lincRNAs were found from the 21p11.2 chromosomal locus (CH507-513H4.3-6) and several other lncRNAs from different categories (lincRNAs, antisense RNAs, rRNAs, miRNA precursors, etc.), the majority of which are completely unannotated (their details are included in Table VI). Regarding deregulated circRNAs, significant deregulation was observed in a circRNA processed from KDM1A histone demethylase gene, a circRNA from MED21 (transcriptional regulator of RNA polymerase II), two circRNAs processed from zing-finger transcription factors (ZNF208 and ZNF608), and circRNAs backspliced from transmembrane proteins (SLC8A1/3, PLXNC1), etc. Among TEs, only L1 element L1PBB_5 showed significant down-regulation between RESPs and nRESPs, but some other TEs were moderately altered at raw p<0.05 (Table IV).
The top deregulated transcripts between responders and nonresponders from each RNA category defined by differential expression analysis.
Differentially expressed lncRNAs between azacitidine responders and nonresponders (|logFC|>1, FDR<0.05, all of them were down-regulated in RESP patients).
Differentially expressed circRNAs between azacitidine responders and nonresponders (|logFC|>1, FDR<0.05).
Molecular markers for the prediction of patient response to AZA. Based on differences in gene expression levels, this study aimed to predict future responses to AZA treatment, which can help us individualize patient care in routine clinical practice. Although DEA was effective in the identification of individual transcripts stratifying RESPs from nRESPs, machine learning methods also take into consideration relations among molecules and might thus eliminate similarly expressed genes from the final set of predictive markers. A classifier combining several independent molecular markers of diverse characteristics should provide more reliable classification of patients compared to simple information on an individual gene level. In this study, a machine learning approach was used and three different classification algorithms (SVM-RFE, voomNSC, and GLMNET) were applied to identify the most powerful set of molecular markers. Using all datasets together (i.e., PCGs, lncRNAs, circRNAs, and TEs), the highest classification accuracy was clearly reached by GLMNET (Figure 5) and this algorithm was therefore chosen for further investigations. The GLMNET performance graph suggests that the classifier was able to reasonably distinguish between RESPs and nRESPs because the observed AUC areas were significantly higher than the value of 0.5 (corresponding to a random classification), and a significant difference could be observed for a broad scale of signature sizes. The GLMNET classifier also reached the optimum performance for a limited signature size of 100 transcripts (AUC=0.85) and its performance did not further grow with its extension. This observation led us to the conclusion that it was possible to distinguish between the patient classes with a relatively small set of transcripts.
Application of three different machine learning classification algorithms (SVM-RFE, voomNSC, and GLMNET) for classification of response to the azacitidine treatment. The graph suggests that GLMNET outperforms the other methods and that a limited set of transcripts is sufficient for the best split. The larger signatures only lead to overfitting.
In addition to the choice of an appropriate data-classification algorithm, contribution of different transcript categories to the response prediction was tested by evaluating each transcript category separately. Surprisingly, lncRNAs overperformed the other types of molecules, as they provided the best classification power between RESPs from nRESPs (AUC=0.83 defined by GLMNET) (Figure 6A).
The most powerful markers of azacitidine response defined by machine learning (GLMNET). (A) The four classes of transcripts had substantially different stratification power to discriminate azacitidine responders and nonresponders. (B) Heatmap showing the expression of the final set of response markers. (C) Boxplot of expression values of individual transcripts discriminating azacitidine responders and nonresponders. The lncRNAs marked with asterisks are discussed in the text in detail.
After optimization of the methodology, GLMNET was employed to identify a particular small set of the most powerful markers. A compromise was reached between the signature size (the number of used transcripts) and its performance, and the size was kept as small as possible while approaching the best possible AUC of 0.85 observed earlier. The final set contained 19 transcripts, of which 14 were lncRNAs, only 4 PCGs and 1 circRNA, and none TE (Figure 6A and B, Table VII). This distribution reinforces the earlier observation that lncRNA molecules represented the core of the RESP/nRESP classification signature. Furthermore, a formula calculating the prediction score of a patient was defined based on the weights of the 19 transcripts selected by GLMNET. The regression score was calculated as follows:
Prediction score=–12.43+βTx
Final set of the most powerful markers of future patient response to azacitidine defined by GLMNET. The lncRNAs marked with asterisks are discussed in the text in detail.
The score uses expression values (x) and regression coefficients (β) for these 19 transcripts (T) (β coefficients can be found in Table VII). When simplified and binarized, a score >0 predicts that the tested patient will respond to AZA, whereas a score <0 suggests future AZA failure.
Pathway analysis of transcripts associated with the AZA response. To explore cellular processes associated with patient capability to benefit from AZA treatment, GSEA was performed on PCGs (Figure 7). Many of the processes associated with the AZA response were also generally deregulated in myelodysplasia (described above and summarized in Figure 2). In particular, epigenetic regulation of gene expression and DNA repair, which were generally suppressed in patients (compared to healthy donors), were suppressed specifically in AZA nRESPs and activated in AZA RESPs. Based on these observations, two GO gene sets were selected, recombinational repair (NES=1.95; p=0.006; core enrichment genes included BRCA1/2, H2AFX, POLQ, RAD51, MCM8, and XRCC2; Figure 8A and B) and negative epigenetic regulation of gene expression (NES=1.91; p=0.014; core enrichment genes included multiple HIST1H transcripts and other epigenetics-related factors, e.g., DNMT1/3B and EZH2; Figure 9A and B), for further investigation.
Pathway analysis of differentially expressed genes in azacitidine responders and nonresponders. Plot of gene set enrichment analysis results shows p-values (logarithmically scaled) and normalized enrichment scores (NES) of the top 10 most enriched pathways in each group of samples sorted according to descended NES values. Two significant gene sets (marked by asterisks; recombinational repair and negative regulation of gene expression - epigenetic) were selected for further evaluation. Heatmap shows the top 50 protein coding genes differentially regulated between the two groups of samples.
Association of the recombinational repair pathway with the azacitidine response. (A) gene set enrichment analysis enrichment plot, (B) expression heatmap of core enriched genes, and (C) coexpression network of transcripts associated with this process. Blue – down-regulation in nonresponders, red – up-regulation in nonresponders, circle-noncoding transcript, square - protein coding gene.
Association of negative epigenetic regulation of gene expression with azacitidine response. (A) gene set enrichment analysis enrichment plot, (B) expression heatmap of core enriched genes, and (C) coexpressional network of transcripts associated with this process. Blue – down-regulation in nonresponders, red – up-regulation in nonresponders, circle - noncoding transcript, square - protein coding gene.
This study further aimed to explore specific noncoding transcripts related to these two processes (epigenetic regulation and DNA repair) in the context of future AZA responses. Computational modeling was applied and coexpressional networks were constructed based on correlation rates among all transcripts (PCGs, lncRNAs, circRNAs, and TEs) deregulated between AZA RESPs and nRESPs. Interaction graphs for both GO terms displaying 20 PCGs the most representative given GO term (the most deregulated leading-edge genes) and 12 noncoding transcripts (the most correlated with these PCGs) are shown in Figure 8C and Figure 9C. Although no TE was identified as related to any of the two processes, four circRNAs (hsa_circ_0001580 from KDM1B gene, hsa_circ_0008737 from CAMTA1 gene, chr2:136,620,176-136,622,733 from MCM6 gene, and chr17:59,853,761-59,878,835 from BRIP6 gene) were associated with recombinational repair pathway, and three circRNAs (chr2:61,722,589-61,727,029 from XPO1 gene, chr5: 134,681,712-134,688,690 from H2AFY gene, and chr14:97,299,803-97,342,457 from VRK1 gene) with epigenetic regulation processes. Regarding lncRNAs, several transcripts with unknown functions were included in the networks, associating them with a cellular process for the first time (e.g., epigenetic regulation: CTC-298B17.1, CTB-83J4.2, RP11-449P15.1; DNA repair: AC004159.2, RP5894D12.5, C20orf166; others can be found in Figure 8C and Figure 9C). Interestingly, three of the lncRNAs defined by GLMNET as the strongest predictors of AZA response were defined by computational modeling as significantly associated with recombinational repair (CTC-482H14.5, RP11-419K12.2, and RP11-736I24.4) and/or negative epigenetic regulation of gene expression (CTC-482H15.5). CTC-482H14.5 (up-regulated in nRESPs) is an antisense RNA of the KDM4B gene, RP11-419K12.2 (down-regulated in nRESPs) belongs to the lincRNA category, and RP11-736I24.4 (down-regulated in nRESPs) is a transcribed unprocessed pseudogene of the OTUD7A gene. Additional descriptions of these lncRNAs are included in Table VIII.
LncRNAs predicting future azacitidine response (assessed by GLMNET as the best predictors) and significantly (FDR<0.05) associated with recombinational repair and/or negative epigenetic regulation of gene expression GO terms.
Discussion
Over the last several years, the understanding of pathogenesis of MDS and AML-MRC has expanded. Multiple publications have described alterations of PCGs involved in chromatin modification, DNA methylation, RNA splicing, cohesin complex, transcription factors, cell signaling, and DNA damage (40). In this study, suppression of a multitude of PCGs associated with epigenetic mechanisms, DNA repair, and RNA splicing was documented, confirming the previous findings and proving validity for our current data. However, the expression of various classes of noncoding transcripts and their involvement in the abovementioned MDS/AML-associated pathways were focused in relation to their contribution to the response to AZA treatment. Classic types of lncRNAs ordinarily identified by classic RNA-seq data analyses were examined but special attention was paid to two additional classes of noncoding transcripts: circRNAs and TEs. The deregulation of lncRNAs in MDS has previously been addressed (34), but circRNAs and TEs have only marginally been studied in the context of this disease.
The expression profiles of the four types of transcripts (PCGs, lncRNAs, circRNAs, and TEs) were compared, and it was noticed that although PCGs and lncRNAs were able to stratify MDS/AML-MRC patients from healthy controls, the transcription of circRNAs and TEs was more heterogeneous, with substantial interindividual variability among patients. Therefore, PCGs and lncRNAs seem to be more relevant molecules in disease pathogenesis. However, circRNAs and TEs are much less numerous groups of transcripts than PCGs and lncRNAs, which might worsen the clustering output. Additionally, some bioinformatical bias might affect the results. For example, low numbers of reads covering circRNA junctions in comparison to the detection of linear transcripts or repetitiveness in TE sequences affecting read alignment might influence data processing. As information on circRNAs and TEs has only recently started being uncovered, the possibility of false discoveries and some differences in the number of detected molecules must be expected. For data annotation, only one bioinformatical tool for each type of transcripts was chosen (i.e., CIRI2 for circRNAs and SalmonTE for TEs). Although these tools are widely recommended, multiple annotation tools with different algorithm strategies should ideally be combined in more detailed, ongoing studies to achieve more reliable results. However, the comparison of multiple bioinformatics tools was out of the scope of this early study. The major focus of this study was an early characterization of expression changes of circRNAs and TEs and their comparison with standard gene expression data in relation to AZA treatment prediction.
CircRNAs are circularized transcripts produced by the so-called back-splicing process. They can act as regulators of gene expression but only limited information about the precise functions of individual circRNAs is available thus far. Due to their unusual stability and expression specificity, they appear to be important candidates for clinical biomarker research (41). In MDS/AML-MRC, several circRNAs processed from genes encoding known regulators of hematopoiesis and/or oncogenesis were deregulated, e.g., down-regulation of circRNAs derived from METRNL (42), EBF1 (43), SKI (44), ABCC1 (45), and AKAP12 (46), and up-regulation of circRNAs derived from ZEB1 (47, 48) and FNDC3B (49, 50). For example, the SKI protooncogene is an important regulator of hematopoietic stem cell activity and its overexpression leads to myeloproliferative disease (44). In this study, the expression of SKI-derived hsa_circ_0007120 was down-regulated in MDS/AML-MRC. This circRNA was found to be strongly down-regulated during the differentiation of mesenchymal stem cells into cardiomyocyte-like cells, suggesting its involvement in cell proliferation and differentiation-related processes (51). Furthermore, the protein encoded by the AKAP12 gene acts as a tumor suppressor regulating cell cycle progression and inhibiting Src-mediated oncogenic signaling (52). The AKAP12 promoter is unmethylated and expressed in normal HSPCs but exhibits profound methylation and transcriptional silencing in leukemic myeloblasts (46). In MDS/AML-MRC, down-regulation of both, linear and circular (hsa_circ_0078299) forms of AKAP12 was detected, suggesting common regulation within the locus. Interestingly, hsa_circ_0078299 was identified as a candidate biomarker of human stroke with a potential role in gene expression control and inflammation, and SRSF2 splicing factor was predicted as its target gene (53). Another potentially relevant circRNA deregulated in MDS/AML-MRC patients is hsa_circ_0000228 transcribed from the ZEB1 gene. The product of the ZEB1 gene, a zinc finger transcription factor, acts as a transcriptional regulator in hematopoiesis, coordinating HSC self-renewal, apoptosis, and multilineage differentiation fates. Additionally, it is required to suppress leukemic potential in AML (48). Several ZEB1-derived circRNAs have been shown to function as oncogenes in breast cancer (54), gastric cancer (55, 56), and hepatocellular carcinoma (57) by sponging various complementary miRNAs. Although the expression of the linear form of ZEB1 was not changed in the MDS/AML-MRC patient samples, an increased rate of circularization of ZEB1 transcripts (i.e., an increase in hsa_circ_0000228) might exert specific roles in regulatory processes in CD34+ cells.
TEs represent detrimental units of DNA whose harmfulness lies in their inherent mobile nature. Their expression can lead to insertional mutagenesis, chromosomal rearrangements, and genomic instability (58). The significance of TEs has previously been addressed in various subtypes of primary AML (59-62), and it has been shown that TE expression was able to predict prognosis independent of mutation-based and coding gene expression-based risk stratification (59). In addition to AML, Colombo et al. (60) also compared the activity of TEs in high-risk and low-risk MDS on a limited number of samples (six high-risk and six low-risk MDS patients) and noticed that TE expression is significantly down-regulated in high-risk disease. They proposed that the induction of TEs in low-risk cases may be a potential mechanism for immune-mediated clearance of cancer cells via the viral recognition pathway, whereas TE suppression in high-risk cases may enable immune escape of cells (60). In this study, CD34+ BM cells from only high-risk MDS patients were examined. Therefore, it was not possible to compare TE activity in low-risk disease. However, the rate of TE deregulation in our high-risk MDS samples was much lower than that of PCGs and other types of noncoding transcripts, negating the fundamental roles of TEs in high-risk disease pathogenesis.
Prediction of response to AZA therapy is an important challenge of applied research in MDS/AML. Current research activities are usually focused on the identification of an ideal, single biomarker whose alterations (such as a mutational change or transcription deregulation) would be predictive of treatment response with high sensitivity and specificity. The patient capability to respond to AZA and the identification of prediction markers have been repeatedly addressed (11, 12, 18). However, given the genomic heterogeneity of the disease, no universal biomarker has been identified thus far. This failure is likely triggered by complex interactions among multiple genomic alterations, affecting the resulting phenotype. Therefore, a combined approach studying different types of molecules and their relations is required to yield an understanding of how these alterations define responsiveness to therapy. In addition to alterations in PCGs, the efficiency of AZA therapy can be influenced by the expression of various noncoding transcripts, contributing together to the clinical effects of the drug. At the level of individual genes, several noncoding RNAs significantly deregulated between AZA RESPs and nRESPs were identified. However, the pathological importance of most of them was somewhat questionable given that they are completely unannotated in current databases. This was especially the case for the most deregulated lncRNAs (e.g., down-regulation of four lincRNAs CH507-513H4.3-6 in AZA RESPs) because the only information available about them is often their genomic location. CircRNAs appeared to be a more relevant category of noncoding transcripts, as altered circularization was observed of several genes potentially relevant for AZA responsiveness (e.g., KDM1A histone demethylase or MED21 regulating transcription with RNA polymerase II). Regarding TEs, only one TE was significantly deregulated between AZA RESPs and nRESPs (down-regulation L1 element L1PBB_5 in AZA RESPs), with merely slight changes in the expression profiles of other TEs. The absence of altered TE transcription doubts the hypothesis suggesting that the clinical effects of AZA might be caused, at least partly, by demethylation of TE regions leading to increased production of TE transcripts, which might in turn activate the innate immune system through the viral defense pathway (60, 63).
Considering the above-discussed data on alterations of individual transcripts predisposing to AZA response or failure, a multibiomarker approach combining several features of different natures seems more appropriate to predict the response with high robustness. By machine learning methods, we compared the four categories of transcripts, and unexpectedly, lncRNAs substantially overperformed the prediction potential of PCGs. However, neither circRNAs nor TEs provided significant potential to stratify AZA RESPs and nRESPs. The combined set of the best predictors of AZA response finally included 19 genes, and surprisingly, 14 of them were of lncRNA nature, whereas only four PCGs, one circRNA, and no TE were added to this final list of the most powerful stratifying molecules.
To better understand the contribution of deregulated noncoding transcripts to the molecular biology of MDS/AML treatment, important cellular pathways associated with AZA were first defined. Previously, the quiescent phenotype of HSCs was proposed as a key factor in the prediction of treatment outcome in MDS (64). In a preceding publication, it was identified that the BM of AZA RESPs contains actively cycling cells poised for erythromyeloid differentiation, whereas in nRESPs, BM HSPCs display a signature of a quiescent state (repression of genes related to active cell cycling, replication, and DNA damage response) (18). These metabolic pathways seem to be closely interconnected with epigenetic modifications; the noncycling status has been associated with decreased expression of the HIST1 cluster and therefore with the transcription-nonpermissive chromatin configuration maintaining cell quiescence (18). In this study, it was confirmed that AZA nRESPs display suppressed expression of epigenetic regulators and many chromatin related genes (especially those from the HIST1H and HIST2H families) in addition to genes involved in recombinational repair. In the process of epigenetic modulation, DNMT1 (a maintenance DNA methyltransferase) was defined as one of the core PCGs related to the AZA response. Previously, it was demonstrated that DNMT1 serves as a pharmacologic target of HMAs (65), and down-regulation of DNMT1 expression in AZA nRESPs has been attributed to noncycling status and reduced replication of HSPCs in these patients (18). Another core PCG of epigenetic mechanisms involved in the response to AZA was the EZH2 gene, which was also specifically down-regulated in our cohort of AZA nRESPs. EZH2, a histone methyltransferase gene, is recurrently mutated in MDS/AML patients, and its mutations are associated with a poor prognosis of these patients (66). EZH2 is significantly suppressed in MDS (67), and its CpG methylation level is reduced by AZA (68). Regarding the DNA damage response, another pathway described in the context of the AZA response (18), we identified RAD51, BRCA1/2, and several members of the MCM gene family as core PCGs down-regulated in the patients with future AZA failure. Alterations in RAD51 (36, 69) and BRCA (70) genes have previously been reported in MDS. The expression of MCM genes, which are involved in the initiation and elongation phases of DNA replication, was reduced following epigenetic therapy (71). Similar to downmodulation of epigenetic machinery, attenuation of DNA repair and response pathways has also been described as a characteristic feature of HSC quiescence, which underlies DNA damage accumulation in these cells (72). In summary, suppression of epigenetic regulation and recombinational repair seem to be crucial for the cell behavior of MDS/AML-MRC HSPCs with respect to the AZA response. Therefore, ncRNAs related to these two processes were explored to identify noncoding transcripts necessary for AZA-triggered mechanisms.
Although many deregulated noncoding transcripts are being described in cancer, their functional characterization is difficult. Construction of co-expression networks of PCGs with other RNAs, which are regulated similarly to PCGs, enables functional analysis of transcripts with unknown functions. This “guilt-by-association” strategy works with the principle that genes with related functions tend to have similar expression profiles. PCGs in a coexpressional module are associated with signaling pathways, attributing the same functions to unknown RNAs in this network. In this study, coexpressional networks were constructed for genes associated with epigenetic regulation and recombinational repair pathways. Unsurprisingly, no TE was related to these processes, repeatedly suggesting little significance of these elements in high-risk MDS. In contrast, several lncRNAs and circRNAs co-regulated with core PCGs of the two processes were defined. Interestingly, the majority of the circRNAs found by co-expressional networks were processed from genes already known to be involved in MDS-related processes (recombinational repair - MCM6 and BRIP6; epigenetic regulation - KDM1B and H2AFY; transcription, cell cycle and proliferation regulators - CAMTA1, XPO1 and VRK1). In particular, KDM genes seem to be closely related to responsiveness to AZA. KDMs encode demethylases of histone lysine residues (KDM1 mediates demethylation of H3K4 and H3K9, KDM4 can act on H3K9, H3K36, and H1K26), yet they can also demethylate a number of nonhistone proteins such as p53, E2F1, DNMT1, or STAT3. KDMs are often differentially expressed in leukemia, and cooperate with transcription factors to activate or repress gene expression (73). In this study, deregulation of several noncoding transcripts was observed related to KDMs – CTC-482H14.5 (asRNA to KDM4B gene, down-regulated in AZA RESPs), hsa_circ_0003889 (circRNA processed from KDM1A, down-regulated in AZA RESPs), and hsa_circ_0001580 (circRNA processed from KDM1B, up-regulated in AZA RESPs). Unlike these noncoding transcripts, the levels of KDM PCGs (namely KDM1A/1B/4B) were not altered between AZA RESPs and nRESPs, suggesting that these ncRNAs play specific roles independent of their host genes.
Finally, the best predictor genes defined by GLMNET were linked to cellular processes related to MDS. Among these predictors, three ncRNAs potentially associated with epigenetic regulation and/or recombinational repair were identified. First, CTC-482H14.5 is an antisense RNA of the KDM4B gene and is coexpressed with several MDS-related PCGs, e.g., MCMs (recombination repair) and DNMT1 (epigenetic regulation). Second, RP11-419K12.2 and RP11-736I24.4 link to recombinational repair through coregulation with several RAD genes, XRCC2/3, BRCA2, or FANCB. RP11-419K12.2 is a lincRNA with an unknown function located in proximity to GSDMC (gasdermin C, melanoma-derived leucine zipper-containing extranuclear factor), and RP11-736I24.4 is a transcribed unprocessed pseudogene of OTUD7A (OTU deubiquitinase 7A). OTUD7A is a possible tumor suppressor gene controlling NF-κB expression through TRAF6 (74), a critical process implicated in the pathogenesis of MDS (75). Therefore, novel functions in epigenetic regulation and recombinational repair processes were predicted for these three unexplored lncRNAs, pointing to possible mechanisms by which their deregulation can be associated with AZA responsiveness.
In conclusion, this study addressed the expression and roles of various types of noncoding transcripts in AZA treatment of high-risk MDS and AML-MRC patients, suggesting particular ncRNAs associated with AZA responsiveness. Importantly, it was shown that lncRNAs can be utilized as novel predictive markers and that the combination of several markers of different natures is a more robust approach than monitoring the expression of a single gene. To the authors’ knowledge, this is the first study identifying particular members of circRNA and TE transcript families deregulated in high-risk MDS and AML-MRC. In particular, deregulation of several circRNAs produced from hematopoiesis- and oncology-relevant genes was identified, whereas altered transcription of TEs seemed of less significance in this disease. However, the data in this pilot study have to be validated in independent cohorts of patients, different bioinformatical approaches analyzing circRNA and TE levels should be used, and functional studies proving the predicted roles in cellular processes have to be performed before rendering final conclusions about the roles of these noncoding transcripts in myelodysplasia.
Acknowledgements
This work was supported by grants from AZV CR (No. 17-31398A and NU20-03-00412), GA CR (No. 20-19162S), and MH CZ-DRO (UHKT, 00023736). Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum provided under the program “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042) was greatly appreciated.
Footnotes
Authors’ Contributions
Conceptualization, Michaela Dostalova Merkerova; Funding acquisition, Michaela Dostalova Merkerova; Investigation, Katarina Szikszai, Andrea Hrustincova, Iva Trsova; Methodology, Zdenek Krejcik; Bioinformatics, Jiri Klema, David Kundrat, Anh Vu Le; Resources, Jaroslav Cermak, Anna Jonasova; Supervision, Monika Belickova, Jaroslav Cermak, Michaela Dostalova Merkerova; Writing – original draft, Michaela Dostalova Merkerova; Writing – review & editing, Zdenek Krejcik, Monika Belickova and Jaroslav Cermak.
Data Availability
The RNA-seq data have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under accession number PRJNA679200.
Conflicts of Interest
The Authors declare that there are no competing interests.
- Received November 9, 2021.
- Revision received January 7, 2022.
- Accepted January 14, 2022.
- Copyright © 2022, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved
References
In this issue
Jump to section
Related Articles
Cited By...
- No citing articles found.