Machine Learning Sequence Modeling Identifies Gene Regulatory Responses to Bone Marrow Stromal Interactions in Multiple Myeloma

Abstract

Introduction Multiple myeloma (MM) is a hematological malignancy of plasma cells which affects different organs including the bone marrow. The interaction between MM and the bone marrow stromal cells (BMSC) has been shown to affect the disease progression and response to treatment (Kumar et al. Multiple myeloma. Nat. Rev. Dis. Primers, 2017). An in vitro system to isolate the signaling mechanisms between BMSC and MM cells has been developed by (Dziadowicz et al., Cancers, 2022), who co-cultured MM cell lines in transwell media (TSW) exposed to BMSC secreted factors and compared their ATAC-seq profiles to MM cell lines (MONO) in isolation without the exposure to BMSC. While some of the gene expression changes induced by the BMSC-MM crosstalk were reported in (Dziadowicz et al., Cancers, 2022), here we model the entire chromatin accessibility profiles with a DNA-sequence based machine learning model to detect the altered Transcription Factor (TF) and enhancer activity responsible for driving the gene expression changes. We used sequence-based machine learning models (Beer et al., 2020, Annu. Rev. Genom. Hum. Genet.) and applied a systematic approach to model the chromatin accessibility changes induced by the interaction of BMSC and MM cells (MM-TSW) vs. MM cell lines cultured in isolation (MM-MONO). Methods We trained gapped-kmer SVM (gkm-SVM) machine learning models on MM1S and RPMI8226 MM cell line ATAC-seq data generated by Dziadowicz et al. as well as 1270 DNase-seq chromatin accessibility profiles of various human primary tissues and cell lines (including MM1S and RPMI8226) from the ENCODE consortium. We trained gkm-SVM on the top 2000 differentially-accessible distal ATAC peaks (>2000bp from transcription start site, e.g. enhancers) to find the differential TF activity between MM-TSW and MM-MONO (Figure 1) and found the TF binding site motifs explaining the gkmSVM output kmer weight distribution (Figure 2). We also identified the set of differentially-expressed TFs and performed gene set enrichment using MSigDB. Results We found that RUNX and Ebox-binding (5‘-CACCTG-3‘) TF family members are the most active TFs in MM1S-MONO and RPMI8226-MONO enhancers (gkmSVM AUROCs~0.92). The contribution of the Ebox TFs to distal chromatin accessibility in the MM cell lines is greater than in any other ENCODE sample (n=1270). Principal Component Analysis (PCA) on distal ATAC peaks compared to promoters showed that the MM transcriptional response to BMSC is more clearly explained by the differences in enhancer activity rather than promoters. This response is driven by dysregulated activity of a set of enhancer regions in MM-TSW compared to MM-MONO (Figure 1) and our machine learning method can identify these regulatory changes (AUROC~0.90). In addition to the altered AP1 and JAK/STAT activity reported by Dziadowicz et al., we discovered that RUNX and Ebox-binding transcriptional activity are downregulated in MM cells following the interaction with BMSC (Figure 2). Although our RNA-seq analysis found differential expression of some Ebox-binding TFs such as TCF3 and ZEB1 and RUNX members, the identification of the exact TF member(s) driving these RUNX and Ebox/E2A downregulation in MM-TSW is still challenging. Furthermore, PCA analysis of the gene expression profiles of MM1S and RPMI8226 cells showed that the targets of Bone Morphogenetic Protein 2 (BMP2) are upregulated in MM-TSW (FDR < 10 -28). Conclusion While some of the gene expression changes induced by the BMSC-MM interaction were reported in (Dziadowicz et al., Cancers, 2022), our modeling of the entire chromatin landscape discovered the previously unreported downregulation of RUNX and Ebox/E2A TFs in MM as the dominant response of MM interacting with BMSCs, accompanied by an upregulation of AP1 and JAK/STAT. RUNX and E2A-binding families of TFs include RUNX1/2/3, ZEB1/2, SNAI1/2 and TCF3, some of which we previously found to be strongly involved in epithelial to mesenchymal transition in gastric cancer and its invasiveness (Razavi-Mohseni et al., submitted, 2023). Bone marrow microenvironment interactions with MM are known to affect the disease progression and treatment, therefore, identifying the set of TFs driving the MM response to stroma may provide better understanding of the disease and better therapeutic opportunities.