The TIMER database (https://cistrome.shinyapps.io/timer/) includes 10,897 samples across 32 cancer types from The Cancer Genome Atlas (TCGA) to estimate the STAT family gene expression patterns (31, 32). We analyzed the differential expression of any gene of interest in tumors and adjacent normal tissues using the Diff Exp module provided by TCGA. In box plots, gene expression levels were captured, in Wilcoxon tests being used to determine if there were significant differences between levels. Based on the expression status of the STATs family, we performed a survival analysis using the Sangerbox database (http://vip.sangerbox.com/home.html). A median gene RNA-seq by Expectation Maximization (RSEM) value was used to classify tumor samples into high and low groups. The overall survival (OS) outcomes were analyzed by calculating the Hazard Ratio (HR) and 95% CI, along with a log-rank p-value.

cBioPortal (http://www.cbioportal.org/) was used to analyze genomic changes to STAT genes in pan-cancer (33). All types of Pan-Cancer Atlas available (32 studies and 10,967 samples) were selected for calculation. According to standard software parameters, altered frequencies of each gene were determined and the fraction of genome altered and survival curves were calculated.

By using GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/), the Genotype-Tissue Expression (GTEx) database was analyzed for expression levels of STAT family genes in each tissue (34). The STAT family genes were investigated for single-nucleotide variations (SNVs), copy number variations (CNVs), methylation, pathway function, and drug sensitivity. In selected cancer types, the SNV module displays the frequency and variant types of each gene set. The statistics of heterozygous and homozygous CNV of each cancer type were displayed as pie charts for gene sets on the CNV module. For selected cancer types, the Methylation module examined differences in methylation between tumors and matching normal tissue, the relationship between methylation and expression, and the impact of methylation on OS. The Pathway Activity module displayed differences in gene expression across pathway activity categories (activation and inhibition), which were determined by pathway scores. Drug sensitivity and gene expression profiling data from cancer cell lines in GDSC were combined in the Drug Sensitivity module. Spearman’s correlation analysis was used to compare the expression of each gene in the gene set to the sensitivity to small molecules and drugs [50% inhibiting concentration (IC50)].

With Xena Shiny (https://shiny.hiplot.com.cn/ucsc-xena-shiny/), we could analyze immune cell infiltration across 32 cancer types from TCGA. We used the module TCGA: Associations between Molecular Profile and Immune Signature module of Xena Shiny to evaluate the interrelationships between STAT expression in 20 immune cells and their infiltration by using “CIBERSORT.” Pearson’s correlation coefficient of gene and immune infiltration score in each tumor was also generated.

The stromal and immune scores are positively related to the stromal and immune components in the tumor microenvironment (TME). In light of estimate scores, a composite score combining stromal scores and immune scores can determine the relative proportion of stromal and immune components in the TME. Sangerbox tool was used to calculate stromal, immune, and estimate scores. Spearman’s correlation coefficient of gene and immune infiltration score in each tumor was also generated by the Sangerbox tool. The Spearman correlation test was used to determine the association between tumor stemness and each STAT expression. An RNA stemness score (RNAss) based on mRNA expression and a DNA stemness score (DNAss) based on DNA methylation patterns were used in this pan-cancer study to measure tumor stemness.

GeneMANIA was used in bioinformatics methods to predict functions for genes or gene lists, as well as to construct PPI networks (35).

To perform enrichment analyses, related genes of STAT families were discovered using GEPIA (http://gepia.cancer-pku.cn/index.html) (36). Metascape (http://metascape.org) was used to analyze the function of STAT members and related genes (37). Gene Ontology (GO) analysis in Metascape could identify STATs and the similar genes by three categories, including molecular function (MF), cellular component (CC), and biological process (BP). By analyzing the Kyoto Encyclopedia of Genes and Genomes (KEGG), we identified the signaling pathways associated with STAT factors, as well as related genes.

The Spearman’s correlation test or the Pearson’s correlation test was used to analyze correlations. Cox proportional hazards models were calculated to determine survivorship risk and HR. Each variable was analyzed using survival plots and compared with log-rank tests. Once two sets of data were detected, a p-value of 0.05 was declared significant.

We examined the expression levels of STAT family genes in all 33 cancer types available in TCGA pan-cancer data. Compared to other adjacent normal tissues, STATs tended to be extensively changed, suggesting a statistically significant difference in nearly half of all pan-cancers. STAT1 was highly expressed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal adenocarcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC) and poorly expressed in kidney chromophobe (KICH), and the difference was statistically significant ( Figure 2A ). STAT2 was significantly upregulated in BLCA, CHOL, ESCA, HNSC, kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), LIHC, LUAD, LUSC, prostate adenocarcinoma (PRAD), and THCA but significantly downregulated in CHOL ( Figure 2B ). STAT3 expression remained high in CHOL, ESCA, HNSC, and STAD and low in BLCA, BRCA, COAD, KICH, KIRC, KIRP, LUAD, LUSC, and PRAD, and the difference was statistically significant ( Figure 2C ). STAT4 was significantly upregulated in ESCA, HNSC, KIRC, KIRP, STAD, and THCA but significantly downregulated in BRCA, COAD, KICH, LUSC, and rectum adenocarcinoma (READ) ( Figure 2D ). STAT5A was highly expressed in CHOL, ESCA, HNSC, LIHC, STAD, and THCA and poorly expressed in BLCA, BRCA, KICH, LUAD, LUSC, PRAD, and UCEC, and the difference was statistically significant ( Figure 2E ). STAT5B expression remained high in CHOL and LIHC and low in BLCA, BRCA, KICH, KIRP, LUAD, LUSC, PRAD, and THCA, and the difference was statistically significant ( Figure 2F ). STAT6 was significantly upregulated in CHOL, ESCA, HNSC, KIRC, KIRP, LIHC, STAD, and THCA but significantly downregulated in BLCA, BRCA, COAD, LUAD, LUSC, PRAD, and UCEC ( Figure 2G ).

To investigate the prognostic value of the seven STAT factors in pan-cancer, survival analysis was performed. High expression of STAT1 was significantly linked with the shortened OS in glioma (GBMLGG) [HR = 1.72 (1.54, 1.92), p = 9.9e-22], brain lower-grade glioma (LGG) [HR = 1.84 (1.56, 2.17), p = 2.1e-13], LUAD [HR = 1.17 (1.01, 1.35), p = 0.04], KIRP [HR = 1.54 (1.18, 2.01), p = 1.6e-3], pan-kidney cohort (KICH+KIRC+KIRP) (KIPAN) [HR = 1.25 (1.12, 1.38), p = 5.0e-5], uveal melanoma (UVM) [HR = 1.37 (1.09, 1.70), p = 4.4e-3], acute myeloid leukemia (LAML) [HR = 1.22 (1.06, 1.41), p = 5.3e-3], and adrenocortical carcinoma (ACC) [HR = 1.79 (1.24, 2.59), p = 2.0e-3], while lower STAT1 expression was significantly associated with lower OS rates in sarcoma (SARC) [HR = 0.82 (0.68, 0.98), p = 0.03], skin cutaneous melanoma (SKCM) [HR = 0.78 (0.71, 0.85), p = 1.1e-7], and ovarian serous cystadenocarcinoma (OV) [HR = 0.86 (0.79, 0.94), p = 1.3e-3] ( Figure 3A ). Higher STAT2 expression was significantly associated with poorer OS in GBMLGG [HR = 1.71 (1.45, 2.01), p = 2.5e-10], LGG [HR = 1.95 (1.55, 2.46), p = 3.3e-8], KIPAN [HR = 1.69 (1.44, 1.98), p = 3.1e-10], KIRC [HR = 1.62 (1.32, 1.99), p = 4.9e-6], LAML [HR = 1.22 (1.05, 1.42), p = 7.5e-3], and ACC [HR = 1.81 (1.08, 3.04), p = 0.02] ( Figure 3B ). GBMLGG, LGG, KIPAN, glioblastoma multiforme (GBM), and UVM patients with high STAT3 expression had worse OS [HR = 2.50 (1.98, 3.15), p = 1.4e-14; HR = 2.56 (1.84, 3.56), p = 4.0e-8; HR = 1.20 (1.03, 1.41), p = 0.02; HR = 1.49 (1.08, 2.06), p = 0.01; HR = 1.88 (1.13, 3.14), p = 0.01, respectively] than those with low STAT3 expression, while colon adenocarcinoma/rectum adenocarcinoma esophageal carcinoma (COADREAD) and SKCM patients with low STAT3 expression had worse OS [HR = 0.71 (0.51, 0.99), p = 0.04; HR = 0.79 (0.67, 0.93), p = 6.1e-3), respectively] than those with high STAT3 expression ( Figure 3C ). High expression of STAT4 was significantly linked with the shortened OS in KIRP [HR = 1.50 (1.15, 1.96), p = 3.2e-3], KIPAN [HR = 1.30 (1.17, 1.45), p = 1.4e-6], and GBM [HR = 1.28 (1.07, 1.53), p = 7.4e-3], while lower STAT4 expression was significantly associated with lower OS rates in BRCA [HR = 0.88 (0.79, 0.98), p = 0.02], SARC [HR = 0.85 (0.74, 0.98), p = 0.02], SKCM [HR = 0.81 (0.75, 0.88), p = 1.3e-7], OV [HR = 0.88 (0.80, 0.96), p = 6.0e-3], and pancreatic adenocarcinoma (PAAD) [HR = 0.82 (0.69, 0.98), p = 0.03] ( Figure 3D ). Higher STAT5A expression was significantly associated with poorer OS in GBMLGG [HR = 1.86 (1.61, 2.15), p = 4.7e-17], LGG [HR = 1.79 (1.46, 2.18), p = 6.3e-9], and testicular germ cell tumors (TGCTs) [HR = 7.02 (1.07, 46.31), p = 0.04], while low expression of STAT5A was significantly associated with lower OS rates in SARC [HR = 0.64 (0.49, 0.82), p = 4.9e-4], KIRP [HR = 0.63 (0.46, 0.87), p = 5.8e-3], HNSC [HR = 0.86 (0.75, 0.98), p = 0.03], SKCM [HR = 0.86 (0.75, 0.99), p = 0.04], mesothelioma (MESO) [HR = 0.54 (0.33, 0.89), p = 0.01], and UVM [HR = 0.52 (0.28, 0.99), p = 0.05] ( Figure 3E ). GBMLGG, KIRC, SKCM, and PAAD patients with low STAT5B expression had worse OS [HR = 0.51 (0.41,0.64), p = 5.4e-9; HR = 0.67 (0.56, 0.80), p = 1.2e-5; HR = 0.75 (0.61, 0.91), p = 3.6e-3; HR = 0.77 (0.60, 0.97), p = 0.03, respectively] than those with high STAT5B expression ( Figure 3F ). High expression of STAT6 was significantly linked with the shortened OS in GBMLGG [HR = 1.99 (1.68, 2.37), p = 4.8e-15], LGG [HR = 1.97 (1.53, 2.52), p = 1.2e-7], UVM [HR = 2.83 (1.34, 5.97), p = 5.5e-3], and LAML [HR = 1.39 (1.19, 1.63), p = 3.2e-5], while lower STAT6 expression was significantly associated with lower OS rates in SARC [HR = 0.67 (0.54, 0.84), p = 5.3e-4] and BLCA [HR = 0.83 (0.70, 0.98), p = 0.03) ( Figure 3G ).

