{'fraction_non_alphanumeric': 0.032061266874350985, 'fraction_numerical': 0.023148148148148147, 'mean_word_length': 5.2215343203230145, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identifying the potential regulators of neutrophils recruitment in hepatocellular carcinoma using bioinformatics method

Background Neutrophils play a crucial role in the development and progression of hepatocellular carcinoma (HCC); however, the mechanism underlying neutrophil recruitment is not fully understood. Therefore, we aimed to explore the potential genes or pathways related to neutrophil recruitment in the cancer microenvironment. Methods We downloaded TCGA HCC gene expression profiles, the abundance of 22 different immune cells in HCC patients, and patient survival information. We used Kaplan-Meier survival analysis to determine if neutrophils were related to survival. Next, we screened different expression genes (DEGs) between patients with high and low level of neutrophils. We then identified the transcription factor and its targets in the fence of DEGs. Then, we carried out enrichment analysis and gene set variation analysis (GSVA) for targets. Finally, we explored the potential mechanism of targets via calculating correlation scores. Results Our survival analysis results showed that neutrophils were significantly associated with patient survival. A total of 736 DEGs were screened. Next, we identified transcription factor larger E26 transformation-specific (ETS) homologous factor (EHF) and 702 targets of EHF from 736 DEGs. Among these targets, the level of FGD6 expression had the highest correlation with the level of EHF expression. Enrichment and GSVA analysis for FGD6 showed that the level of GO:0043547 had a positive regulatory effect on GTPase activity and the GO:0007010 cytoskeleton organization was significantly difference between the high and low neutrophils counts. By calculating the correlation between FGD6 and genes in GO:0043547 and GO:0007010, we identified RIC8B and SIPA1L3. Conclusions These findings demonstrated that transcription factor EHF can influence recruitment of neutrophils by mediating the transcription of FGD6. Further investigations are needed to shed new light on EHF and its target FGD6.


Introduction
Currently, hepatocellular carcinoma (HCC) is the most frequent liver malignancy worldwide, and the morbidity and mortality rates are high (1). Within 5 years after curative resection, metastasis and recurrence of HCC are up to 70% (2). The tumor immune microenvironment plays a crucial role in initiation, progression, metastasis, and recurrence of HCC. As major components in the tumor immune microenvironment, neutrophils have been increasingly recognized as important factors in cancer behavior and progression (3). Recent evidence has demonstrated that neutrophils contribute to the metastatic process of tumor cells via cell-to-cell interactions (4). Evidence suggests that the neutrophil-to-lymphocyte ratio can predict the prognosis of HCC (5). Therefore, a deep understanding of the mechanisms leading to the interaction between neutrophils and cancer cells will provide better strategies to orchestrate the immune system to eradicate cancers.
Recently, a number of studies have focused on the roles of neutrophils in the tumor microenvironment. Neutrophils have been confirmed to induce the stem cell characteristics in HCC cells by secreting BMP2 and TGFB2 (6). Xiao et al. (7) reported that neurotensin increases the neutrophil count of the local HCC microenvironment by inducing the expression of IL-8. Haider et al. (8) concluded that the CXCL5-dependent attraction of neutrophils depends on TGF-β and Axl signaling. Despite substantial progress in research involving the roles of neutrophils in the tumor microenvironment, the precise molecular mechanisms underlying neutrophil recruitment are only partially understood.
In this study we showed that the transcription factor, EHF, and its target, FGD6, had significantly different levels of expression in patients with high-and lowlevel neutrophil counts. The complete nomenclature of transcription factor EHF is larger E26 transformationspecific (ETS) homologous factor. Increasing evidence over recent years has revealed that EHF has a critical function in tumorigenesis. It has been demonstrated that EHF activates the TGF-β signaling pathway in colorectal carcinoma cells (9). Increased expression of EHF leads to poor survival in gastric cancer by activating HER family signaling (10). Liu et al. (11) reported that the level of EHF expression has a close relationship with the efficacy of anti-PD-1 therapy in pancreatic ductal adenocarcinoma. Our outcome showed that the expression level of FGD6 had significant correlation with the expression level of RIC8B and SIPA1L3. Previous evidence demonstrates that FGD6 can interact with CDC42 which contributes to cell adhesion, cell polarity, and membrane recycling during bone degradation (12). Taken together, our data suggest that the transcription factor, EHF, influenced recruitment of neutrophils by mediating the transcription of FGD6. EHF-FGD6-CDC42 plays an important role in neutrophil recruitment.

HCC data sets downloading and preprocessing
RNA-seq data and clinical data of HCC patient samples were extracted from the National Cancer Institute Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/ pancanatlas). We performed log2 normalization for RNAseq data. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Assess and analyze the infiltration of immune cells in HCC
Infiltrations of immune cells in HCC were assessed using the CIBERSORT algorithm (13) with default parameters. Next, we carried out Kaplan-Meier survival analysis to screen the immune cells that were significantly correlated with patient survival. Patients were divided into high and low neutrophil groups based on the median neutrophil count. Kaplan-Meier survival analysis was implemented for the neutrophil fraction using the survival package in R (14). The P values for the Kaplan-Meier survival analysis were adjusted using the false discovery rate (FDR) method. We screened the different expression genes (DEGs) in the two groups utilizing the limma package in R (15) [filter criteria: fold change (FC) ≥1.5; adjusted P value <0.01]. The P values were further adjusted using the Benjamini-Hochberg (BH) multiple testing correction method.

