The Oncomine database was used to compare the expression levels of the three RUNX members across a variety of cancer types. The p-value was set as 1E−4, and the fold change was set as 2.

The RUNX expression profile was analyzed using the Gene Expression Profiling Interactive Analysis (GEPIA) tool, which was developed based on The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases (Tang et al., 2017). The RUNX expression profiles were compared between different breast cancer stages.

The UALCAN portal (Chandrashekar et al., 2017) was used to analyze how the clinical characteristics of breast cancer patients were related to RUNX expression profiles and RUNX promoter methylation status. The pan-cancer analysis of RUNX factor expression was also carried out using the UALCAN web.

The relationship between RUNX expression levels and prognosis of breast cancer patients, was analyzed using the Kaplan–Meier plotter (KM plotter; Gyorffy et al., 2010). A total of 2,032 patients were analyzed, and they were split according to the median expression level of RUNX transcription factors. The two groups were compared to determine OS, distant metastasis-free survival (DMFS), relapse-free survival (RFS), and progression-free survival (PFS).

The alterations of RUNX members in breast cancer subtypes were analyzed using the Breast Invasive Carcinoma (TCGA, PanCancer Atlas) dataset, which includes information from 1,084 samples. The alteration types include mutation, fusion, amplification, deep deletion, and multiple alterations.

LinkedOmics is a multi-omics database containing information on 32 cancer types from TCGA (Vasaikar et al., 2018). The significantly differentially expressed genes that were correlated with RUNX members were analyzed using the TCGA_BRCA cancer cohort (HiSeq RNA platform) in LinkedOmics. Kyoto Encyclopedia Genes and Genomes (KEGG) pathway analysis was performed using gene set enrichment analysis (GSEA). The genes were also classified using gene ontology (GO) according to: biological processes, cellular components, and molecular functions.

The TIMER web server (Li B. et al., 2016; Li et al., 2017) was used to analyze the infiltration of six types of immune cells in breast cancer, including B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells. The gene module was used to evaluate the relationship between target gene expression and immune cell infiltration, whereas the mutation module was used to analyze the gene mutation with an abundance of immune infiltrates. The SCNA module was used to analyze the correlation between RUNX gene somatic copy number alterations and tumor infiltration levels.

The Gene Set Cancer Analysis (GSCA) database (Liu C.J. et al., 2018) was used to analyze the relationship between RUNX factor methylation levels and the infiltration of six immune cell types: B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells in breast cancer.

TISIDB is an integrated web portal for the analysis of tumor-immune system interactions (Ru et al., 2019). The correlations between RUNX gene methylation and the expression of TGF-β1 and TGFBR1, two of the critical immunomodulators in breast cancer, were analyzed using TISIDB.

To study the role of RUNX family members in cancer, we performed a pan-cancer analysis using Oncomine. Transcription levels of RUNX1, RUNX2, and RUNX3 were all increased in most types of cancers, including breast, esophageal, head and neck, and pancreatic cancer (Figure 1A). Expression levels of RUNX factors were then compared across TCGA tumors. RUNX1 and RUNX2 expression in breast cancer was relatively high, ranking second after acute myeloid leukemia (LAML) for all the analyzed tumor types, indicating their potential role in breast cancer (Figure 1B).

The expression of RUNX members in breast cancer cells was analyzed using GEPIA. We found that the expression level of RUNX1 was significantly associated with the stage of breast cancer (p = 0.0001), while there was no significant difference in RUNX2 and RUNX3 expression (Figure 2). To further dissect the relationship between RUNX expression levels and breast cancer patient characteristics, we performed UALCAN analysis. The transcriptional level of RUNX1 in breast cancer was found to be associated with the TP53 mutation status. Breast cancer patients with mutant TP53 had significantly lower RUNX1 levels than normal controls, whereas patients with non-mutant TP53 had significantly higher RUNX1 levels (Figure 3A). There was a decreasing trend in RUNX1 transcriptional levels from stage one to stage four in breast cancer patients (Figure 3B). There was also a significant difference in RUNX1 levels between different breast cancer subclasses, with luminal type having the highest level and triple-negative showing the lowest level (Figure 3C). Furthermore, RUNX1 transcriptional levels were associated with different histological types (Figure 3D). Similarly, we also used the UALCAN database to further analyze the expression levels of RUNX2 and RUNX3 in breast cancer. We found that RUNX2 expression significantly varied in different stages, subclasses, and histological subtypes of breast cancer (Supplementary Figure 1). RUNX3 was down-regulated in luminal subclass but up-regulated in triple-negative subclass (Supplementary Figure 2). To validate these findings, we compared the levels of RUNX transcripts in different breast cancer subtypes to those of normal controls (Table 1). RUNX transcriptional levels have been shown to be significantly increased in various subtypes of breast cancer.

