Comprehensive analysis of the RSK gene family in acute myeloid leukemia determines a prognostic signature for the prediction of clinical prognosis and treatment responses

ABSTRACT

Objective

The prognosis of acute myeloid leukemia (AML) remains poor although the basic and translational research has been highly productive in understanding the genetics and pathopoiesis of AML and a plethora of targeted therapies have been developed. Consequently, it is crucial to deepen the knowledge of molecular pathogenesis underlying AML for the advancement of new treatment options.

Method

A RSK gene family-related signature was constructed to investigate whether RSK gene family members were useful in predicting the prognosis of AML patients. The relationship between the RSK gene family-related signature and the infiltration of immune cells was further assessed using the CIBERSORT algorithm. The ‘oncoPredict’ package was used to analyze relationships between the RSK gene family-related signature and the sensitivity to drugs or small molecules.

Results

Patients were classified into two groups using the RSK gene family-related signature following the median risk score. Overall survival (OS) was significantly longer in patients with low-risk scores than that in patients with high-risk scores as showed by both training and validation datasets. Moreover, the signature was helpful in predicting 1-year, 3-year, and 5-year OS in training and validation datasets. In addition, it was identified that low-risk patients exhibited greater sensitivity to 20 drugs or small molecules and that high-risk patients had higher sensitivity to 38 drugs or small molecules.

Conclusion

RSK gene family members, particularly RPS6KA1 and RPS6KA4, may help to predict prognosis for AML patients. Furthermore, RPS6KA1 may serve as a novel drug target for AML.

KEYWORDS:

• Acute myeloid leukemia
• RSK gene family
• prognostic signature
• prognosis
• treatment response
• RPS6KA1
• RPS6KA4

Background

Acute myeloid leukemia (AML) is one of the most frequently occurred hematopoietic cancers in adults [1]. During the last 15 years, substantial effects have been yielded through translational and basic studies in understanding the molecular and genetic mechanism behind AML pathogenesis. Consequently, numerous novel and efficacious targeted treatments have been established, like protein kinase small-molecule inhibitors, such as gilteritinib [2], enasidenib [3], and ivosidenib [4], immune checkpoint antibodies, such as gemtuzumab [5], and mitochondrial inhibitors, such as venetoclax [6,7]. AML still has poor outcomes in spite of these advancements, with about 50% of 5-year overall survival (OS) rate in patients before 40 years and 10% of 5-year OS rate in patients older than 65 years [8,9]. Therefore, more knowledge about the molecular mechanism of AML pathogenesis is urgently needed to develop novel therapies [10,11].

The RSK family comprises four human isoforms (RPS6KA1, RPS6KA2, RPS6KA3, and RPS6KA6) and two structurally related homologues (RSP6KA4 and RSP6KA5). RSK gene family members share a high degree of sequence homology (75–80% amino-acid identity) and are uniquely characterized by the presence of two distinct functional kinase domains [12,13].

Arsenic trioxide (As₂O₃) is one of the most efficacious agents in the treatment of acute promyelocytic leukemia (APL) [14,15] but is ineffective against other forms of AML [16]. Inhibition of RPS6KA1 (RSK1) enhances As₂O₃-dependent apoptosis of AML cells, suggesting the potential of RPS6KA1 targeting in combination with As₂O₃ as a novel treatment option for AML [17].

A previous study documented that RPS6KA1 represents a promising target for novel therapeutic strategies targeting therapy-resistant internal tandem duplication mutations in FLT3 juxtamembrane domain (FLT3-ITDs)-positive AML [18]. However, there is no research on the possibility of using RPS6KA1 to predict the OS of AML patients. Likewise, the drug sensitivity of patients with different RPS6KA1 expression has not yet been explored.

On this basis, this study systematically probed into the expression of RSK gene family members and their prognostic values and constructed a prognostic signature in 130 samples with AML from The Cancer Genome Atlas (TCGA) cohorts. Then, two microarray datasets were used to validate the RSK gene family-related prognostic signature, followed by analyses of the association of the signature with clinical outcomes, immune heterogeneity, and drug response features.

Materials and methods

AML dataset source

