Prognostic value and immune characteristics of RUNX gene family in human cancers: a pan-cancer analysis

TCGA data and processing

The TCGA database (http://cancergenome.nih.gov) comprises the results of a groundbreaking cancer genomics study that has molecularly described over 20,000 primary cancer and normal samples encompassing 33 cancer types up to 2021. The TCGA data included 33 different tumor types: adrenocortical carcinoma (ACC); bladder urothelial carcinoma (BLCA); breast cancer (BRCA); cervical squamous cell carcinoma and endocervical adenocarcinoma; cholangiocarcinoma (CHOL); colon adenocarcinoma (COAD); lymphoid neoplasm diffuse large B-cell lymphoma; esophageal carcinoma (ESCA); glioblastoma multiforme (GBM); head and neck squamous carcinoma (HNSC); kidney chromophobe (KICH); kidney renal clear cell carcinoma (KIRC); kidney renal papillary cell carcinoma (KIRP); acute myeloid leukemia (LAML); brain lower grade glioma (LGG); liver hepatocellular carcinoma (LIHC); lung adenocarcinoma (LUAD); lung squamous cell carcinoma (LUSC); mesothelioma (MESO); ovarian serous cystadenocarcinoma (OV); pancreatic adenocarcinoma (PAAD); pheochromocytoma and paraganglioma (PCPG); prostate adenocarcinoma; rectum adenocarcinoma (READ); sarcoma (SARC); skin cutaneous melanoma (SKCM); stomach adenocarcinoma (STAD); testicular germ cell tumors; thyroid carcinoma (THCA); thymoma (THYM); uterine corpus endometrial carcinoma (UCEC); uterine carcinosarcoma; and uveal melanoma (UVM). The UCSC Xena (https://xenabrowser.net/) website was used to collect RUNX gene family data from TCGA, including RNA sequencing results, clinical data, stemness scores, DNA methylation, and immune subtype data.

Gene expression data with immunotherapy

To investigate the predictive value of the RUNX gene family, two separate datasets, GSE78220 and IMvigor210, were obtained and examined. GSE78220 (platform GPL11154 Illumina HiSeq 2000) gene expression profiles and clinical data were acquired from GEO (https://www.ncbi.nlm.nih.gov/geo/) databases, which was chosen for this study as the anti-PD-1 immunotherapy cohort. The IMvigor210 dataset was obtained under the Creative Commons 3.0 license from a freely available, fully documented software and data package available, which was chosen for this study as the anti-PD-L1 immunotherapy cohorts.

RUNX gene family expression analysis

In the TCGA, we utilized the “ggpubr” R package to compare the expression of RUNX gene family members in cancers and matching normal tissues. Gene expression was standardized in the pan-tumor study. Transcripts per kilobase million values were created by converting fragments per kilobase million values. Using the “pheatmap” R package, the findings were shown as heatmaps.

The Oncomine database (http://www.oncomine.org) has 715 datasets and 86,733 samples for genome-wide expression analysis [15]. We used Oncomine to compare the expression levels of RUNX gene family members in cancer tissues and neighboring normal tissues. In our research, we employed the student’s t-test to compare the expression of RUNX gene family members in different forms of cancer; a p-value of 0.05, a fold change of 2, and a gene rank in the top 10% were set as the significance thresholds.

Survival prognosis analysis

The link of the RUNX gene family with overall survival (OS), disease-specific survival (DSS), progression-free interval (PFI), and disease-free interval (DFI) across all TCGA cancers was investigated using Kaplan–Meier curves. With 95 percent confidence intervals, we calculated log-rank p-values and hazard ratios (HR) (95 percent CI). Forest plots (created with the “forestplot” R package) and survival curves were used to visualize the data.

Mutation profiles

The cBioPortal (https://www.cbioportal.org) website is an interactive tool for exploring, visualizing, and analyzing multidimensional cancer genomics data [16]. On a pan-cancer basis, we investigated the copy number changes and genetic modification features of the RUNX gene family using the “Cancer Types Summary” module. The “Cancer Types Summary” module of cBioPortal was used to acquire the alteration frequency, mutation type, and copy number alteration of the RUNX gene family across all tumors in TCGA. To query the genetic modification features of the RUNX gene family, we set “TCGA Pan Cancer Atlas Studies” in the “Quick pick” area.

Tumor microenvironment and tumor stemness analysis

The “estimate” R package was used to measure tumor purity in 33 human malignancies from the TCGA. The ESTIMATE score is made up of the immunological and stromal scores, which indicate the amount of immune and stromal components in the body, respectively. Lower tumor purity is associated with higher ESTIMATE scores. C1 (wound healing), C2 [interferon (IFN)-r dominant], C3 (inflammatory), C4 (lymphocyte depletion), C5 (immunologically quiet), and C6 [tumor growth factor b (TGFb) dominant] are the six immune subtypes used to quantify immune infiltration in the tumor environment. The correlations between the expression of RUNX gene family members and the ESTIMATE immune and stromal scores were discovered using Spearman’s correlation analysis. We also used epigenetic and transcriptome data to calculate tumor RNA stemness scores (RNAss) and DNA stemness scores (DNAss). Specifically, RNAss is calculated using RNA sequencing data, whereas DNAss is calculated using DNA methylation data. The associations between the RUNX gene family and tumor stemness data were investigated using Spearman’s correlation analysis. The relationships between TGFBI expression and immunoinhibitory and immunostimulatory gene subsets, tumor mutation burden (TMB), and microsatellite instability were investigated using Pearson’s correlation analysis (MSI). Using R’s “pheatmap” package, the findings were visualized as heatmaps.

Drug responses

The RUNX gene family medication response information were gathered from the National Cancer Institute (NCI)-60 database through the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do). The NCI-60 database comprises molecular and pharmacological information for 60 different human cancer cell lines. Pearson’s correlation analysis was used to assess the relationship between RUNX gene family expression (as measured by transcript levels) and drug response (GI50). All of the medications used in the correlation study were identified by the Food and Drug Administration or through drug clinical studies.

Statistical analyses

Oncomine was used to perform correlation analysis for the expression of RUNX gene family members between cancer and neighboring tissues, including p-values, fold changes, and gene rankings. R version 4.0.4, 64-bit (https://www.r-project.org/) and its associated packages were utilized in all studies. The outcomes of survival are displayed with HR, 95% CI, and log-rank p-values. For all statistical analyses, p<0.05 was considered statistically significant.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: TCGA database (http://cancergenome.nih.gov), UCSC Xena (https://xenabrowser.net/), Oncomine database (http://www.oncomine.org), the National Cancer Institute (NCI)-60 database via the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do), and GeneMANIA (http://www.genemania.org).

Gene expression analysis of RUNX

The mRNA levels of RUNX gene family members in malignancies and matched normal tissues were evaluated using the Oncomine database to compare expression levels between tumor and normal tissues. RUNX1 and RUNX2 expression was found to be quite high in head and neck cancer, kidney cancer, leukemia, and pancreatic cancer, according to our findings. RUNX1 was also shown to be greater in malignancies of the brain and central nervous system, colorectal cancer, and sarcoma. RUNX3 expression was shown to be greater in esophageal cancer, head and neck cancers, kidney cancer, lymphoma, and sarcoma than in other malignancies (Figure 1A).

Using the R program, we profiled and compared the expression of RUNX gene family members across all TCGA cancers to further analyze their differential expression. RUNX1 expression was considerably increased in BLCA, CHOL, COAD, ESCA, GBM, HNSC, LIHC, LUAD, READ, KIRC, KIRP, STAD, THCA, and UCEC, according to our findings. In prostate adenocarcinoma (PARD), however, RUNX1 expression was shown to be decreased (Figure 1B). RUNX2 expression was shown to be greater in BLCA, BRCA, CHOL, ESCA, GBM, HNSC, KIRC, KIRP, LUAD, LUSC, STAD, and THCA, but lower in PARD (Figure 1C). RUNX3 expression was shown to be substantially higher in ESCA, GBM, LUSC, HNSC, KIRC, KIRP, STAD, and UCEC than in BRCA, COAD, LIHC, LUAD, and THCA (Figure 1D).

We also investigated at the pan-cancer and inter-tumor heterogeneity of RUNX gene family expression, and found that RUNX1 and RUNX3 were highly expressed at the pan-cancer level, whereas RUNX2 was modestly expressed (Figure 2A). We found considerable variation in gene expression of RUNX gene family members across various tumor types in TCGA, as illustrated in Figure 2B. Then we looked at the relationships between members of the RUNX gene family; RUNX1 and RUNX2 had the strongest positive association (correlation coefficient = 0.66, Figure 2C).

Prognostic value of RUNX

We used Cox analysis to look at the relationship between RUNX gene family expression and patient prognosis across different tumor types. According to our results, all RUNX gene family members showed different expressions associated with the prognosis of patients. Specifically, RUNX1 had a detrimental role in UVM, LGG, MESO, KIRC, PAAD, GBM, KIRP, and OV (HR>1, p<0.05); however, it had a protective role in SKCM, LUAD, BRCA, ESCA, and THYM (HR<1, p<0.05). RUNX2 was detrimental in UVM, KICH, ACC, LGG, KIRC, BLCA, PAAD, MESO, GBM, and SARC (HR>1, p<0.05) but protective in SKCM (HR<1, p<0.05). Furthermore, greater RUNX3 expression was linked to worse survival outcomes in LGG and COAD (HR>1, p<0.05, Figure 3). Detailed results of the Cox analysis of the RUNX gene family across the TCGA database are summarized in Supplementary Table 1.

We utilized Kaplan–Meier survival curves to assess patient OS based on clinical data retrieved from TCGA in order to better understand the RUNX gene family’s predictive significance. Of the three family members, RUNX1 predicted the best prognosis in patients with BRCA, ESCA, and SKCM (all p<0.05, Figure 4A–4C). By contrast, in patients with KIRC, KIRP, LGG, MESO, OV, and UVM, RUNX1 indicated a poor prognosis (all p<0.05, Figure 4D–4I). In PCPG and SKCM, RUNX2 played a protective effect (all p<0.05, Figure 4J, 4K). On the other hand, RUNX2 was found to be harmful in seven cancer types: BLCA, KIRC, KICH, LGG, MESO, SARC, and UVM (all p<0.05, Figure 4L–4R). RUNX3 was shown to be beneficial in BRCA and ESCA (all p<0.05, Figure 4S, 4T) but a detrimental role in LGG (all p<0.05, Figure 4U).

We also examined the correlation between RUNX gene expression and DSS in pan-cancer. Higher RUNX1 expression was linked with poor DSS in GBM, LGG, KIRC, OV, MESO, and UVM, whereas increased RUNX1 expression predicted favorable DSS in SKCM, THYM, and BRCA (all p<0.05, Figure 5A–5I). RUNX2 played a negative influence in KICH, KIRC, LGG, MESO, and UVM (all p<0.05, Figure 5J–5N). RUNX3 was shown to be harmful in LGG (p<0.05, Figure 5O).

DFI in 33 TCGA tumors is analyzed using the same way. RUNX1 played a protective function in BRCA (p<0.05, Figure 6A). On the other hand, RUNX1 had a negative impact in PAAD and CESC (all p<0.05, Figure 6B, 6C). In LIHC, PARD, and TGCT (all p<0.05, Figure 6D–6F), high RUNX2 expression was a protective factor, while it was all a risk factor in the LUAD (all p<0.05, Figure 6G). In LIHC and UCS, RUNX3 played a protective function (p<0.05, Figure 6H, 6I).

Finally, the PFI was examined in 33 TCGA cancers. RUNX1 expression had a detrimental role in CESC, COAD, GBM, KIRC, LGG, and UVM (all p<0.05, Figure 7A–7F), but was all a protective role in the BRCA and SKCM (all p<0.05, Figure 7G, 7H). The high expression of RUNX2 was a protective factor in LIHC (all p<0.05, Figure 7I), but was all a risk factor in the KICH, LGG, and UVM (all p<0.05, Figure 7J–7L). RUNX3 had a detrimental role in GBM, LGG, PARD, and THYM (p<0.05, Figure 7M–7P).

We also examined the correlation between the expression of the RUNX gene family and the clinicopathologic stage. The expression levels of the RUNX gene family in multiple tumor tissues were considerably varied in different clinical phases, according to our findings. RUNX1 expression rose in the UVM as tumor grade increased, while RUNX2 expression increased in the STAD, ESCA, and KICH. In contrast, when tumor grade grew, RUNX1 expression in BLCA and BRCA declined, and RUNX3 expression in TGCT decreased (Supplementary Figure 1A–1G).

Correlation of RUNX gene family expression with tumor microenvironment and tumor stemness

The tumor microenvironment plays a crucial part in the initiation, progression metastasis, and drug resistance of cancer [17, 18]. We used the ESTIMATE method to generate immune and stromal scores, as well as tumor purity, in a variety of cancer types from the TCGA database to investigate the relationship between the RUNX gene family and the tumor microenvironment. We discovered a favorable pan-cancer correlation between the RUNX gene family and ESTIMATE scores after a series of investigations (Figure 8A–8C). This shows that the RUNX gene family’s three members have comparable roles in the tumor microenvironment.

In addition, we assessed the pan-cancer correlations between the RUNX gene family and tumor stemness. Using the TCGA tumor stemness database, RNAss and DNAss were analyzed. Our results showed that the RUNX gene family had a negative correlation with RNAss (Figure 8D). By contrast, it had a positive correlation with DNAss in CHOL, KIRP, LGG, MESO, THCA, and UVM (Figure 8E).

We also analyzed the relationship between TMB or MSI and the RUNX gene family expression in various malignancies. RUNX1 expression was strongly related with TMB in MESO, LUSC, LIHC, LGG, LAML, KIRP, HNSC, ESCA, CESC, BRCA, and BLCA, among the three family members (Supplementary Figure 2A). In LUSC, LUAD, KIRC, HNSC, DLBC, and CHOL, we discovered that gene expression was substantially linked with MSI (Supplementary Figure 2B). RUNX2 expression was shown to be connected with TMB in many malignancies, including LIHC, KIRP, KICH, HNSC, ESCA, CHOL, and BRCA (Supplementary Figure 2C), whereas RUNX2 expression was found to be correlated with MSI in LIHC, HNSC, DLBC, and COAD (Supplementary Figure 2D) (Supplementary Figure 2D). In PAAD, LUAD, LIHC, LGG, KIRP, ESCA, and BRCA, RUNX3 expression was shown to be substantially linked with TMB (Supplementary Figure 2E). Gene expression was also shown to be strongly linked with MSI in OV, LGG, KIRP, ESCA, DLBC, COAD, and BRCA in OV, LGG, KIRP, ESCA, DLBC, COAD, and BRCA (Supplementary Figure 2F).

Correlation of RUNX gene family expression with tumor immune subtype and immune checkpoints

Tumor-infiltrating immune cells are critical elements of the tumor microenvironment and assist malignancy elude the immune system. We investigated the possible association between the expression of the RUNX gene family and tumor immune subtypes to better understand the molecular signature of the RUNX gene family in tumor-infiltrating immune cells. On a pan-cancer basis, members of the RUNX gene family were all linked to immunological subtypes. RUNX1 expression was greater in C1, C2, C3, C4, and C6. RUNX2 levels were found to be higher in C1, C2, and C6. RUNX3 was found to be more expressed in C1, C2, C3, and C6 cells (Figure 9A). We also looked at the relationships between the RUNX gene family and 47 frequent immune checkpoint genes in the TCGA database across all cancer types. RUNX1 expression was substantially linked with immunoinhibitory genes in most malignancies, except CESC, CHOL, MESO, LAML, and UCS, as shown in Figure 9B. RUNX2 expression was also highly linked with immunostimulatory genes, with the exception of CESC, MESO, and UCS (Figure 9C). Except for UCS, RUNX3 expression was substantially linked with immunostimulatory genes (Figure 9D). This showed a possible RUNX gene family synergy with known immunological checkpoints.

Genetic alteration analysis of RUNX gene family

Using the cBioPortal database (10967 samples from 32 studies), we analyzed the genetic alterations of the RUNX gene family in different cancer types. Our results show that LAML had a relatively high mutation level, with RUNX1 alterations exceeding 13% with “mutation” as the primary type; the hot spot mutation of RUNX1 was D96Gfs*15/Gfs*11/Mfs*10 in the Runt domain, which occurred in nine cancers (LAML and BRCA) in nine patients and resulted in a truncated protein (Figure 10A–10C). The highest alteration frequency of RUNX2 (>6%) appeared for patients with esophageal adenocarcinoma with primary type “amplification” (Figure 10D, 10E). The RUNX2 hotspot mutation S31Pfs*9 resulted in a shortened protein and was found in five malignancies (STAD and COAD) in five individuals (Figure 10F). RUNX3 has a modest amount of genetic change when compared to RUNX1 and RUNX2. Patients with SKCM had the highest prevalence of RUNX3 mutations, with “mutation” being the most common form (Figure 10G–10I).

Drug response analysis of RUNX gene family

We used a correlation analysis to screen out the association between the RUNX gene family and drug response in various human cancer types from the CellMiner database to study the potential relationship between the RUNX gene family and drug response in various human cancer types. The Z-score was employed to assess medication response in this investigation. RUNX1 expression was found to be linked to a good response to batracylin, tricirbine phosphate, nelarabine, chelerythrine, irinotecan, bendamustine, clofarabine, and asparaginase (Figure 11A–11H). The expression of RUNX2 was related to a good response to staurosporine, bleomycin, and dasatinib (Figure 11I–11K). Hypothemycin, selumetinib, PD-98059, dabrafenib, and dasatinib were all linked to RUNX3 expression (Figure 11L–11P).

Functional enrichment analysis of the RUNX gene family

We utilized GeneMANIA to screen out other genes or transcription factors linked with the RUNX gene family in order to do a full functional and pathway analysis of the RUNX gene family. Several transcription factors (Figure 12) were identified as possible targets of RUNX gene family members. RUNX members and related transcription factors were also found to be involved in the following pathways or functions: chromatin, regulatory region DNA binding, regulatory region nucleic acid binding, histone deacetylase binding, transcriptional repressor complex binding, and activating transcription factor binding, according to the study.

RUNX gene family in immunotherapy

We also explored the connection between RUNX gene family expression and therapeutic benefit in patients who had anti-PD-1 and anti-PD-L1 immunotherapy. In the GSE78220 cohort, the low expression of RUNX2 group had a greater objective response rate than the high expression of RUNX2 group (p = 0.025; Figure 12A). The IMvigor cohort, which received anti-PD-L1 treatment, had a similar result (p = 0.075; Figure 12B). In the IMvigor210 cohort, we also discovered that decreased RUNX1 expression is associated to objective response to anti-PD-L1 treatment (p = 0.015; Figure 12C). Overall, these findings suggest that ICI scores may be linked to immunotherapy responsiveness.