To study the relationship between RUNX transcription factor expression levels and the outcomes of patients with breast cancer, we performed KM plotter analysis. The results showed that patients with higher RUNX levels had longer DMFS, OS, and RFS (Figure 4). Specifically, higher RUNX1 expression was associated with longer DMFS, OS, and RFS. Higher RUNX2 expression was associated with longer RFS. Higher RUNX3 expression was associated with longer OS and RFS. These results indicate that RUNX1 and RUNX3 may play a more important role in the OS and RFS of breast cancer patients, and RUNX2 only affects the RFS of breast cancer patients.

The alterations in RUNX members in breast cancer subtypes were analyzed using cBio Cancer Genomics Portal (cBioPortal). Alteration types include mutation, fusion, amplification, deep deletion, and multiple alterations. Interestingly, the three RUNX members seemed to be altered in different subtypes of breast cancer with varying tendencies. The subtypes with the highest frequency of all alternate types were breast invasive lobular carcinoma, breast invasive carcinoma NOS, and breast invasive ductal carcinoma for RUNX1, RUNX2, and RUNX3, respectively (Figure 5A). For RUNX1, mutation ranked first in all alternate types in breast invasive lobular carcinoma (8.02%), followed by breast invasive carcinoma NOS (4.11%), and breast invasive ductal carcinoma (3.10%). RUNX1 was mainly distributed in the Runt domain and the linker region (Supplementary Figure 3). For RUNX2, the alteration with the highest frequency was amplification, with a percentage of 2.74% in breast invasive carcinoma NOS and 1.62% in breast invasive ductal carcinoma. For RUNX3, alterations were mainly found to be amplified in breast invasive ductal carcinoma (0.13%). RUNX1 showed the highest number of alterations among the three members (Figure 5B).

We analyzed the differentially expressed genes that were correlated with RUNX1 in breast cancer using LinkedOmics (Figure 6A). The top 50 positively and 50 negatively correlated genes were visualized in heatmaps (Figures 6B,C). GSEA was performed to identify the enriched KEGG pathways for the top significantly differentially expressed genes (Figure 6D). Among them, ECM-receptor interaction and focal adhesion were significantly upregulated, whereas cell cycle, RNA transport, pyrimidine metabolism, and DNA replication pathways were significantly downregulated. These genes were also classified using GO (Figure 6E). The top three enriched biological process terms were biological regulation, metabolism, and responses to stimulus. The top three enriched cellular component terms were the membrane, nucleus, and membrane-enclosed lumen. The top three enriched molecular function terms were protein binding, ion binding, and nucleic acid binding. Similarly, differentially expressed genes correlated with RUNX2 (Supplementary Figure 4) and RUNX3 (Supplementary Figure 5) were also analyzed and enriched.

Promoter methylation is an important regulator of gene expression. A previous study has shown that promoter methylation contributes to RUNX3 inactivation in breast cancer (Lau et al., 2006). We performed UALCAN analysis to evaluate RUNX1 promoter methylation levels and its relationship with breast cancer patient characteristics. Interestingly, we found that the RUNX1 promoter methylation level was significantly lower in breast cancer tissues than in normal tissues, regardless of stage, subclass, or histological type (Figure 7). This was consistent with the increased RUNX1 expression in breast cancer. Meanwhile, we also found that the RUNX2 promoter methylation level was significantly lower in breast cancer tissues than normal tissues, regardless of stage, subclass, or histological type (Figure 8). In contrast, RUNX3 promoter methylation level was significantly higher in breast cancer tissues than normal tissues, regardless of stage, subclass, or histological type (Figure 9). Moreover, we used the MEXPRESS database to further verify the methylation level of the RUNX family in breast cancer. The results of the MEXPRESS database analysis were consistent with the UALCAN database analysis. These results indicate that different changes in the methylation levels of RUNX factors in breast cancer may account for the distinct effects of RUNX proteins on patient outcome.