Construct prognosis prediction model using least absolute shrinkage and selection operator (LASSO) regression
In this step, we constructed a prognostic predictive module using LASSO regression. We screened survival-related genes using univariate Cox proportional risk regression model in survival R package (14). P<0.001 was selected as the threshold. We screened the most powerful prognostic markers from the survival-related genes using LASSO regression in the glmnet R package (16). The workflow of the prognostic predictive module was as follows: (I) patients were randomly divided into training (70% patients) and test sets (30% patients); (II) the training set was trained on the prognostic model and the performance of the model on the test set was assessed; (III) and the concordance index (c index) of the prognostic model in the training and test sets were calculated utilizing the Hmisc R package (17). We used a univariate Cox proportional risk regression model for age, gender, stage, stage T, stage N, stage M, cigarette smoking, and subtypes. The survival-related clinical features were selected for multivariate Cox proportional risk regression for risk scores and abundance of immune cells.

Identifying transcription factor and the targets
The detailed information of human transcription factor was downloaded from Ensembl (18). We screened all transcription factors from DEGs based on information about transcription factors downloaded from Ensembl. Next, we selected the transcription factor with the maximum logFC for further analysis. The ChIP-seq outcome of transcription factor was downloaded from Cistrome Data Browser (19). Yan et al. (20) contributed to the ChIP-seq outcome. We selected the genes which were at the intersection of the ChIP-seq outcome and DEGs as the target of transcription factor for further analysis. We uploaded genes we selected to DAVID (21) and assessed the gene ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Next, we calculated the correlation between the transcription factor and the differentially-expressed targets. The target that had the greatest correlation was selected for further analysis.

Function analysis of transcription factor targets
We selected the target that had the greatest correlation with the transcription factor. Then, we submitted the target to DAVID and obtained the GO functions and KEGG pathways. The GSVA package in R (22) was used to calculate the GO function and KEGG pathway scores in all patients. The limma package in R to screen GO functions and KEGG pathways were shown to be different between the high and low neutrophil group (filter criteria: log2FC ≥0.1; adjusted P value <0.01). To further evaluate the possible mechanisms underlying the target with the greatest correlation, we calculated the correlation between the target with the greatest correlation and the screened genes with different GO functions and KEGG pathways. The functions and pathways of significantly correlated genes were annotated by DAVID.

