Comprehensive Analysis of Prognostic and immune infiltrates for FOXPs Transcription Factors in Human Breast Cancer

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.


Methods
Oncomine database analysis. The expression level of the FOXPs in various types of cancers was identified in the Oncomine database (https:// www. oncom ine. org/ resou rce/ login. html) 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 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 (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 (www. kmplot. com) 21,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 cBio-Portal (www. cbiop ortal. org) 23,24 . The genomic map contains putative copy-number alterations (CNA) from GIS-TIC, mRNA expression z-scores and Protein expression z-scores mutations. STRINGS analysis. STRINGS (www. string-db. org) is an online tool for analysis of all publicly data of protein-protein interaction (PPI) 25 . In this study, we used a PPI network analysis on FOXPs to explore their functions in human breast cancer.
GeneMANIA analysis. GeneMANIA (www. genem ania. org) is an online tool for analysis of gene functions 26 . In this study, we performed GeneMANIA to select the 50 most important genes to construct genegene interaction network for FOXPs.
Drug-gene interactions. DGidb 27 (https:// dgidb. genome. wustl. edu/), a web server for discovering druggene 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 (www.r-proje ct. org). gene expression data and clinical information were visualizes by R package "ggplot2" R package 28  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. 1A,B). The results showed that the transcriptional levels of FOXP3/4 were up-regulated in BRCA patients, but FOXP2 were downregulated. 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 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. 1C

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. 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. 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. 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. 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. 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. 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. 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. 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. 2).
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. 3B)  www.nature.com/scientificreports/  Table 2. Of note, the FOXP3 alteration showed a statistically significant co-occurrence rather than mutual exclusivity with extensive immune checkpoints, such as 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. 4B). However, FOXP1/2/4 was weakly associated with PDCD1, CD274, CTLA4 (Fig. 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. 5. FOXP1 expressions was positively correlated with the infiltration of CD4+ T cells, CD8+ T cells, neutrophils and macrophages (Fig. 5A). FOXP2 expressions was positively correlated with CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells (Fig. 5B). FOXP3 expression was positively correlated with infiltration of B cells, CD4+ T cells, CD4+ T cells, neutrophils, macrophages and dendritic cells (Fig. 5C). FOXP4 expressions was positively correlated with the infiltration of B cells, CD4+ T cells, neutrophils and macrophages, while it was negatively correlated with   www.nature.com/scientificreports/ FOXPs expression in BRCA to have a more intuitive understanding of the immune infiltration (Fig. 5E). We can see that there were most immune cell infiltration in FOXP2 and FOXP3.  www.nature.com/scientificreports/

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. 6A), the PPI network diagram contains four FOXPs proteins and 10 proteins that are closely related to FOXPs. A GGI network of four FOXPs was constructed, and their functions were analyzed using the GeneMANIA database (Fig. 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 S2) are mainly distributed in the biological process (5 items) (Fig. 6C), the cell component (5 items) (Fig. 6C) and the molecular function (5 items) (Fig. 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. 6C (Table 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. 6D) and KEGG network map to show the connection between pathways (Fig. 6E).
Drug-gene interactions. The result showed expressions of FOXP3 had significant correlations with druggene interaction. A total of 3 drugs were explored using DGIDB that might have potential to treat affected patient, Epirubicin, Tacrolimus and Methotrexate (Fig. 6F).

Discussion
The previous studies have revealed that the dysregulation of FOXPs is significantly related to the carcinogenesis and progression of many tumors [13][14][15]31,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.  www.nature.com/scientificreports/ 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 33 , particularly of breast cancer with BRCA2 mutations 34 . Fox SB et al. 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 36 . Similarly,a report confirmed that the increased FOXP1 protein expression could predict a good effect to tamoxifen in breast carcinoma patients 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 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 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. 39 found that downregulation of FOXP2 strengthened tumor initiation in breast carcinoma and promoted cancer stem cell metastasis. Furthermore, Chen et al. 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 41 . FOXP3 also inhibited growth and induced the cell death of a breast cancer cell line MCF-7 42 . In addition, Some studies have demonstrated that FOXP3 is an important tumor suppressor of oncogenes in breast cancer with good prognosis [42][43][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 45 , or as an oncogene in in patients with hepatocellular carcinoma 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. 47 results and its expression levels were not related to the tumor stage, histological type and pathologic stage. In addition, Ma et al. 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 cooccurring 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 48,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 50,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. 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) 53 , cytooxic T lymphocyte associated antigen 4 (CTLA-4) and CD25 54 , TGFβ 55 , nuclear factor of activated T cells (NFAT) 55 , Runt-related transcription factor 1 (RUNX1) 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 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.

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.