Using the cBioPortal tool on 10,967 samples in 32 studies, we calculated the mutation frequency of the seven STAT genes. The gene of STAT genes was altered in 1,003 (9%) samples; STAT1, STAT2, STAT3, STAT4, STAT5A, STAT5B, andSTAT6 were altered in 2.3%, 1.7%, 2.0%, 2.4%, 1.6%, 1.9%, and 1.9% of the demanded pan-cancer samples ( Figure 4A ). We found that mutations of STAT family genes more frequently occur in melanoma, mature B-cell neoplasms, UCEC, esophagogastric adenocarcinoma, BLCA, NSCLC, COADREAD, CHOL, and cervical squamous cell carcinoma (>10%) ( Figure 4B ). An analysis of SNVs found the highest mutation levels in STATs in tumor tissues in UCEC and SKCM. These mutations increased expression levels in genes. In the greatest number of tumors, STAT4 was correlated with tumor growth, then STAT1, STAT3, STAT6, STAT5B, STAT2, and STAT5A ( Figure 4C ). The alterations were not related to the survival of STAT family gene mutations, and the OS and disease-free survival (DFS) of patients with alterations were not shortened as opposed to those without alterations (all p > 0.05; Figures 4D, E ).

According to the CNV pie chart, most types of CNV were heterozygous amplifications and deletions ( Figure 5A ). ACC and TGCT showed heterozygous amplification of STAT6; KIRP displayed STAT5A, STAT5B, and STAT3; and ACC and TGCT showed greater than 50% amplification of STAT2 (p < 0.05, Figure 5B ). OV and KICH exhibited heterozygous deletions of STAT5B, STAT5A, and STAT3; STAT4 and STAT1 in KICH were both greater than 50% deleted (p < 0.05, Figure 4B ). As a result of CNV, expression levels of every STAT factor increased. This interrelationship has been found in the largest number of tumors for STAT5B, followed by STAT3, STAT6, STAT2, STAT5A, STAT1, and STAT4 ( Figure 5C ).

