{'fraction_non_alphanumeric': 0.031512183372807194, 'fraction_numerical': 0.01830516639235679, 'mean_word_length': 5.12690605017216, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development and Validation of a Novel 8 Immune Gene Prognostic Signature Based on the Immune Expression Profile for Hepatocellular Carcinoma

Background The immune microenvironment plays a vital role in the development of hepatocellular carcinoma (HCC). This study explored novel immune-related biomarkers to predict the prognosis of patients with HCC. Methods RNA-Seq data were downloaded from The Cancer Genome Atlas (TCGA). Univariate Cox regression was used to identify prognosis-related genes; the Lasso method was used to construct the prognosis risk model. Validation was performed on the International Cancer Genome Consortium (ICGC) cohort, and the C-index was calculated to evaluate its overall predictive performance. Western blots were conducted to evaluate the expression of genes. Results There were 320 immune-related genes, 40 of which were significantly related to prognosis. Eight immune gene signatures (CKLF, IL12A, CCL20, PRELID1, GLMN, ACVR2A, CD7, and FYN) were established by Lasso Cox regression analysis. This immune signature performed well in different cohorts and can be an independent risk factor for prognosis. In addition, the overall predictive performance of this model was higher than the other models reported previously. Conclusion The predictive immune model will enable patients with HCC to be more accurately managed in immunotherapy.


Introduction
Hepatocellular carcinoma (HCC) is the second leading cause of cancer-related mortality and the fifth most common malignancy worldwide. 1 In 2012, it was estimated that approximately 700,000 people globally die of hepatocellular carcinoma each year. 2 The incidence of HCC varies by geographic region, gender, age, and risk factors associated with cancer development. 3 Despite advances in treatments such as liver transplantation, surgical resection, and radiofrequency ablation, the prognosis in patients with HCC is still very poor, with a 5-year survival rate of less than 30%. 4 In addition, tumor heterogeneity and immune microenvironment play key roles in the progression, treatment, and prognosis of HCC. Therefore, an in-depth exploration of immune phenotypic changes and their impact on prognosis is necessary to improve patient survival.
HCC is mainly caused by chronic inflammation of the liver and is considered a typical immunogenic cancer. 5 The immune response plays an important role in the development of HCC. 6 There are many antigen-presenting cells in the liver. When presenting antigens to T cells in the liver, they usually express immunosuppressive surface molecules or secrete immunosuppressive cytokines, thereby mediating liver immune tolerance. 7 These antigen-presenting cells play an important role in HCC and are involved in the immune escape process of HCC cells. 8 In addition, in the microenvironment of HCC, infiltrating immune cells are impaired by the expression of inhibitory signal molecules such as PD-1 on the surface, leading to tumor immunosuppression. These inhibitory signal molecules are immune checkpoints. 9 In recent years, immune checkpoint inhibitors have become an effective means of treating tumors. 10,11 Therefore, it is of vital importance to further explore the clinical relevance of HCC immunophenotypes and their relationship with prognosis.
Several previous studies have explored the role of tumor microenvironment in HCC. One investigation detected the expression of 49 immune genes in 61 HCC tissues by immunohistochemistry and PCR, indicating that these genes were correlated with the prognosis of HCC. 12 This research clarified the potential effect of the immune microenvironment on prognosis and initially confirmed that the expression level of immune system genes had great prognostic value in HCC. In addition, tumor-infiltrating lymphocytes (TILs), as an important part of the immune microenvironment, have been widely confirmed to be correlated with the prognosis of HCC. 13,14 For example, the low expression of IL-36αwas associated with to the poor prognosis of HCC, the specific mechanism of which was related to T cells. 15 Several studies have also confirmed the role of immune checkpoint-related genes in predicting prognosis. For example, studies have demonstrated that overexpression of PD-L1 and PD-L2 was associated with poor prognosis. 16,17 However, as this research has thus far been limited, there is still no systematic understanding of the immunophenotype in the HCC microenvironment and its relationship to prognosis.
In this study, we integrated multiple gene expression cohorts to develop an HCC prognostic signature based on immune-related genes. The stability and reliability of the model are proven through the external independent cohort and laboratory experiments.

Inclusion Criteria
RNA-seq expression profiles were retrospectively collected from two public cohorts: The Cancer Genome Atlas hepatocellular carcinoma cohort (TCGA-LIHC), and the International Cancer Genome Consortium cohort (ICGC-LIRI-JP). All patients had undergone primary tumor removal. Only those patients with available followup time and gene expression profiles were included. The primary outcome was overall survival, and staging was evaluated according to the AJCC system.
For the TCGA-LIHC cohort, the GDC (https://portal. gdc.cancer.gov/) API was used to download the RNA-seq expression profile and clinical follow-up information. The cohort contained 423 samples, which were further converted into gene symbols by the hgu133plus2.db. R package. We downloaded the ICGC-LIRI-JP cohort as the external validation cohort, which has 212 samples.

Processing of Gene Expression Profiles
All data were mapped to Entrez id. If multiple probes corresponded to the same Entrez ID, the probe with the highest average signal was selected as the expression level of its corresponding gene. The details of sample quality control are described in Table 1. Finally, 583 cases of HCC samples including TCGA-LIHC and ICGC-LIRI-JP cohorts were included in this study.

Identification of Immune-Related Genes
The immune-related genes were downloaded from the Molecular Signatures Database v4.0 (immune system process, M13664; immune response, M19817; immune effector process, M14818; and immune system development, M3457). Finally, 321 immune-related genes were included in the subsequent analysis.

Establishment of Immune Prognostic Signature of HCC (IPSHCC)
First, the TCGA training cohort was used to identify genes related to survival. There were a total of 321 genes in the immune-related gene set, and a total of 320 genes in this gene set were detected in the TCGA-LIHC training cohort.
A Cox proportional hazards model was used to exclude genes with p-values > 0.05. Because there might be interactions between genes, the Lasso method was used for further screening of genes. Lasso was more suitable for high-dimensional, strong correlation, and small sample survival data. The basic goal of using Lasso was to minimize the sum of squared residuals under the constraint that the sum of the absolute values of the regression coefficients was less than a constant, so that some regression coefficients were strictly equal to 0, thereby obtaining an interpretable model. minβ∑i=1n(yi−∑j=1pβjxij)2, subject to∑j=1p||βj||≤λ minβ∑i=1n(yi−∑j=1pβjxij)2, subject to ∑j=1p|βj|≤λ As l increased, the terms ∑j = 1p||βi||∑j = 1p | βi | would decrease, and the coefficients of some independent variables would be gradually compressed to 0, to achieve dimensionality reduction of high-dimensional data. The dimension reduction of the Lasso method was achieved by penalizing the number of regression coefficients.
Cross-validation, such as k-fold cross-validation, was a standard method to estimate the adjustment parameter l. In the k-fold cross-validation, all data observations were roughly divided into k equal parts, and the data were fitted in turn with all possible k-1 parts as the training cohort, and the remaining one was used as the testing cohort. The ratio of the number of observations of the testing cohort to the training cohort was approximately 1: (k-1). A total of k calculations were performed to obtain k indicators (such as the error rate or other indicators) when fitting the testing cohort and then averaged. This performance was repeated for each model; the one with the smallest error rate was finally selected.
The R package glmnet was used to perform Lasso regression using 10-fold cross-validation. The coefficients of each category were further determined by multivariate Cox regression using the reduced-dimensional genes: β i E i β i was the regression coefficient of the i-th gene, and E i was the expression level of the i-th gene.

Validation of IPSHCC
We used the median risk score to divide patients into highrisk and low-risk groups. After the IPSHCC prognostic risk score was determined in the training cohort, it was calculated in the validation cohort.

Pathway Enrichment Analysis
To further understand the gene function in IPSHCC, single-sample GSEA analysis was performed using the GSVA. R package, based on the KEGG database, and the enrichment scores of each sample in each pathway were calculated. Pathways significantly related to IPSHCC were further screened based on the correlation between the enrichment score and IPSHCC.

Comparison with Previous Prognostic Signatures
To assess the survival classification and predictive performance of IPSHCC, five published prognostic signatures were retrospectively collected for comparison, ranging from 3 to 9 genes (Table 2). [18][19][20][21][22] A continuous prognostic risk score was calculated for each signature. The p-value and the overall concordance index (C-index) of the univariate Cox model were compared among the other five signatures.

Validation of Expression of the Immune Genes
The Human Protein Atlas (HPA) provides information on the tissue and cell distribution of 26,000 human proteins. It primarily uses specific antibodies to study protein expression in cell lines, normal tissues, and tumor tissues. We explored the expression of 8 genes (ACVR2A, CCL20, CD7, CKLF, FYN, GLMN, IL12A, PRELID1) in normal and tumor tissues. We explored the expression of 8 genes in HCC and normal tissues in GSE6764 and GSE14520, and draw boxplots on the expression of genes.
Genetic Alterations of the 8 Immune Genes

Sample Collection
HCC and adjacent tissues were collected from 3 patients, immediately placed in liquid nitrogen, and preserved at −80°C. Patients and their families were fully informed, and informed consent was obtained from the participants. The Ethics Committee of Department of hepatobiliary and pancreatic surgery, Hainan General Hospital, approved this study.

Western Blotting
Western blotting was carried out according to standard protocols. We used primary antibodies raised against GAPDH (Santa-Cruz Biotechnology, CA, USA),

Statistical Analysis
The Cox proportional hazard model was used to assess the relationship between clinical features and overall survival. Kaplan-Meier (KM) survival curves were plotted, and Log rank tests were used to detect differences among subgroups. The survival and survRM2. R packages were used to estimate the C-index.

Immune Pathways Were Significantly Suppressed in HCC
We performed a single-sample Gene Set Enrichment Analysis (ssGSEA) analysis in TCGA-LIHC cohort and calculated the enrichment scores of each sample in the four immune pathways ( Figure 1A). We observed that the scores in tumor samples were significantly lower ( Figure 1B), indicating that immunity was significantly suppressed.

Development and Definition of 8-IPSHCC Signature
The TCGA cohort was randomized into the training cohort (N = 186) and the internal testing cohort (N = 185). In the training cohort, a univariable analysis was performed on each immune gene to identify a total of 40 prognostic genes. Lasso regression was used for dimensionality reduction analysis, and 10-fold cross-validation was selected. The minimum error rate when λ = 0.07388702 is shown in Figure 2A and B. The 8-IPSHCC signature was established using multivariable Cox regression analysis (Table 3).

Predictive Performance of the Signature in Different Cohorts
The model was used to predict risk scores of the training cohort, testing cohort, and external validation cohort, and samples were classified according to their median risk score. The average area under the curve (AUC) of the 8-IPSHCC signature in the training cohort was 0.76 ( Figure 3A). The prognoses of the highrisk and the low-risk groups were significantly different (Hazard ratio, HR = 2.72, 95% CI: 2.07-3.57, Figure 3B). The average AUC in the testing cohort was 0.65 ( Figure 3C). The prognosis of high-risk samples was significantly worse than low-risk samples (HR = 1.71, 95%  Figure 3D), and the average AUC in the external validation cohort was 0.73 ( Figure 3E). The overall survival of high-risk samples was significantly worse than the low-risk samples with p = 0.00016, which was consistent with the training cohort (HR = 73.75, 95% CI: 11.59-469.4, Figure 3F). In short, 8-IPSHCC signature provided effective prognostic classification in multiple cohorts.

Univariate and Multivariate Analysis of IPSHCC
To identify the independence of the IPSHCC in clinical applications, we used Cox regression to analyze the clinical information, including age, stage, gender, and risk score (Table 4). It showed that the risk score was independent prognostic risk factors in multivariate Cox regression analysis (HR = 2.383, 95% CI = 1.7849-3.1825, p = 3.94E-09). We plotted survival curves for the age and stage ( Figure 4A-D).

Pathway Enrichment Analysis of IPSHCC
The risk score of each sample was calculated based on the immune gene signature in the training cohort, and Kyoto

DovePress
Encyclopedia of Genes and Genome (KEGG) pathway analysis was performed using ssGSEA. The Pearson correlation between the immune gene signature and sample pathway enrichment score was calculated. A total of 50 pathways had false discovery rate (FDR) <0.01. Among these pathways, there was a total of 12 significant positive correlations and 38 significant negative correlations ( Figure 5A). There were pathways related to metabolism, cancer, the cell cycle,

8132
DNA repair, and other functions. Most of these pathways were significantly related to tumors: the six most relevant pathways included "FATTY ACID METABOLISM", "BETA-ALANINE METABOLISM", "SPLICEOSOME", "CALCIUM SIGNALING PATHWAY", "CELL CYCLE", and the "ADIPOCYTOKINE SIGNALING PATHWAY". These pathways had significant different enrichment scores between the high and low-risk groups ( Figure 5B).

Comparison of Prognostic and Clinical Characteristics with Other Signatures
To assess the survival classification and predictive performance of IPSHCC, the accuracy of the IPSHCC model in predicting HCC prognosis was compared with 5 published prognostic signatures of HCC (Table 2). Continuous prognostic scores were calculated based on each feature, and cohorts were compared. The C-indices of prognostic signatures in the training cohort and the external validation cohort are shown in Figure 6A and B. Among these, IPSHCC had the highest (0.725) average C-index in both cohorts, demonstrating that IPSHCC had a good predictive performance.

External Validation of the IPSHCC Expression
The protein expressions of the presumptive 8 predictive genes were analyzed using the HPA database. Among them, CCL20 and CKLF were not collected in HPA database, CD7 was negative in tumor and normal tissue. There were no significant differences in the expression of FYN and IL12A in tumor and normal tissue, and ACVR2A, GLMN, PRELID1 was relatively highly expressed in tumors (Figure 7). The differential expression of the 8 genes in GSE6764 and GSE14520 cohort was analyzed, among which expression of GLMN, CKLF, and CCL20 was high in tumors, whereas the expression of FYN, IL12A in tumors was low ( Figure 8).

Genetic Alterations of 8 Genes in IPSHCC
The mutations of the 8 genes were explored in the cBioportal database. Among them, the gene with the highest mutation rate was CD7, accounting for 4%, and the primary type of mutation was amplification. ACVR2A mutations accounted for 3%, and the mutation form is mainly deep deletion (Figure 9).

Clinical Validation of Protein and mRNA Levels of 8 Genes
We analyzed 3 pairs of HCC tissues and para-cancerous controls to validate the protein and mRNA levels of 8 genes. The results showed that the expressions of CD7, IL12A, and FYN in cancer and normal tissue was not significant. The expressions of GLMN, CKLF, CCL20, PRELID1, and ACVR2A was relatively highly expressed A B

DovePress
OncoTargets and Therapy 2020:13 DovePress in tumors. The trend of our experimental results was almost consistent with data analysis (Figure 10).

Discussion
HCC is a highly aggressive malignancy and is currently the second leading cause of tumor-related death. 1 Despite plenty of advances in comprehensive treatment strategies for hepatocellular carcinoma, effective molecular targeted therapies are still lacking. 23 The complex immune microenvironment promotes the development of HCC and mediates the poor prognosis of patients. 24,25 Therefore, it is critical to integrate immune-related genes and conduct comprehensive research on HCC, to improve prognoses. In this study, we conducted the first ssGSEA analysis of HCC samples. The results showed that immunity was significantly suppressed, suggesting the important role of the immune microenvironment in the progression of HCC. Furthermore, an 8-IPSHCC signature was constructed by Lasso and multivariable analysis, and it showed stable and consistent predictive performance in both the internal and external validation cohort. More importantly, it was an independent prognostic risk factor in HCC patients.
In IPSHCC signature, CKLF, IL12A, CCL20, PRELID1, GLMN, ACVR2A, and CD7, were risk factors, and FYN was protective. Previous studies have confirmed that CKLF is related to the prognosis of patients with colon cancer, and the specific mechanism was related to leukocyte infiltration. 26 In addition, the gene polymorphism of IL12A was correlated with the risk and prognosis of various tumors, including lung cancer, 27 HCC, 28 cervical cancer, 29 and gastric cancer. 30 As a chemokine, CCL20 plays an important role in the progression of tumors and is closely related to the prognosis of patients with solid tumors such as colon cancer, 31 lung cancer, 32 and pancreatic cancer. 33 PRELI was originally thought to be expressed on germinal center B lymphocytes. 34 Recent studies have confirmed that it is associated with the prognosis of patients with cervical squamous cell carcinoma, pancreatic adenocarcinoma, and acute myeloid leukemia. 35 ACVR2A has been shown to correlate with the prognosis of patients with colon cancer. 36,37 CD7 is a common antigen in lymphoid hematopoietic tissues, often expressed in T cells and NK cells, and used in the diagnosis of T-cell lymphoma and leukemia. 38,39 The primary type of mutation of gene CD7 was amplification, accounting for 4%. FYN is a non-receptor tyrosine kinase that belongs to the Src family of kinases. 40 By regulating cell growth, apoptosis, morphogenetic transformation, and motor capacity, it promotes the development and progress of tumors and is closely related to the prognosis of patients with various tumors. 41,42 GLMN is not associated with the prognosis of tumor patients, and it was first discovered in this study as a prognostic marker for HCC.
Finally, the results of GSEA enrichment analysis revealed potential regulatory pathways of IPSHCC. The liver is vital in lipid metabolism, 43 and the occurrence and development of liver cancer are often accompanied by complex metabolic abnormalities. 44,45 The 8-IPSHCC signature may promote the metabolic reprogramming of liver cancer cells by regulating the metabolic pathway, such as the "FATTY ACID METABOLISM," "BETA-ALANINE METABOLISM," and "ADIPOCYTOKINE SIGNALING PATHWAY", thereby exerting a cancer-promoting effect. The "SPLICEOSOME" pathway catalyzes the process of alternative splicing, which is a post-transcriptional modification mechanism to ensure high transcript and protein diversity. 46,47 We suggest that the IPSHCC may regulate fundamental processes in HCC cells by alternative splicing via the "SPLICEOSOME" pathway. Ca 2+ signaling mediates fundamental cellular processes such as gene transcription, motility, cell proliferation, and apoptosis. 48,49 The "CALCIUM SIGNALING PATHWAY" was enriched, suggesting that the IPSHCC may function via regulation of Ca 2+ ion transporters and channels. In addition, an infinitely limited growth and cell cycle disorder is the main characteristic of tumors. 50 Notably, the occurrence of the vast majority of HCC is related to abnormalities in the regulatory mechanism of the G1 phase of the cell cycle. 51 So, the presence of the "CELL CYCLE" pathway in the GSEA results suggests the 8 immune gene signature may be related to the proliferation ability of HCC cells. Therefore, the enrichment analysis pathway in this study requires further research to elucidate the mechanism of HCC occurrence and progression.
Moreover, the results of GSEA enrichment analysis also showed that the IPSHCC signature was significantly related to pathways such as cell cycle, and DNA repair. Many previous studies have confirmed the significant effects of metabolic pathways and cell cycle changes on the biological behavior of HCC cells and the prognosis of clinical patients, 52,53 which was consistent with our results. Therefore, the pathway of enrichment analysis in this study needed further research to deepen the understanding of the mechanism of pancreatic cancer occurrence and progression.
Many previous studies have attempted to screen and construct prognostic marker models for HCC, including one that utilized 6 genes based on the TCGA cohort, 18 and another the utilized 3 genes (UPB1, SOCS2, and RTN3) from the same cohort. 19 A prognostic signature based on the GEO and TCGA cohorts constructed for patients with HCC contained 4 genes (SPINK1, TXNRD1, LCAT, and PZP). 20 A genetic marker combining hypoxia and lymphocyte infiltration was used to construct an 8-gene prognostic signature for patients with multiple solid tumors, including HCC. 21 A 9-gene prognostic signature has also been constructed. 22 To confirm the advantages of our signature based on immunophenotypes, we analyzed the above 5 genetic models simultaneously. The results showed that the 8-IPSHCC signature in our study was better than these models in predicting the prognosis. Further C-index analysis confirmed our model's overall superior performance. Therefore, this model could help clinicians predict individual risk of patients and guide patient assessment and treatment decisions.
Finally, we used 3 pairs of HCC tissues and normal controls to verify the protein and mRNA levels of 8 genes. The results showed that the expressions of CD7, IL12A, and FYN in the tumor and control groups were relatively low. GLMN, CKLF, CCL20, PRELID1, and ACVR2A were highly expressed in tumor tissues. Our experimental results are generally consistent with the results from the HPA and GEO databases, which confirms our model's reliability.

Conclusions
In summary, we integrated multiple gene expression cohorts of 583 samples to develop and validate an 8-IPSHCC signature. It performed well in different cohorts and could improve prediction accuracy of overall survival in HCC patients.

Abbreviations
HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; IPSHCC, immune prognostic signature of HCC; AUC, area under the curve; FDR, false discovery rate; HR, hazard ratio; KEGG, Kyoto Encyclopedia of Genes and Genomes; KM, Kaplan Meier; OS, overall survival; ROC, receiver operating characteristic; ssGSEA, single-sample Gene Set Enrichment Analysis; C-index, concordance index; RMS, Restricted mean survival.

Data Sharing Statement
The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Publish your work in this journal
OncoTargets and Therapy is an international, peer-reviewed, open access journal focusing on the pathological basis of all cancers, potential targets for therapy and treatment protocols employed to improve the management of cancer patients. The journal also focuses on the impact of management programs and new therapeutic agents and protocols on patient perspectives such as quality of life, adherence and satisfaction. The manuscript management system is completely online and includes a very quick and fair peer-review system, which is all easy to use. Visit http://www.dovepress.com/ testimonials.php to read real quotes from published authors.