Statistical analysis
Wilcoxon test was utilized to compare continuous variables. Correlations were calculated with Spearman correlation. All the tests were two-sided, and P<0.05 was regarded as indicating significance, unless otherwise stated. Survival analysis was performed using the Kaplan-Meier method, and the survival of the clusters was compared using the log-rank test. The FDR correction was utilized in multiple tests to decrease false positive rates. All of the analyses were performed with R software (version 3.4.2, http://www. R-project.org).

HCC data sets downloading and pre-processing
RNA-seq data and clinical data of 323 HCC patients were extracted from the National Cancer Institute Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/ pancanatlas). We performed log2 normalization for RNAseq data.

Assessment and analysis of immune cell infiltration in HCC
Utilizing the CIBERSORT algorithm, we assessed infiltration of 22 immune cells in HCC ( Figure 1A). Next, we divided patients into high and low immune cell groups based on the median of each immune cell count. We carried out Kaplan-Meier survival analysis for 22 immune cells. The P values of Kaplan-Meier survival analysis were adjusted using the FDR method ( Table 1). Our outcome showed that neutrophil infiltration was significantly associated with patient survival ( Figure 1B). Next, we screened DEGs between the high (n=147) and low neutrophil groups (n=206) using the limma R package. A total of 736 DEGs were screened between the high and low neutrophil groups ( Figure 2). Five hundred sixty-three DEGs had significantly high expression, whereas 173 DEGs were downregulated in the high neutrophil group.

Construction of a prognosis predictive model using LASSO regression
A total of 323 patients were selected to fit the prognosis predictive model. We screened 1,616 survival-related genes using a univariate Cox proportional risk regression model with a P<0.001 threshold. Next, we performed LASSO regression to screen survival-related markers. A total of three survivalrelated markers were screened. The risk score for each patient was calculated by combining the expression of the prognosisrelated markers with the corresponding LASSO coefficients. The formula of the risk score is listed as follows: Risk score RAP GAP To assess the performance of our model, we plotted cross-validated time-dependent receiver operating characteristic (ROC) curves. The 12-month area under the curve (AUC) was 0.823 ( Figure 3A). The 36-and 60-month AUCs were both >0.79. The AUCs indicated that our model had good performance in prognosis prediction both   in the short-and long-term. We calculated the optimal cut-off value according to the optimum sensitivity and specificity of the 5-year survival ROC curve. Patients with a risk score ≥−0.027 were assigned to the high-risk group, and the remaining patients comprised the low-risk group. The survival of patients in the low-risk group was significantly better than patients in the high-risk group ( Figure 3B). The number of patients in high risk group was less than the number of patients in low risk group ( Figure 4A). The highrisk group had more deaths than the low-risk group ( Figure  4B). As presented in Figure 4C, RAP1GAP and SLC39A1 in patients in the high-risk group exhibited high expression levels while ZNF23 showed the opposite tendency. Finally, we carried out multivariate Cox regression analysis and the outcome revealed that our prognosis predictive model independently predicted prognosis ( Table 2).

Identifying the transcription factor and the targets
The information of human transcription factor was downloaded from Ensembl. Through comparison, we identified 18 transcription factors in DEGs ( Figure 5A). Transcription factor EHF with maximum logFC (logFC =1.1717) was selected for further analysis ( Figure 5B). We downloaded the ChIP-seq outcome of EHF from Cistrome Data Browser. Based on the ChIP-seq outcome, we screened 702 targets of transcription factor EHF in DEGs. Among these targets, the expression level of FGD6 was most correlated with transcription factor EHF ( Figure 5C). Next, we submitted the targets of EHF to DAVID and carried out pathways and GO function analysis. The outcome revealed that the targets of EHF mainly enriched in GO:0016021 integral component of membrane, GO:0006954 inflammatory response, GO:0005887 integral component of plasma membrane, GO:0005886 plasma membrane and GO:0005615 extracellular space ( Figure 6).

