1 INTRODUCTION

The incidence of liver cancer has increased rapidly worldwide in recent years. Study shows that its morbidity has elevated to 6th in cancer, while it also ranks the third death cause of cancer.1 Hepatocellular carcinoma (HCC) represents approximately 90% of primary liver cancer.2 In the tides of immune therapy, patients with many other kinds of cancers are benefited from the novel therapy.3 However, the prognosis of patients with HCC is not optimistic.4-6 By deeper understanding the role of the immune system in hepatocellular carcinoma, we may discover novel prognostic biomarkers and advanced treatment.

The complement system is a fundamental branch of innate immunity of human body, containing more than 50 inherent components and membrane binding receptors and regulators. Three main pathways of complement activation depending on C37 have been described: the classical pathway (CP), the alternative pathway (AP), and the lectin pathway. Other factors, such as proteases, kallikrein, plasmin, and factor XIIa, can generate complement activation products without C3.8 Thrombin, a member of the coagulation pathway, can locally generate C5a in C3-deficient mice which are unable to generate the conventional C5 convertase.9

It is reported that the complement system has a complicated relationship with human cancer.10 On one hand, the complement activation may act as antitumor defense mechanism because of its participation in immune surveillance and complement-dependent cytotoxicity. Rituximab, an anti-CD20 antibody against malignant B cells, acts its antitumor effect by activating the classical complement pathway in vitro and in vivo studies.11, 12 On the other hand, it is now widely believed that the activation of the complement system may promote tumorigenesis and cancer progression by causing inflammation in the tumor microenvironment.13 In this situation, complement activation may not strengthen immune surveillance but, instead, cause immunosuppression by upregulating cytokines such as IL-10, PDL-1, CD46 in tumor cells as reported.14 Also, complement receptors, C3aR and C5aR, are considered to promote cancer cell proliferation in cancers, including, ovarian cancer,15 colon cancer,16 lung cancer,17-19 and so on. The effects of complement on HCC have also been explored but remain unclear. It is reported that Aristolochic acid I, which is considered to be hepatotoxic, promotes the invasion and migration of HCC cells by activating the C3a/C3aR complement system.20 Chen et al.21 reported downregulation of C3aR/C5aR inhibits cell proliferation and epithelial–mesenchymal transition (EMT) in HCC. Complement factor H, a negative regulator of the alternative pathway of the complement system, has been reported that its deficiency in mice caused spontaneous hepatic tumors.22 Recent studies show that enhanced CFH or CFHR (CFHR1 and CFHR3) levels either by gene therapies or by CFH reconstitution might lower tumor burden in HCC.23 However, previous studies have focused on particular complement members' effect on HCC. The role of other complement members remained unknown in the development and progression of HCC.

To address this, we thoroughly explored the expressions of the complement-related genes and their correlation with prognosis and immune infiltration in patients with HCC via The Cancer Genome Atlas (TCGA) public and GEO database for the first time.

2 MATERIALS AND METHODS 2.1 Study patients

TCGA is a comprehensive database containing information on more than 30 kinds of human cancers, including mRNA sequencing, DNA sequencing, pathological data, and clinical pathological information.24 In this study, we downloaded clinic pathological information and mRNA expressions of 364 HCC patients from TCGA. Table 1 showed clinical data of all the 364 patients, including gender, age, weight, AFP, Child-Pugh stage, adjacent tissue inflammation, cirrhosis, histologic grade, pathologic stage, and so on.

Table 1. Basic characteristics of 364 HCC patients Variables HCC patients (N = 364) Gender (Male/female) 246/118 Age (years, mean ± SD) 59.67 ± 13.37 Weight (kg, median) 69 (40–172) PLT (10e9/L, median) 211 (4–499000) Albumin (g/L, median) 4 (0.2–5200) Creatinine (mg/dl, median) 0.9 (0.4–12.4) PT (s, median) 1.1 (0.8–36.4) TB (μmol/L, median) 1.2 (0.2–2.1) AFP (ng/ml, median) 15 (1–2035400) Child-Pugh stage A N = 126 B N = 21 C N = 1 Adjacent tissue inflammation Non N = 117 Mild N = 97 Severe N = 17 Cirrhosis Non-cirrhosis N = 74 Cirrhosis N = 134 Histologic grade 1 N = 55 2 N = 174 3 N = 118 4 N = 12 Pathologic stage 1 N = 170 2 N = 83 3 N = 83 4 N = 4 Abbreviations: AFP, alpha fetoprotein; HCC, hepatocellular carcinoma; PLT, platelet; PT, prothrombin time; SD, standard deviation; TB, total bilirubin.

