HKCI-2, 4, 9, 10, 11, C1, C2, and C3 cell lines were maintained and cultured as described previously. HKCI-2 and 10 cell lines were derived from HCC patients of NASH etiology. HKCI-4,9 and 11 were derived from HCC patients of HBV etiology. HKCI-C1, C2, and C3 were derived from HCC patients of HCV etiology (Gho et al., 2008).

The subcellular fractionation methodology was established previously by Djebali et al. in a study of the transcriptional landscape of ENCODE reference human cell lines (Djebali et al., 2012). We used a similar stepwise lysis protocol to generate both cytoplasmic and nuclear fractions from each sample. Briefly, the tissues were initially processed by disrupting the outer cellular membrane to release the cytoplasmic contents. The lysates were subjected to high-speed centrifugation to pellet the intact nuclei (nuclear fraction) from the supernatant (cytoplasmic fraction). A RNeasy MiniElute Cleanup kit (Qiagen) was used to extract the total RNA from the respective fractions. To validate the fractionation protocol, qPCR was conducted to quantitate the relative enrichments of cytoplasmic marker genes (RPS14, GAPDH) and nuclear marker genes (MALAT1) in a control fractionation experiment. qPCR primers were designed to cover exon-exon junctions and to capture only spliced transcripts ( Supplementary Table 1 ). Results confirmed >10-fold enrichments of cytoplasmic and nuclear markers in respective fractions, suggesting effective separation of the two subcellular fractions ( Supplementary Figure 1 ).

Total RNA samples were converted into strand-specific, rRNA-depleted libraries, and pair-end sequenced on an Illumina HiSeq 2500 platform by Macrogen Co. (Seoul, South Korea). A total of 7.04 billion reads were generated for eight in-house HCC cell lines. The numbers of clean reads after trimming obtained from cytoplasmic RNA ranged from 59.8 to 99.2 million across eight HCC cell lines and the effective sequencing depth was estimated between 21.3x and 31.6x. The final number of clean reads obtained from nuclear RNA ranged from 95.9 to 171.2 million with the effective sequencing depths between 5.5x and 12.84x. To confirm the outcomes of subcellular fractionation, the expression levels of cytoplasmic marker genes (RPS14, GAPDH) and nuclear marker genes (MALAT1, NEAT1, PVT1) (Chan and Tay, 2018; Sun et al., 2018; Yu et al., 2018) were quantitated from RNA-seq data using featureCounts (version 1.6.3) (Liao et al., 2014). The RNA-seq expression fold changes of RPS14, GAPDH, and MALAT1 were concordant with the qPCR results from the control fractionation experiment. In addition, the RNA-seq expression of NEAT1 and PVT1 were enriched in the nuclear fraction by >30-fold and >10-fold respectively. The observation indicated effective separation of the subcellular fractions among all HCC cell line samples ( Supplementary Tables 2A–C ).

RNA-seq datasets (in FASTQ format) and corresponding clinical data from patients in the TCGA-LIHC cohort were obtained from TCGA following authorization (dbGaP controlled dataset phs000178.v10.p8) (Weinstein et al., 2013). A total of 371 tumor tissue RNA-seq datasets and 50 paired (tumor tissue and tumor-adjacent tissue) RNA-seq datasets were included in subsequent analyses. Among the 50 patients with paired RNA-seq datasets, 7, 5, and 3 of them are with HBV, HCV, and NASH etiology, respectively.

RNA-seq datasets from the eight HCC cell line samples and the 50 paired TCGA-LIHC samples were first trimmed using trimmomatic (Bolger et al., 2014) (version 0.36). For each individual sample:

Finally, the transcript merge mode of StringTie was used to unify and merge the sample-specific transcriptomes into a high-consensus transcriptome. The high-consensus transcriptome was composed of transcripts assembled in both HCC cell line samples and TCGA samples.

The assembly support of a transcript was defined as the number of sample-specific transcriptome derived from the 50 paired TCGA samples that included the transcript. The assembly support of transcripts in tumor and tumor-adjacent tissues was counted separately.

The expression level of each transcript was quantified using Kallisto (Bray et al., 2016) (version 0.43.0). EBSeq (Leng et al., 2013), DESeq2 (Love et al., 2014), and edgeR (McCarthy et al., 2012) were used to identify transcripts that were differentially expressed between TCGA LIHC tumor and tumor-adjacent tissue samples. In DESeq2 and edgeR, differentially expressed transcripts were defined as those that satisfied 2 criteria: |log₂(fold-change)| > 1 and p < 0.01 after the Benjamini–Hochberg correction. In EBSeq, differentially expressed transcripts were defined as those with a PPDE (posterior probability that a transcript is differentially expressed) > 0.99 and |log₂(fold-change)| > 1. Dysregulated transcripts related to HCC progression were defined as the consensus subset of differentially expressed transcripts among the three methods.

