3.1 Gene expression analysis of RUNX
To compare expression levels of RUNX family members between tumor and normal tissues, the mRNA levels of RUNX family members in tumors and matched normal tissues were analyzed via the Oncomine database. Our results showed relatively high expression of RUNX1 and RUNX2 in breast cancer, head and neck cancer, kidney cancer, leukemia, and pancreatic cancer. In addition, RUNX1 was also relatively higher in brain and central nervous system cancers, colorectal cancer, and sarcoma. The expression of RUNX3 was relatively higher in esophageal cancer, head and neck cancers, kidney cancer, lymphoma, and sarcoma (Figure 1A).
To further evaluate the differential expression of RUNX family members, we profiled and compared their expression across all TCGA tumors using the R software. Specifically, our results showed that RUNX1 expression was significantly elevated in BLCA, CHOL, COAD, ESCA, GBM, HNSC, KIRC, KIRP, LIHC, LUAD, READ, STAD, THCA, and UCEC. By contrast, lower RUNX1 expression was found in PARD (Figure 1B). RUNX2 showed higher expression in BLCA, BRCA, CHOL, ESCA, GBM, HNSC, KIRC, KIRP, LUAD, LUSC, STAD, and THCA but decreased expression in PARD (Figure 1C). RUNX3 expression was significantly elevated in ESCA, GBM, HNSC, KIRC, KIRP, LUAD, LUSC, STAD, and UCEC, but lower RUNX3 expression was found in BRCA, COAD, LIHC, LUAD, and THCA (Figure 1D).
We also analyzed the pan-cancer and inter-tumor heterogeneity of RUNX family expression; the results showed that RUNX1 and RUNX3 were highly expressed whereas RUNX2 was moderately expressed at the pan-cancer level (Figure 2A). As shown in Figure 2B, we observed significant heterogeneity in gene expression of RUNX gene family members across different tumor types in TCGA. Then, we analyzed the correlations among RUNX family members; RUNX1 and RUNX2 were the two genes with the most significant positive correlation (correlation coefficient = 0.66, Figure 2C).
3.2 Prognostic value of RUNX
In order to investigate the correlation of RUNX family expression with prognosis of patients across different tumor types, we performed Cox analysis. According to our results, all RUNX family members showed different expression associated with 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). In addition, higher expression of RUNX3 was associated with poor survival outcomes in LGG and COAD (HR>1, p<0.05, Figure 3). Detailed results of the Cox analysis of the RUNX family across the TCGA database are summarized in Supplementary Table S1.
In order to understand the prognostic value of the RUNX gene family, we used Kaplan–Meier survival curves to determine patient OS based on the clinical information downloaded from TCGA. Of the three family members, RUNX1 predicted the best prognosis in patients with BRCA, ESCA, and SKCM (all p<0.05, Figure 4A-C). By contrast, RUNX1 predicted poor prognosis in patients with KIRC, KIRP, LGG, MESO, OV, and UVM (all p<0.05, Figure 4D-I). RUNX2 had a protective role in PCPG and SKCM (all p<0.05, Figure 4J-K). On the other hand, RUNX2 had a detrimental role in seven cancer types: BLCA, KIRC, KICH, LGG, MESO, SARC, and UVM (all p<0.05, Figure 4L-R). RUNX3 had a protective role in BRCA and ESCA (all p<0.05, Figure 4S-T) but a detrimental role in LGG (all p<0.05, Figure 4U).
In addition, we analyzed the relationship of the RUNX gene family expression and DSS in pan-cancer. The results suggested that increased expression of RUNX1 was associated with poor DSS in GBM, LGG, KIRC, OV, MESO, and UVM, whereas, increased RUNX1 expression predicted good DSS in SKCM, THYM, and BRCA (all p<0.05, Figure 5A-I). RUNX2 had a detrimental role in KICH, KIRC, LGG, MESO, and UVM (all p<0.05, Figure 5J-N). RUNX3 had a detrimental role in LGG (p<0.05, Figure 5O).
The same method is used to analyze DFI in 33 TCGA tumors. RUNX1 had a protective role in BRCA (p<0.05, Figure 6A). On the other hand, RUNX1 had a detrimental role in PAAD and CESC (all p<0.05, Figure 6B-C). The high expression of RUNX2 was a protective factor in LIHC, PARD, and TGCT (all p<0.05, Figure 6D-F), but was all a risk factor in the LUAD (all p<0.05, Figure 6G). RUNX3 had a protective role in LIHC and UCS (p<0.05, Figure 6H-I).
Finally, we also analyzed the PFI in 33 TCGA tumors. RUNX1 expression had a detrimental role in CESC, COAD, GBM, KIRC, LGG, and UVM (all p<0.05, Figure 7A-F), but was all a protective role in the BRCA and SKCM (all p<0.05, Figure 7G-H). 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-L). RUNX3 had a detrimental role in GBM, LGG, PARD, and THYM (p<0.05, Figure 7M-P).
We also explored the relationship of RUNX gene family expression and clinicopathologic stage. Our results demonstrated that the RUNX gene family expression levels in several tumor tissues were significantly different in different clinical stages. With the increase of tumor grade, the expression of RUNX1 increased in the UVM, and the expression of RUNX2 increased in the STAD, ESCA,and KICH. In contrast, as the tumor grade increased, the expression of RUNX1 decreased in BLCA and BRCA, and the expression of RUNX3 decreased in TGCT (Supplementary Figure 1).
3.3 Correlation of RUNX 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[18, 19]. In this study, we explored the association between the RUNX gene family and the tumor microenvironment by using the ESTIMATE algorithm to calculate immune and stromal scores and tumor purity in diverse cancer types from the TCGA database. After a series of analyses, we observed a positive pan-cancer association between the RUNX gene family and ESTIMATE scores (Figure 8A-C). This suggests that the three members of the RUNX gene family show similar role 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).
Furthermore, we analyzed the correlation between TMB or MSI and the RUNX gene family expression in different types of cancers. Of the three family members, we found that RUNX1 expression was significantly associated with TMB in MESO, LUSC, LIHC, LGG, LAML, KIRP, HNSC, ESCA, CESC, BRCA, and BLCA (Supplementary Figure 2A). We found that gene expression was significantly correlated with MSI in LUSC, LUAD, KIRC, HNSC, DLBC, and CHOL (Supplementary Figure 2B). A significant correlation was found between RUNX2 expression and TMB in various cancers, such as LIHC, KIRP, KICH, HNSC, ESCA, CHOL, and BRCA (Supplementary Figure 2C), while, RUNX2 expression was correlated with MSI in LIHC, HNSC, DLBC, and COAD (Supplementary Figure 2D). We found that RUNX3 expression was significantly associated with TMB in PAAD, LUAD, LIHC, LGG, KIRP, ESCA, and BRCA (Supplementary Figure 2E). We also found that gene expression was significantly correlated with MSI in OV, LGG, KIRP, ESCA, DLBC, COAD, and BRCA (Supplementary Figure 2F).
3.4 Correlation of RUNX family expression with tumor immune subtype and immune checkpoints
Tumor-infiltrating immune cells are crucial components of the tumor microenvironment and promote the immune escape of cancer. To clarify the molecular significance of the RUNX gene family in tumor-infiltrating immune cells, we analyzed the potential relationship between expression of the RUNX family and tumor immune subtypes. The members of the RUNX family were all related to immune subtypes on a pan-cancer basis. RUNX1 had higher expression in C1, C2, C3, C4, and C6. Higher levels of RUNX2 were associated with C1, C2, and C6. RUNX3 showed higher expression in C1, C2, C3, and C6 (Figure 9A). We also analyzed the correlations between RUNX gene family expression and 47 common immune checkpoint genes, across all cancer types in the TCGA database. As shown in Figure 9B, in most cancers, except CESC, CHOL, MESO, LAML, and UCS, RUNX1 expression was significantly correlated with immunoinhibitory genes. Additionally, except CESC, MESO, and UCS, RUNX2 expression was significantly correlated with immunostimulatory genes (Figure 9C). RUNX3 expression was significantly correlated with immunostimulatory genes, except UCS (Figure 9D). This suggested a potential synergy of RUNX gene family with known immune checkpoints.
3.5 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-B, G). The highest alteration frequency of RUNX2 (>6%) appeared for patients with esophageal adenocarcinoma with primary type “amplification” (Figure 10C-D). The hot spot mutation S31Pfs*9 of RUNX2 occurred in five cancers (STAD and COAD) in five patients and also encoded a truncated protein (Figure 10H). Compared with RUNX1 and RUNX2, the genetic alteration level of RUNX3 was relatively low. Specifically, SKCM patients had the highest alteration frequency of RUNX3, with “mutation” as the primary type (Figure 10E-F, I).
3.6 Drug response analysis of RUNX gene family
In order to investigate the potential relationship between the RUNX gene family and drug response in various human cancer types, we performed a correlation analysis to screen out the association between RUNX gene family and drug response from the CellMiner database. In this study, Z-score was used to measure drug response. Specifically, RUNX1 expression was positively associated with drug response to batracylin, tricirbine phosphate, nelarabine, chelerythrine, irinotecan, bendamustine, clofarabine, and asparaginase (Figure 11B-D, F-H, N, O). RUNX2 expression was positively associated with response to staurosporine, bleomycin, and dasatinib (Figure 11L, M, P). RUNX3 expression was positively associated with hypothemycin, selumetinib, PD-98059, dabrafenib, and dasatinib (Figure 11A, E, I-K).
3.7 Functional enrichment analysis of the RUNX gene family
To perform a comprehensive functional and pathway analysis of the RUNX gene family, we used GeneMANIA to screen out other genes or transcription factors associated with the RUNX gene family. The results showed that several transcription factors (Figure 12) were potential targets of RUNX family members. In addition, according to the analysis, RUNX members and related transcription factors were involved in the following pathways or functions: chromatin, regulatory region DNA binding, regulatory region nucleic acid binding, histone deacetylase binding, transcriptional repressor complex, and activating transcription factor binding.
3.7 RUNX gene family in UVM
We also studied the relationship between the expression level of the RUNX gene family and the ESTIMATE score and immune subtype. The expression of RUNX gene family members were all related to immune subtypes in UVM. Notably, only C3, C4, and C5 immune subtypes were observed in the RUNX gene family (Figure 13A). We also conducted a correlation analysis, which showed that expression of RUNX2 was negatively associated with RNAss in UVM. The expression levels of RUNX1 and RUNX2 were positively associated with DNAss in UVM. Besides, for all three members of the RUNX gene family, expression levels were positively associated with stromal score, immune score, and ESTIMATE score (Figure 13B). Using the data from HPA cohort, we also analyzed the differential protein expression of RUNX family. Our result shown that the aberrant expression of RUNX1 was detected in metastatic melanoma tissues, when compared normal tissues. However, there was no significantly expression different of RUNX2 and RUNX3 between metastatic melanoma tissues and normal tissues (Figure 13C).