Landscape of associations between long non‐coding RNAs and infiltrating immune cells in liver hepatocellular carcinoma

Abstract Liver hepatocellular carcinoma (LIHC) can be detected by the immune system; however, it acquires features for evasion of immune surveillance during its origin and development. Long non‐coding RNAs (lncRNAs) are critical as immune regulators in cancers; nevertheless, the biological functions and mechanisms of lncRNAs in evasion of immune system by LIHC remain unclear. In this study, an integrated and computational approach was developed to identify immune‐related lncRNAs and to divide LIHC patients into diverse immune‐related risk groups based on the expression profiles of coding genes and lncRNAs. LIHC‐specific genes and lncRNAs in 17 immune cell populations were identified and analysed. Gene and lncRNA co‐expressing networks for diverse immune cell populations were constructed and analysed. Some imported immune‐related lncRNAs, such as MIR9‐3HG, were also identified. The LIHC patients comprised three different groups based on immune‐related risk. LIHC patients possessing a greater diversity of immune cell populations had better survival prognosis. The collective data are evidence of a credible computational model that can prioritize novel immune‐related lncRNAs and depict the atlas of immune‐related lncRNAs in LIHC. These findings will further the understanding of lncRNA function and advance the identification of immunotherapy targets in LIHC.

resistance and disease relapse after treatment are still key problems that hinder LIHC therapy. There is an urgent need to improve LIHC diagnosis and treatment.
Long non-coding RNAs (lncRNAs) are a subclass of functional non-coding RNAs. LncRNAs are over 200 nucleotides in size and do not encode protein. 7,8 However, lncRNAs are reportedly essential in many diseases, including cancers. 9 As well, lncRNAs have been associated with the regulation of biological functions of LIHC. For example, the up-regulation of lncRNA DLX6-AS1 can promote proliferation and tumour formation. 10  The tumour microenvironment is considered a key risk factor for cancer development and progression. Immune system imbalance and disorders are major causes of cancer. 13 Tumours can be detected by the immune system, and they acquire features to evade immune surveillance throughout the origin and development of cancer. Immunotherapy has emerged as a promising cancer treatment strategy. Dysregulation of immune genes could directly influence the immune system. Recent increasing evidence has indicated the participation of lncRNAs in the immune system. 14,15 For example, inducible degradation of the Sros1 ln-cRNA promotes interferon-gamma-mediated activation of innate immune responses by stabilizing Stat1 mRNA. 16 The HOTTIP ln-cRNA enhances interleukin 6 expression to potentiate immune escape of ovarian cancer cells by up-regulating the expression of programmed death-ligand 1 in neutrophils. 17 However, data concerning the diversity and function of immune-related lncRNAs in cancer are insufficient. Further studies of lncRNAs and their roles in immune regulation are essential to identify immunotherapy targets in LIHC.
In the present study, an integrated computational pipeline was developed based on gene and lncRNA expression profiles to identify immune-related lncRNAs in LIHC. LIHC patients were divided into three immune-related risk groups. We constructed 17 sets of immune cell populations and identified differentially expressed genes in these sets. Diverse immune cell-related genes and lncRNA co-expressed networks were constructed and analysed. Some imported immune-related lncRNAs, such as MIR9-3HG, were also identified. LIHC patients harbouring more diverse immune cell populations displayed a better survival prognosis. The data will provide a valuable resource for exploration of lncRNA function and identification of immunotherapy targets in LIHC.
Clinical follow-up information was retained for further analysis.
Genes and lncRNAs with expression values of 0 in the samples were removed. All expression values were log 2 transformed to satisfy normal distribution. The summary of characters and clinical followup information for LIHC patients and normal control was listed in Table S1.

| Collection of multiple immune cell populations-related (ICPR) genes
Seventeen immune cell populations were selected. They included B cells, eosinophils, macrophages, mast cells, natural killer (NK)

CD56bright cells, NK CD56dim cells, neutrophils, T helper cells, Tcm cells, Tem cells, Tfh cells, activated dendritic cells (aDCs), induced dendritic cells (iDCs), activated CD8 T cells, gamma delta T cells,
regulatory T cells (Tregs) and cytotoxic cells. 18 Information on multiple ICPR genes was also obtained from previous studies. 19,20 All the ICPR genes were analysed.