From the merged assembly assemblies, transcripts in non-repeat-masked genomic regions with ≥2 exons, a length >200 nucleotides, available strand information, and an expression level of TPM > 1 were selected. This subset of transcripts was compared with the GENCODE (v27) annotation, and gffcompare was used to assign an annotated or unannotated status. Annotated transcripts classified as “lincRNA” or “antisense” in GENCODE (v27) were considered lncRNAs. Unannotated transcripts were considered lncRNAs if both CPC2 (Kang et al., 2017) and COME (Hu et al., 2017) predicted their coding potential as “noncoding RNA”.

TCGA LIHC patients with clinical survival data were classified into the high-risk or low-risk group based on lncRNA transcript expression. A Kaplan-Meier survival analysis and log-rank test were used to estimate differences in the overall survival times between patients in the 2 groups. All analyses were conducted on the R-3.4.1 framework.

Processed H3K4me3 ChIP-seq peak calling data from 28 cell lines in 4 categories were obtained from the Human Epigenome Atlas repository (Roadmap Epigenomics Consortium et al., 2015) (Release 9).

LncRNA transcripts and genes were considered associated with an active promoter mark if an H3K4me3 ChIP-seq peak was present at ±1000 bp of the 5´ end.

Single-cell RNA-seq datasets (in FASTQ format) from human preimplantation embryonic cells in seven stages (oocyte, zygote, 2-cell, 4-cell, 8-cell, morula, and late blastocyst) were downloaded from NCBI SRA (accession: SRP011546) (Yan et al., 2013). Transcript assembly support was calculated as described above. Transcript expression levels were quantified using Kallisto (Bray et al., 2016). LncRNAs were filtered away if no assembly support was detected for the lncRNA transcript and its isoforms in all preimplantation embryonic cell datasets. LncRNAs were considered to be expressed in a preimplantation embryonic stage if a TPM measurement >1 was detected in any cell in ≥2 embryo samples. LncRNAs were further designated as possibly lineage-specific if a TPM >1 was observed in ≥2 cells but <67% of all cells in an embryo sample.

Differentially expressed, tumor-enriched protein-coding genes with a H3K4me3 active promoter histone mark in ≥4 (of 8) ESC samples were selected and tested for enrichment against a background of all annotated protein-coding genes using GOATOOLS (Klopfenstein et al., 2018). P values were corrected for multiple testing using the Benjamini–Hochberg procedure at an alpha = 0.05.

The FPKM values of each gene in cytoplasmic and nuclear fractions from the eight HCC cell lines were computed using featureCounts (version 1.6.3) (Liao et al., 2014). DESeq2 was used to conduct the in-sample normalization of FPKM values (Love et al., 2014). For each cell line sample, only genes with an expression of ≥0.1 FPKM in both the cytoplasmic and nuclear fractions were retained when computing the log₂(FPKM_{cytoplasmic fraction}/FPKM_{nuclear fraction}) metric [abbreviated as log₂(FPKM_cyto/FPKM_nuc)]. A positive log2(FPKM_cyto/FPKM_nuc) value indicated a transcript abundance bias towards the cytoplasmic fraction, while a negative value indicated a bias towards the nuclear fraction.

The list of human housekeeping genes was obtained from a previous study (Eisenberg and Levanon, 2013). The 5^th and 95^th percentiles of the mean log₂(FPKM_cyto/FPKM_nuc) values of these housekeeping genes were selected as the lower and upper reference limits, respectively, for non-fraction-specific genes. Genes with log₂(FPKM_cyto/FPKM_nuc) measurements below the lower limit in ≥4 cell lines were considered nuclear-enriched, while those with measurements exceeding the upper limit in ≥4 cell lines were considered cytoplasmic-enriched.

The mean log₂(FPKM_cyto/FPKM_nuc) values for the cytoplasmic marker genes GAPDH and RPS14 were 3.860 and 3.306, and mean value for the nuclear marker gene MALAT1 was -4.132. Our proposed metric cutoff scheme suggested cytoplasmic-enrichment status of GAPDH and RPS14 and the nuclear-enrichment status of MALAT1. The findings from fractionation-then-sequencing data agreed with the results of prior qPCR experiments.

