Immune Implication of ASF1B Gene in Hepatocellular Carcinoma

Background: Anti-silencing function 1B (ASF1B) has been demonstrated to contribute 16 to tumorigenesis. However, its carcinogenic and immune effects in hepatocellular 17 carcinoma (HCC) have not been reported. This study aimed to identify immune role of 18 ASF1B in HCC. 19 Methods: HCC datasets obtained from The Cancer Genome Atlas (TCGA) database 20 were used to investigate the role of ASF1B gene in HCC, followed by validation using 21 Gene Expression Omnibus (GEO) datasets and Gene Expression Profiling Interactive 22 Analysis (GEPIA) website. CIBERSORT analysis was performed to evaluate immune 23 cell infiltration levels. The TISIDB and cBioPortal network tool were used to seek 24 ASF1B-associated immunomodulators and its co-expressed genes. TCGA cohort was 25 divided into train set and test set according to the ratio of 7:3. Cox regression was used 26 to identify ASF1B-associated prognostic immunomodulators in train set, followed by 27 internal validation using the test set. Based on the median risk-score, HCC patients were 28 divided into high- and low-risk group for the further survival curves and receiver 29 operating characteristic (ROC) analysis, as well as nomogram and calibration curves 30 analysis. Finally, the dataset collected from the GEO was adopted for external 31 validation. 32 Results: ASF1B was over-expressed in TCGA HCC cohort and contributed poor 33 prognosis, which was verified in two GEO datasets (GSE14520 and GSE6764) and 34 GEPIA, as well as Kaplan Meier Plotter network tool. The immune cell infiltration 35 levels were found to be associated with the ASF1B copy numbers and mRNA 36 expression. A total of 78 ASF1B-associated genes were screened out, including 7 37 immunoinhibitors, 21 immunostimulators and 50 tightly co-expressed genes. Finally, 5 38 ASF1B-associated genes (TNFSF4, TNFRSF4, KDR, MICB and CST7) were 39 identified to be strongly related to HCC survival. study. Conclusion: Our findings showed that the ASF1B was associated with HCC immunity. 45 The selected ASF1B-asociated immune markers could be promising biomarkers for the 46 prognosis of HCC.


49
Liver cancer is the sixth most common tumor and the fourth most common cause of 50 cancer-related death worldwide (1). An increased risk of HCC is associated with the 51 fundamental chronic liver disease such as hepatitis B and C virus infection (2, 3). 52 Hepatocellular carcinoma (HCC) constitutes 75%-85% of primary liver cancers with 53 an unclear mechanism and poor prognosis (4, 5). Although immunotherapy has become 54 a promising alternative therapy for patients with HCC (6, 7), just a proportion of HCC 55 patients benefit from immunotherapy (8), which may result from the limited effective 56 target. Prognostic immune biomarkers could help to recognize immunotherapy-57 responsive subgroups. Some evidence has indicated that tumor infiltrating leukocytes 58 are related to the prognosis of cancers (9). Therefore, the molecular characteristics of 59 the immune micro-environment within HCC need to be further explored. It is essential 60 to fully understand HCC immunology and the molecular regulatory mechanisms to 61 ensure the success of immunotherapy.