To identify the epigenetic regulation in TCGA database, the methylation characteristics of STAT genes were examined. In different tumors, methylation of STAT genes was heterogeneous: in PRAD, LUSC, BRCA, UCEC, and COAD, more highly hypermethylated genes could be found, while in THCA, KIRC, and LIHC, more highly hypomethylated genes can be found ( Figure 6A ). A correlation analysis of the methylation levels and mRNA expression levels revealed that a large proportion of genes, especially STAT1, STAT5A, and STAT6, correlated negatively with the methylation levels ( Figure 6B ). A survival study of various types of tumors showed that hypermethylation of STAT5A and STAT5B was mainly linked to poorer survival, while hypomethylation of STAT1, STAT2, and STAT4 was generally linked to shorter survival (p < 0.05, Figure 6C ).

The linked pathway network indicated that STAT family genes acted significantly in nine famous cancer-related signaling pathways. For STAT1, the principal inactivated pathways were apoptosis, cell cycle, and EMT. STAT2 was mostly involved in the activation of apoptosis and EMT, while the key pathways inactivated were cell cycle, DNA damage response, and hormone Androgen Receptor (AR). STAT3 was mostly involved in the activation of EMT, hormone Estrogen Receptor (ER), rat sarcoma/mitogen activated protien kinase RAS/MARK, and Receptor tyrosine kinase (RTK), while the major inactivated pathways were cell cycle and DNA damage response. STAT4 was generally involved in the inhibition of DNA damage response and hormone AR, while the main activated pathways were apoptosis, EMT, and hormone ER. STAT5A was mainly involved in the activation of EMT and hormone ER, while the main inactivated pathways were cell cycle and PI3K/AKT. In STAT5B, the primary inactivated pathways were apoptosis and cell cycle. STAT6 was mostly involved in RAS/MARK activation, while the main inactivated pathways were apoptosis, cell cycle, DNA damage response, and EMT (all p < 0.05, Figure 7 ).