Runt-related transcription factors can regulate the development and maturation of immune cells, but the relationship between the expression levels of RUNX factors and immune cell infiltration in breast cancer is unclear. Therefore, we used the TIMER database to analyze the expression levels of RUNX factors and their relationship with the infiltration of six types of immune cells (B cells, CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells) in breast cancer. The results indicated that the expression of RUNX1 was positively correlated with the infiltration of CD8⁺T cells (COR = 0.183, P = 8.66E−09), CD4⁺T cells (COR = 0.13, P = 5.24E−05), neutrophils (COR = 0.103, P = 1.43E−03), and macrophages (COR = 0.271, P = 5.77E−18), whereas it was negatively correlated with B cell infiltration (COR = −0.084, P = 8.26E−03) (Figure 10A). The expression level of RUNX2 was positively correlated with the infiltration of all types of immune cells that were analyzed: B cells (COR = 0.066, P = 3.87E−02), CD8⁺T cells (COR = 0.294, P = 5.57E−21), CD4⁺T cells (COR = 0.196, P = 8.40E−10), macrophages (COR = 0.473, P = 7.75E−56), neutrophils (COR = 0.303, P = 1.25E−21), and dendritic cells (COR = 0.298, P = 6.53E−21) (Figure 10B). The expression of RUNX3 was also positively correlated with the infiltration of B cells (COR = 0.434, P = 3.38E−46), CD8⁺T cells (COR = 0.495, P = 1.81E−61), CD4⁺T cells (COR = 0.582, P = 4.04E−88), macrophages (COR = 0.14, P = 1.1E−05), neutrophils (COR = 0.635, P = 1.13E−108), and dendritic cells (COR = 0.667, P = 2.99E−123) (Figure 10C).

Furthermore, we analyzed the relationship between the expression of RUNX1 and the infiltration of six types of immune cells in BRCA-basal, BRCA-HER2, and BRCA-luminal breast cancer subtypes. We found that RUNX1 expression was only positively correlated with CD4⁺T cell invasion (COR = 0.227, P = 1.21E−02) in the BRCA-basal type. For the BRCA-HER2 subtype, RUNX1 expression was only positively correlated with macrophage infiltration (COR = 0.29, P = 2.72e−02). However, the expression level of RUNX1 in the BRCA-luminal subtype was positively correlated with the infiltration of CD8⁺T cells (COR = 0.161, P = 1.84E−04), CD4⁺T cells (COR = 0.208, P = 1.15E−06), macrophages (COR = 0.212, P = 6.26E−07), neutrophils (COR = 0.192, P = 7.81E−06), and dendritic cells (COR = 0.131, P = 2.39E−03) (Supplementary Figure 6). We also analyzed the relationship between the expression levels of RUNX2/RUNX3 and the infiltration of immune cells in BRCA-basal, BRCA-HER2, and BRCA-luminal breast cancer subtypes (Supplementary Figures 7, 8).

Various genes are mutated during breast cancer pathogenesis. Through analysis of the TIMER database, we found that RUNX1 had highly frequent mutations in breast cancer, while RUNX2 and RUNX3 had a lower frequency of mutation rate. Therefore, we analyzed the relationship between RUNX1 mutations and the infiltration of immune cells in breast cancer. The results showed that RUNX1 mutations significantly increased the levels of CD8⁺T cells, CD4⁺T cells, and macrophage infiltration in breast cancer (Supplementary Figure 9). This suggests that RUNX1 mutation may be used as a potential marker for detecting immune cell infiltration in breast cancer.

Runt-related transcription factors not only affect breast cancer cells, but also modulate the function of immune cells. Through TIMER analysis, we found a positive correlation between RUNX1 gene copy number and the infiltration of B cells, CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells in breast cancer (Figure 11A). The copy number change of the RUNX2 gene was only positively correlated with the infiltration of four types of immune cells in breast cancer: CD8⁺T cells, CD4⁺T cells, macrophages, and neutrophils (Figure 11B). In addition, the change in RUNX3 gene copy number was significantly correlated with the infiltration of CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells (Figure 11C). These results suggest that changes in the copy number of RUNX genes may reflect the infiltration of multiple immune cells in breast cancer.