Gene expression data of tumor specimens and clinical data were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Clinical information, such as survival, gender, and age of patients, was also acquired from the TCGA website. After samples without completed clinical information were removed, this research included 130 tumor tissues.

Microarray datasets GSE37642 [19–22] and GSE12417 [23] were harvested from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). In GSE37642, gene expression data of 136 samples were screened from the platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. In GSE12417, gene expression data of 79 samples were retrieved from the GPL570 platform. This was because only the data retrieved from the GPL570 platform contained the expression of the RSK family.

Data preprocessing

For GEO database-derived microarray datasets, the Robust Multi-array Average approach available through the ‘affy’ package was utilized to obtain background-adjusted, probe-level, and quantile-normalized data-summarized values for probe sets. For the TCGA-LAML dataset-derived high-throughput sequencing, transcripts per million (TPM) values were generated by converting row read count values because TPM values were more comparable between samples and more alike to values produced by microarrays [24].

The association between RSK family genes and clinical features

UALCAN (Birmingham, AL, USA) publicly available at http://ualcan.path.uab.edu [25] was used to evaluate the association between RSK family genes and clinical features including age, gender, FLT3 mutation, PML/RAR-fusion, RAS activation, and race. The 2017 European Leukemia Net (ELN) recommendations for genetic risk stratification of ALM have been widely adopted [26]. The expression of RSK family genes in each group of the ELN risk system was also explored.

Identifying an RSK family-related prognostic signature

Univariate Cox proportional hazards regression analyses were firstly implemented on expression matrix of genes in the RSK family to evaluate the relationship of genes in the RSK family with prognosis (OS) in the TCGA-LAML cohort. Genes in the RSK family with p < 0.1 were regarded as the potential prognosis-associated genes.

In the training step, an optimal prognostic signature with the fewest RSK genes was constructed using the least absolute shrinkage and selection operator (LASSO) penalty in the discovery cohorts. The optimal value of the penalty parameter l was tuned using ten-fold cross-validation, which yielded the minimal partial likelihood deviance. Finally, a prognostic signature related to two RSK genes with non-zero coefficients was identified.

The risk score was computed for the signature in each sample using the equation as follows: $\mathrm{Risk}\mathrm{Score}={\sum }_{(i=1)}^n{Coef}_i\times E_i$, where the Coef_i marks the coefficient and the E_i indicates the normalized expression of each chosen RSK gene with log₂ and Z-score transformations. The median risk score was employed as cutoff to categorize patients into high-risk and low-risk groups. Kaplan-Meier curves were plotted for analyzing survival rates, with the difference measured using the log-rank test. The time-dependent receiver operating characteristic (ROC) curve was computed using the timeROC package.

Differentially expressed genes (DEGs)

DEGs between the predicted high-risk and low-risk groups were determined by DESeq2 packages [27]. Significant DEGs were selected with the criteria of absolute log2FC > 1 and BH-adjusted p value < 0.05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed with clusterProfiler packages [28].

Analysis of drug sensitivity in cancer