According to the drug sensitivity analysis, low levels of STAT5B and STAT5A showed resistance to 56 and 42 drugs, respectively. Drug resistance toward vorinostat, tubastatin A, NPK76-II-72-1, I-BET-762, TPCA-1, TL-1-85, NG-25, navitoclax, and methotrexate negatively correlated with STAT5B expression (correlation coefficient >-0.20). Similarly, drug resistance toward BHG712, TPCA-1, TL-1-85, and NG-25 negatively related to the STAT5A expression (correlation coefficient >-0.20) ( Figure 8 ).

Scatterplot analysis revealed Spearman’s correlation coefficient between each STAT family gene and the tumor infiltration of each cell type analyzed in different cancer types ( Figures 9A-G ). STAT1 expression positively correlated with the infiltration of T regulatory cells (Tregs), T follicular helper cells (TFHs), CD8+ and activated memory CD4+ T cells, M1 macrophages, resting dendritic cells (DCs), and naive B cells in most tumors and negatively related to the infiltration of naive CD4+ T, plasma, resting NK, monocyte, resting and activated mast cells, M2 and M0 macrophages, activated DCs, and memory B cells in the majority of tumors ( Figure 9A ). Infiltration of resting memory CD4+ T cells and M1 macrophages positively associated with STAT2 expression; however, STAT2 negatively correlated with the infiltration of activated mast and memory B cells in most cancer types ( Figure 9B ). STAT3 expression positively correlated with the infiltration of resting memory CD4+ T and naive B cells in most tumors and negatively associated with the infiltration of TFHs, CD8+ T, activated NK, and memory B cells in different cancer types ( Figure 9C ). Infiltration of Tregs, TFHs, CD8+, resting and activated memory CD4+ T cells, M1 macrophages, resting DCs, and naive B cells positively correlated with STAT4 expression; however, STAT4 negatively correlated with the infiltration of naive CD4+ T, resting NK and resting mast cells, M0 macrophages, activated DCs, and memory B cells in most tumor types ( Figure 9D ). STAT5A expression positively correlated with the infiltration of Tregs, CD8+, resting and activated memory CD4+ T cells, M1 macrophages, resting DCs, and naive B cells in most tumors and negatively correlated with the infiltration of naive CD4+ T cells, M0 macrophages, activated DCs, and memory B cells in the majority of tumors ( Figure 9E ). The levels of STAT5B positively associated with the infiltration of resting memory CD4+ T, resting mast, and naive B cells in most tumors; however, STAT5B negatively correlated with the infiltration of Tregs, gammadelta T cells, TFHs, activated NK cells, and memory B cells in different cancer types ( Figure 9F ). STAT6 expression positively associated with the infiltration of resting memory CD4+ T, monocyte, and resting mast cells in most tumors and negatively related to the infiltration of gammadelta T cells, activated memory CD4+ T cells, and M0 macrophages in the majority of tumors ( Figure 9G ).

