Glucocorticoid Resistance in Multiple Myeloma: a Dataset Re-analysis

Posted on by

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. 

  1. 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).

Figure 1.  PCA plot showing the clustering of the samples according to their phenotype (resistant or GC-sensitive) and treatment (WA or UT). The samples form distinct groups in the plot: WA-treated GC-sensitive cells (purple), WA-treated GC-resistant cells (blue), untreated GC-sensitive cells (orange) and untreated GC-resistant cells (green). 

Figure 2. (A) Clustered heatmap showing gene expression of the samples. Red indicates overexpression, blue indicates underexpression. Two distinct main clades separating GC-resistant (MM1R) and GC-sensitive (MM1S) samples are present. Furthermore, within each of the main clades, samples are separated based on exposure to WA or DMSO (UT). The four main gene clusters are indicated on the left hand side by S1-S4. (B) Annotation of the four main gene clusters in the samples, based on the heatmap in Figure 2. Cluster S1 distinguishes samples treated with WA from untreated samples both in GC-sensitive and GC-resistant cells.

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). 

Figure 3. MA plot of the pairwise comparison of gene expression between untreated GC-resistant and GC-sensitive samples. Y-axis shows the effect size (logFC), X-axis the average expression

Table 1. Top 20 dysregulated genes obtained from the Omics Playground Platform. Genes in red are also present in the top 20 list from the original article.

Table 2. Top 20 dysregulated genes from Lugie et al (2021) with p-values and logFC from the platform.

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). 

Figure 4. Gene expression across groups, expressed in logCPM. UTR= Untreated GC-resistant cell line, UTS= Untreated GC-sensitive cell line, WAR=WA-treated GC-resistant cell line, WAS= WA-treated GC-sensitive cell line.

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. 

Figure 5. MA plot of the pairwise comparison of gene expression between WA treated and untreated GC-resistant samples. Y-axis shows the effect size (logFC), X-axis the average expression.

Table 3. Significantly dysregulated heat shock proteins and chemokine-related proteins after WA treatment of GC-resistant MM cells.

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).

Figure 6. 24 most positively correlated drug expression profiles to WA treatment in the GC-resistant cell line by q-value. Matches with two public WA datasets (highlighted in orange) are present.

Figure 7. 24 most positively correlated drug expression profiles to WA treatment in the GC-sensitive cell line by q-value. Matches with two public WA datasets (highlighted in orange) are present.

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.

Figure 8. (A) MOA plot showing the top most frequent drug class having similar enrichment compared to the WA treated vs. untreated GC-resistant cell line samples. (B) MOA plot showing the top most frequent drug class having similar enrichment compared to the WA treated vs. untreated GC-sensitive cell line samples

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.

Figure 9. GSEA enrichment plot for the untreated GC-resistant vs GC-sensitive cell lines expression profile against the strophanthidin perturbation gene set. The normalised enrichment score (NES, green line), indicates a statistically significant (q<0.05) negative correlation.

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.

Figure 10. (A) One of the possible decision trees for classification based on the top most important genes by expression level. UTR= Untreated GC-resistant cell line, UTS= Untreated GC-sensitive cell line, WAR= WA-treated GC-resistant cell line, WAS= WA-treated GC-sensitive cell line. (B) Expression profiles for gene DNAJA4 across the four sample groups. UTR= Untreated GC-resistant cell line, UTS= Untreated GC-sensitive cell line, WAR= WA-treated GC-resistant cell line, WAS= WA-treated GC-sensitive cell line.

CONCLUSION

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.