The co-expression analysis was performed using the Weighted Correlation Network Analysis (WGCNA) R package (Langfelder and Horvath, 2008) and was based on expression data from all protein-coding and lncRNA genes (including unannotated lncRNAs assembled in this study) quantified from 421 tumor datasets in the TCGA-LIHC cohort using Kallisto (Bray et al., 2016). Genes expressed in <20% of datasets were excluded from the analysis. Subsequently, pairwise Pearson’s correlation coefficients (PCCs) were calculated between 26 nuclear-enriched HCC-associated lncRNA genes and all other genes. Genes exhibiting significant co-expression with nuclear-enriched HCC-associated lncRNA genes (PPC ≥ 0.6, p < 0.05) were selected for a hierarchical clustering analysis using Python.

Initially, we conducted a transcriptome assembly to detect lncRNAs in HCC samples ( Figure 1A ). The resolution of RNA-seq data from TCGA is not sufficient to decode clear strands for assembled transcripts because of a lack of strand information. We overcame this limitation by integrating a set of strand-specific RNA-seq data from in-house HCC cell lines into our pipeline. A total of 216 unique transcriptome assemblies were produced from 108 datasets (8 from in-house HCC cell lines, 100 from 50 pairs of TCGA-LIHC tumor/tumor-adjacent samples). The reference-based and ab initio transcriptome assemblies were incorporated using StringTie (Pertea et al., 2015) and Trinity (Grabherr et al., 2011), respectively. To generate a high-consensus transcriptome, 13, 835 multi-exonic transcripts (including 10, 591 genes) with definitive strand information identifiable by both StringTie and Trinity in both types of datasets were extracted for downstream analyses. These datasets included 506 lncRNA transcripts (444 genes) and 11,084 protein-coding transcripts (9, 038 genes) annotated by GENCODE (v27). Using the reference lncRNA annotation from GENCODE, we identified 450 high-quality unannotated lncRNA transcripts (411 genes) with support from CPC2 (Kang et al., 2017) and COME (Hu et al., 2017) ( Figure 1B ). These novel lncRNA transcripts were first divided into 4 categories based on the “transfrag class codes” defined by gffcompare ( Supplementary Table 3 ). Specifically, > 50% of lncRNAs were assigned a class code “j” and were identified as potential isoforms that shared at least 1 splice junction with reference transcripts. Furthermore, 24.4% and 10.2% of lncRNAs were assigned class codes “u” and “i” and were classified as intergenic and intronic transcripts, respectively. The remaining 15% of lncRNAs were assigned class codes of “x”, “o”, and “c,” indicating lncRNAs that overlapped with currently annotated transcripts in either strand. Through a comparison with other noncoding annotations, we observed that approximately half of our newly assembled lncRNAs were not previously reported. For example, the MiTranscriptome (Iyer et al., 2015) and NONCODE (v5) (Fang et al., 2018) annotations yielded overlap rates of 55.5% and 41%, respectively, with the novel lncRNAs in our study ( Supplementary Figure 2 ). Accordingly, our RNA-seq data seemed to have a sufficient sequencing depth for capturing novel transcripts.

Using the assembled transcripts from our annotation, we revisited the basic characteristics of lncRNAs and protein-coding transcripts. Previous studies mentioned that lncRNAs have a shorter length, lower expression, and smaller exon number than protein-coding transcripts. To address whether the newly assembled lncRNAs also exhibited these features, we compared 450 unannotated lncRNAs with the protein-coding transcripts and annotated lncRNAs in our final transcriptome annotation. We observed that our newly assembled lncRNA transcripts were significantly shorter in length than protein-coding transcripts (p < 0.001, Student’s t-test), but were comparable to the lncRNAs annotated in GENCODE ( Figure 1C ). We observed similar results in comparisons of exon numbers and expression ( Figures 1D–F ). Taken together, these results suggest that the sequencing data and lncRNA identification scheme applied in this study allowed us to assemble highly credible lncRNAs with considerable lengths and expression levels.

To identify lncRNAs with broad involvement in the prognosis of HCC, we conducted differential expression tests of tumor and tumor-adjacent samples in TCGA using the EBSeq (Leng et al., 2013), DESeq2 (Love et al., 2014) and edgeR (McCarthy et al., 2012) methods. Notably, these methods yielded comparable results. Overall, 87 and 117 deregulated lncRNAs predicted by all 3 methods were identified in the annotated and unannotated lncRNA sets, respectively ( Figure 2A ). To validate the outcomes of a differential expression analysis, all 421 RNA-seq datasets from the 371 patients in the TCGA-LIHC cohort were clustered according to the expression profiles of all dysregulated lncRNAs. We observed that the upregulated and downregulated lncRNAs could be used effectively to separate tumor and nontumor samples (p < 0.001, chi-square test) ( Figure 2B ), indicating that the differentially expressed lncRNA set detected in this study was reliable.