| Identification of LIHC-specific ICPR genes and lncRNAs based on expression profiles
In order to identify LIHC-specific ICPR genes and lncRNAs, the t test was used for differential expression analyses of the expression profiles of LIHC and normal samples. Significantly expressed ICPR genes (P < .05) were considered LIHC-specific ICPR genes.
Significant lncRNAs (P < .01) were considered LIHC-specific lncR-NAs. Only the LIHC-specific ICPR genes and lncRNAs were used for subsequent analyses. with absolute values of difference >0.5 were used for subsequent analyses and network construction. The significantly changed coexpressed networks of LIHC-specific ICPR genes and lncRNAs were constructed using Cytoscape 3.3.0 (https://cytos cape.org/). Degree analysis was also performed using Cytoscape 3.3.0.

| Evaluating immune-associated risks for LIHC patients based on co-expressed networks
An integrated and computational pipeline was developed to evaluate immune-associated risks for each LIHC patient based on co-expressed networks ( Figure 1). Firstly, the expression profiles of genes and lncRNAs in co-expressed network for each immune F I G U R E 1 Workflow of evaluating the immune risk for LIHC patients based on ICPR genes and lncRNAs.
Step 1: Identification of LIHCspecific ICPR genes and lncRNAs based on expression profiles. Step 2: Construction of LIHC-specific ICPR genes and lncRNAs co-expression networks for diverse kinds of immune cell populations.
Step 3: Evaluating immune risk for LIHC patients and dividing the patients into groups cell population were extracted from the whole expression profile.
The differentially expressed genes and lncRNAs were allocated to up-regulated and down-regulated groups. Secondly, the expression values were ranked from high-to-low for all LIHC patients.
For up-regulated genes and lncRNAs, the genes or lncRNA of a specific sample ranked last 70% were considered risk genes or lncRNAs for the LIHC patient. For down-regulated genes and lncRNAs, the genes or lncRNA of a specific sample ranked before 30% were considered risk genes or lncRNAs for the LIHC patient.
Thirdly, for a specific LIHC patient, the risk genes in each kind of ICPR co-expression network were counted. If there were more than 50% risk genes or lncRNAs for an immune cell population,

| Survival analyses for multiple kinds of immune risk-related groups
In order to evaluate the prognosis of the three immune risk-related groups, Kaplan-Meier (K-M) survival analyses were performed.
Statistical significance assessed using the log-rank test. All analyses were performed within the R 2.6.6 framework.

| Some ICPR genes and lncRNAs were specific in LIHC
We obtained ICPR gene information for the 17  T cells, and 19 up-and nine down-regulated ICPR genes in neutrophils. Differentially expressed lncRNAs, considered as LIHCspecific lncRNAs, were also identified (56.76%) ( Figure 2C). Similar to LIHC-specific ICPR genes, most of the differentially expressed lncRNAs were down-regulated ( Figure 2D). These results revealed that immune-related genes and some lncRNAs may participate in the development of LIHC.

| Significant change in co-expression between LIHC-specific ICPR genes and lncRNAs in LIHC compared with normal samples
We used the absolute difference of PCCs for ICPR genes and lncR-

| Particular features of co-expressed networks of LIHC-specific ICPR genes and lncRNAs
In order to better describe the role of immune changes in LIHC, more significant changes (absolute difference value >0.7) for LIHCspecific ICPR gene and lncRNA PCCs were used to construct networks. The numbers of LIHC-specific ICPR gene and lncRNA pairs in each immune cell population were determined ( Figure 4A). iDCs and B cells displayed the highest number of LIHC-specific ICPR gene and lncRNA pairs. A comprehensive network of LIHC-specific ICPR gene and lncRNA pairs was constructed by integrating the multiple immune cell populations ( Figure 4B). The integrated network revealed obvious differences between the diverse LIHC-specific ICPR gene and lncRNA pairs. The same ICPR gene could also interact with different lncRNAs with diverse strengths and changing patterns. Only a few ICPR genes were prominent. These may be important in LIHC.
A similar network was also constructed in single immune cell population. For example, a co-expressed network of LIHC-specific ICPR genes and lncRNAs was constructed in B cells ( Figure 4C). This network contained 141 pairs of LIHC-specific ICPR genes and lncRNAs, 18 ICPR genes and 138 lncRNAs. The integrated degree distributions of the network displayed a scale-free network power-law distribution (R 2 = .83, Figure 4D). This scale-free distribution indicated that this network was a meaningful biological network.

| Identification of important ICPR genes and corresponding lncRNAs in LIHC
Some ICPR genes and corresponding lncRNAs were more frequently present in the aforementioned co-expression networks of LIHCspecific ICPR genes and lncRNAs. Thus, we inferred that at least some of these frequent ICPR genes and corresponding lncRNAs may play important roles in LIHC. These ICPR genes included ABCB4, CD101, CD1A, CD1B and FABP4 ( Figure 5A). ABCB4 is an immunerelated gene as well as a phospholipid translocator at the canalicular membrane of the hepatocyte, which "flops" phosphatidylcholine into bile. 21 Figure 5F). These collective results revealed essential roles of some lncRNAs in the regulation of the immune system of LIHC.