Function analysis of transcription factor targets
We submitted FGD6 to DAVID and obtained 16 GO function terms and zero KEGG pathways. GSVA and limma analysis demonstrated that the high and low neutrophil groups had significant differences in 2 GO terms [GO:0043547 (positive regulation of GTPase activity) and GO:0007010 (cytoskeleton organization); Figure 7]. To further evaluate the possible mechanisms underlying FGD6 in GO:0043547 and GO:0007010, we calculated the correlation between FGD6 and genes in the two GO terms. The outcome showed that the RIC8B and SIPA1L3 genes had the greatest correlation ( Figure 8). The KEGG pathways analysis utilizing DAVID showed that SIPA1L3 was involved in the hsa04015: Rap1 signaling pathway.

Discussion
Our findings demonstrated that the neutrophil count was significantly correlated with patient survival. The transcription factor, EHF, and its target, FGD6, had significantly different levels of expression between patients with high and low neutrophil counts. GO:0043547 exerts positive regulation of GTPase activity and GO:0007010 facilitates cytoskeleton organization; FGD6 participates in both functions and there were significant differences in patients with high and low neutrophil counts. Correlation analysis showed that FGD6 had the greatest correlation with RIC8B and SIPA1L3 in the positive regulation of GTPase activity and cytoskeleton organization. SIPA1L3 is a member of the hsa04015: Rap1 signaling pathway. EHF is a member of the highly diverse ETS superfamily. Accumulation evidence over recent years has demonstrated that EHF has a great effect on tumorigenesis and the tumor microarray environment. Studies have revealed that EHF inhibits tumor invasion, metastasis, and the mesenchymal phenotype (23)(24)(25). Liu et al. (11) reported that deficiency of EHF induces T reg cell and myeloid derived suppressor cell (MDSC) accumulation and causes a decrease in CD8 + T cell infiltration. Wang et al. (9) found that EHF promotes colorectal carcinoma progression via activation pf TGF-β1 transcription and canonical TGF-β signaling. Liu et al. (11) revealed that overexpression of EHF can significantly lead to better response to anti-PD-1 therapy. High expression of EHF promotes oral squamous cell carcinoma metastasis, cancer stemness, and drug resistance (26). Taken together, EHF has emerged as a promising target in cancer research; however, the contributions of EHF to neutrophils in HCC are not fully understood. We found that the level of EHF expression was different between patients with high and low neutrophil counts. Our finding suggested that EHF plays an To understand the mechanism underlying EHF, we screened its downstream target, FGD6, based on the outcome of ChIP-seq. FGD6 is a member of the faciogenital dysplasia family. FGD6 includes six Rho guanine nucleotide exchange factors (GEFs) that activate Rho proteins (12). FGD6, which is highly homologous to FGD1, has two pleckstrin homology (PH) domains and one Dbl homology (DH) domain. Previous studies have shown that mutations in the FGD1 gene are linked to the inherited disease, faciogenital dysplasia, and Aarskog-Scott syndrome (27,28). Recent evidence demonstrated that FGD6 plays a crucial role in cell polarity and membrane recycling by regulating the assembly of different actin-based protein networks and activating CDC42 at different locations (12). Our data showed that FGD6 is involved in GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. In addition, the level of FGD6 expression was highly correlated with the level of RIC8B and SIPA1L3 expression in GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. This finding indicated that FGD6, RIC8B, and SIPA1L3 may work together in the positive regulation of GTPase activity and cytoskeleton organization.
SIPA1L3, which has N-terminal RapGAP (a Rap-GTPase activating domain), PDZ, and C-terminal leucine zipper domains, is a member of the signal-induced proliferation-associated (SIPA) proteins. A previous study revealed that loss of SIPA1L3 interferes with cell polarity and cytoskeletal organization (29). Previous evidence showed that RIC8B is GEFs for the α subunits of heterotrimeric guanine nucleotide-binding proteins (30). In addition, it has been suggested that RIC8B contributes to the steady state of G protein α subunits (31). Some studies have showed that overexpression of RIC8B enhances activation of adenylyl cyclase (AC) via Gs and Golf (32,33). We discovered that the level of FGD6 expression had a high correlation with the level of RIC8B and SIPA1L3 expression. Therefore, it is essential to explore the relationship and potential mechanisms underlying FGD6, RIC8B, and SIPA1L3.
Our data showed that SIPA1L3 participates in the has04015: Rap1 signaling pathway. Accumulating evidence demonstrates that the Rap1 signaling pathway contributes a major role in the endothelial barrier regulation, vasculogenesis, and angiogenesis (34,35). A recent study involving the Rap1 signaling pathway has focused on integrin-mediated cell adhesion (36). Based on the diagram of the has04015: Rap1 signaling pathway on the KEGG website (https://www.genome.jp/kegg-bin/ show_pathway?map=hsa04015&show_description=show) (Figure 9), we found that the G protein-coupled receptors (GPCRs) signaling pathway activates the Rap1 signaling pathway via EPAC. EPAC, which includes GEFs, can exchange bound GDP of Rap1 for free GTP. In contrast, SIPA1L3, which contains the RapGAP domain, can activate the Rap-GTPase inhibiting Rap1 signaling pathway. Active Rap1 can activate CDC42 via Src and FRG. CDC42 will activate the process of cell adhesion, migration, and polarity. Collectively, these processes can be characterized in three sections: (I) the process of GPCRs regulating Rap1; (II) the process of regulation between Rap1-GDP and Rap1-GTP; We showed that FGD6, RIC8B, and SIPA1L3 are crucial regulatory factors in these three sections. RIC8B enhances activation of AC via Gs. This process will trigger EPAC to promote Rap1-GDP conversion into Rap1-GTP, which will activate the Rap1 signaling pathway. SIPA1L3, which has an N-terminal RapGAP domain, promotes the transformation of Rap1-GTP to Rap1-GDP. This process inhibits the Rap1 signaling pathway. FGD6, which is regulated by the transcription factor, EHF, can activate CDC42. This process will positively regulate the Rap1 signaling pathway and contribute to cell adhesion, migration, and polarity.
We inferred that the transcription factor, EHF, mediates GO:0043547 positive regulation of GTPase activity, GO:0007010 cytoskeleton organization, and the has04015: Rap1 signaling pathway by regulating the transcription of FGD6. These processes will mediate local adhesion. Increasing local adhesion may contribute to the recruitment and persistence of neutrophils in tumor lesions. Also, tumor cells might be protected and transported by the surrounding neutrophils following cell-cell interactions (3). As the key step in tumor metastasis, increasing adhesion will raise the ability of metastasis and colonization of tumor cells. This feature can contribute to the attachment of tumor cells to matrix collagen fibers of vascular endothelial cells and secondary lesions (37). These reasons suggest that patients with a high neutrophil count had worse survival. Increased EHF upregulates FGD6. Increased FGD6 and RIC8B or decreased SIPA1L3 will contribute to activation of the Rap1 signaling pathway. Activation of the Rap1 signaling pathway will increase cell adhesion, which is advantageous to neutrophil recruitment and metastasis, and colonization of tumor cells.

Conclusions
These findings demonstrated that transcription factor EHF can influence recruitment of neutrophils by mediating the transcription of FGD6. Further investigations are needed to shed new light on EHF and its target FGD6 which may indicated with dashed gray lines. We carried out log-rank statistics to identify the optimal cut-off for transforming the continuous variable of gene expression into categorical high-and low-expression groups in a survfit model. The test score at each candidate cut-off across the logtransformed gene expression values was plotted. The highest test score (indicated with a blue arrow) was applied for best separating patients into 4 different risk groups (using solid blue lines; named groups I to IV).

Figure 9
The diagram of has04015: Rap1 signaling pathway. The pertinent proteins have been marked by red color on the plot. RIC8B can enhance activation of adenylyl cyclase (AC) via Gs. This process will trigger EPAC to promote Rap1-GDP converting into Rap1-GTP.
SIPA1L3, which has an N-terminal Rap-GAP domain (Rap-GTPase activating domain), can promote the transforming of Rap1-GTP to Rap1-GDP. This process will inhibit Rap1 signaling pathway. FGD6 which includes six Rho GEFs activating Rho proteins can activate CDC42.
contribute to immunotherapy in HCC.
the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.