Transcript data from 378 HCC patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). After excluding cases without complete information, clinical profiles of 240 patients were downloaded and are shown in Supplementary Table S1. DNA methylation information was obtained from the UCSC Cancer Browser (https://xena.ucsc.edu/). 430 DNA methylation data and annotation files were downloaded from the Illumina Human Methylation 450 platform (https://xenabrowser.net/datapages/?dataset=TCGA-LIHC.methylation450.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) for further study. The threshold for rejection of the CpG sites is given below: The site was missing in more than 70% of the samples; it was located in the sex chromosomes and single-nucleotide polymorphisms; CpGs above 2 kb upstream to 0.5 kb downstream; and cross-reactive genome CpG sites (Chen et al., 2013). Ultimately, we employed 226 HCC samples with both expression matrix and methylation data for further analysis. We obtained fragments per thousand base million (FPKM) of HCC patients from the TCGA database and converted the FPKM value to the transcript per million (TPM) value. The LIRI-JP dataset including 260 HCC samples from the ICGC database was employed as the external validation group. The corresponding expression profiling information and the clinical data (Supplementary Table S2) were downloaded from the ICGC database (https://dcc.icgc.org). Four categories of somatic mutation data of HCC samples were obtained from The Cancer Genome Atlas (TCGA) portal. We singled out the mutation files that were obtained through the “SomaticSniper variant aggregation and masking” platform for subsequent analysis. We prepared the Mutation Annotation Format (MAF) of somatic variants and implemented the “maftools” (Mayakonda et al., 2018) R package.

In order to classify prognostic-related HCC molecular subtypes, the methylation level of each CpG site, age, gender, clinicopathological stage, tumor status (T, M, and N), and survival data were employed to construct univariate Cox proportional risk regression models. Supplementary Table S3 presents 610 significant CpG sites (p < 0.05). Then, these significant sites were introduced into multivariate Cox proportional risk regression models. Finally, we selected the CpG sites that were still significant as the classification features (p < 0.05, Supplementary Table S4) that can significantly affect the prognosis of HCC patients.

In order to define the HCC subgroups based on the most variable methylated CpG sites (Supplementary Table S5), consensus clustering was conducted using the ConcensusClusterPlus R package (Wilkerson and Hayes, 2010). Each case we collated was partitioned into k groups on the basis of a user-specified clustering algorithm (k-means). This process was repeated for a user-specified number of repetitions, offering a method of establishing consensus values and estimating the stability of the established clustering. Pairwise consensus values, referring to “the proportion of clustering runs in which two subjects gathered together,” were calculated and stored in a consensus matrix for each k. The final results from agglomerative hierarchical consensus clustering based on one-Pearson correlation distances were also divided into k groups. The exported graphs contained the consensus matrices, consensus cumulative distribution function (CDF) plot, delta area plot, and tracking plot. We determined the k value if there were a low relative change in the area under the CDF curve, relatively high conformity in the clusters, and a low variation coefficient. The pheatmap package in R was used to create the heatmap (Diao et al., 2018).

To explore the overall survival of HCC patients whose subgroups were defined from DNA methylation information, we employed the “survival” package (Huang et al., 2017; Li et al., 2019; Xiong et al., 2019), and the results are illustrated in Kaplan–Meier plots. The log-rank test was used to identify the significant differences among different clusters. To determine the distinction in categorical profiles as the clinical variables among different subtypes, the chi-squared test was analyzed. All tests were two-sided with a significance level of 0.05.

To further reveal the CpG site–related gene expression pattern in the standpoint of fundamental biology, Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) enrichment analyses on CpG site–related genes significantly affected prognosis. We employed gene set enrichment analysis (GSEA) to explore underlying mechanisms correlated with the prognostic signature. The gene sets of “c2. cp.kegg.v7.2. symbols.gmt [Curated]” from the Molecular Signatures Database were analyzed through gene set enrichment analysis (GSEA) (Subramanian et al., 2005) with a Java program (http://software.broadinstitute.org/gsea/index.jsp).

Two hundred twenty-six HCC samples were randomly divided into the discovery and validation groups at the ratio of 1:1 for integrated analysis utilizing the “caret” package. Both discovery and validation groups were required to meet the following criteria (Forner et al., 2018): samples were randomly assigned to discovery and validation cohorts (Bray et al., 2018) and the clinical characteristics of subjects in these groups were similar (Supplementary Table S6). The validation group with 112 cases and the whole cohort were used to validate results obtained from the discovery set. The TCGA cohort was randomly classified into discovery and validation groups again to demonstrate the repeatability of the risk model. The detailed clinical variables of samples are recorded in Supplementary Table S7. The “survival” package with coxph function was employed to create a Cox proportional hazards model based on the combination of expression profiles for remarkably methylated CpG sites corresponding to genes in cluster 1 (Supplementary Table S8) and follow-up data (Zhang, 2016; Bhattacharjee et al., 2020). The candidate methylation-related genes significantly affect prognosis (p < 0.05), which were identified by employing a univariate proportional hazards model of the expression level of seven methylation-related genes. Next, the risk coefficient of each gene was obtained by performing the LASSO regression algorithm with the “glment” package after the deletion of highly correlated genes. Finally, two genes were selected and introduced into a prognostic predictive signature in HCC. The risk score of each patient was calculated by the following equation: risk score = sum of risk coefficients * methylation-related gene expression level.

Each patient obtained the risk score according to the above formula together with their clinical data. Based on the median risk score, patients were classified into low/high‐risk subgroups for further study. Kaplan–Meier survival curves were analyzed to compare prognosis of different subgroups. We then performed the time-dependent receiver-operating characteristic (ROC) curve analysis to assess the prognostic value, which was achieved by comparing the specificity and sensitivity in predicting prognosis on the basis of the risk score. Furthermore, multivariate Cox regression was employed to demonstrate the risk score was an independent prognostic indicator for HCC patients. The predictive performance of the as-constructed risk signature was then confirmed in the validation group (n = 112) and combined cohort. Additionally, the ICGC cohort was employed as an external validation group to facilitate extensive application. Each test was two-sided, and p < 0.05 was deemed statistically significant.

To comprehensively assess the prognosis predictive ability of risk signature, stage, gender, age, and WHO grade for one/two/three-year OS, time-dependent receiver-operating characteristic (ROC) curves were performed to calculate the area under the curve (AUC) values (Blanche et al., 2013). To contribute to a quantitative manner predicting the overall survival of patients with HCC, we established a nomogram containing the risk score and other clinical variables to estimate one‐, two‐, and three‐year overall survival possibility. Subsequently, we analyzed the calibration curve that shows the prognostic value of the as-constructed nomogram.

To better identify the difference of TIME characteristics between different subgroups, we performed several analyses as follows. Firstly, the “GSEABase” R package with regard to 29 immunity-related signatures was used to further reveal distinction of TIME characterization between different subgroups. Then, the R package “CIBERSORT” was performed to calculate 22 immune cell subpopulations in HCC. Subsequently, immune infiltration information consists of every specimen’s immune cell fraction (i.e., B cells, CD4⁺ T cells, CD8⁺ T cells, dendritic cells, macrophages, and neutrophils) downloaded from tumor immune estimation resource (TIMER) (https://cistrome.shinyapps.io/timer/).

Referring to existing studies, it was found that the expression level of immune checkpoint blockade–related key genes might be correlated with the clinical outcome of immune checkpoint inhibitor blockade treatment⁴². Herein, we employed six key genes of immune checkpoint blockade therapy: programmed death ligand 1 (PD‐L1, also known as CD274), programmed death 1 (PD‐1, also known as PDCD1), programmed death ligand 2 (PD‐L2, also known as PDCD1LG2), T‐cell immunoglobulin domain and mucin domain–containing molecule‐3 (TIM‐3, also known as HAVCR2), indoleamine 2,3‐dioxygenase 1 (IDO1), and cytotoxic T‐lymphocyte antigen 4 (CTLA‐4) in HCC (43–45). To elucidate the potential player of the as-constructed risk signature in ICB treatment of HCC, we correlated the prognostic signature and the expression level of six immune checkpoint blockade key genes (i.e., IDO1, CTLA-4, HAVCR2, CD274, PDCD1, and PDCD1LG2). Subsequently, we systematically determined the expression value of 47 immune checkpoint blockade–related genes (i.e., PDCD1) between patients from different subgroups to investigate the potential role of the risk score in ICB treatment.

To further explore the role of LRRC41 in TIME, we employed single-cell transcriptome sequencing data from GSE140228 (Zhang et al., 2019), which are the transcriptome data of CD45⁺ immune cells made by the Zemin Zhang team for HCC patients. The researchers uploaded the hepatic carcinoma single-cell RNA sequencing data of the study to an interactive website (http://cancer-pku.cn:3838/HCC/) to facilitate the researcher’s in-depth exploration of related fields. Herein, we use 10× Genomics sequencing data to analyze the expression of LRRC41 in tumor, hepatic lymph node, adjacent liver, ascites, and blood and compare the abundance of LRRC41 in immune cell subsets in tumor tissues.

TMB was defined as the number of somatic, coding, base replacement, and insertion–deletion mutations per megabase of the genome examined using non-synonymous and code-shifting indels under a 5% detection limit. TMB scores for each sample were calculated through dividing the number of somatic mutations by the total length of exons (38 million). The R package “maftools” (Mayakonda et al., 2018) was used to calculate the total number of somatic non-synonymous point mutations within each sample.

L02 (human hepatic cell line) and two human HCC cell lines (MHCC-97H cells and HCC-LM3 cells) were purchased from the Cell Bank of the Type Culture Collection of the Chinese Academy of Sciences, Shanghai Institute of Biochemistry and Cell Biology. The cell lines were all cultured in Dulbecco’s minimum essential media (DMEM) plus 10% fetal bovine serum (FBS; Invitrogen, Carlsbad, CA, United States). All cell lines were grown without antibiotics in a humidified atmosphere of 5% CO₂ and 99% relative humidity at 37°C. Three different cell lines were subjected to quantitative real-time polymerase chain reaction (qRT-PCR).

Total RNA was extracted from cells using TRIzol (Invitrogen, Carlsbad, CA, United States) according to provided instructions. RNA concentration and purity were measured in triplicates utilizing the NanoDrop 2000 spectrophotometer (Thermo Scientific Inc., Waltham, MA, United States). Then, total RNA was reverse transcribed to cDNA using the cDNA Reverse Transcription Kit (Vazyme, Nanjing, China). To determine the expression of LRRC41, cDNAs were subjected to qRT-PCR using SYBR Green Real-Time PCR Master Mix (Takara) in Applied Biosystems 7500/7500 Fast Real-Time PCR System (Thermo Fisher Scientific). All samples were analyzed in triplicates. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) levels were used as the endogenous control, and the relative expression of LRRC41 was calculated using the 2-ΔΔCt method. The sequences of primers used for PCR were as follows: LRRC41, 5′- TGG​CTG​GCG​AGA​AGG​AGG​ATG -3′ (forward) and 5′- CAA​GGT​GGA​GAT​GCT​GCG​GAA​TC -3′ (reverse), and GAPDH, 5′-CAG​GAG​GCA​TTG​CTG​ATG​AT-3′ (forward) and 5′-GAA​GGC​TGG​GGC​TCA​TTT-3′ (reverse).

The correlation between data with a normal distribution was analyzed with Pearson’s correlation analysis. Comparisons of two groups were made with the t-test; if there were more than two groups, the comparison was made with the analysis of variance (ANOVA) test. The Wilcox test was used to compare the methylation levels between the seven clusters. Samples with a CIBERSORT output value of p < 0.05 were screened for further analysis (Chen et al., 2018). All the statistical analysis results described below were obtained with the R software. The significance level was set to p < 0.05.

After methylation data were obtained from the TCGA cohort and pretreated as described above, we screened 19392 CpG sites among which 610 methylation sites were determined as candidate DNA methylation indicators for prognostic prediction in HCC patients via univariate Cox regression analysis (Supplementary Table S3). A subsequent multivariate Cox regression analysis was employed to identify these CpG sites closely correlated with prognosis. And 136 sites were identified and deemed potential prognosis-related CpG sites (Supplementary Table S4).

Taking advantage of consensus clustering, we found that when k = 7, there is high consistency in these clusters with a relatively low coefficient of change but no appreciable variation in the area under the CDF curve and delta area (Supplementary Figures S1A,B). The heatmap of the seven clusters with a good polymerization effect was mostly diagonal (Supplementary Figure S1C). Kaplan–Meier analysis indicated there was a remarkable distinction among different clusters (p = 1.221e−15), suggesting that clustering results classified the patients into seven categories with distinct overall survival (Figure 1A). The distribution of clinical characteristics (age, gender, historical staging, T status, N status, and M status) in each subtype is presented in Figure 1B. Subsequently, the expression level of the CpG site–related genes was determined in the subtypes. Figure 1C presents the heatmap of gene expression values. To better comprehend the methylation site expression from the viewpoint of biology level, KEGG and GO enrichment analyses on CpG site–related genes were employed. Functional annotation results indicated that overexpressed methylation-related genes were primarily enriched in autophagy, endosome membrane, protein C−terminus binding, and p53 signaling pathway (Figures 1D–G, corresponding GO and KEGG, respectively). Subsequently, we analyzed intra-cluster fraction for the seven clusters based on clinical characteristics (Supplementary Figures S2A–G, age, gender, stage, and TNM category, respectively). Tendencies for correlations between clinical parameters and categories were as follows: clusters 1, 4, 6, and 7 with higher clinical grade; clusters 1 and 5 with advanced stage; clusters 4, 5, and 7 with higher T status; clusters 5, 6, and 7 with higher N status; and clusters 1 and 6 with higher M status. These results suggest that each clinical variable correlated with a distinct intra-cluster distribution.

One hundred thirty-six prognostic methylation sites’ methylation levels are shown in Supplementary Table S5. A heatmap presents distinction of methylation levels between selected clusters and the other clusters (Figure 2A). The heatmap result shows that the significantly changed CpG sites were mostly distributed in cluster 1. Combined with survival rate analysis (Figure 1A), indicating that cluster 1 had a middle survival time among these clusters, we picked the differential CpG sites’ methylation data according to cluster 1 to develop boxplot. A consistent result was obtained in methylation level analysis, which shows that cluster 1 had a middle methylation level among different clusters (Figure 2B). The differential methylated CpG sites’ corresponding genes in cluster 1 were introduced as candidate genes for further study.

To further reveal the potential interaction between these candidate methylation genes, correlation was employed to comprehensively visualize the methylation interaction network (Figure 2C; Supplementary Figure S3A, p < 0.05). Notably, LRRC41 was positively and significantly correlated with five hub genes, suggesting its indispensable role in methylation regulation. In order to investigate the prognosis predictive value of differential methylated CpG sites’ corresponding genes from cluster 1 (Supplementary Table S8), univariate Cox regression on candidate genes’ expression profiles was performed, finding two out of seven candidate genes were significantly correlated with overall survival (Supplementary Table S9, p < 0.05). Next, the LASSO algorithm was employed to identify methylation-related genes with the most robust prognosis prediction capability (Supplementary Figures S3B,C; Table S10). Finally, the two candidate genes LRRC41 and KIAA1429 were introduced into a methylation-based risk prognostic signature in HCC. The relationship of each methylation-associated gene with overall survival is presented in Supplementary Table S11.

The risk score was obtained as follows: risk score = (0.1682 ∗ expression value of LRRC41) + (0.0981 ∗ expression value of KIAA1429). Then, each HCC patient together with the corresponding risk score was classified into low/high-risk groups based on the median value.

Patients in the TCGA-LIHC cohort were randomly divided into the discovery set and validation set (Supplementary Tables S11, S12). The distributions of two methylation-related genes’ expression levels with patients and corresponding groups are displayed in Supplementary Figure S4A. Supplementary Figures S4D,G show that distributions of the risk score and survival status are plotted in, indicating that low-risk HCC patients had better prognosis. The Kaplan–Meier curve further demonstrated that low-risk patients had significantly longer overall survival relative to patients with high risk (p = 3.473e-01; Figure 3A). To estimate the specificity and sensitivity of risk signature for differentiating low-risk patients from high-risk patients, we analyzed ROC curve analysis. We observed that the area under curves of prognostic signature at one‐year overall survival was up to 0.706, indicating encouraging efficiency of prognostic performance (Figure 3D). Furthermore, univariate Cox analysis showed that the risk score significantly correlated with overall survival [hazard ratio (HR): 1.473, 95% CI: 1.268–1.711, p < 0.001; Figure 3G], and multivariate Cox analysis further presented that the risk score was an independent prognostic factor in HCC (HR: 1.448, 95% CI: 1.225–1.711, p < 0.001; Figure 3J).

In order to better assess its prognosis predictive accuracy, we examined above results in the validation group and combination cohort. The figures show the distributions of methylation-related genes’ expression levels, survival, and risk scores of patients in the testing set (Supplementary Figures S4B,E,H) and the combination cohort (Supplementary Figures S4C,F,I). The survivorship curve shows that patients with low risk had higher survival probability than high-risk patients in both the validation group (Figure 3B, p = 3.725e-02) and the combination cohort (Figure 3C, p = 1.852e-03). The value of area under the ROC curve (AUC) was up to 0.713 in the testing set (Figure 3E) and 0.712 in the combined cohort (Figure 3F), indicating strong prognosis predictive power in different groups. Likewise, the risk signature was an independent prognostic indicator significantly correlated with overall survival in both univariable and multivariable regression analyses of the testing set as well as the combined cohort (Figures 3H,I,K,L). Besides, the signature was employed in the LIRI-JP cohort to validate the external prognosis predictive performance. Figures 4A–C show the distributions of genes’ expression levels, risk scores, and survival time in the ICGC validation cohort. Survival analysis (p = 2.095e-03; Figure 4D) and ROC curve analysis (AUC = 0.675; Figure 4E) also indicated that this risk score model had an excellent prognosis predictive performance in the LIRI-JP group. Additionally, consistent results were obtained and great prognostic value of the risk model was corroborated in another random classification (Supplementary Figures S5, S6, corresponding training set and testing set, respectively).

The results of correlation of the risk score with clinical variables presented that the advanced clinical stage (most p < 0.05, Figure 4F) and risk score remarkably increased. However, there was no significant difference of risk score in age subgroups and gender subgroups (Figures 4G,H).

Subsequently, we plotted an ROC analysis curve, and the observed AUC value at one-, two-, and three-year overall survival was 0.693, 0.620, and 0.628, respectively, suggesting excellent prognostic prediction performance (Figure 5A). In order to compare the prognostic predictive efficiency of risk signature with other traditional clinical parameters (age, gender, stage, and grade), we gathered above clinical characteristics and then performed the ROC curve analysis for one-, two-, and three-year survival time and demonstrated that the value of AUC of risk score was the highest (Figures 5B–D). The risk score and clinicopathological stage were introduced into plotting a nomogram to quantitatively predict the overall survival of HCC patients (Figure 5E). Age, gender, and grade were eliminated owing to AUCs lower than 0.55. Calibrated curves demonstrated good prognosis prediction accuracy of one-, two-, and three-year overall survival in the as-constructed nomogram (Figures 5F–H).

Furthermore, to confirm whether prognostic signature remained a robust prognosis prediction validity indicator in patients subdivided into different subtypes according to clinicopathological variables, stratification analysis was performed. Relative to patients with low risk, patients in the high-risk group presented lower survival probability in both the young (<=65) and old (>65) subgroups (Supplementary Figures S7A,B). And the prognostic signature presented a powerful prognostic prediction value in both patients gendered male and female (Supplementary Figures S7C,D) and patients in the early stage (Supplementary Figure S7E), whereas there was no significant difference of overall survival between low/high-risk patients in the advanced stage (Supplementary Figure S7F). These results indicated that the prognostic signature can be a novel and powerful overall survival predictor in HCC.

In order to explore whether the risk score could be a potential indicator of TIME, we analyzed the relationship of the risk score with ssGSEA signatures and TIC proportion and level (assessed using the CIBERSORT algorithm). Firstly, the CIBERSORT results showed that activated memory CD4 T cells were significantly more abundant in high-risk patients, whereas patients with low risk presented higher infiltration of resting memory CD4 T cells (Figure 6A), indicating the phenotype of memory CD4 T cells contributes to the formation of molecular risk to further predict prognosis. Then, the population of aDCs, macrophages, and Tregs and expression of MHC class I were more in the high-risk group; meanwhile, the infiltration of B cells, neutrophils, NK cells, and pDCs, cytolytic activity, T-cell costimulation, and IFN response were higher in the low-risk group (Figure 6B). Furthermore, we analyzed whether the prognostic signature was correlated with immune infiltration. We observed that the risk score had significantly positive correlation with infiltrating B cells (r = 0.239; p = 3.995e−06), infiltrating CD4 T cells (r = 0.226; p = 1.252e−05), infiltrating CD8 T cells (r = 0.195; p = 1.711e−04), infiltrating dendritic cells (r = 0.341; p = 2.051e−11), infiltrating macrophages (r = 0.357; p = 2.133e−12), and infiltrating neutrophils (r = 0.354; p = 3.099e−12; Figures 6C–H).

The above results provided robust evidence to support that the methylation sites’ prognostic signature may act as a vital player to elucidate the context of TIME and further predict clinical immunotherapy efficiency for HCC patients.

Subsequently, we correlated six key immune checkpoint inhibitor genes (PDCD1, CD274, PDCD1LG2, CTLA‐4, HAVCR2, and IDO1) (Kim et al., 2017; Nishino et al., 2017; Zhai et al., 2018). And we analyzed the correlation between the prognostic signature and ICB vital genes to explore its potential player in immunotherapy (Figure 6I). We observed that the risk score had significantly positive correlation with CD274 (r = 0.2; P = 1e−04; Figure 6J), CTLA4 (r = 0.17; p = 0.0014; Figure 6K), HAVCR2 (r = 0.24; p = 2.4e−06; Figure 6L), IDO1 (r = 0.13; p = 0.014; Figure 6M), and PDCD1LG2 (r = 0.11; p = 0.031; Figure 6N), suggesting the risk prognostic signature might play a crucial role in the monitoring of the ICB therapy outcome in patients with HCC. According to correlation analysis, we observed that 15 of 47 ICB-related genes’ (i.e., CD274) expression levels were remarkably higher in patients with high risk relative to low-risk ones (Supplementary Figure S8A). These results demonstrated that the methylation sites’ risk signature may provide a novel approach to predict immunotherapeutic efficiency in HCC.

To better understand the potential player of the methylation sites’ prognostic signature mediated in the underlying mechanism of HCC, we performed GSEA analysis in not only the high‐risk group but also the low‐risk group. GSEA enrichment analysis results presented that the high-risk score mainly enriched in pathways, including the ERBB signaling pathway, Wnt signaling pathway, mTOR signaling pathway, and VEGF signaling pathway (Supplementary Figure S8B).

LRRC41 whose expression level was upregulated was the methylation site–related gene and deemed the negative indicator. Thus, the potential player of LRRC41 in HCC was investigated in further validation experiments. Firstly, we determined the expression level of LRRC41 between paracancerous samples and cancer tissues according to the TCGA database. Compared with normal samples, the expression level of LRRC41 was higher in tumor samples (Figure 7A). By using qRT-PCR, we compared the expression level of LRRC41 between two different tumor cell lines and normal liver cell line. In accordance with previous results, LRRC41 was overexpressed in HCC cells compared with hepatic cells (Figure 7B). In order to better evaluate the prognosis predictive ability of LRRC41, survival curve analysis was performed and plotted between LRRC41 high- and low-expressed patients. We found that lower LRRC41 expression significantly suggested higher overall survival probability (Figure 7C, p = 6.539e−04). The expression level analysis among major clinical stages showed that LRRC41 was expressed significantly different among distinct clinicopathological stages (Figure 7D, F = 3.68 and p = 0.0124). We observed that the higher the N status, the higher the LRRC41 expression level (Figure 7E). Furthermore, the expression level of LRRC41 was significantly and negatively correlated with the methylation level of LRRC41 (Figure 7F, r = −0.129, p = 8.735e–03).

To better research the association between the LRRC41 expression level and immune infiltration, we explored the correlation between the expression level of LRRC41 and immune infiltration via TIMER. The results indicated that the LRRC41 expression was significantly correlated with B cells (r = 0.333; p = 2.37e−10), CD8⁺ T cells (r = 0.298; p = 1.91e−08), CD4⁺ T cells (r = 0.369; p = 1.67e−12), macrophages (r = 0.448; p = 3.02e−18), neutrophils (r = 0.5; p = 3.23e−23), and dendritic cells (r = 0.485; p = 1.97e−21; Figure 8A). Furthermore, distinct mutational types of LRRC41 were correlated with infiltration of neutrophils and CD8⁺ T cells (Figure 8B).

Next, we analyzed the correlation between the LRRC41 expression level and ICB key genes’ expression levels adjusted by tumor purity by TIMER to explore the potential role of LRRC41 in ICB treatment. TIMER results showed that LRRC41 presented significantly positive correlation with CD274 (r = 0.537; p = 3.24e−27), CTLA4 (r = 0.173; p = 1.24e−03), HAVCR2 (r = 0.472; p = 1.66e−20), IDO1 (r = 0.228; p = 1.99e−05), PDCD1 (r = 0.202; p = 1.64e−04), and PDCD1LG2 (r = 0.325; p = 6.31e−10; Figures 8C–H), suggesting LRRC41 has a vital role in ICB therapy.

In order to further reveal the role of LRRC41 in the formation of characteristics in TIME of HCC, we performed correlation analysis of expression value of LRRC41 with ssGSEA enrichment (by the GSEABase method) and immune infiltration fraction and level (using the CIBERSORT algorithm) and further conducted single-cell transcriptome sequencing data analysis. Patients with HCC were divided into low/high-LRRC41 subtypes based on the median value of LRRC41 expression level. The CIBERSORT results presented that the LRRC41 expression level was positively correlated with the abundance of macrophages M0 and M1 whereas negatively correlated with monocytes and macrophage M2 (Figure 9A). According to ssGSEA results, inflammation promotion and parainflammation were activated, CCR and HLA were upregulated, and infiltration of macrophages and T helper cells was remarkably elevated with the LRRC41 expression being increased (Figure 9B). According to the results of single-cell transcriptome sequencing data analysis, we observed that the expression value of LRRC41 is more in tumor tissues than paracancerous tissues (Figure 9C). And LRRC41 is mainly expressed in Mast-c1-IL7R cells and Mast-c2-CPA3 cells in HCC tumor tissues (Figure 9D). Figures 9E,F present the distribution of LRRC41 in immune cells of tumor. Based on the previous findings, mast cell proteases play the crucial role in promoting tumor angiogenesis (de Souza Junior et al., 2015), suggesting LRRC41 may act as a positive player in HCC progression.

To further investigate the potential role of the LRRC41 methylation level in TIME characterization, HCC samples were grouped into hypermethylation and hypomethylation subtypes based on the median value of LRRC41 methylation level. The CIBERSORT results pointed out that the infiltration of memory B cells was elevated, but plasma cells were downregulated in the LRRC41 hypermethylation group (Figure 10A).

Additionally, APC costimulation, T-cell costimulation, and parainflammation were activated, Checkpoint, CCR and HLA were upregulated and infiltration of aDCs, CD8⁺ T cells, DCs, macrophages, neutrophils, pDCs, T helper cells, Tfh, TIL, and Treg was remarkably elevated with LRRC41 methylation being decreased (Figure 10B). Collectively, these results suggested that both methylation and expression levels of LRRC41 might be pivotal players in the context of TIME and immune response of HCC.

To further reveal the biological role of KIAA1429 in immune cell infiltration, the correlation of the expression value of KIAA1429 with immune cell infiltration was analyzed via TIMER. The results indicated that the KIAA1429 expression was significantly correlated with B cells (r = 0.279; p = 1.39e−07), CD8⁺ T cells (r = 0.121; p = 2.54e−02), CD4⁺ T cells (r = 0.269; p = 4.20e−07), macrophages (r = 0.232; p = 1.53e−05), neutrophils (r = 0.274; p = 2.46e−07), and dendritic cells (r = 0.278; p = 1.85e−07; Supplementary Figure S9A).

Next, the expression value of KIAA1429 in normal tissues and tumor samples was analyzed based on the TCGA database. Relative to normal tissues, KIAA1429 was upregulated in cancer samples (Supplementary Figure S9B). To assess the prognostic value of KIAA1429, the survival curve was analyzed between KIAA1429 high- and low-expressed groups. The result showed that patients with low LRRC41 presented significant advantage of overall survival time (Supplementary Figure S9C, p = 0.0039).

Subsequently, the correlation of the KIAA1429 expression level and immunotherapy key genes’ expression levels adjusted by tumor purity by TIMER was analyzed to uncover the potential player of KIAA1429 in immunological treatment. TIMER results showed that KIAA1429 presented significantly positive correlation with CD274 (r = 0.226; p = 2.23e−05), CTLA4 (r = 0.213; p = 6.45e−05), HAVCR2 (r = 0.264; p = 6.29e−07), IDO1 (r = 0.117; p = 2.99e−02), PDCD1 (r = 0.15; p = 5.17e−03), and PDCD1 (r = 0.141; p = 8.90e−03; Supplementary Figures S9D–H), suggesting KIAA1429 has a vital role in immunotherapy.

To further explore the biological mechanism of methylation regulation, correlation networks based on LRRC41 and KIAA1429 were constructed, respectively. A total of 47 interactors were identified in the LRRC41-based network (Supplementary Figure S10A), while 186 interactors were determined in the KIAA1429-based network (Supplementary Figure S10B). We reviewed the literature correlated with these interactors in HCC, CUL5 (Ma et al., 2013), RNF7 (Yu et al., 2018), and SOCS1 (Yang et al., 2020), which are potential targets of LRRC41 in methylation regulation of HCC. Additionally, KIAA1429 might interact with EGFR (Ye et al., 2016), HSPA8 (Khosla et al., 2019), and HSP90AA1 (Shi et al., 2020) to modulate methylation in HCC. As such, these underlying targets exhibited promising potential to act as critical regulators involving in the DNA methylation in HCC and further mediated tumorigenicity and progression.

As summarized in the waterfall map, 327 out of 364 HCC patients had the somatic mutation altered, accounting for 89.84%. And we observed that TP53, CTNNB1, and TTN mutations are the top three mutated genes in HCC samples, and the frequency was 30, 25, and 24%, respectively (Supplementary Figure S11A). Missense mutations occupied an absolute position in the total mutation classification (Supplementary Figure S11Ba), and single-nucleotide polymorphism (SNP) accounted for more proportion than deletion (DEL) or insertion (INS, Supplementary Figures S11Bb,Be). Meanwhile, C > T had the highest frequency, 13933 times, in variant types of single-nucleotide variant (SNV) (Supplementary Figures S11Bc,D). Supplementary Figure S11Bd presents that the number of variants per sample and the median value of mutation variants were 71. Besides, the top 10 genetical variated genes were TP53, TTN, CTNNB1, MUC16, ALB, PCLO, MUC4, APOB, RYR2, and ABCA13 (Supplementary Figure S11Bf). The rainfall plot of the sample TCGA−UB−A7MB−01A−11D−A33Q−10 is presented in Supplementary Figure S11C. Each dot represents the SNV mutation type with corresponding color. To further elucidate the intrinsic connection between these genetic altered genes, the exclusive and co-occurrence correlations are presented in Supplementary Figure S11E. To further reveal the intrinsic connection of gene mutation status with DNA methylation, the top three mutated genes (TP53, TTN, and CTNNB1) were fetched for correlation analysis. The results showed that the TP53 mutation and CTNNB1 mutation were significantly higher in the KIAA1429 hypomethylation group (Figures 10G,H), but not TTN mutation (Figure 10F). However, there was no significant correlation of the LRRC41 methylation level with gene mutation (Figures 10C–E).