Gene expression profiles of HCC along with FPKM and count number expression profile data and clinical information of normal samples were downloaded from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/) (TheLegacy. Cell (20, 2018). The expression profile was divided into lncRNA and mRNA by GENCODE gene annotation file. To retain the biologically significant genes, genes with expression values of 0 in both disease and normal samples were excluded. The work flow chart is shown in Figure 1.

We used the immune function collection of the ImmPort database (Bhattacharya et al., 2014), which contains a large number of human immune-related genes. We obtained a total of 17 immune function-related pathways, namely, antigen processing and presentation, antimicrobials, BCR signaling, chemokines, chemokine receptors, cytokines, cytokine receptors, interferons, interferon receptors, interleukins, interleukin receptors, natural killer cell cytotoxicity, TCR signaling pathway, TGFb family member, TGFb family members receptors, TNF family member, and TNF family member receptors. These 17 immune pathways contained a total of 1,811 related protein-coding genes available for download from Gene Lists (https://www.immport.org/shared/genelists).

To investigate the relationship between gene co-expression in HCC samples and normal samples, the correlation between mRNA and lncRNA expression in the samples was calculated using Pearson’s correlation method. Logarithmic transformation of the two expression spectra was performed to ensure expression values close to normal distribution. Next, the Cor.test function in R language (Dou et al., 2020) was used to obtain the Pearce correlation coefficient and the significant p value of all mRNAs and lncRNAs in disease and normal samples.

The enrichment between lncRNA and 17 immune function-related pathways was studied using gene set enrichment analysis (GSEA) (Subramanian et al., 2007) to investigate the relationship between lncRNAs and immune function. The R-values and p-values of the correlation coefficients between the expression of each pair of mRNAs and lncRNAs in liver cancer and normal samples represented the two score values for genetic correlations with the following formula: Score = − log 10 P × sign ( R )

For each lncRNA, all mRNAs were ranked from small to large based on their correlation scores, and the significance of enrichment results and lncRES score (Li et al., 2020) was calculated for each pair of lncRNA and immune function-related pathway using the GSEA method. Under the condition of false discovery rate (FDR) < 0.05 and |lncRES| > 0.995, lncRNA regulators of immune function in HCC samples and normal samples were obtained.

To investigate the differential expression of lncRNAs in HCC, an analysis was performed using the edgeR in R package (Robinson et al., 2010). The DGEList object was first developed through the DGEList method (DeLisi et al., 2018) with count number expression profiles of 373 disease samples and 50 normal samples. Next, genes with CPM >1 in at least two samples were retained. estimateDisp (Cheng et al., 2020) was used to estimate data dispersion, and the negative binomial generalized log-linear model (exactTest) (Shirazi et al., 2016) was used to analyze differential expression of the samples. Finally, significantly differentially expressed lncRNAs in HCC were identified with FDR < 0.05.

To further verify the immune role of dysregulated immune lncRNA in liver cancer, marker genes of 24 immune cells (CD4_naïve, CD8_naïve, Cytotoxic, Exhausted, Tr1, nTreg, iTreg, Th1, Th2, Th17, Tfh, Central_memory, Effector_memory, NKT, MAIT, DC, Bcell, Monocyte, Macrophage, NK, Neutrophil, Gamma_delta, CD4_T and CD8_T) were obtained from a previous study (ImmuCellAI) (Miao et al., 2020).

Next, dysregulated immune lncRNAs showing a significant correlation with marker genes of 24 immune cells were extracted, and a hypergeometric enrichment analysis method was used to identify the significant enrichment relationship between candidate dysregulated immune lncRNAs and immune cells. The hypergeometric calculation formula was as follows: p i j = C M k C N − M n − k C N n

Specifically, p _ij represents the enrichment significance of a marker gene to an lncRNA I and immune cell J; N represents the number of mRNA showing significant correlation with lncRNA I in HCC samples; M represents the number of marker genes in immune cell J; K represents the number of significant correlation between lncRNA I and marker gene of immune cell J; and N represents the total number of mRNAs. In addition, when K was less than 3, the test results were not significant. Based on this method, we calculated the expression correlation and enrichment significance between each candidate lncRNA and 24 immune cells and identified a significant enrichment relationship between all dysregulated immune candidate lncRNA regulators and immune cells under the threshold of p < 0.05. Finally, we screened dysregulated immune lncRNAs specific to HCC and significantly enriched in at least 10 immune cells.

According to the expression values of the dysregulated immune lncRNAs in cancer samples, the R-packs were used in unsupervised clustering of HCC samples.

Survival analysis was conducted on samples according to survival time and survival status of different sample categories. Samples with survival shorter than 30 days were excluded to ensure the accuracy and validity of analysis results. Survival rate and median survival time were evaluated by the Kaplan–Meier method, and the subclasses of different samples were compared among each group based on the log-rank test.

Differences of biological characteristics among different types of liver cancer samples were analyzed based on immunological characteristics of TCGA HCC samples (macrophage regulation, intratumor heterogeneity, proliferation, wound healing, IFN-gamma response, and TGF-beta response) from a previously published article (Thorsson et al., 2018).Wilcoxon rank sum test was used to compare biological characteristics and characterization of subclasses.

HCC cell lines HepG2, HuH7, Li7, SNU387, and SNU182 were kindly provided by the Stem Cell Bank, Chinese Academy of Sciences, and MHCC97H cells were obtained from Beyotime Company (Shanghai, China). The human liver cells THLE3 were obtained from ATCC. HepG2, HuH7, and MHCC97H cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. Li7, SNU387, and SNU182 cells were cultured in Roswell Park Memorial Institute (RPMI) 1640 medium supplemented with 10% FBS and 1% penicillin/streptomycin. Then, THLE3 cells were cultured in Endothelial Cell Medium (ECM) supplemented with 5% FBS, 1% endothelial cell growth supplement (ECGS), and 1% penicillin/streptomycin. All cells were maintained in a humidified atmosphere containing 5% CO₂ at 37°C.

Total RNA was extracted from the cells by using RNeasy Mini Kit (Qiagen, Germany). Complementary DNA (cDNA) synthesis was performed using the PrimeScript™ RT Master Mix (Takara). The primers were synthesized by Takara (Tokyo, Japan). The β-actin gene was used as an endogenous control. Quantitative polymerase chain reaction (PCR) was performed using the TB Green® Premix Ex Taq™ II (Takara) under standard conditions according to the manufacturer’s instructions and conducted with the LightCycler 480 System (Roche, Switzerland). The data were analysed using the 2^−△△Ct method. The primers are listed in Supplementary Table S1.

In Figure 1A a total of 374 samples were included, including 373 patients with follow-up information, so a survival analysis was performed using 373 samples. The gene expression profiles of 373 liver cancer samples with 19,406 mRNA and 13,507 lncRNA and the gene expression profiles of 50 normal samples with 18,905 mRNAs and 11,770 lncRNA were obtained.

We identified 13,810 significant lncRNA–immune function pathway pairs in HCC samples, and there were 3,732 lncRNA regulators enriched to different immune function sets. Meanwhile, 14,129 significant lncRNA–immune function pathway pairs were detected in the normal samples; here, 5,322 lncRNAs were found to be enriched to different immune functions. Furthermore, a comparison of the immune function of lncRNA regulators in normal and HCC samples showed the enrichment relationships of 3,134 lncRNA–immune pathways in both diseased and normal samples; specifically, HCC and normal samples account for 22.7% and 22.04% of the total lncRNA–immune pathway pairs, respectively (Figure 2). In addition, 1,748 lncRNAs showed immunoregulatory functions in both cancer and normal samples, accounting for 46.84% and 32.84% of all the immunoregulatory functions, respectively, and a total of 1,984 immune-regulating lncRNAs were only found in HCC (Figure 2). These results indicated that the immune function of lncRNAs was significantly changed during cancer progression; therefore, we further focused on analyzing those lncRNAs with immune functions only in HCC samples.

After comparing the quantity of lncRNAs enriched in different immune pathways in normal and cancer samples, we found that more lncRNA regulators were enriched in cytokines, cytokine receptors, interleukins, and interleukins receptors in HCC samples (Figure 3). In addition, a higher proportion of HCC-specific immune function lncRNAs such as cytokines and interferons were involved in the regulation of immune-related molecular pathways (Figure 3).

As most genes with a key role in cancer development show dysregulated expression, we performed a differential expression analysis on the count expression profile of liver cancer from TCGA. Under the threshold of FDR < 0.05, we identified 5,244 significantly differentially expressed lncRNAs (4,218 highly expressed lncRNAs and 1,026 significantly lowly expressed) in HCC, accounting for 39.3% of all lncRNAs. Next, the intersection of these 5,244 significantly dysregulated lncRNAs and 1,984 lncRNAs specific to immune function of HCC showed a total of 498 lncRNA regulators with both characteristics, that is, significant expression differences in liver cancer samples and immune function specific to HCC. These 498 candidate lncRNA regulators included 300 significantly up-regulated lncRNAs and 198 significantly down-regulated lncRNAs in HCC (Figure 4A). Six of these lncRNAs were significantly dysregulated in liver cancer (Figure 4B).

We found that 498 differentially expressed immune function lncRNAs regulated a total of 17 immune molecular pathways, involving 1,482 lncRNA–immune function pathway relationship pairs, with an average of 93 candidate lncRNAs enriched in each pathway. Among them, more lncRNAs were enriched in cytokines, cytokine receptors, and TCR signaling pathway. Compared with liver cancer-specific immune function lncRNAs, dysregulated immune lncRNAs were more involved in cytokines, cytokine receptors, TGFb family member, TNF family member, and regulation of TNF family member receptor pathway (Figure 5). Cytokines mediate the interaction between cells and have a variety of biological functions, such as regulating cell growth, differentiation, and maturation; function maintenance; regulating immune response; participating in inflammatory response and wound healing, and a wide range of anti-tumor activities. TGFb family member enhances tumor migration and inhibits immune cell function by influencing the tumor microenvironment. Therefore, these candidate immunocompromised lncRNAs may play an important role in tumor invasion and metastasis.

To further examine immune regulation of candidate lncRNAs in liver cancer, we collected marker genes of 24 immune cells. Among them, effector memory cells contained 32 marker genes, with an average of 19 marker genes in each immune cell. We obtained a total of 538 significant lncRNA–immune cell relationship pairs, of which lncRNA LINC02723 and LINC02416 showed 22 related immune cells. Furthermore, 18 lncRNAs with a significant enrichment were considered as HCC-specific deregulated immune lncRNA collections and had at least 10 immune cells; moreover, these 18 lncRNAs were significantly enriched in 24 immune cells (Figure 6A). We found that 18 dysregulated immune lncRNAs were involved in lncRNA–immune function pathways in a total of 120 HCC samples and that T cell antigen receptor (TCR) signaling pathway, B cell antigen receptor (BCR) signaling pathway, cytokines, and cytokine receptors as well as other pathways had more lncRNA regulators. LncRNA HOXb-AS1 and LINC02550 were highly expressed in HCC and were enriched in the cytokines and cytokine_sensors molecular pathways (Figure 6B). Among the 18 lncRNAs, 12 lncRNAs were highly expressed in HCC samples, and six lncRNAs were lowly expressed in HCC samples (Figure 6C).

The classification of molecular subclasses of cancer samples is of great significance for personalized treatment of HCC; in this study, the samples were classified based on the expression of 18 HCC-specific dysregulated immune lncRNAs. Using the R package ConsensusClusterPlus, HCC samples in TCGA were well classified into two categories (group 1 and group 2), which contained 188 and 185 samples, respectively (Figure 7A). In addition, we integrated follow-up prognostic data, including 373 patients for prognostic analysis, and we found a significant survival difference between the two groups (p = 0.014), in which HCC samples in group 1 showed a higher survival rate and a longer median survival time (Figure 7B).

To study the differences in immunological characteristics of these two types of samples, we obtained the scores of important immune characteristics such as intratumor heterogeneity, proliferation, wound healing, and IFN-gamma response of TCGA HCC samples. The data demonstrated that patients in group 1 had significantly higher TGF-beta response, IFN-gamma response, and wound healing scores than patients in group 2 (Figure 8A). Moreover, group 1 patients showed significantly higher macrophage and monocyte scores than group 2 patients (Figure 8B). Compared with group 1 patients, group 2 patients had significantly higher proliferation levels, neutrophil distribution, and intratumoral heterogeneity (Figure 8C). In addition, the distribution of initial B cell, leukocyte ratio, Th1, Th17, and stromal ratios were all higher in HCC patients in group 1 (Supplementary Figure S1). These results indicated that group 1 patients had a better cytokine response and a higher percentage of immune cells, which may contribute to a better prognosis. As the samples in group 2 had higher proliferation and intratumoral heterogeneity and more neutrophil-promoting cancer metastasis, they tended to develop a poor survival.

To verify the differences of lncRNA in different subtypes, we screened the six most relevant lncRNAs from 18 immune-related lncRNAs for experimental verification. Specifically, the classification performance of each lncRNA on two molecular subtypes was calculated for 18 lncRNAs. The lncRNA with an AUC greater than 0.6 was selected, and the optimal AIC was selected by stepwise regression analysis to obtain the optimal model. Finally, the most critical lncRNAs of the six subtypes were identified (AC002480.1, MIR503HG, AC012368.1, LINC01357, HOXB-AS1, and LINC02416). A differential expression analysis showed that AC002480.1, AC012368.1, HOXB-AS1, and LINC02416 were highly expressed, while MIR503HG and LINC01357 were lowly expressed in tumors (Supplementary Figure S2A). The expressions of six lncRNAs in HCC cell lines were analyzed by RT-qPCR. As expected, AC002480.1, AC012368.1, HOXB-AS1, and LINC02416 were highly expressed, while MIR503HG and LINC01357 were lowly expressed in HCC cell lines (Supplementary Figure S2B). These results supported the conclusion of data analysis.