Higher stromal or immune scores indicated that stromal or immune elements have infiltrated deeper into the TME. Among all STAT genes, but not STAT6, expression was positively related to stromal scores. Specifically, STAT1 was strongly correlated with READ, LAML, and COAD. STAT2 was strongly correlated with COAD, LAML, and READ. STAT3 was strongly correlated with TGCT, PAAD, and LAML. STAT4 was strongly correlated with ACC, BLCA, BRCA, COAD, ESCA, HNSC, KICH, LUSC, OV, PAAD, PRAD, READ, SARC, SKCM, THCA, and UCEC. STAT5A was strongly correlated with CHOL, LGG, LUAD, PAAD, PRAD, and TGCT. STAT5B was strongly correlated with PAAD, READ, and STAD, and STAT6 was strongly correlated with GBM and TGCT. The expressions of STAT1, STAT2, STAT3, STAT4, and STAT5A were mostly positively related to immune scores, and STAT5B and STAT6 expressions were partly positively related to immune scores. To be specific, STAT1 was closely related to BLCA, CESC, COAD, DLBC, HNSC, KIRC, READ, SKCM, TGCT, THCA, and UVM. STAT2 was closely related to COAD, DLBC, and READ. STAT3 was closely related to COAD, DLBC, and READ. STAT4 was closely related to ACC, BLCA, BRCA, CESC, COAD, DLBC, ESCA, HNSC, KICH, KIRC, LUAD, LUSC, MESO, OV, PAAD, PRAD, READ, SARC, SKCM, STAD, TGCT, THCA, UCEC, and UCS. STAT5A was closely related to CHOL, DLBC, HNSC, KIRC, LUAD, LUSC, PAAD, PRAD, and TGCT. STAT5B was positively correlated with DLBC but negatively related to GBM and SARC. All STAT genes except STAT6 exhibited a positive correlation with estimate scores. In particular, STAT1 was highly associated with BLCA, COAD, DLBC, KIRC, READ, SKCM, TGCT, THCA, and UVM. STAT2 was highly associated with COAD and READ. STAT3 was highly associated with COAD, DLBC, LGG, PAAD, and READ. STAT4 was highly associated with ACC, BLCA, BRCA, CESC, COAD, DLBC, ESCA, HNSC, KICH, KIRC, LUAD, LUSC, MESO, OV, PAAD, PRAD, READ, SARC, SKCM, STAD, TGCT, THCA, UCEC, and UCS. STAT5A was highly associated with CHOL, LUAD, LUSC, PAAD, PRAD, and TGCT. STAT5B was highly associated with READ and DLBC ( Table 1 ; Figure 10 ).

In DNAss, the expression of STAT1 was positively related to GBMLGG, LGG, PRAD, THYM, THCA, and PCPG, while STAT1 expression was negatively related to CESC, ESCA, Stomach and Esophageal carcinoma (STES), SARC, KIPAN, STAD, UCEC, HNSC, KIRC, LUSC, LIHC, UCS, and BLCA ( Figure 11A ). The expression of STAT2 was positively related to GBMLGG, LGG, KIRP, PRAD, THYM, THCA, and PCPG, while STAT2 expression was negatively related to COAD, COADREAD, BRCA, ESCA, STES, STAD, HNSC, MESO, and BLCA ( Figure 11B ). The expression of STAT3 was positively related to GBMLGG, LGG, THYM, THCA, ACC, and CHOL, while STAT3 expression was negatively related to COAD, COADREAD, LAML, BRCA, SARC, STAD, PAAD, TGCT, and BLCA ( Figure 11C ). The expression of STAT4 was positively related to LGG, PARD, THCA, and PCPG, while STAT4 expression was negatively related to CESC, LUAD, COAD, ESCA, COADREAD, BRCA, STES, SARC, HNSC, KIPAN, STAD, UCEC, KIRC, LUSC, LIHC, PAAD, UCS, BLCA, KICH, and DLBC ( Figure 11D ). The expression of STAT5A was positively related to GBMLGG, LGG, PARD, THYM, THCA, PCPG, and ACC, while STAT5A expression was negatively related to GBM, LUAD, COAD, UCEC, COADREAD, BRCA, STES, SARC, KIPAN, STAD, HNSC, LUSC, LIHC, OV, TGCT, UCS, BLCA, and DLBC ( Figure 11E ). The expression of STAT5B was positively related to LIHC and PCPG, while STAT5B expression was negatively related to GBMLGG, LGG, COAD, LAML, COADREAD, BRCA, STES, SARC, STAD, UCEC, HNSC, LUSC, PAAD, and BLCA ( Figure 11F ). The expression of STAT6 was positively related to GBMLGG, LGG, LUSC, THYM, THCA, PCPG, and KICH, while STAT6 expression was negatively correlated with LUAD, BRCA, SARC, PAAD, and TGCT ( Figure 11G ).

