Abstract Forkhead-box-P family include FOXP1/2/3/4 and its clinical significance still remains unclear in breast cancer (BRCA). We analysed the expressions of FOXPs in BRCA patients to determine diagnostic and prognostic values. Our results indicated that the transcriptional levels of FOXP3/4 were up-regulated in BRCA patients, but FOXP2 were down-regulated. No statistically significant correlation were found between the expression levels of FOXPs in Pathologic stage. FOXP2/3 had a significantly high AUC value in the detection of breast cancer, with 96.8% or 95.7% in accuracy respectively. Our study also suggested that BRCA patients with high transcription levels of FOXP1/2/4 were significantly associated with longer Overall Survival (OS). In contrast, BRCA patients with high transcription level of FOXP3 was not statistically related with OS. Our work revealed that FOXPs were closely related to the alteration of extensive immune checkpoints in breast invasive carcinoma. Additionally, FOXP3 has a significant positive correlation with PDCD1, CD274, CTLA4 and TMB in breast cancer, and FOXP3 expression showed a statistically significant correlation with infiltration of immune cells. Finally, we found that FOXP3 expression predicted the breast cancer cells response to anticancer drugs. Altogether, our work strongly suggested that FOXPs could serve as a biomarker for tumor detection, therapeutic design and prognosis. Subject terms: Breast cancer, Immunology Introduction Surpassing lung cancer, female breast cancer has now been in the leading position of the world cancer incidence in 2020, with about 2.3 million new cases, accounting for 11.7% of all cancer cases. It is the fifth cause of cancer mortality in the world, with 685,000 deaths^[38]1. Despite the recent advanced treatments made in breast cancer treatments include earlier detection, surgery, chemotherapy, immunotherapy and targetting drug therapy^[39]2,[40]3. In 5–10% of the breast cancer patients, tumors have already expanded to the advanced stage with extensive lymph invasion or distant metastasis when diagnosed, there is still poor survival rate of these metastatic breast cancer patients, of these patients only one-fifth survive 5 years^[41]4. Due to individual difference and tumor heterogeneity, the current diagnostic and prognostic biomarkers for breast cancer have some limitations^[42]5, It is therefore imperative demand to investigate more effective diagnostic and prognostic biomarkers for optimizing the management of breast cancer. Cancer progression is involved in epigenetic and genetic alterations including transcription factors, growth factors, cytokines, and proteases under tumor microenvironment^[43]6. As one of transcription factors, Forkhead box P (FOXP) family include FOXP1, FOXP2, FOXP3 and FOXP4 with similar 110 amino acid DNA-binding domain termed forkhead domain^[44]7. FOXP proteins can regulate gene transcription in connection with carcinogenesis^[45]8, immune function^[46]9, invasion and metastasis of carcinoma^[47]10, differentiation^[48]11 and angiogenesis^[49]12. Accumulating evidence shows that FOXP family proteins have dual biologic functions as an oncogene or a tumor suppressor. Diffuse large B-cell lymphoma with overexpressed FOXP1 has poor prognosis^[50]13, while FOXP1 acts as a tumor suppressor in breast and lung carcinoma^[51]14,[52]15. CD4+/CD25+/FOXP3+ Treg cells work in gastric cancers through immunosuppression as oncogenes, while FOXP3 overexpression in patients with breast cancers indicates good prognosis as a tumor suppressor^[53]16,[54]17. It remains unclear that clinical significance of FOXP family proteins act as an entirety in human breast cancer^[55]18. Bioinformatics analysis has been applied to survey the role of transcription factors in breast cancer. In the present study, according to the analyses of gene expressions or variations in published online, we analysis the expressions and mutations of different FOXP factors in patients with breast cancer in detail to determine diagnostic and prognostic values of FOXP in breast cancer. Methods Oncomine database analysis The expression level of the FOXPs in various types of cancers was identified in the Oncomine database ([56]https://www.oncomine.org/resource/login.html)^[57]19. The threshold was determined according to the following values: P value of 0.01, fold change of 1.0, and gene ranking of all. Tumor immune estimation resource (TIMER) database TIMER is a comprehensive online resource for systematic analysis of immune infiltrates across various cancer types^[58]20. In this study, we observed the expression difference of FOXPs between tumor and adjacent normal tissues for the BRCA of the TCGA project. Meanwhile, we performed TIMER to determine the relationship between FOXPs expression in BRCA and 6 immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells). RNA-sequencing data of FOXPs in human BRCA The RNA-Seq expression data of FOXPs in BRCA was downloaded from TCGA ([59]https://portal.gdc.cancer.gov/). Therefore, 113 adjacent normal tissues and 1109 BRCA data were retained. The samples selected contained FOXPs gene expression data and associated clinical information, including age, gender, Pathological stage, Race, Histological type. Immunohistochemistry Clinical samples were obtained from breast cancer patients who were surgically treated at Hunan Provincial People’s Hospital/The First Affiliated Hospital of Hunan Normal University. Tumor tissue and its adjacent normal tissues were prepared into 4 mm paraffin sections and incubated with primary rabbit monoclonal antibodies of FOXP1, FOXP2, FOXP3, FOXP4 (1:200 dilution; Santa Cruz Biotechnology, Santa Cruz, USA) at 4° overnight. The sections were coupled with goat anti-rabbit antibody labeled with horseradish peroxidase (1:400, Abcam, USA) at room temperature for 60 min, then each section was stained with 3,3-diaminobenzidine (DAB) reagent, and finall weakly counterstained with hematoxylin. The Kaplan–Meier plotter analysis The Kaplan–Meier plotter ([60]www.kmplot.com)^[61]21,[62]22 was used to assess the prognostic value of FOXPs mRNA expression in BRCA patients and analyzed the overall survival (OS), progression-free survival (PFS), post-progression survival (PPS) and distant metastasis-free survival (DMFS) of patients with BRCA. The patients divided into high expression groups and low expression groups according to the median values of FOXPs mRNA expression. The cBioPortal analysis We selected a Breast Invasive Carcinoma dataset (TCGA, Firehose Legacy) that contained 1109 pathological reports to analyzed the expression of FOXPs and immune checkpoints using cBioPortal ([63]www.cbioportal.org)^[64]23,[65]24. The genomic map contains putative copy-number alterations (CNA) from GISTIC, mRNA expression z-scores and Protein expression z-scores mutations. STRINGS analysis STRINGS ([66]www.string-db.org) is an online tool for analysis of all publicly data of protein–protein interaction (PPI)^[67]25. In this study, we used a PPI network analysis on FOXPs to explore their functions in human breast cancer. GeneMANIA analysis GeneMANIA ([68]www.genemania.org) is an online tool for analysis of gene functions^[69]26. In this study, we performed GeneMANIA to select the 50 most important genes to construct gene–gene interaction network for FOXPs. Drug–gene interactions DGidb^[70]27 ([71]https://dgidb.genome.wustl.edu/), a web server for discovering drug–gene interactions or potentially available drug categories, was used to explore the potential druggable genes and drugs of FOXPs in patients. Statistical analyses All statistical analyses were implemented with R ([72]www.r-project.org). gene expression data and clinical information were visualizes by R package “ggplot2” R package^[73]28. ROC curve was performed to detect the cutoff value of FOXPs, BRCA1, BRCA2 and ERBB2. Lollipop chart used Spearman’s correlation analysis to describe the correlation between 24 immune cell^[74]29,[75]30. TMB Score: We used Spearman’s correlation analysis to describe the correlation between quantitative variables without a normal distribution. A P value of less than 0.05 was considered statistically significant. Ethics statement All experiments involving clinical samples (Informed consent was obtained from all subjects and/or their legal guardian(s)) were approved by the Ethics Committee of Hunan Provincial People’s Hospital/The First Affiliated Hospital of Hunan Normal University (document NO.202193). All the experiments were conducted in accordance with ethical guidelines, including the tenets of the Declaration of Helsinki. Significance Our study strongly suggests that FOXPs could serve as a biomarker for tumor detection, therapeutic design and prognosis. Results Transcriptional levels of FOXPs in BRCA patients Four FOXPs are generally found in mammalian cells, but are expressed abnormally in different tumor tissues. We performed the Oncomine and TIME to compare the mRNA expression of FOXPs in different cancer and normal tissue samples (Fig. [76]1A,B). The results showed that the transcriptional levels of FOXP3/4 were up-regulated in BRCA patients, but FOXP2 were down-regulated. However, FOXP1 in Oncomine is up-regulated in cancer tissues and FOXP1 in TIME is up-regulated in normal tissues. In addition, the significant changes of FOXPs expression in transcription level between breast cancer and normal breast tissues showed in Table [77]1 (Oncomine database). Unpaired data analysis also showed that the mRNA expression levels of FOXP3/4 in BRCA tissues (n = 1109) were significantly higher than those in adjacent normal tissues (n = 113), and FOXP1/2 in BRCA tissues (n = 1109) were significantly lower than those in adjacent normal tissues (n = 113). (Fig. [78]1C, FOXP1 N: 4.851 ± 0.457 vs. T: 4.534 ± 0.804; FOXP2 N: 1.465 ± 0.481 vs. T: 0.539 ± 0.462; FOXP3 N: 1.428 ± 0.646 vs. T: 2.818 ± 1.02; FOXP4 N: 5.319 ± 0.503 vs. T: 5.819 ± 0.834, Mann–Whitney U-test, P < 0.001) (ns, no significance, *P < 0.05, **P < 0.01, ***P < 0.001) . Figure 1. [79]Figure 1 [80]Open in a new tab FOXPs expression levels, clinical characteristics and immunohistochemistry in human cancer. (A) The transcription levels of FOXPs in different types of cancers (Oncomine). (B) Human FOXPs expression levels in human BRCA from TCGA database were determined (TIMER). (C) The mRNA expression levels of FOXPs in 1109 BRCA samples and 113 normal samples from TCGA. (D) FOXP2,3,4 mRNA expression were significantly higher in patients (> = 60) than in patients (< 60), however the result was the opposite in FOXP1. (E) No statistically significant correlation were found between the expression levels of FOXPs in Pathologic stage. (F) The expression of FOXPs was different among different races. (G) FOXP1,2 mRNA expression were significantly higher in Infiltrating Lobular Carcinoma than Infiltrating Ductal Carcinoma, however no statistically significant correlation were found between FOXP3,4. (H) The Expression of FOXPs in BRCA (IHC) (ns, no significance, *P < 0.05, **P < 0.01, ***P < 0.001). Table 1. The significant changes of FOXP expression in transcription level between different types of breast cancer and normal breast tissues (Oncomine database). Type of breast cancer versus normal breast tissue Fold change P value t test Source FOXP1 Benign breast neoplasm 1.101 8.27E−06 5.374 Curtis breast statistics Ductal breast carcinoma in situ epithelia 1.18 0.001 3.805 Curtis breast statistics Invasive ductal and invasive lobular breast carcinoma 1.199 1.30E−12 7.497 Curtis breast statistics Invasive breast carcinoma 1.216 7.96E−04 3.573 Curtis breast statistics Invasive lobular breast carcinoma 1.208 1.49E−13 7.676 Curtis breast statistics Mucinous breast carcinoma 1.195 7.76E−06 4.69 Curtis breast statistics Invasive ductal breast carcinoma 1.195 5.00E−23 11.181 Curtis breast statistics Tubular breast carcinoma 1.084 5.31E−05 3.958 Curtis breast statistics Invasive breast carcinoma stroma 4.861 1.01E−19 18.498 Finak breast statistics Invasive ductal breast carcinoma stroma 1.584 5.45E−04 4.541 Karnoub breast statistics FOXP2 Ductal breast carcinoma 1.273 1.73E−05 4.81 Richardson breast 2 statistics Invasive ductal breast carcinoma 1.017 0.003 2.732 TCGA breast 2 Invasive breast carcinoma stroma 1.803 8.79E−05 5.936 Finak breast statistics FOXP3 Mucinous breast carcinoma 1.309 0.002 3.839 TCGA Invasive breast carcinoma 1.136 0.007 2.483 TCGA Breast carcinoma 1.148 0.002 3.384 Curtis breast statistics Invasive breast carcinoma 1.075 0.002 3.046 Curtis breast statistics Invasive ductal and invasive lobular breast carcinoma 1.073 2.99E−05 4.115 Curtis breast statistics Invasive lobular breast carcinoma 1.042 0.001 2.998 Curtis breast statistics Invasive ductal breast carcinoma 1.053 9.01E−07 4.943 Curtis breast statistics Invasive breast carcinoma stroma 1.592 6.06E−07 10.59 Finak breast statistics FOXP4 Mucinous breast carcinoma 3.061 2.28E−05 10.005 TCGA Invasive ductal and lobular carcinoma 1.544 2.23E−06 9.435 TCGA Invasive ductal breast carcinoma 1.681 2.36E−24 12.939 TCGA Invasive breast carcinoma 1.66 8.40E−13 7.866 TCGA Mixed lobular and ductal breast carcinoma 1.224 0.002 3.233 TCGA Invasive lobular breast carcinoma 1.343 6.07E−05 4.132 TCGA Benign breast neoplasm 1.07 0.005 4.380 Curtis breast statistics Medullary breast carcinoma 1.227 0.003 2.883 Curtis breast statistics Invasive ductal breast carcinoma 1.113 1.30E−07 5.337 Curtis breast statistics Ductal breast carcinoma 2.122 3.40E−06 5.272 Richardson breast 2 statistics Ductal breast carcinoma in situ epithelia 1.46 0.004 3.282 Ma breast 4 statistics Invasive breast carcinoma stroma 1.308 1.80E−10 9.534 Finak breast statistics Invasive ductal breast carcinoma 1.041 2.28E−08 5.534 TCGA breast 2 [81]Open in a new tab Relationships between FOXPs mRNA levels and clinical characteristics of BRCA patients To evaluate the association between the mRNA expression of FOXPs and clinical pathological characteristics of BRCA samples, we performed Mann–Whitney U-test analysis. As shown in Fig. [82]1D, FOXP2,3,4 mRNA expression were significantly higher in patients (< 60) than in patients (> = 60) (p.adj < 0.001, p.adj = 0.006, p.adj < 0.001), however the result was the opposite in FOXP1(P = 0.009). No statistically significant correlation were found between the expression levels of FOXPs in Pathologic stage (p.adj > 0.05) (Fig. [83]1E). We used the Bonferroni method to correct the multiple hypothesis test (Dunn's test) of significance level. The results of Race showed that FOXP1 expression was lower in black or African American than in Asian, and the difference was statistically significant (p.adj = 0.001); FOXP2 expression was lower in black or African American than in Asian, and the difference was not statistically significant (p.adj = 1); No statistically significant correlation were found between the expression levels of FOXPs and Race (p.adj > 0.05); FOXP4 expression was lower in black or African American than Asian, and the difference was not statistically significant (p.adj = 1) (Fig. [84]1F). The results of Histological type showed that FOXP1,2 mRNA expression were significantly higher in Infiltrating Lobular Carcinoma than Infiltrating Ductal Carcinoma (p.adj < 0.001), however no statistically significant correlation were found between FOXP3,4 and Histological type (p.adj > 0.05) (Fig. [85]1G). (ns, no significance, *P < 0.05, **P < 0.01, ***P < 0.001). We used immunohistochemistry (IHC) to detect the protein expression of FOXPs in BRCA and its paired adjacent tissues. The results showed that FOXPs protein express in the nucleus. FOXPs protein express in the nucleus. We found that the protein levels of FOXP3 and FOXP4 were higher in BRCA tissues than in the adjacent tissues, however the result was the opposite in FOXP1 and FOXP2 (Fig. [86]1H). The prognostic value of FOXPs in BRCA patients To evaluate the value of FOXPs at different transcription levels in the progression of BRCA, we evaluated the correlation between FOXPs at different transcription levels and clinical outcome using Kaplan–Meier plotter analysis. The OS curve is shown in Fig. [87]2. BRCA patients with high transcription levels of FOXP1/2/4 were significantly associated with longer OS. In contrast, BRCA patients with transcription levels of FOXP3 was not statistically related with OS. In addition, the studies showed that BRCA with high expression of FOXP1/2/3/4 was significantly associated with longer PFS (Fig. [88]2). While transcription levels of FOXP1/2 were not related to PPS in BRCA patients. In contrast, BRCA patients with low transcription levels of FOXP3/4 were significantly associated with longer PPS (Fig. [89]2). The value of FOXPs at different transcription levels in DMFS of BRCA patients was also evaluated. Only BRCA patients with high mRNA expression of FOXP1 were significantly related with longer DMFS (Fig. [90]2). Figure 2. [91]Figure 2 [92]Open in a new tab The prognostic value of mRNA Level of FOXPs in BRCA patients (Kaplan–Meier plotter). Overall survival (OS), progression-free survival (PFS), post-progression survival (PPS) and distant metastasis-free survival (DMFS). Co-expression of FOXPs in BRCA patients Evaluation of the mutual exclusion between the four FOXPs genes in the TCGA BRCA cohort showed that there were co-expressed relationships between FOXP1 and FOXP2/FOXP3/FOXP4/ERBB2/BRCA1, FOXP2 and FOXP3/BRCA2, FOXP3 and FOXP4/ERBB2/BRCA2, FOXP4 and ERBB2/BRCA1/BRCA2, ERBB2 and BRCA1/BRCA2, BRCA1 and BRCA2 (P < 0.05). However, FOXP1 with BRCA2, FOXP2 with FOXP4/ERBB2/BRCA1 and FOXP3 with BRCA1 didn’t have co-expressed relationships (P > 0.05) (Fig. [93]3A). Figure 3. [94]Figure 3 [95]Open in a new tab Relevance of different FOXPs and important genes, diagnostic value, TMB score and mutation analysis in BRCA. (A) Correction between different FOXPs and Important Genes in BRCA. (B) ROC curve analysis was performed to evaluate the diagnostic power of the FOXPs and Important Genes. (C) Correlation analysis of FOXPs expression and TMB in BRCA. The horizontal axis in the figure represents the expression distribution of the gene, and the ordinate is the expression distribution of the TMB score. The density curve on the right represents the distribution trend of the TMB score; The upper density curve represents the distribution trend of the gene; The top side: the value represents the correlation p value, correlation coefficient and correlation calculation method. (D) Summary of alterations in different expressed FOXPs and immune checkpoint in BRCA (cBioPortal). Differential RNA-Seq levels of FOXPs as prospective biomarker to distinguish BRCA samples from normal samples and the TMB score of FOXPs in BRCA To investigate the value for FOXPs to distinguish Breast Invasive Carcinoma samples from normal smples, we performed a ROC curve analysis using ERBB2, BRCA1 and BRCA2 as controls. As showed in (Fig. [96]3B) (Table [97]S1), the ROC curve analysis showed FOXP1 (AUC: 0.530, accuracy: 0.862), FOXP2 (AUC: 0.927, accuracy: 0.968), FOXP3 (AUC: 0.932, accuracy: 0.957), FOXP4 (AUC: 0.764, accuracy: 0.928), ERBB2 (AUC: 0.801, accuracy: 0.900), BRCA1 (AUC: 0.861, accuracy: 0.973), BRCA2 (AUC: 0.849, accuracy: 0.957). These fingerings indicated that FOXP2 and FOXP3 could be a promising biomarker to differentiate Breast Invasive Carcinoma tissues from normal tissues. In addition, we made a TMB scoring model (Fig. [98]3C). Genetic alteration of FOXPs and immune checkpoint in BRCA patients We used the cBioPortal online tool to analyze changes and correlations of FOXPs in Breast Invasive Carcinoma. Among 1108 Breast Invasive Carcinoma patients, FOXPs were changed in 258 samples (23.28%) (Fig. [99]3D). The genomic investigation revealed that FOXPs was actually involved in the alteration of immune checkpoints in Breast Invasive Carcinoma. The general landscape of FOXPs and immune checkpoint alteration in Breast Invasive Carcinoma was compactly visualized, including structural variant, mRNA, amplification, deep deletion, truncating, splice and missense mutations (Fig. [100]4A). The detailed relationship between FOXPs and each representative immune checkpoint was individually presented as indicated in Table [101]2. Of note, the FOXP3 alteration showed a statistically significant co-occurrence rather than mutual exclusivity with extensive immune checkpoints, such as CTLA4, CD48, PDCD1, CD70, PDCD1LG2, CD86, CD80, CD274, ICOSLG. These findings strongly indicate that FOXP3 is a potential coregulator of immune checkpoints in Breast Invasive Carcinoma. Figure 4. [102]Figure 4 [103]Open in a new tab FOXPs and immune checkpoints mutation analysis and the relationship between the expression of FOXPs and immune checkpoints genes (PDCD1, CD274, CTLA4) in BRCA. (A) Landscape of FOXPs and immune checkpoint alteration in Breast Invasive Carcinoma (cBioPortal). (B) FOXP3 has a significant positive correlation with PDCD1, CD274, CTLA4, but FOXP1, 2, 4 didn’t. Table 2. Mutual-exclusivity analysis between FOXPs and multiple-immune checkpoints in BRCA. A B Neither A Not B B Not A Both Log2 odds ratio p value q-Value Tendency Significant FOXP1 PVR 651 41 34 12 2.486 < 0.001 < 0.001 Co-occurrence Yes FOXP1 NECTIN2 645 41 40 12 2.239 < 0.001 < 0.001 Co-occurrence Yes FOXP1 CD70 650 45 35 8 1.723 0.008 0.027 Co-occurrence Yes FOXP1 TNFSF4 604 41 81 12 1.126 0.025 0.064 Co-occurrence No FOXP1 VTCN1 637 45 48 8 1.238 0.039 0.092 Co-occurrence No FOXP1 CD86 646 46 39 7 1.334 0.039 0.092 Co-occurrence No FOXP1 HHLA2 661 48 24 5 1.521 0.05 0.108 Co-occurrence No FOXP1 TNFRSF14 643 46 42 7 1.22 0.053 0.11 Co-occurrence No FOXP1 CD200 648 47 37 6 1.161 0.079 0.15 Co-occurrence No FOXP1 ICOSLG 633 46 52 7 0.889 0.12 0.205 Co-occurrence No FOXP1 PDCD1 642 47 43 6 0.931 0.13 0.218 Co-occurrence No FOXP1 CD48 576 41 109 12 0.629 0.14 0.229 Co-occurrence No FOXP1 CD80 658 49 27 4 0.992 0.176 0.272 Co-occurrence No FOXP1 PDCD1LG2 636 47 49 6 0.729 0.194 0.286 Co-occurrence No FOXP1 VSIR 646 48 39 5 0.787 0.201 0.293 Co-occurrence No FOXP1 TNFSF18 604 45 81 8 0.407 0.301 0.393 Co-occurrence No FOXP1 CD274 631 48 54 5 0.284 0.42 0.495 Co-occurrence No FOXP1 LGALS9 629 48 56 5 0.227 0.45 0.513 Co-occurrence No FOXP1 CTLA4 648 51 37 2 − 0.542 0.457 0.517 Mutual exclusivity No FOXP2 HHLA2 658 51 21 8 2.297 0.001 0.006 Co-occurrence Yes FOXP2 VTCN1 633 49 46 10 1.49 0.009 0.029 Co-occurrence Yes FOXP2 CD200 644 51 35 8 1.529 0.016 0.045 Co-occurrence Yes FOXP2 NECTIN2 636 50 43 9 1.413 0.017 0.047 Co-occurrence Yes FOXP2 CTLA4 647 52 32 7 1.445 0.029 0.072 Co-occurrence No FOXP2 CD70 643 52 36 7 1.266 0.047 0.104 Co-occurrence No FOXP2 VSIR 642 52 37 7 1.224 0.053 0.11 Co-occurrence No FOXP2 PDCD1LG2 632 51 47 8 1.077 0.062 0.126 Co-occurrence No FOXP2 CD86 640 52 39 7 1.143 0.065 0.131 Co-occurrence No FOXP2 TNFRSF14 637 52 42 7 1.03 0.086 0.159 Co-occurrence No FOXP2 CD80 653 54 26 5 1.218 0.093 0.167 Co-occurrence No FOXP2 CD48 571 46 108 13 0.579 0.15 0.242 Co-occurrence No FOXP2 CD274 627 52 52 7 0.699 0.182 0.275 Co-occurrence No FOXP2 ICOSLG 627 52 52 7 0.699 0.182 0.275 Co-occurrence No FOXP2 TNFSF4 596 49 83 10 0.551 0.195 0.286 Co-occurrence No FOXP2 PVR 638 54 41 5 0.527 0.303 0.393 Co-occurrence No FOXP2 LGALS9 624 53 55 6 0.361 0.36 0.444 Co-occurrence No FOXP2 PDCD1 633 56 46 3 − 0.44 0.436 0.504 Mutual exclusivity No FOXP2 TNFSF18 597 52 82 7 − 0.029 0.581 0.6 Mutual exclusivity No FOXP3 CTLA4 676 23 29 10 > 3 < 0.001 < 0.001 Co-occurrence Yes FOXP3 CD48 600 17 105 16 2.427 < 0.001 < 0.001 Co-occurrence Yes FOXP3 PDCD1 666 23 39 10 2.892 < 0.001 < 0.001 Co-occurrence Yes FOXP3 CD70 671 24 34 9 2.888 < 0.001 < 0.001 Co-occurrence Yes FOXP3 PDCD1LG2 660 23 45 10 2.673 < 0.001 < 0.001 Co-occurrence Yes FOXP3 CD86 668 24 37 9 2.759 < 0.001 < 0.001 Co-occurrence Yes FOXP3 CD80 681 26 24 7 2.933 < 0.001 0.001 Co-occurrence Yes FOXP3 CD274 655 24 50 9 2.296 < 0.001 0.004 Co-occurrence Yes FOXP3 ICOSLG 655 24 50 9 2.296 < 0.001 0.004 Co-occurrence Yes FOXP3 CD200 667 28 38 5 1.648 0.037 0.089 Co-occurrence No FOXP3 LGALS9 650 27 55 6 1.393 0.047 0.104 Co-occurrence No FOXP3 VTCN1 654 28 51 5 1.195 0.097 0.172 Co-occurrence No FOXP3 TNFRSF14 660 29 45 4 1.016 0.169 0.263 Co-occurrence No FOXP3 TNFSF4 618 27 87 6 0.659 0.227 0.32 Co-occurrence No FOXP3 PVR 662 30 43 3 0.622 0.339 0.425 Co-occurrence No FOXP3 TNFSF18 621 28 84 5 0.401 0.366 0.449 Co-occurrence No FOXP3 NECTIN2 656 30 49 3 0.421 0.415 0.494 Co-occurrence No FOXP3 VSIR 663 31 42 2 0.026 0.599 0.609 Co-occurrence No FOXP3 HHLA2 677 32 28 1 − 0.404 0.624 0.627 Mutual exclusivity No FOXP4 NECTIN2 643 43 38 14 2.462 < 0.001 < 0.001 Co-occurrence Yes FOXP4 ICOSLG 634 45 47 12 1.847 < 0.001 0.005 Co-occurrence Yes FOXP4 PDCD1LG2 636 47 45 10 1.588 0.006 0.022 Co-occurrence Yes FOXP4 PDCD1 640 49 41 8 1.35 0.028 0.07 Co-occurrence No FOXP4 TNFRSF14 640 49 41 8 1.35 0.028 0.07 Co-occurrence No FOXP4 CD70 645 50 36 7 1.327 0.04 0.093 Co-occurrence No FOXP4 CD48 574 43 107 14 0.805 0.066 0.131 Co-occurrence No FOXP4 VSIR 638 56 43 1 − 1.916 0.128 0.216 Mutual exclusivity No FOXP4 TNFSF18 602 47 79 10 0.697 0.134 0.222 Co-occurrence No FOXP4 CD200 643 52 38 5 0.702 0.232 0.324 Co-occurrence No FOXP4 CD80 651 56 30 1 − 1.368 0.291 0.388 Mutual exclusivity No FOXP4 CD274 628 51 53 6 0.479 0.3 0.393 Co-occurrence No FOXP4 CTLA4 644 55 37 2 − 0.66 0.405 0.491 Mutual exclusivity No FOXP4 TNFSF4 596 49 85 8 0.195 0.431 0.502 Co-occurrence No FOXP4 VTCN1 630 52 51 5 0.248 0.439 0.505 Co-occurrence No FOXP4 CD86 638 54 43 3 − 0.279 0.517 0.558 Mutual exclusivity No FOXP4 PVR 638 54 43 3 − 0.279 0.517 0.558 Mutual exclusivity No FOXP4 LGALS9 625 52 56 5 0.102 0.518 0.558 Co-occurrence No FOXP4 HHLA2 654 55 27 2 − 0.183 0.609 0.616 Mutual exclusivity No [104]Open in a new tab The relationship between the expression of FOXPs with PDCD1, CD274, CTLA4 in BRCA The two-gene correlation map was realized by the R software package “ggstatsplot”, and data from TCGA.We found that FOXP3 had a significant positive correlation with PDCD1, CD274, CTLA4 (Fig. [105]4B). However, FOXP1/2/4 was weakly associated with PDCD1, CD274, CTLA4 (Fig. [106]4B). The relationship between FOXPs expression levels and immune infiltration levels in BRCA TIMER online analysis tool is used to evaluate the relationship between the transcription level of FOXPs and the level of immune infiltration in BRCA. It was found that FOXPs are involved in inflammatory response and immune cell infiltration. The analysis results are shown in Fig. [107]5. FOXP1 expressions was positively correlated with the infiltration of CD4+ T cells, CD8+ T cells, neutrophils and macrophages (Fig. [108]5A). FOXP2 expressions was positively correlated with CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells (Fig. [109]5B). FOXP3 expression was positively correlated with infiltration of B cells, CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells (Fig. [110]5C). FOXP4 expressions was positively correlated with the infiltration of B cells, CD4+ T cells, neutrophils and macrophages, while it was negatively correlated with infiltration of Fig. [111]5D. These studies indicated that the level of FOXPs expression was associated to the level of immune infiltration in BRCA. At the same time, we made four Lollipop charts 24 immune cells associated with FOXPs expression in BRCA to have a more intuitive understanding of the immune infiltration (Fig. [112]5E). We can see that there were most immune cell infiltration in FOXP2 and FOXP3. Figure 5. [113]Figure 5 [114]Open in a new tab The relationship between FOXPs expression levels and immune infiltration levels in BRCA. The correlation between the abundance of immune cell and the expression of (A) FOXP1, (B) FOXP2, (C) FOXP3, (D) FOXP4 in BRCA (TIMER). (E) Lollipop chart about correlation between 24 immune cell and the expression of FOXPs in BRCA. PPI and neighbor gene network of FOXPs in BRCA patients We performed a PPI network analysis of FOXPs at different transcription levels using STRING to study the potential interactions between them. As shown in (Fig. [115]6A), the PPI network diagram contains four FOXPs proteins and 10 proteins that are closely related to FOXPs. Figure 6. [116]Figure 6 [117]Open in a new tab PPI network, neighbor gene network, interaction analyses of FOXPs, the functions of FOXPs, genes significantly associated with FOXPs alterations in BRCA and drug–gene interaction. (A) Protein–protein interaction network of different expressed FOXPs (STRING). (B) Gene–gene interaction network of different expressed FOXPs (GeneMANIA). (C) GO enrichment analysis predicted the functional roles of target host genes based on three aspects, including biological processes, cellular components, and molecular functions. The functions of FOXPs and genes significantly associated with FOXPs alterations were predicted by the analysis of KEGG. (D,E) Network of GO and KEGG enriched terms. (F) FOXP3 has significant correlations with drug-gene interaction. A GGI network of four FOXPs was constructed, and their functions were analyzed using the GeneMANIA database (Fig. [118]6B). The four central nodes of FOXPs are surrounded by 50 nodes representing genes that are strongly associated with FOXPs in shared protein domains, physical interactions, colocalization, co-expression, prediction, genetic interactions, and pathway. The top 50 genes most associated with FOXPs are NFATC2, CTLA4, RORC, MKI67, SIVA1, HDAC9, GATAD2B, IKZF3, SFTPC, TBR1 and so on. Among them, NFATC2 is related to FOXP1/2/3/4 in terms of physical interaction and genetic interaction. NFATC2 and CTLA4 have a pathway relationship with FOXP3. CTLA4 is related to FOXP3 in terms of co-localization. However, RORC, MKI67, SIVA1, HDAC9, GATAD2B, IKZF3, SFTPC, TBR1 and FOXP1/2/3/4 all have physical interaction. Further functional analysis showed that these genes indicated the greatest correlation with lymphocyte differentiation (FDR = 1.54E-4). In addition, these genes were correlated with regulatory T cell differentiation, T cell differentiation, regulation of lymphocyte differentiation, regulation of leukocyte differentiation, regulation of hemopoiesis and regulation of T cell differentiation. Functional enrichment analysis of FOXPs in BRCA patients In this study, we used “ClusterProfiler” R package to perform functional annotation and pathway enrichment analysis of FOXPs from 50 nodes representing genes. The first 15 items of GO enrichment (Table [119]S2) are mainly distributed in the biological process (5 items) (Fig. [120]6C), the cell component (5 items) (Fig. [121]6C) and the molecular function (5 items) (Fig. [122]6C). Three of the first five projects are in the T cell function, which are regulation of regulatory T cell differentiation, regulatory T cell differentiation, lymphocyte differentiation, and the other two are regulation of leukocyte differentiation and DNA-binding transcription activator activity, RNA polymerase II-specific. The first 5 KEGG pathways of FOXPs are illustrated in Fig. [123]6C (Table [124]S2). Among them, Th17 cell differentiation, FoxO signaling pathway, T cell receptor signaling pathway, Inflammatory bowel disease and Th1 and Th2 cell differentiation are significantly associated with the occurrence and development of various tumors and are also involved in the tumorigenesis of BRCA. At the same time, we made more intuitive GO network map (Fig. [125]6D) and KEGG network map to show the connection between pathways (Fig. [126]6E). Drug–gene interactions The result showed expressions of FOXP3 had significant correlations with drug-gene interaction. A total of 3 drugs were explored using DGIDB that might have potential to treat affected patient, Epirubicin, Tacrolimus and Methotrexate (Fig. [127]6F). Discussion The previous studies have revealed that the dysregulation of FOXPs is significantly related to the carcinogenesis and progression of many tumors^[128]13–[129]15,[130]31,[131]32. In this study, database analysis showed that the transcription levels of FOXPs in many human tumors were frequently altered. Although the role of FOXPS in the carcinogenesis, development and prognosis of some cancers has been partially elucidated, there have been no further bioinformatics studys of different FOXPs expression and function in breast cancer. This study is the first time to investigate the mRNA expression, gene variation, molecular mechanism, and biological function of different FOXP factors in breast cancer and its influence on the prognosis and immune infiltration in patients with breast cancer through bioinformatics analysis. Among all members of FOXPs, FOXP1 gene has been mapped to chromosome 3p14.1, a region that it has been detected widespread loss of heterozygosity in breast cancer^[132]33, particularly of breast cancer with BRCA2 mutations^[133]34. Fox SB et al.^[134]35 found that nuclear protein expression of FOXP1 was significantly positively related to estrogen receptor status but not associated with tumor size, age, lymph node status, or grade. In addition, FOXP1 co-expression with estrogen receptor significantly improved relapse-free survival,it suggests FOXP1 may function as a tumour suppressor in breast cancer. FOXP1 could modulate cell proliferation in breast cancer cells and improve 5-year recurrence-free survival of patients with tamoxifen-treated breast cancer from Shigekawa T et al.’s study^[135]36. Similarly,a report confirmed that the increased FOXP1 protein expression could predict a good effect to tamoxifen in breast carcinoma patients^[136]37. A systematic review and meta-analysis revealed that decreased FOXP1 protein expression was significantly associated with an unfavorable relapse-free survival (RFS) in breast cancer patients^[137]38. In this study, database analysis showed that the transcription levels of FOXP1 in human breast cancer were lower than in normal tissues, and immunohistochemical staining from our breast carcinoma specimen also demonstrated this result. But the expression of FOXP1 in patients with breast cancer were not associated with the tumor stage. In addition, Our research was similar to Jian X et al. ’s study^[138]38, low FOXP1 transcription levels were associated with poor OS, PFS and DMFS in patients with breast cancer. Recently, some roles of FOXP2 have been verified in cancer development as a tumor suppressor, though its mutations could cause language disorders. Also, Cuiffo et al.^[139]39 found that downregulation of FOXP2 strengthened tumor initiation in breast carcinoma and promoted cancer stem cell metastasis. Furthermore, Chen et al.^[140]40 reported that the transcription level of FOXP2 in breast cancer tissue was also markedly lower than in normal breast tissue and these patients also had poor RFS rate. Similarly, In this study, database analysis found that the transcription levels of FOXP2 in human breast cancer were lower than in normal tissues, and similar result also was found in our specimen by immunohistochemical staining. But the expression of FOXP2 in breast cancer patients has nothing to do with the tumor histological type. It was also found that the low transcription levels of FOXP2 in breast cancer patients correlated with poor OS, PFS. FOXP3 plays an important role in regulating Treg cells development and functions for immune response against cancer^[141]41. FOXP3 also inhibited growth and induced the cell death of a breast cancer cell line MCF-7^[142]42. In addition, Some studies have demonstrated that FOXP3 is an important tumor suppressor of oncogenes in breast cancer with good prognosis^[143]42–[144]44. The database analysis and IHC in this study testified that the expression levels of FOXP3 in human breast cancer were higher than in normal tissues, and its expression levels were not related to the tumor stage, histological type and race. The survival analysis found that high transcription levels of FOXP3 in breast cancer patients resulted in worse PPS and had better PFS. Previous studies indicated that FOXP4 had dual biologic function as a tumor suppressor in patients with kidney cancer^[145]45, or as an oncogene in in patients with hepatocellular carcinoma^[146]46. In present study, the database analysis and IHC revealed that the expression levels of FOXP4 in human breast cancer were higher than in normal tissues, which was consistent with Ma et al.^[147]47 results and its expression levels were not related to the tumor stage, histological type and pathologic stage. In addition, Ma et al.^[148]47 results showed that high expression of FOXP4 predicted a poor OS in breast canccer patients, contrastly, in this study, low FOXP4 expression levels were associated with poor OS and PFS in patients with breast cancer, interestingly, except PPS. Furthermore, a high gene alteration rate of FOXPs was foundin breast cancer patients, and there were difference gene alteration rate in different histological type of breast cancer. Moreover, a mutually exclusive or co-occurring connection between FOXPs or between FOXPs and BRCA or ERBB2 was different, suggesting that these gene play an different role in development of breast cancer. Previous reports showed that FOXP2 and FOXP3 might be a potential biomarker for breast cancer^[149]48,[150]49. Consistantly, in this study, we performed ROC curve analysis. Our results showed that FOXP2 or FOXP3 had a significantly high AUC value in the detection of breast cancer, with 96.8% or 95.7% in accuracy respectively. On the basis of these findings, we conclude that FOXP2 and FOXP3 might act as a potential diagnostic biomarker to differentiate breast cancer from normal normal tissues. Accumulating evidence demonstrated that FOXPs proteins play important roles in the regulation of immune function^[151]50,[152]51. In this work, genomic analysis revealed that FOXPs was closely related to the alteration of extensive immune checkpoints in breast invasive carcinoma. Importantly, the connection between FOXP3 alteration with extensive immune checkpoints was co-occurrence but not mutual exclusivity. Additionally, we found that FOXP3 had a significant positive correlation with PDCD1, CD274, CTLA4 and TMB in breast cancer. This study yet demonstrated that FOXPs were involved in inflammatory response and immune cell infiltration, of note, FOXP3 expression showed a statistically significant correlation with infiltration of B cells, CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells. Consistantly, West et al.^[153]52 reported that the breast cancer patients with FOXP3+ TILs had better survival. These findings strongly indicate that FOXP3 is a potential regulator of immune in breast invasive carcinoma. Previous studies showed that FOXPs were associated with a great deal of genes or proteins, such as TNF receptor family-related gene (GITR)^[154]53, cytooxic T lymphocyte associated antigen 4 (CTLA-4) and CD25^[155]54, TGF-β^[156]55, nuclear factor of activated T cells (NFAT)^[157]55, Runt-related transcription factor 1 (RUNX1)^[158]56. In this study, we also analysed relation of FOXPs and its neighboring genes or proteins, and found that the main 50 genes were associated with FOXPs. Further analysis showed that the functions of these proteins are mainly related to lymphocyte differentiation and regulation of lymphocyte function. Another significant result of this study revealed that FOXP3 expression predicted the breast cancer cells’ response to anticancer drugs, whereas FOXP1, FOXP2 and FOXP4 did not predict. Consistantly, Ladoire et al.’s^[159]47 report showed that FOXP3 expression in breast cancer was independently related to improved OS in patients treated with anthracycline-based adjuvant chemotherapy. Conclusions In conclusion, Our results suggested that BRCA patients with high transcription levels of FOXP1/2/4 had better prognosis and FOXPs was closely related to the alteration of extensive immune checkpoints in breast invasive carcinoma. FOXP3 expression showed a statistically significant correlation with infiltration of B cells, CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells and predicted the breast cancer cellsʼ s response to anticancer drugs, the main 50 genes were involved in FOXPs, our study suggested that FOXPs could serve as a biomarker for tumor detection, therapeutic design and prognosis. Supplementary Information [160]Supplementary Table S1.^ (11KB, xlsx) [161]Supplementary Table S2.^ (10.6KB, xlsx) Abbreviations FOXP Forkhead box P BRCA Breast cancer PPI network Protein protein interaction network Author contributions J.N.Y, S.Y.T, Y.J.Z, L.H.Z and L.Y.L helped with data collection; J.N.Y wrote the manuscript; J.N.Y, J.Z, P.Z.F and C.J.Z designed and oversaw the study. Data availability These data are drawn from the public domain. The datasets used and/or analysed during the current study available from the corresponding author on reasonable request. Competing interests The authors declare no competing interests. Footnotes Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Luyao Liu, Email: liuluyaodoctor@126.com. Peizhi Fan, Email: Fanpzh64@163.com. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-022-12954-3. References