| Different survival outcomes in LIHC for the diverse immune risk groups
Evaluating the immune risk of each LIHC patient could be beneficial in clarifying patient status and for personalized treatment. We inferred that LIHC patients harbouring more diverse immune cell populations could be considered higher immune risk patients. The LIHC patients differed in their diversity of immune cell populations ( Figure 6A). For example, 81.57% and 68.29% of the LIHC patients were associated with Treg and NK CD56dim cells, respectively. However, only 9.21% of the LIHC patients were related to neutrophils. We also analysed the numbers of immune cell populations for each LIHC patient. Forty-seven LIHC patients harboured two kinds of immune cells populations and 43 harboured one kind ( Figure 6B).  Figure 6C). The data were used to divided the LIHC patients into three immune risk groups of multi-immune (n = 64), mediaimmune (n = 110) and min-immune groups (n = 189; Figure 6D).
Survival increased as the immune risk decreased ( Figure 6E). The survival analysis also revealed that diverse immune risk groups were significantly associated with survival (P = .003, Figure 6F). More immune risk types for a LIHC patient implied stronger immune participation in the specific sample. Thus, a specific LIHC patient would have a better prognosis with the alteration of multiple immune cell populations in LIHC. Liver hepatocellular carcinoma is highly heterogeneous immunologically, and immune-related genes have been identified and studied. 24,25 Most of these data involve coding genes. Some important non-coding RNAs, including lncRNAs, have been ignored.

| D ISCUSS I ON
LncRNAs are emerging as critical regulators of gene expression and play fundamental roles in immune regulation for many kinds of cancers. 15 In the present study, we integrated gene and lncRNA expression in order to globally describe immune characteristics for LIHC patients. We discovered that some of the identified immune-related genes and lncRNAs were associated with immune functions and cancer in previous experiments (Table S2). Some key immune-related lncRNAs, such as MIR9-3HG, SYP-AS1 and LINC02542, were identified. LINC01564 could interact with 12 ICPR genes and was considered a cancer-related lncRNA. 26 More significantly, LINC01564 was also successfully detected in patient urine samples and was up-regulated when compared with urine from normal (healthy) individuals. These immune-related lncRNAs were co-expressed with many ICPR genes in LIHC. They may regulate ICPR gene expression or may work with ICPR genes together in LIHC. Presently, a co-expression network was used to identify immune-related lncRNAs. Weighted correlation network analysis (WGCNA) is also an effective method to determine relationships between genes and lncRNAs. WGCNA retains the continuity of network node connectivity and can identify core modules.
However, different data preprocessing and analysis parameters will also lead to different results for WGCNA. Thus, appreciable differences can be evident for selected core modules. With a small number of known immune genes, comparing to WGCNA the co-expressed network may contain more global information.
In future studies, we intend to investigate the differences in the co-expression network and WGCNA approaches in the identification of immune-related lncRNAs.
The tumour microenvironment can be classified into immunologically active "inflamed" tumours and inactive "non-inflamed" tumours based on the infiltration of cytotoxic immune cells. 24  The present study provides a standardized computational pipeline to identify immune-related lncRNAs in LIHC. It also provides a data resource concerning candidate immune-related lncRNAs for further analyses in LIHC. Future studies will focus on investigating more LIHC samples to validate the accuracy and stability of the method presented here. Experimental validation of our results is also necessary.
In summary, lncRNAs are an important layer in the immune complexity in LIHC. Dividing LIHC patients into diverse immune risk groups would aid in evaluating prognosis. Our study provides an insight for the identification of novel immune-related biomarkers and immunotherapies for LIHC.

ACK N OWLED G EM ENTS
Not applicable.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

E TH I C S A PPROVA L A N D CO N S E NT TO PA RTI CI PATE
Not applicable.

PATI E NT CO N S E NT FO R PU B LI C ATI O N
Not applicable.

CO N S E NT FO R PU B LI C ATI O N
Not applicable.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.