We first collected gene expression data and clinical data from 424 patients with HCC in TCGA (https://portal.gdc.cancer.gov). In addition, the RNA sequencing data (GSE14520) were downloaded from the Gene Expression Omnibus (GEO) database. HTSeq-FPKM of level 3 format was converted into transcripts per million (TPM). Screening was performed to exclude patients with incomplete information, and the TPM data from 374 patients were used in subsequent analyses (Supplementary Table 1). The evolution process used the “ggplot2” R package.

Normal human liver cells (WRL68) were purchased from AcceGen (Fairfield, USA), and HCC cell lines (BEL7402, SK-hep1, HCCLM3, HepG2 and Huh7) were purchased from the Chinese Academy of Science (Shanghai, China). WRL68 cells were cultured in RPMI-1640 medium, and other cell lines were cultured in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin. All cells were incubated in a 37°C incubator with 5% CO₂.

We collected 30 pairs of HCC tissues and adjacent liver tissues at the First Affiliated Hospital of Harbin Medical University from 2006 to 2013 after obtaining informed consent from patients. The research project was conducted under the supervision of the Ethics Committee of the First Affiliated Hospital of Harbin Medical University.

Quantitative real-time PCR was performed on the samples as described previously (5). The following primers were used: CFHR4-F, 5’-TGCGGTTTAAGCTCCATGACA -3’; CFHR4-R, 5’-CCCATCTTCACCACACACTATG-3’; GAPDH-F, 5’ -TGACTTCAACAGCGACACCCA-3’ and GAPDH-R, 5’-CACCCTGTTGCTGTAGCCAAA-3’. GAPDH was used as a control to determine changes in mRNA levels using the 2^-ΔΔCT method.

The differentially expressed genes (DEGs) between high CFHR4 expression and low CFHR4 expression samples from TCGA database were analyzed using the DEseq2 (1.26.0) R package (22) with Student’s t test. Differences were considered statistically significant for an adjusted p value < 0.05 and absolute log2-fold change > 1.5. Moreover, volcano plots and heatmaps were constructed to visualize the DEGs.

Pathway enrichment analyses were performed with the “clusterProfiler” R package (23, 24). The c2.cp.v7.2.symbols.gmt curated gene sets were retrieved from the Molecular Signatures Database (MSigDB). Each analytical technique was conducted repeatedly a thousand times. An FDR-corrected q value < 0.25 and adjusted p value< 0.05 were considered statistically significant.

We analyzed the levels of infiltration of 24 types of immune cells in HCC using the ssGSEA method with the GSVA package in R. We then quantified the enrichment score for each immune cell by performing gene expression profiling of each HCC sample based on the signature of immune cells (25, 26).

The univariate Cox regression analysis of the correlation between CFHR4 expression and the values multiple clinical prognostic parameters in patients with HCC was performed using R software with the “survival” package. Using the RMS package (version 6.2-0) and survival package (version 3.2-10), nomograms including important clinical features and calibration plots were constructed. The 45° line represents the best-predicted value, and calibration curves were graphically evaluated by mapping the nomogram-predicted probability against observed occurrences. The consistency index (C-index) was used to measure the discriminative capability of the nomogram and to compare the predictive accuracy of nomograms and individual prognostic indicators. This process was calculated using the bootstrap method and repeated 1000 times. In the present study, one-way analysis of variance (ANOVA) and two-tailed Student’s t test were used to analyze the data. A P value < 0.05 was considered statistically significant.

The TargetScan (http://www.targetscan.org), DIANA-microT (http://diana.imis.athena-innovation.gr/DianaTools/index) and RNAinter (http://www.rnainter.org) online sites were used together to predict and analyze the target miRNAs of CFHR4, compare the correlations between the expression of CFHR4 and target miRNAs and screen miRNAs that were more compatible with ceRNA networks. The target lncRNAs of the screened miRNAs were predicted and analyzed using miRNet2.0 (www.mirnet.ca/miRNet/home.xhtml) and starBase3.0 (www.starbase.sysu.edu.cn), and the correlation between the two was further analyzed to screen for additional eligible ceRNAs. A comprehensive analysis of negatively correlated miRNA–mRNA and miRNA-lncRNA expression levels was performed to establish an HCC-related lncRNA-miRNA–mRNA (CFHR4) ceRNA network.

The R package (version 3.6.3) was used for statistical analyses and plotting. CFHR4 expression in unpaired and paired samples was analyzed using the Wilcoxon rank sum test and Wilcoxon signed rank test, respectively, with the pROC (1.17.0.1) package for ROC analysis. In addition, the Kruskal–Wallis test and univariate Cox analysis were applied to investigate whether CFHR4 expression was associated with clinicopathological factors. Using the KM method and log-rank test, we compared the differences in 10-year OS, DSS and PFI between patients with high CFHR4 expression and those with low CFHR4 expression in TCGA. In all studies, a P value < 0.05 was defined as statistically significant.

By analyzing GTEx and TCGA datasets, we investigated the CFHR4 mRNA levels across cancer types using the Wilcoxon rank sum test, including adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), uterine corpus endometrial carcinoma (UCEC) and uterine carcinosarcoma (UCS). We found that CFHR4 expression was significantly decreased in LIHC and CHOL compared with normal tissues (Figure 1A). We obtained similar results from the Timer and GEPIA databases (Supplementary Figures 1A, B). According to the expression of CFHR4 in 374 HCC tissues and 50 normal liver tissues, we confirmed that the CFHR4 expression level was also noticeably decreased in HCC tissues (P<0.001) (Figure 1B). Furthermore, CFHR4 was underexpressed in the GSE14520 HCC cohort (P<0.001) (Figure 1C). Similar results were obtained for adjacent HCC tissues among the 50 matched HCC tissues and adjacent HCC tissues (P<0.05) (Figure 1D). We extracted protein from human normal hepatic cells (WRL68) and HCC cells (BEL7402, SK-hep1, HCCLM3, HepG2 and Huh7) and confirmed the low expression of CFHR4 in HCC cells using Western blot (Figure 1E). Subsequently, 30 pairs of HCC samples were validated, and similar conclusions were reached (Figure 1F). CFHR4 mRNA expression levels were further validated using quantitative real-time PCR analyses (P<0.001) (Figures 1G, H). In addition, we constructed the receiver operating characteristic (ROC) curve. The area under the curve (AUC) for CFHR4 was 0.698, and it has a significant diagnostic value for HCC (Figure 1I).

According to the CFHR4 expression level, we divided the data from patients with HCC into high and low CFHR4 expression groups for comparison. The DESeq2 package was used to infer CFHR4-associated genes and analyze the DEGs between the high and low expression groups. An adjusted p value < 0.05 and absolute log2-fold change > 1.5 were considered statistically significant. A total of 721 significant DEGs were identified. 113 DEGs were associated with the high CFHR4 expression group, and 608 DEGs were associated with the low CFHR4 expression group (Figure 1J and Supplementary Table 2). The top 10 DEGs were identified, further analyzed using HTSeq-Counts and sorted by relative expression (Figure 1K).

GO and KEGG enrichment analyses were performed using the “clusterProfiler” R package to further analyze the potential biological functions of CFHR4-related DEGs. The GO analysis indicated that CFHR4-related DEGs may be involved in gated channel activity, regulation of signal release, regulation of ion transmembrane transport, metal ion transmembrane transporter activity, synaptic membrane, transmembrane transporter complex and passive transmembrane transporter activity (Figures 2A, B; Supplementary Table 3). In the KEGG enrichment analysis, CFHR4-related DEGs were mainly involved in chemical carcinogenesis, retinol metabolism, the calcium signaling pathway, the PPAR signaling pathway, bile secretion, insulin secretion and gastric acid secretion (Figures 2C, D).

GSEA was conducted between the high and low CFHR4 expression groups to further reveal CFHR4-related signaling pathways in HCC. The following pathways were significantly enriched in patients with low CFHR4 expression: FCGR-activated reaction, PLK1 pathway, reactant FCERI-mediated MAPK activation, ATR pathway, MCM pathway, cascade reaction of PI3K and FGFR1, reactant-mediated MAPK activation and FOXM1 pathway (Figures 2E–J; Supplementary Table 4).

We explored the association between 721 DEGs in the HCC group using the STRING database by setting the interaction threshold to 0.70 and constructed a PPI network to further investigate the underlying mechanisms (Figure 3A; Supplementary Table 5). Subsequently, 301 proteins and 420 edges were screened, and five central gene clusters were identified using a total score ≥5000 (Figures 3B–F). In addition, the top 7 central genes were screened, including CENPA, CDC20, UBE2C, CEP55, BIRC5, FAM64A and TRIP13 (Figure 3G). By analyzing the GeneMANIA and STRING online datasets, potential CFHR4-interacting target genes were identified (Supplementary Figures 2A, B). CFHR4-related genes were selected by performing a crossover analysis, including C3, CRP, CFHR1, CFHR3 and CFHR5 (Supplementary Figure 2C). We subsequently analyzed the association between CFHR4 and the 5 intersecting genes (Supplementary Figures 2D–H).

Based on the ssGSEA algorithm, we confirmed and quantified the correlations between CFHR4 expression and the immune cell infiltration levels (Figure 4A). The expression of CFHR4 was negatively correlated with aDCs, TFH cells, NK CD56bright cells and Th2 cells, and it has positive correlations with Th17 cells, DCs, neutrophils, mast cells, Tgd cells, Tcm cells, cytotoxic cells, Tregs, NK cells, pDCs, eosinophils, iDCs, B cells, T cells, CD8 T cells, Tems, NK CD56dim cells, T helper cells, macrophages and Th1 cells (Figures 4B–H). We further confirmed the correlation between CFHR4 expression with immunomarker of various immune cells in HCC. The results showed that CFHR4 expression was significantly correlated with the immunomarkers IRF5 and INOS of M1 macrophages in HCC (Table 1). It indicated that CFHR4 may induce macrophages to M1 polarization in HCC. This analysis of immune markers of different functions T cells showed that CFHR4 expression was highly correlated with the most immunomarkers (CD8B, CD3D, STAT1, IFN-γ, STAT5A, IL21, TGFβ, PD-1, CTLA4, LAG3 and TIM-3) of T cells in HCC (Table 1). It turns out that CFHR4 may perform an indispensable role in the T cells’ immune response to HCC. Especially for T cells exhaustion, consistent results with the GISTIC analysis were obtained. The somatic copy number alteration (SCNA) module demonstrated that the arm-level deletion of CFHR4 was markedly associated with immune cell infiltration levels in HCC (Figure 4I). In addition, the results also showed a correlation between CFHR4 expression and the immunomarkers of TAMs, neutrophils and dendritic cells (Table 1). Subsequently, according to the expression level of CFHR4, HCC samples were dichotomized into CFHR4-high and low expression groups, we aimed to reveal whether different expression groups of CFHR4 differ in the tumor immune microenvironment of HCC (Figure 4J). We found that cytotoxic cells, DCs, iDCs, mast cells, neutrophils, NK cells, pDCs, Tcm cells, Tem cells, Tgd cells, Th17 cells and Tregs were increased in the CFHR4 high expression group (P < 0.05), while the NK CD56bright cells, TFH cells and Th2 cells decreased (P < 0.05). These findings confirmed that reduced expression of CFHR4 in HCC was closely associated with immune cell infiltration.

The clinical data from patients with HCC in TCGA database were obtained to investigate the clinical characteristics of patients with different CFHR4 expression levels. After removing patients with incomplete clinical data, 374 patients remained for further analysis; the average age was 61.5 years (49.25 to 70.00 years), and 67% were male. Table 2 provides a detailed description of the clinical data. We evaluated the differences in clinicopathological variables after stratifying patients based on CFHR4 expression using the Kruskal–Wallis test, and the level of CFHR4 was strongly correlated with age, sex, race, TNM stage, histologic grade, pathological stage, tumor status, residual tumor, vascular invasion and AFP level (Figures 5A–L). Notably, CFHR4 was expressed at higher levels in the older age group (age>60 years) than in the younger age group (age ≤ 60 years) (P<0.05). Significant differences in CFHR4 expression levels were also noted in different races (P<0.001). Moreover, a higher histological grade, TNM grade, pathological stage and tumor status were also significantly associated with low CFHR4 expression. Subsequently, we further confirmed the lower CFHR4 expression level in the group with a high AFP level (>400 ng/mL) (P<0.001). Based on these results, patients with HCC presenting lower CFHR4 expression seemed to have a more advanced tumor stage.

Afterward, we performed a series of studies to determine the association of CFHR4 expression levels with the prognosis of patients with HCC. The Kaplan–Meier Plotter analysis revealed an association between low CFHR4 expression and a poor prognosis (Figures 6A–C). Moreover, we performed subgroup analyses of OS, DSS and PFI. Patients with high CFHR4 expression had a correspondingly better prognosis for OS, DSS and PFI in the Asian group (Figures 6D–F). However, OS, DSS and PFI in the white and black or African–American subgroups were not significantly different (Supplementary Figures 3A–C). In addition, patients with HCC presenting high CFHR4 expression who were aged ≤ 60 years experienced longer OS and DSS but had a worse prognosis in terms of PFI (Figure 6G–I). However, no significant differences were observed in the younger age subgroups for OS, DSS and PFI (age ≤ 60 years) (Supplementary Figures 3D–F). We further confirmed that the T3 and T4 subgroups and the stage III and stage IV subgroups experienced poorer OS (Supplementary Figures 3G, H).

A univariate Cox regression analysis was performed with TNM stage, pathological grade, tumor status and CFHR4 expression levels to further identify factors associated with different prognoses (Supplementary Table 6). The forest plot illustrated that low expression of CFHR4 was a risk factor for the OS (Figure 6J; Supplementary Table 7), DSS (Figure 6K; Supplementary Table 6) and PFI (Figure 6L; Supplementary Table 8) of patients with HCC. According to the results of the univariate Cox regression analysis, CFHR4 expression and other independent clinicopathological factors were used to construct the point scale of the nomogram. Each variable was scored with reference to the scale of the nomogram, and the total scores were dispatched to the outcome line and predicted the prognosis of patients with at 1, 3 and 5 years. The C-index of the nomogram was 0.706 (95% confidence interval: 0.671-0.741). This result suggested that the prognostic nomogram of CFHR4 had good discriminatory power (Figure 6M). The deviation correction line in the calibration analysis approached the ideal curve, indicating that the predicted values were consistent with the observed values (Figure 6N). Consistent results were obtained with the univariate Cox regression analysis.

As reported in previous studies, m6A RNA methylation exerts an important effect on the development of HCC (27–29). The correlations between CFHR4 expression and the expression of 23 m6A-related genes were analyzed in TCGA (Figure 7A). The correlation analysis showed significant negative correlations between the expression of CFHR4 (P < 0.05) and 15 m6A-related genes in HCC (Figures 7B–P). Furthermore, groups were established based on the median CFHR4 expression, and 211 patients were assigned to the high expression group and 210 patients were assigned to in the low expression group. We determined the relationship between the CFHR4 expression level and m6A modification level in HCC by analyzing the differential expression of 23 m6A-related genes in different expression groups (Figure 7Q). The expression of YTHDC1, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, HNRNPA2B1, LRPPRC, HNRNPC, RBMX, METTL16, METTL3, RBM15, RBM15B, VIRMA, WTAP and ALKBH5 was reduced in the high CFHR4 expression group (P < 0.05). In summary, a strong correlation was observed between m6A RNA methylation in HCC and the CFHR4 expression level.

Accumulating evidence highlights the regulatory role of lncRNA–miRNA–mRNA ceRNA networks in cancers. Therefore, we analyzed and constructed a ceRNA regulatory network for CFHR4 in HCC. Through TargetScan, DIANA-microT and RNAinter database predictions, the following 11 miRNAs were jointly predicted: hsa-miR-32-3p, hsa-miR-142-5p, hsa-miR-146a-5p, hsa-miR-302c-5p, hsa-miR-361-5p, hsa-miR-4775, hsa-miR-4786-5p, hsa-miR-4795-3p, hsa-miR-5590-3p, hsa-miR-580-3p and hsa-miR-590-3p (Figure 8A). Based on the regulatory relationship in the ceRNA network, a negative correlation was observed between mRNAs and miRNAs. Four miRNAs negatively correlated with CFHR4 expression were identified and screened by performing a correlation analysis. The scatter plots showed the correlation between CFHR4 expression and the target miRNAs, and the TargetScan database was used to predict the potential binding sites in CFHR4 for target miRNAs (Figures 8B–E). Subsequently, the lncRNAs that may interact with the target miRNAs (hsa-miR-146a-5p, hsa-miR-361-5p and hsa-miR-580-3p) were further predicted using the miRNet and starBase databases (Figures 8F–H). This interaction is due to the negative correlation between the expression of lncRNAs and miRNAs. Consequently, using the starBase database, we further screened and confirmed the lncRNAs in HCC that were negatively correlated with the three target miRNAs. Based on these results, the following 10 ceRNA regulatory networks that play a role in HCC were constructed: TMEM161B-AS1-hsa-miR-146a-5p-CFHR4, CCDC183-AS1-hsa-miR-146a-5p-CFHR4, NEAT1-hsa-miR-146a-5p-CFHR4, MALAT1-hsa-miR-146a-5p-CFHR4, XIST-hsa-miR-146a-5p-CFHR4, DNAAF4-CCPG1-hsa-miR-361-5p-CFHR4, NEAT1-hsa-miR-580-3p-CFHR4, LINC00641-hsa-miR-580-3p-CFHR4, DNAAF4-CCPG1-hsa-miR-580-3p-CFHR4 and DSCAM-AS1-hsa-miR-580-3p-CFHR4 (Figure 8I).