Of the 204 deregulated lncRNAs, 105 unannotated and 76 annotated transcripts were upregulated in tumor samples. Nevertheless, low but detectable levels of expression in the tumor-adjacent samples ( Figures 2C, D ) suggest that the biological functions of some lncRNAs might not be specific to cancer-related mechanisms. Alternatively, as TCGA datasets are non-strand-specific, this observation might also be attributable to artifacts from transcript quantification. To evaluate the tumor specificity of lncRNA accurately, we revisited our transcriptome assembly and calculated the assembly support for each within the 50 pairs of TCGA-LIHC samples. As a successful assembly requires moderate and uniform sequence coverage across the entire transcript region, assembly support may be a better and more stringent indicator of the presence of transcripts. Here, 81/45 annotated/unannotated and tumor-upregulated lncRNAs were found to have no assembly support from tumor-adjacent samples. These lncRNAs were supported by up to 19 TCGA-LIHC tumor samples ( Figure 2E ), coherent with the general observation that lncRNA expression is highly heterogeneous across tumor samples. Results of differential expression analysis and the assembly support of these lncRNAs were listed in Supplementary Table 4 .

Next, the gene expression landscapes in the nuclear and cytoplasmic fractions were investigated using fractionation-then-sequencing data from in-house HCC cell lines. To evaluate the relative transcript abundances between the two fractions, the normalized gene expression (in FPKM) in the cytoplasmic fraction were divided by the expression in the nuclear fraction to obtain a log₂(FPKM_cyto/FPKM_nuc) value for each gene. A positive log₂(FPKM_cyto/FPKM_nuc) value indicated a transcript abundance bias toward the cytoplasmic fraction, while a negative value indicated a bias towards the nuclear fraction.

The overall mRNA population was not apparently biased toward either fraction, with a median log₂(FPKM_cyto/FPKM_nuc) close to 0. In contrast, the overall lncRNA population was significantly more biased toward the nuclear fraction (p < 0.001, Mann-Whitney U test), with a negative median log₂(FPKM_cyto/FPKM_nuc). Moreover, the HCC-associated lncRNAs identified in this study also resembled the characteristics of nuclear bias from the overall lncRNA population ( Figure 3A ). These observations were also consistent when the measurements of the eight HCC cell lines were averaged ( Figure 3B ) and were consistent with findings from previous fractionation-then-sequencing studies of other cell lines, such as HepG2 (Benoit Bouvrette et al., 2018).

Nevertheless, fraction bias in individual transcripts might be due to temporal fluctuations in cellular transcript levels (Blake et al., 2003; Kaern et al., 2005; Raj and van Oudenaarden, 2008; Eldar and Elowitz, 2010) or due to mechanisms such as burst transcription (Dar et al., 2012; Bahar Halpern et al., 2015a; Bahar Halpern et al., 2015b), rather than genuine fraction-enrichment events. To select the most confident fraction-enriched subset of genes, a list of 3,804 RNA-seq-derived human housekeeping genes was downloaded (Eisenberg and Levanon, 2013). Of these, 3,543 (93.1%) are mRNA genes and consistently expressed (FPKM ≥ 0.1) in both the nuclear and cytoplasmic fractions of all eight HCC cell lines. The mean log₂(FPKM_cyto/FPKM_nuc) values of these housekeeping genes were then assumed a reference range for normal variations in fraction bias. This implied that only genes with a substantially negative or positive log₂(FPKM_cyto/FPKM_nuc) values would be considered fraction-enriched. Based on this assumption, the 5^th and 95^th percentiles of the reference range (-1.278–3.056) were selected as the cutoffs for detecting nuclear-enriched and cytoplasmic-enriched genes. Asymmetry in the cutoff values was consistent with the general knowledge that mRNAs are predominantly enriched in the cytoplasmic fraction, which is also the site of protein translation ( Figure 3B ).