The R package ‘oncoPredict’ [29] was used to clarify the difference in drug sensitivity between the high-risk group and the low-risk group. As a public database, the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) [28] can be utilized for exploring cancer molecular therapy and mutation. The gene expression data of cell lines and half-maximal inhibitory concentration (IC50) values of variable anticancer drugs were collected from a website (https://osf.io/c6tfx/files/osfstorage). A ridge regression model was utilized to analyze the difference in the sensitivity to different anticancer drugs between the high-risk group and the low-risk group [29,30].

Protein–protein interaction (PPI)

To explore the function of RSK family genes, we performed PPI analysis by using STRING database [31].

Results

Six RSK genes were adopted to evaluate the relationship of gene expression with prognosis in the TCGA-LAML cohort with a univariate Cox proportional hazards regression model. Only two genes (RPS6KA1 and RPS6KA4) were determined as potential prognosis-associated genes (p < 0.1, Table 1). In addition, the association of these two genes with clinical features such as age, gender, FLT3 mutation, PML/RAR-fusion, RAS activation, the ELN 2017 risk system, and race was analyzed (Figure 1). We found that the expression of RPS6KA1 in adverse and intermediate groups of the ELN 2017 risk system was higher than that in the favorable group.

Figure 1. The association of RPS6KA1/RPS6KA4 expression with age, gender, FLT3 mutation, PML/RAR-fusion, RAS activation, race and the ELN 2017 risk system.

Table 1. The relationship of RSK gene family members with OS in AML patients of TCGA cohorts.

Gene symbol	Hazard ratio	p value
RPS6KA1	1.53 (1.19 - 1.98)	0.0011
RPS6KA2	0.77 (0.57 - 1.06)	0.11
RPS6KA3	1.06 (0.83 - 1.35)	0.65
RPS6KA4	1.41 (1.11 - 1.79)	0.0050
RPS6KA5	0.92 (0.72 - 1.17)	0.49
RPS6KA6	0.96 (0.78 - 1.20)	0.74

Morphologically, AML was classified into eight groups in the French-American-British classification system (FAB M0-M7) [32]. The expression of these two genes was further tested in different stages. Both genes were highly expressed in AM5, without significant differences (Figure 2).

Figure 2. The gene expression of RPS6KA1 and RPS6KA4 in different stages of AML.

To explore the function of RPS6KA1 and RPS6KA4, we performed PPI analysis (Figure 3). The result showed that RPS6KA1 interacted with S100B, TSC2, MAPK1, and CREB1, while RPS6KA4 interacted with CREB1, MAPK1, and MAPK14.

Figure 3. The PPI network displayed the interaction between RSK family genes (RPS6KA1 and RPS6KA4) and other genes.

The LASSO analysis was carried out to construct an optimal prognostic signature containing the fewest RSK genes in the TCGA-LAML cohort. Two genes were observed to have non-zero coefficients. The LASSO-selected features were utilized to construct the RSK family-related signature (Figure 4A). The following equation was adopted to compute the corresponding risk score in all training and validation datasets: Risk Score = 0.222 × RPS6KA1 + 0.103 × RPS6KA4.

Figure 4. Construction of the RSK gene family-related prognostic model. (A) The corresponding hazard ratios of two genes included in the multivariate Cox regression model. (B) Risk score distribution of AML patients in TCGA cohorts, and the heatmap of the survival status of patients and the expression of characteristic genes. (C) The Kaplan-Meier curve showing the difference detected by the log-rank test for patients in TCGA cohorts. (D) The time-dependent ROC curve evaluating the accuracy of the prognostic model in TCGA cohorts.

Based on the calculated risk score, the patients in the three datasets were arranged into high-risk and low-risk groups with their respective median risk scores as cutoff (Figure 4B). As risk scores in both groups increased, the number of living patients decreased (Figure 4B). The expression levels of both RPS6KA1 and RPS6KA4 were higher in patients with high-risk scores than in patients with low-risk scores (Figure 4B). To determine the predictive performance of the signature, the Kaplan-Meier analysis was performed to determine whether OS was different in patients with low-risk scores and patients with high-risk scores. The ROC curve was utilized to assess the predictive accuracy of the signature. In the TCGA-LAML cohort, patients with high-risk scores displayed prominently shorter OS than patients with low-risk scores (p = 0.0072, Figure 4C). Additionally, the area under the ROC curve (AUC) of the signature was 0.70 for 1-year OS and 0.76 for 5-year OS (Figure 4D). The performance of the RSK family-related signature was slightly better than that of the ELN 2017 risk system especially for 1-year OS (Figure S1). In the GSE12417 dataset, OS was conspicuously longer in patients with low-risk scores than in patients with high-risk scores (p = 0.028, Figure 5A), with the AUC of the signature for 1-year OS of 0.66 (Figure 5B). In the GSE37642 dataset, the OS time also evidently varied in different groups (p = 0.0036, Figure 5C), with the AUC of the signature for 1-year OS of 0.66 (Figure 5D).

Figure 5. Evaluation of the accuracy of the prognostic model in validation datasets. (A, C) Kaplan-Meier curves showing the difference detected by the log-rank test for patients in GSE12417 (A) and GSE37642 (C) datasets. (B, D) Time-dependent ROC curves evaluating the accuracy of the prognostic model in GSE12417 (B) and GSE37642 (D) datasets.

Furthermore, univariate and multivariate analyses were performed to assess whether the RSK family-related signature had prognostic significance (Table 2). The results of univariate analysis showed the association of OS with age (HR, 1.03; 95% CI: 1.02-1.05; p value = 2.21 × 10⁻⁵), gender (HR, 1.01; 95% CI: 0.64-1.57; p value = 0.98), the RSK family-related signature (HR, 4.63; 95% CI: 1.99-10.79; p value = 3.88 × 10⁻⁴), and the ELN 2017 risk system (HR, 1.83; 95% CI: 1.30-2.58; p value = 5.21 × 10⁻⁴). The multivariate analyses showed that age (HR, 1.03; 95% CI: 1.01-1.05; p value = 2.13 × 10⁻⁴), the RSK family-related signature (HR, 4.25; 95% CI: 1.80-10.00; p value = 9.30 × 10⁻⁴), and ELN 2017 risk system (HR, 1.49; 95% CI: 1.03-2.16; p value = 0.035) remained independent prognostic factors for OS.

Table 2. Uni – and multi-variate Cox regression analysis for overall survival of TCGA-LAML cohort.

Factors	Univariate		Multivariate
	HR (95% CI)	p value	HR (95% CI)	p value
Age	1.03 (1.02, 1.05)	2.21 × 10^−5	1.03 (1.01, 1.05)	2.13 × 10^−4
Gender	1.01 (0.64, 1.57)	0.98
risk score	4.63 (1.99, 10.79)	3.88 × 10^−4	4.25 (1.80, 10.00)	9.30 × 10^−4
ELN 2017 risk	1.83 (1.30, 2.58)	5.21 × 10^−4	1.49 (1.03, 2.16)	0.035

A nomogram plot that integrated risk scores, age, gender, and the ELN risk system was constructed (Figure 6A) to more accurately assess the signature. The calibration curves showed that the signature could accurately predict the utility of 1-year, 3-year, and 5-year OS in AML patients (Figure 6B).

Figure 6. Nomogram (A) and the calibration curve of nomogram (B) showing clinicopathological characteristics, such as age, gender, and risk scores and manifesting favorable predictive accuracy for 1-, 2-, and 3-year survival rates.

The DESeq2 package was adopted to screen out DEGs between the predicted high-risk and low-risk groups. With the criteria of absolute log2FC > 1 and BH-adjusted p value < 0.05, 1,291 DEGs were identified including 706 significantly upregulated genes and 585 significantly downregulated genes (Figure 7A). Afterwards, DEGs were mainly enriched in cytokine-cytokine receptor interaction and chemokine signaling pathway, as revealed by KEGG functional analyses (Figure 7B). The result of GO analyses manifested that DEGs were mostly enriched in neutrophil degranulation and neutrophil activation involved in immune response (Figure 7C).

Figure 7. Functional analysis. (A) The volcano plot of the 1,291 significant DEGs between the high-risk group and the low-risk group. The enrichment analysis showed the enriched KEGG pathways (B) and enriched GO terms (C).

To clarify the relationship between the RSK gene family-related signature and the sensitivity of drugs or small molecules, the gene expression data of 809 cell lines and the IC50 value of 198 anticancer drugs or small molecules were downloaded from the GDSC database. The IC50 values of each patient in TCGA cohorts were computed in the light of the gene expression of cell lines to 198 anticancer drugs or small molecules. Afterwards, the difference in IC50 values between patients with low-risk scores and patients with high-risk scores was analyzed (Figure 8A, B). The total IC50 values of 20 anticancer drugs or small molecules in patients with high-risk scores were considerably higher than those in patients with low-risk scores subsequent to Bonferroni correction (p < 2.5 × 10⁻⁴, 0.05/198, Figure 8A). Furthermore, patients in the high-risk group were sensitive to 38 anticancer drugs or small molecules (Figure 8B).

Figure 8. Evaluation of the sensitivity to different drugs or small molecules with risk scores of the RSK gene family according to the GDSC database. (A) Drugs or small molecules that patients in the low-risk group were more sensitive to. (B) Drugs or small molecules that patients in the high-risk group were more sensitive to.

Discussion

In this study, the RSK gene family-related signature was constructed to explore whether the RSK gene family assisted in predicting the prognosis of AML patients. Based on the median risk score, the patients were assigned to low-risk and high-risk groups. Low-risk patients had markedly longer OS than high-risk patients in both validation and training datasets. The results manifested that the signature was useful in predicting 1-year, 3-year, and 5-year OS in validation and training datasets. The high-risk score illustrated that there existed a higher percentage of immune cell infiltration.

The obtained RSK gene family-related signature was employed to identify the difference in OS between the low-risk group and the high-risk group. The result was evaluated in training and validation datasets (Figure 4C, Figure 5A, C), which illustrated the usefulness of the signature for predicting the prognosis of AML patients. The performance of the RSK family-related signature was slightly better than that of the ELN 2017 risk system especially for 1-year OS. The multivariate analysis results also indicated that the RSK family-related signature could be used to further stratify AML patients on the basis of the ELN 2017 risk system. Thereafter, the nomogram plot also proved the accuracy of the RSK gene family-related signature combined with clinical risk factors and the ELN 2017 risk system to predict 1-year, 3-year, and 5-year OS (Figure 6).

RPS6KA1 and RPS6KA4 expression was positively linked to risk scores, and the survival rate was poor in the high-risk group, demonstrating RPS6KA1 and RPS6KA4 as unfavorable factors. This result was consistent with a previous finding that RPS6KA1 (RSK1) was able to enhance the survival and proliferation of AML cells [18]. FLT3-ITD is the most frequent tyrosine kinase mutation in AML associated with poor prognosis. FLT3-ITD activates RPS6KA1 through activation of the MEK/ERK pathway [18]. Phosphorylated RPS6KA1 and PIM can synergistically activate the mTORC1/S6 K/4EBP1 pathway and eIF4B [18]. RPS6KA1 downregulation increased the sensitivity to BH3 mimetics and induced the activation of Bax which led to AML cell apoptosis [33,34]. However, the function of RPS6KA4 in the survival of AML cells was unknown. A previous study showed that the overexpression of element-binding protein (CREB1) increased AML cell growth by driving the expression of key regulators of apoptosis [35]. The PPI results suggested that RPS6KA4 may enhance the survival and proliferation of AML cells by interacting with and phosphorylating CREB1.

Additionally, the current study identified 20 drugs or small molecules that low-risk patients exhibited greater sensitivity to and 38 drugs or small molecules that high-risk patients demonstrated more sensitivity to. The drugs or small molecules that high-risk patients had greater sensitivity to may be beneficial to patients with high-risk scores. Several of them have been reported to inhibit the survival of AML cells. For example, previous studies showed that Aurora kinase A inhibitor Alisertib and Aurora kinase B inhibitor Barasertib strongly inhibited the growth and proliferation of AML cells. ADC2014, a mTORC1/2 dual inhibitor, reduced pH in the acidic organelles, including lysosomes, and synergistically enhanced the cytotoxic effect of Gemtuzumab ozogamicin in primary leukemia cells [36]. The combine of PIM and AKT kinase inhibitors (AZD5363) showed synergistic cytotoxicity in AML [37]. Other drugs (like camptothecin, dactinomycin and dactolisib) were also verified to be useful for treating AML [38,39]. Further research is warranted to delve into the mechanism and availability of other drugs or small molecules in the clinical treatment of AML. These drugs may provide potential clues for precise treatment for AML patients with different risk scores.

In summary, the current study obtained two RSK genes linked to prognosis in AML patients. The constructed signature was helpful in predicting 1-year, 3-year, and 5-year OS. Of note, it was also confirmed that RPS6KA1 might be a novel drug target for the treatment of AML. The identified drugs or small molecules might provide new ideas for the therapy of AML patients with high expression of RPS6KA1.

Acknowledgements

Thanks to all the contributors.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.

Additional information

Funding

This research was supported by Science and Technology Foundation of Guizhou Provincial Health Commission in 2023 (No. gzwkj2023-412).