As there are only 50 normal liver samples in TCGA database, we also collected HCC data set GSE25097 from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). This data set contains 268 HCC samples and 243 normal controlled samples from HCC patients.25 The data set was estimated thoroughly through the full text.

2.2 Inclusion and exclusion criteria

The inclusion criterion was as follows: (1) patients in the TCGA and GSE25097 databases were diagnosed with HCC. (2) There were clear criteria for the diagnosis and staging of HCC. (3) The study data provided OR (odds ratio)/HR (hazard ratio) and its 95% confidence interval (CI), or could be converted to OR and its 95% CI. The exclusion criterion was as follows: (1) Repeat recorded cases. (2) Data was incomplete.

2.3 Differential expression analysis

First, differences in mRNA expressions of the whole complement system between cancer samples and normal controlled samples were analyzed by TIMER2.TIMER2 (http://timer.cistrome.org/) is an online analysis tool that can not only analyze gene expression profiles in different kinds of cancers, but also provide an advanced algorithm to evaluate the abundance of tumor-infiltrating immune cells by gene expression profiles based on TCGA database.26

Then, we used the pair-controlled samples for HCC samples in GEPIA database and datasets GSE25097 to further verify the different mRNA expressions between HCC samples and their normal counterparts. Gene Expression Profiling Interactive Analysis (GEPIA, www.gepia.cancer-pku.cn) is a web server analyzing RNA sequencing expression data of tumors and normal samples from both TCGA and Genotypic Tissue Expression (GTEx) projects using standard processing pipelines.27

Next, UALCAN database was used to examine the relationship between mRNA expressions of the complement system and patients' clinicopathologic information, such as tumor stage and grade. UALCAN (http://ualcan.path.uab.edu), a web tool based on level 3 RNA-seq and clinical data of 31 cancer types from TCGA database, can be used to explore the association of the transcriptional expression of a specific gene with cancer patients' clinicopathologic parameters.28

2.4 Survival analysis

In this study, we evaluated the prognostic value of mRNA expressions of the complement system in patients with HCC by calculating overall survival (OS), disease-free survival (DFS), and progress-free survival (PFS) using the Kaplan–Meier plotter (www.kmplot.com), which is a web tool for survival analysis by product-limit method based on gene expression data and survival information of patients with different kinds of cancers based on TCGA database.29, 30 Information can be found on the Kaplan–Meier plotter web page, including number of patients, best cutoff values of mRNA expression, 95% CI, HR, and p-value.

2.5 Mutation analysis

We used cBioPortal database to investigate the mutation rate of the whole complement system in HCC patients and to confirm its relation with the survival of HCC patients. cBioPortal, a database that includes data for putative copy-number alterations (CNAs), mRNA expression z-scores, and mutations, allows researchers to explore multidimensional cancer genomics data.31

2.6 PPI and enrichment analysis

We performed a protein–protein interaction (PPI) network analysis on the prognostic genes of the complement system to explore their interactions by STRINGS (www.string-db.org), an online analysis tool that collects, scores and integrates all publicly available sources to predict PPI between various proteins.32

LinkedOmics was used to evaluate both positively and negatively co-expressed genes of the prognostic complement genes of the complement system in HCC patients. LinkedOmics database predicted co-expressed genes of a specific gene by Spearman's correlation coefficient test using multi-omics data for 32 cancer types and a total of 11,158 patients from TCGA project.33

Next, we used the list of co-expressed genes to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis to predict the function of these prognostic genes in HCC by Metascape (http://metascape.org), a free and well-maintained online bioinformatics database for GO and KEGG enrichment analysis.34 In this study, we considered as significant only when p < .01, minimum count > 3, and enrichment factor > 1.5.

2.7 Immune infiltration analysis

We first used TIMER2 to detect the relationship between the mRNA expression of the whole complement system and tumor immune cells infiltration, including CD8+ and CD4+ T cells, macrophages, neutrophils, and dendritic cells (DCs), by spearsmans' correlation coefficients. Then, we examined the relationships between the prognostic complement genes and immune gene marker sets in HCC sample from TCGA using TIMER2. The gene markers of tumor-infiltrating immune cells, including T cells, monocytes, tumor-associated macrophages (TAMs), M1 and M2 macrophages, neutrophils, DCs, Tregs, and T cell exhaustion, were referenced from prior studies.24, 29, 30

2.8 Construction of the prognostic model

A prognostic model was constructed based on our results of the complement expressions and OS prognostic information in the cluster. The method of the least absolute shrinkage and selection operator (LASSO) regression algorithm is applied to find a model that can best explain the data but contains the fewest free parameters.

First, timeROC(v 0.4) analysis was performed to compare the predictive accuracy of each gene and score. Then, LASSO regression algorithm for feature selection, using tenfold cross-validation. The K–M survival analysis with log-rank test was also used to test the constructed model. All analytical methods above were performed using R software version 4.0.3.

2.9 Statistical analysis

The differential expression of complements was explored using the Wilcoxon test in TIMER2, ANOVA in GEPIA, and Students' t-test in UALCAN and GSE25097. A fold change >2 or <0.5, or log2(FC) > 1 or <−1 with p-values < .05 in both Gepia database and GSE25097 were considered statistically significant in mRNA expression in this study. The K–M curve diagrams of patients' survival were displayed as HR and log-rank p-values. In Spearman's correlation coefficients, we defined the strength of correlation as follows: 0.00–0.1, negligible correlation; 0.10–0.39, weak correlation; 0.40–0.69, moderate correlation; 0.70–0.89, strong correlation; and 0.90–1.0, very strong correlation. In this study, p-values < .05 were considered statistically significant.

3 RESULTS 3.1 The differential expression of complement in hepatocellular carcinoma

First, Table S1 showed the complement members evaluated in this study and their classification. Then, mRNA expressions of complement genes were explored via TIMER2, which gave us a general idea of the differences in complement gene expressions between HCC (n = 371) and normal tissues (n = 50) in patients (Table S2). Our results revealed that 27 complement genes, including C1QA, C1QB, C1QC, C1R, C1S, C3, C5, C6, C7, C8A, C8B, C9, CFB, CFD, CFI, CFP, CFHR3, CFHR4, C5AR1, SERPING1, MBL2, MASP1, C4BPA, CLU, CPN1, VSIG4, and CR1, were significantly downregulated in HCC tissues than in normal tissues (p < .05). In contrast, nine complement genes, including C2, C5AR2, C4BPB, CD46, CD59, ITGAM, ITGAX, C1QBP, and C1QR, were significantly upregulated in HCC tissues compared with normal samples (p < .05, Figure 1).

Table 2. The significant changes of complement expression in transcription level between HCC and normal liver tissues Gene Fold change (HCC vs. Normal) p C1QA 0.52 ** C1QB 0.80 ** C1QC 0.55 ** C1R 0.33 ** C1S 0.80 ** C2 1.07 ** C3 0.94 ** C5 0.82 ** C6 0.35 ** C7 0.15 ** C8A 0.61 ** C8B 0.62 ** C9 0.39 ** CFB 0.92 ** CFD 0.62 ** CFI 0.57 ** CFP 0.10 ** CFH 0.86 ** CFHR3 0.45 ** CFHR4 1.11 0.07 C5AR1 0.71 ** SERPING1 0.80 ** C5AR2 1.09 * MASP1 0.28 ** MBL2 0.52 ** C4BPA 0.78 ** C4BPB 0.90 ** CLU 0.89 0.11 CPN1 0.61 ** VSIG4 0.78 ** CR1 0.51 ** ITGAM 1.10 0.22 ITGAX 1.16 0.06 CD93 1.09 0.09 CD46 1.24 ** CD59 0.98 0.44 C1QBP 0.80 ** Abbreviation: HCC, hepatocellular carcinoma. * p < .05 ** p < .01.

Transcriptional expression levels of complements in hepatocellular carcinoma in the TIMER2 database. ***p < .001, **p < .01, *p < .05

To further confirm the differential gene expressions between HCC samples and normal tissues, we used Gepia database and an independent data set (GSE25097), since they contained more paired normal liver samples. According to the result from Gepia database (Figure 2), the mRNA expressions of C1R, C6, C7, C8A, C9, CFP, CFHR3, and MBL2 were still significantly lower in HCC samples (n = 368) than in normal samples (n = 160), while the mRNA expression of C1QA was significantly higher in HCC samples (n = 368) than in normal samples (n = 160). The expression of the other 27 complements showed no difference between HCC samples and normal tissues, with a log2(FC) > 1 or <−1 and p-value < .05. We then used data set (GSE25097) to confirm the mRNA expression of all the complement genes. Our result revealed that the mRNA expressions of C1R, C6, C7, C9, CFP, CFHR3, and MASP1 were significantly downregulated in HCC tissues (n = 268) compared with normal liver tissues (n = 243), with a fold change >2 or <0.5 and p-value < .05 (Table 2), while the expression of the rest complement genes showed no difference between HCC samples and normal tissues.

The mRNA expressions of complements in HCC tissues (GEPIA). HCC, hepatocellular carcinoma. ***p < .001, **p < .01, *p < .05

To sum up, C1R, C6, C7, C9, CFP, and CFHR3 were significantly downregulated in all three HCC data sets. C8A, MASP1, and MBL2 were found downregulated in two HCC data sets. C1QA was found downregulated in TIMER2 database, but was found upregulated in Gepia database. In GSE25097, mRNA expression of C1QA was not significantly altered in HCC samples. So, we considered that it need further validation. Other complement genes, including C1QB, C1QC, C1S, C2, C3, C5, CFB, CFD, CFI, CFHR4, C5AR1, SERPING1, C4BPA, C4BPB, CLU, CPN1, VSIG4, CR1, ITGAM, ITGAX, CD46, CD59, and C1QBP transcription levels were significantly down- or upregulated in HCC samples according to TIMER2, while no difference was found in their transcription level according to Gepia and GSE25097 database. So, we considered the mRNA expression of these genes in HCC tissues and normal liver tissues was not clear and definite.

In general, our result showed that the mRNA expression of six complement genes, including C1R, C6, C7, C9, CFP, and CFHR3, were significantly lower in HCC tissues than their normal counterparts.

3.2 Prognosis analysis of complements in HCC

We analyzed the correlation between transcription levels of complement genes and patients' survival by Kaplan–Meier plotter database as survival analysis might indicate novel potential pathogenesis of HCC and promote further research. The results showed that higher mRNA expressions of 32 genes of the complements, namely C1R, C1S, C2, C3, C5, C6, C7, C8A, C8B, C8G, CFB, CFI, CFP, CFH, CFHR1, CFHR2, CFHR3, CFHR4, CFHR5, SERPING1, MASP2, MASP1, MBL2, C4BPA, C4BPB, VTN, CLU, CD59, CPN1, and CR2 were associated with longer OS in HCC patients while the higher mRNA expression of C5AR2 and ITGAM was associated with shorter OS in HCC patients (Table 3). To further evaluate the prognosis value of complement genes in HCC, we first explored the prognosis value of C1R, C6, C7, CFP, and CFHR3 via DFS and PFS, since their mRNA expressions were significantly downregulated in HCC tissues. We found that higher expressions of C1R, C6, C7, CFP, and CFHR3 were also associated with better DFS and PFS in HCC patients (Figure 3). Therefore, it was conceivable that higher mRNA expressions of C1R, C6, C7, CFP, and CFHR3 were independent biomarkers that predict a better prognosis in patients with HCC.

Table 3. Correlation of complement gene expression and overall survival in HCC patients via Kaplan–Meier plotter Complement gene Medium survival months Hazard ratio p Low expression cohort (months) High expression cohort (months) C1QA 59.7 56.2 0.81 0.24 C1QB 46.6 70.5 0.77 0.15 C1QC 46.6 70.5 0.76 0.14 C1R 30 61.7 0.57 ** C1S 31 70.5 0.55 ** C2 37.8 70.5 0.56 ** C3 33.5 71 0.48 ** C4A 56.2 84.4 0.67 0.06 C5 33.5 70.5 0.49 ** C6 31 61.7 0.54 ** C7 28.3 70.5 0.52 ** C8A 45.7 70.5 0.59 ** C8B 45.7 84.4 0.47 ** C8G 27.9 70.5 0.50 ** C9 56.2 70.5 0.75 0.15 CFB 38.3 61.7 0.57 ** CFD 52 56.5 0.74 0.10 CFI 31 61.7 0.64