Based on our cutoff scheme, approximately 80% of the overall mRNA population was considered non-fraction-specific ( Figure 3C ). In contrast, >50% of the lncRNA population appeared to be fraction-enriched, with a strong bias toward enrichment in the nuclear fraction ( Figure 3D ). Of the 126 proposed HCC-associated lncRNA transcripts, 28 transcripts corresponding to 26 lncRNA genes ( Supplementary Table 5 ) exhibited stable nuclear-enrichment in ≥4 in-house HCC cell lines ( Figure 4A ). Remarkably, the list of genes included PXN-AS1, which was reported to interact with the nuclear-localized splicing factor MBNL3 to promote liver and lung cancer progression (Yuan et al., 2017). Although six of the proposed nuclear-enriched HCC-associated lncRNA genes were not annotated in the GENCODE reference, all had corresponding entries in the MiTranscriptome catalog, and MHCC.16327 and MHCC.17346 were also suggested to be dysregulated in liver cancer by MiTranscriptome study (Iyer et al., 2015) ( Supplementary Table 6 ).

The 26 nuclear-enriched HCC-associated lncRNAs were further characterized in a nonexclusive manner for histone marks on active promoters and correlations with patient survival and embryonic tissue expression. First, we investigated H3K4me3 histone marks on the candidate lncRNAs. H3K4me3 indicates the presence of active promoter structures (Mikkelsen et al., 2007) and is used widely to support lncRNA expression in cell or tissue samples (Guttman et al., 2009; Marques et al., 2013). We focused on H3K4me3 marks associated with stem cells, as some cancer cells may express lncRNAs ectopically to acquire a stem-like phenotype (Jiang et al., 2016). Strikingly, all GENCODE-annotated candidates possessed stem cell-associated H3K4me3 marks and seven, including PXN-AS1, were found to possess H3K4me3 marks in all surveyed ESCs, ESC-derived cells, and iPSCs ( Figure 4 ).

To test our hypothesis, we downloaded and evaluated a set of experimentally supported HCC prognostic lncRNA from the Lnc2Cancer database (Ning et al., 2016). Of the 10 reported HCC-prognostic lncRNAs, the majority possessed stem-cell associated H3K4me3 ( Figure 4 ). Moreover, a gene ontology (GO) analysis of 3, 295 tumor-enriched protein-coding genes with stem-cell associated H3K4me3 marks also revealed the significant enrichment of GO terms related to stem cell functions such as cell cycle, cell division, and cell adhesion ( Supplementary Figure 3 ).

Next, the candidate lncRNA expression patterns in human preimplantation embryonic cells were surveyed based on published single-cell RNA-seq datasets (Yan et al., 2013). This analysis was motivated by the observation that the majority of the 10 reported HCC prognostic lncRNAs were expressed and regulated along the course of preimplantation development ( Figure 4 ). This pattern suggests the possible superimposition of lncRNA functions on developmental regulation and cancer progression (Janic et al., 2010; Planells-Palop et al., 2017; Bruggeman et al., 2018). Excitingly, all nuclear-enriched HCC-associated lncRNAs were expressed in some preimplantation stages, consistent with reports of prognostic lncRNA. Moreover, these transcripts exhibited possible lineage specificity because they were expressed in only a fraction of embryonic cells at different stage(s), typically beginning at the 4-cell or 8-cell stage ( Figure 4 ).

Lastly, we conducted a survival analysis based on the clinical data of TCGA-LIHC patients. Twenty-three of 26 candidates significantly affected the survival outcomes of TCGA HCC patients (p < 0.05) and could serve as potential prognostic markers ( Figure 4 ). Overall, this characterization analysis suggested that the 26 nuclear-enriched HCC-associated lncRNAs identified in our bioinformatics pipeline shared the characteristics of experimentally supported HCC-prognostic lncRNAs and may also possess cancer-driving capabilities.

To connect HCC-associated lncRNAs with their potential biological functions, we conducted a co-expression analysis and identified 33 partner genes that were strongly co-expressed (PCC ≥ 0.6, p < 0.05) with some of the 26 nuclear-enriched HCC-associated lncRNA genes ( Figure 5 ). Of the 33 co-expression partner genes, more than half were suggested cancer prognostic markers (Uhlen et al., 2017), while 13 were identified as possible cancer drivers by either experimental studies or non-TCGA based bioinformatics studies. Subsequently, hierarchical clustering was conducted to group HCC-associated lncRNAs into three clusters based on their co-expression partners ( Figure 5 ):

Overall, our co-expression analysis established connections between some nuclear-enriched, HCC-associated lncRNAs, and reported cancer-driving protein-coding genes. Moreover, our findings enabled the prioritization of more promising lncRNA candidates for further experimental evaluations. The co-expression partners of these candidates provided clues regarding the cellular mechanisms that might involve the candidate lncRNAs.