62
ASF1B is one of the histone H3-H4 chaperone anti-silencing function 1 (ASF1) 63 paralogs (10, 11). Interestingly, studies have reported that high ASF1B expression is 64 linked to diagnosis and prognosis of breast cancer, renal cell carcinoma and cervical 65 cancer (12-14). However, the carcinogenic effect and the immune role of ASF1B in 66 HCC have not been reported yet. 67 In this study, RNA-Seq data from TCGA HCC cohort and GEO databases were (https://portal.gdc.cancer.gov/). Three liver cancer datasets GSE6764, GSE14520 and 76 GSE54236 were obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/). To further 77 process RNA expression data, the limma package was used for R software.

78
Relationship between ASF1B and clinical characteristics 79 The expression pattern of ASF1B was analyzed in the TCGA HCC cohort between 80 tumor and adjacent tissues. GEPIA was used to combine GTEx data and TCGA for 81 validation (http://gepia.cancer-pku.cn/). GSE6764 and GSE14520 were also used for 82 verification. According to ASF1B expression, the tumor samples were divided into high 83 and low expression groups. The survival curve of high and low groups was drawn by 84 the Kaplan Meier method. Kaplan Meier Plotter network tool was adopted to verify the 85 survival difference (https://kmplot.com/analysis/). Moreover, the correlation between 86 ASF1B expression and clinical features were analyzed by logistic regression and Cox 87 regression.

88
Between ASF1B and tumor mutation burden 89 Mutation data from TCGA LIHC was used to determine tumor mutation burden (TMB).

90
The relationship between TMB and ASF1B expression and its impact on the overall 91 survival (OS) were analyzed in HCC dataset.

93
The median expression of ASF1B gene was used as the cutoff value, and all tumor 94 samples in TCGA cohort were divided into two groups with high and low expression.

95
Gene set enrichment analysis was used to analyze the signal pathways that were 96 significantly related to the expression level of ASF1B. Instead of concentrating on only 97 a handful of mostly altered genes, GSEA evaluates the genome-wide expression 98 profiles at the levels of gene sets. A set of genes denotes a set of concordant genes with 99 a comparable biological function, chromosomal location, or regulation (15). were uploaded to cBioPortal (www.cbioportal.org). Based on RNA-seq data from 124 cancer samples, the queried genes were able to seize 50 simultaneously altered genes.

125
The obtained proteins were used to analyze protein interaction by STRING network tool (https://string-db.org/), and performed GO and KEGG analysis 127 (http://consensuspathdb.org/). 129 The TCGA cohort was separated into train set and test set with the ratio of 7:3.

130
Univariate and multivariate COX regression were used in train set to analyze ASF1B-131 related immune genes to screen out prognostic-related immune genes. Risk-score was 132 generated based on coefficient: risk-score = β1x1 + β2x2 +... + βixi. Xi was the expression 133 level of each gene, and βi was the coefficient of each gene obtained from the Cox model.

134
Kaplan Meier curve was used to assess the risk-score of immune-related gene with OS.

135
The time-dependent receiver operating characteristic (ROC) curves were adopted to 136 determine the prognostic accuracy of the risk-score. 138 In cancer prognosis, nomogram is widely used to estimate of the probability of a single 139 event, such as death or recurrence, tailored to a single patient's characteristics. In this 140 study, patients' parameters and risk-scores were integrated to form the nomogram to 141 analyze the prognosis. The nomogram was created by the rms package of R software.

142
To measure the predictive precision of the nomogram, the concordance index (C-index) 143 and calibration plot were used. Calibration curves were used with the application of the   was used to analyze the differential expression of ASF1B, univariate and multivariate 158 COX regression analysis were used to analyze prognosis-related genes. P < 0.05 was 159 considered to be statistically significant.

161
ASF1B was over-expressed in HCC and predicted poor prognosis 162 For the first time, we found that ASF1B expression in tumors was significantly higher 163 than that of adjacent tissues based on differential expression between 374 cases of HCC 164 tissues and 50 cases of adjacent tissues in TCGA RNA-Seq dataset (Figure 1a-b), 165 similar results can also be found in two GEO datasets (GSE6764 and GSE14520)  Table 1-2). 173 We also studied the LIHC mutation data from TCGA, which was divided into two 174 groups with high and low mutations according to the TMB score. The results showed 175 that the survival of patients in the low TMB group was better than that in the high group 176 (Figure 4a). Next, we compared the mutations status in the high and low ASF1B 177 expression groups and found that TMB was significantly associated with ASF1B 178 expression. High expression of ASF1B was associated with high mutation of some 179 genes, with the highest mutation in TP53 (Figure 4b-d). These results suggested a 180 possibility that ASF1B mediated mutation that might involve in the occurrence and 8 development of HCC.

183
To explore the up-regulated ASF1B expression associated pathways in HCC, GSEA 184 analysis was used to study the datasets of tumor samples collected from TCGA. We 185 found that elevated ASF1B participated in several paths, such as cell cycle, p53 186 signaling pathway, mTOR signaling, Fc gamma receptor-mediated phagocytosis, T cell 187 receptor and natural killer cell mediated cytotoxicity, that were involved in cell 188 proliferation and immunity ( Figure 5).

235
A nomogram was created to predict patient survival probability (Figure 14a). The curve 236 of calibration indicated a good fit between nomogram-predicted probability and idea 237 reference line for the 1-year, 3-year and 5-year survival (Figure 14b-d). Additionally, 238 the C-index was 0.69, which suggested a good predictive power.  (Figure 15b). Furthermore, the calibration curves showed a good consistency 253 between predicted survival and actual survival (Figure 15c-d).           Relationship between ASF1B and TP53 mutation. 447  The hazard ratios and coefcient of ASF1B-associated prognostic genes of multivariate COX regression. 464  Table 1        Associations between ASF1B copy numbers in HCC.

Figure 8
Association between ASF1B expression and immune cell in ltration levels via TISIDB network tool. The black asterisks in the correlation heatmap indicated immune cell types signi cantly associated with ASF1B expression levels in TCGA.

Figure 9
Identi cation and analysis of immunomodulators associated with the ASF1B gene.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.