For this post, I decided to re-analyse a dataset from a recent article by Logie et al (2021) using our Omics Playground platform.
Multiple Myeloma (MM) is a malignancy that frequently develops both primary and acquired resistance to anticancer drugs. Resistance is often associated with the constitutive hyperactivation of tyrosine kinase signaling. In this article, the authors tested the activity of two Bruton’s Tyrosine Kinase (BTK) inhibitors on glucocorticoid-resistant MM cells, to better understand their mode of action and also test their antitumor properties. They concluded that the preclinical phytochemical withaferin A (WA) was particularly effective thanks to its more promiscuous tyrosine kinase targeting and also noted its ability to reduce both BTK mRNA and protein expression.
In this re-analysis I focused on a GEO dataset (GSE162475) that contained mRNA profiles of glucocorticoid (GC) sensitive (MM1S) and GC-resistant (MM1R) B lymphoblast cell lines after treatment with WA or DMSO (UT, untreated control).
I re-analysed the dataset first of all by repeating the analysis of the original article and then expanded upon it thanks to the modules available in the Omics Playground platform.
- Clustering Analysis
Both the PCA plot (Fig. 1) and the heatmap based on individual gene expression (Fig. 2A) indicated a clear structure in the dataset. Thus the samples clearly cluster into two main separate groups based on GC-resistance phenotype that were further subdivided based on treatment with WA or exposure to DMSO (UT).
The platform also provides a functional annotation based on different gene sets of the gene clusters used to produce the heatmap. The top 12 Hallmark terms are displayed in Fig. 2B. Cluster S1 distinguishes WA treated cells from untreated cells, independently of their resistance phenotype (Fig. 2A). The annotation includes terms associated with apoptosis, cell cycle regulation (p53 pathway), the IL-6/JAK/STAT3 pathway (which plays a key role in the growth of many cancer types) and the regulation of transcription factors (TNFα signalling via NF-κB). An annotation with GO terms also identified terms associated with the activation and differentiation of B- and T-cells and angiogenesis, consistent with observations in the original article (data not shown).
For the differential expression analysis, non-coding genes of the genome (such as lincRNAs) were also included. Comparing untreated GC-resistant and GC-sensitive cell lines using DESEQ2 (lrt test) with FDR <0.01 and logFC<-1 or logFC>1 produces a smaller number of dysregulated genes compared to the list obtained with the same conditions in the paper (413 vs. 1383) (Fig. 3). Upon closer inspection of the read counts, we found that the difference was due to the number of genes with very low read counts. As the platform employs a standard filtering rule expecting at least 1 cpm for at least two samples, many of the genes considered in the original article (such as HER4 and BLK), which had very low coverage, were discarded. I preferred to keep these parameters intact, due to the unreliability of comparing genes with low expression levels.
The list of the 20 most dysregulated genes by adjusted p-value was different from the list presented in the paper, with a notable increase in the number of downregulated genes (Table 1). BTK was present in both lists. However, all 20 genes from the original list were also identified as significantly dysregulated in the platform output (Table 2).
When I looked more in detail at the BTK gene, I observed that while there was a significant expression difference between GC-resistant and GC-sensitive cell lines, treatment with WA only resulted in a modest (but statistically significant) reduction of expression (Fig. 4).
3) Differential expression analysis: WA treated vs untreated GC-resistant MM cell line
Treatment of GC-resistant MM cells with WA resulted in a total of 362 dysregulated genes, as detected by the DESEQ2 with logFC>1 or logFC<-1 and P<0.01 (Fig. 5).
I observed, in particular, a significant increase in expression for several heat shock proteins (HSPs, e.g. HSPA6, HSPA7, HSPB1, HSP90AA1) and a significant decrease in chemokine ligand and chemokine receptor proteins (e.g. CCL3L1 and CCR2) (Table 3). The upregulation of HSPs is not surprising, as it is part of the cell reaction to WA exposure, while WA also inhibits the expression of various chemokines.
4) Gene Set Enrichment Analysis (GSEA)
Next, I accessed the GSEA module of the platform that can perform KEGG pathway, GO terms analysis, as well as access the drug connectivity map (Drug cMAP) database and over 50’000 gene sets from public databases (Akhmedov et al, 2020). The results were obtained by intersecting the hits from three different statistical methods, gsva, fgsea and fisher’s exact test, with an FDR<0.05 and logFC>1/logFC<-1.
WA treatment in the GC-resistant cell line produced a positive correlation with several anticancer drugs, including WA itself (Fig. 6). The same was also true when GC-sensitive cell lines were treated with WA (Fig. 7) . Among the top 24 most positively correlated drugs were several heat shock protein 90 (HSP90) inhibitors, such as monorden/radicicol, alvespimycin, tanespimycin and geldanamycin (Fig.6 and Fig. 7).
In order to get a better overview of the most common modes of action (MOA) of the statistically significant positive drug profile matches, the platform produces a plot based on the functional annotation (if available) for the public profiles. The plots obtained for the WA treated vs untreated GC-resistant (Fig.8A) and GC-sensitive (Fig.8B) cell line samples are remarkably similar. Heat shock protein (HSP) inhibitors are prominent in both plots, confirming the previous observation on the presence of several HSP90 inhibitors among the top hits. It is worth noticing that the glucocorticoid receptor (GR) in eukaryotic cells forms a heterocomplex with HSP90, which also regulates its function.
Other notable MOA include mTOR, NF-κB pathway and proteasome inhibition. On the other hand, BTK inhibition is absent from the list. This is due to the fact that several genes, whose expression is critical for this MOA (such as BMX and BLK), are not present in the gene dataset due to their insufficient read coverage across samples, as discussed above.
A quick search in the literature reveals that WA’s anticancer properties extend well beyond the inhibition of BTK. Thus, WA is known to interfere with heat shock responses, the mTOR pathway, the NF-κB pathway, the proteasome, protein synthesis and oxidative stress responses (here), as well as displaying inhibitory potential against various epidermal growth factor receptors (EGFR). All these mechanisms are present in the MOA plots.
It is interesting to notice that while the original article focuses exclusively on the anti-BTK activity of WA against MM cancer cells, analysis of the data through the platform neatly reveals in a single plot the multiple modes of action of the drug. This may contribute to WA’s more potent effect on MM GC-resistant cell death than ibrutinib, as noticed in the study.
Expanding beyond the purpose of the original study, I performed a GSE analysis on the gene expression profile of untreated GC-resistant vs. untreated GC-sensitive cell lines against the Drug cMAP L1000 database in order to find drug gene expression profiles that were negatively correlated to it and could thus indicate suitable drugs to overcome GC-resistance in MM. One of the compounds that stood out was strophanthidin (Fig. 9). Strophanthidin is an inhibitor of the MAPK signalling pathway in human cancers and a study demonstrated that inhibition of the MAPK pathway could reverse GC-resistance in leukemia.
5) Biomarkers Selection
To conclude the re-analysis of the data, the “find biomarkers” module was used to identify potential genetic biomarkers that can distinguish the four sample groups available in the study. Several decision trees are possible, but a common trait in all the trees that were generated was the use of the expression of gene DNAJA4 (HSP40) to distinguish WA-treated GC-sensitive samples from the others (Fig. 10A). DNAJA4 shows a drastic increase in expression level in GC-sensitive MM cells upon treatment with WA, but not in GC-resistant cells (Fig. 10B). It is part of a family of molecular chaperone proteins that have various biological functions. DNAJ-like proteins have been shown to interact with the HSP90/GR heterocomplex and potentiate its assembly. It is thus intriguing to see that DNAJA4 upregulation upon exposure to WA is completely abrogated in GC-resistant MM cells.
The re-analysis of the transcriptomic data from Logie et al confirms the main findings from the article, despite the exclusion of genes that were insufficiently expressed according to the standard filters applied in the platform. Furthermore, GSE analysis highlighted that the anti-tumoral properties of WA extended beyond the interference of BTK. This was corroborated in the literature, which shows that WA can interact with various pathways to suppress cancer cells. These included the binding of WA to HSP90, which regulates glucocorticoid receptor function. The MM GC-resistant cell expression profile indicated a possible increased susceptibility to strophanthidin, an anti-tumoral drug. Finally, a biomarker analysis revealed a significantly increased expression of gene DNAJA4 in MM GC-sensitive cells (but not in GC-resistant cells) upon exposure to WA. The lack of upregulation in GC-resistant cells appears to be a consequence of the acquisition of GC-resistance and may be related to the role played by DNAJA4 in the assembly of glucocorticoid receptors.
Read more about Omics Playground and its features here.