In RNAss, the expression of STAT1 was positively associated with LUAD, BRCA, STES, STAD, and PCPG, while STAT1 expression was negatively associated with GBMLGG, LGG, COAD, COADREAD, KIPAN, HNSC, THYM, LIHC, THCA, READ, and PAAD ( Figure 11H ). The expression of STAT2 was almost negatively associated with all tumor types except ACC, SKCM, and PCPG ( Figure 11I ). The expression of STAT3 was negatively associated with all tumor types ( Figure 11J ). The expression of STAT4 was positively associated with PCPG, while STAT4 expression was negatively related to all other tumor types ( Figure 11K ). The expression of STAT5A was positively associated with PCPG, while STAT5A expression was negatively related to all other tumor types ( Figure 11L ). The expression of STAT5B was positively associated with GBMLGG, LGG, and PCPG, while STAT5B expression was negatively related to all other tumor types ( Figure 11M ). The expression of STAT6 was negatively associated with GBM, GBMLGG, LGG, CESC, LUAD, COAD, COADREAD, LAML, BRCA, SARC, KIRP, KIPAN, PRAD, UCEC, KIRC, LUSC, THYM, LIHC, THCA, PAAD, TGCT, PCPG, SKCM, and BLCA ( Figure 11N ).

We constructed a PPI network for each of the STAT genes using the GeneMANIA database. A hub node representing each member of the STAT family was surrounded by 20 nodes corresponding to genes that were significantly correlated with their members ( Figure 12A ).

GEPIA was used to identify the top 20 genes that are similar to those of each STAT member in order to explore possible mechanisms through which STATs can contribute to pan-cancer pathology ( Table 2 ). After removing duplicates, there were 135 genes, including seven STAT family factors and 128 similar genes. To catch the relationship between terms, the Metascape database was used in GO and KEGG enrichment analyses of STATs and their relative genes. Three aspects, including MF, CC, and BP, were considered in GO enrichment analysis for predicting target host gene functions. As displayed in Figure 12B , we found that GO:0045087 (innate immune response), GO:0001819 (positive regulation of cytokine production), GO:0046649 (lymphocyte activation), GO:0034097 (response to cytokine), GO:0032479 (regulation of type I interferon production), GO:0034340 (response to type I interferon), GO:0043122 (regulation of I-kappaB kinase/NF-kappaB signaling), GO:0140289 (protein mono-ADP-ribosylation), GO:0045648 (positive regulation of erythrocyte differentiation), GO:0042803 (protein homodimerization activity), GO:0042393 histone binding, GO:0061629 (RNA polymerase II-specific DNA-binding transcription factor binding), GO:0008047 (enzyme activator activity), GO:0001228 (DNA-binding transcription activator activity, RNA polymerase II-specific), GO:0050778 (positive regulation of immune response), GO:0000209 (protein polyubiquitination), GO:0001779 (NK cell differentiation), GO:0060216 (definitive hemopoiesis), GO:0043124 (negative regulation of I-kappaB kinase/NF-kappaB signaling), and GO:0042113 (B cell activation) were prominently related to STAT factors and their similar genes. Furthermore, six pathways related to the functions of the STAT family were found through KEGG analysis; pathways such as hsa05162: measles, hsa05161: hepatitis B, hsa00310: lysine degradation, hsa05340: primary immunodeficiency, hsa05145: toxoplasmosis, and hsa04664: Fc epsilon RI signaling pathway were involved in the tumorigenesis and pathogenesis of pan-cancer ( Figure 12C ). Additionally, the networks of enriched GO terms and KEGG pathways were displayed in Figures 12D, E .