Furthermore, we analyzed the relationship between RUNX gene copy number and immune cell infiltration in breast cancer subtypes. We found that the infiltration of CD4⁺T cells, macrophages, and neutrophils in BRCA-basal and BRCA-luminal subtypes was affected by changes in the RUNX1 gene copy number. However, the infiltration of immune cells into the BRCA-HER2 subtypes was not affected. In addition, the change in RUNX1 gene copy number in BRCA-luminal cells also affected the infiltration of B cells (Supplementary Figure 10). The change in RUNX2 gene copy number influenced the infiltration of multiple immune cells in the subtypes. In the BRCA-basal type, the RUNX2 gene copy number was associated with the infiltration of CD8⁺T cells, CD4⁺T cells, and neutrophils. The change in RUNX2 gene copy number in the BRCA-HER2 subtype was correlated with the infiltration of B cells, CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells. However, in the BRCA-luminal subtype, the change in RUNX2 gene copy number was only associated with the infiltration of B cells and CD4⁺T cells (Supplementary Figure 11). For RUNX3, changes in gene copy number only affected the infiltration of CD4⁺T cells in BRCA-HER2 and BRCA-luminal subtypes, whereas multiple types of immune cells were affected in the BRCA-basal subtype (Supplementary Figure 12). These results suggest that the change in RUNX1 gene copy number might mainly affect immune cell infiltration in BRCA basal and BRCA-luminal subtypes; the change in RUNX2 gene copy number might mainly affect immune cell infiltration in BRCA-basal and BRCA-HER2 subtypes; whilst the RUNX3 gene copy number change might mainly influence the change in immune cell infiltration in BRCA-basal type.

Through the above analysis, we have shown the methylation level change of RUNX factors in breast cancer, as well as the relationship between RUNX factor expression level and immune cell infiltration in breast cancer. However, it remains unclear whether changes in the methylation of RUNX genes also influence the infiltration of immune cells in breast cancer. Therefore, we further analyzed the association between the methylation levels of RUNX genes and immune cell infiltration (B cells, CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells) in breast cancer using the GSCA database (Figure 12). The results showed that RUNX2 methylation was positively correlated with the infiltration levels of B cells, CD8⁺T cells, and CD4⁺T cells, and negatively correlated with the infiltration levels of macrophages and neutrophils (Figure 12B). In contrast, the methylation level of RUNX3 was negatively correlated with the infiltration of B cells, CD8⁺T cells, and CD4⁺T cells, whereas it was positively correlated with the infiltration levels of macrophages and neutrophils (Figure 12C). This suggests that the methylation levels of RUNX factors may have opposite effects on the same type of immune cell infiltration.

TGF-β is a critical negative immunoregulatory factor in immune balance. TGF-β not only inhibits various targets in the immune system, but also plays a role in tumor immune escape and adverse reactions to tumor immunotherapy (Batlle and Massague, 2019; Derynck et al., 2021). Through analysis of the TIMER database, we found that the expression of TGFB1, TGFBR1, and TGFBR2 are correlated with the infiltration levels of various immune cells in breast cancer (Supplementary Figure 13). We used the TIMER database to further analyze the relationship between the expression level of TGFB1 and immune cell infiltration in different breast cancer subtypes. We found that the expression level of TGFB1 was positively correlated with the infiltration level of various immune cells in different breast cancer subtypes (Supplementary Figure 14). The relationship between the expression levels of TGFBR1 and TGFBR2 and immune cell infiltration in different breast cancer subtypes was also analyzed. We found that the expression level of TGFBR1 was mainly correlated with the level of macrophage infiltration in different breast cancer subtypes (Supplementary Figure 15). The expression level of TGFBR2 was positively correlated with the infiltration level of multiple immune cells in different subtypes of breast cancer (Supplementary Figure 16). These results indicate that the TGF-β signaling pathway is related to the infiltration level of multiple immune cells in different subtypes of breast cancer. However, it remains unclear whether the methylation levels of RUNX genes affect the TGF-β signaling pathway. To further evaluate the potential mechanisms underlying the effects of RUNX gene methylation on immune cell infiltration in breast cancer, we analyzed the relationship between the methylation level of RUNX genes and the TGF-β signaling pathway using the TISIDB database. The results showed that RUNX1 gene methylation level was positively correlated with the expression level of TGF-β1 (Figure 13A), while there was no correlation between RUNX2 methylation level and the expression of TGF-β1 (Figure 13B). However, RUNX3 methylation levels were negatively correlated with the expression levels of TGF-β1 (Figure 13C). The relationship between the methylation level of RUNX genes and the expression of TGFBR1 was also analyzed, but no correlation was identified (Supplementary Figure 17). These results suggest that changes in the methylation level of RUNX genes may affect the expression of TGF-β1, which in turn might influence the activation of the TGF-β signaling pathway and immune cell infiltration in breast cancer.