Whole-genome DNA methylation data from spinal cord ependymal tumors were downloaded from Gene Expression Omnibus version 1 (24) [accession no. GSE65362 (4)]. Specifically, these samples were classified into three classes: i) Spinal ependymoma (SP-EPN); ii) spinal myxopapillary EPN (SP-MPE); and iii) spinal subependymoma (SP-SE). lncRNA annotation data were downloaded from GENCODE (version 28) (25). These lncRNAs were classified into nine groups based on their location and features (http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html) (Fig. 1A). Notably, the To be Experimentally Confirmed (TEC) category has been specifically created for the ENCODE project to indicate regions that could indicate the presence of protein coding genes (PCGs) that require experimental validation (25). Importantly, Ensembl (https://www.ensembl.org/info/genome/genebuild/biotypes.html) and Vega (http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html) databases classify TEC genes as lncRNAs.

The genomic locations of the Infinium Human Methylation 450 (HM450) BeadChip (Illumina, Inc.) probes based on GRCh37 were converted to the genomic locations in GRCh38 (The Genome Reference Consortium; http://www.ncbi.nlm.nih.gov/grc) using the University of California, Santa Cruz (UCSC) Batch Coordinate Conversion liftOver tool version 2 (26). Probes exhibiting single nucleotide polymorphisms (SNPs) >5 bp from their 3′-end and probes with non-unique 3′-subsequences of 30 bases were excluded, as previously described (27). lncRNA promoter regions were defined as 3-kb windows from either side of the transcription start site (TSS), as previously described (22). Subsequently, the methylation probes were mapped to the lncRNA promoter regions using BEDtools version 2.24.0 (28). The same approach was used to map the PCG probes. In addition, methylation probes that simultaneously mapped to lncRNAs and PCGs were excluded.

lncRNA promoter methylation values were calculated as the mean values of all probes located in the corresponding promoter. DNA methylation patterns around TSSs were calculated in 100 bp windows based on the median methylation value across samples. The lncRNA methylation value was used to cluster samples using the pheatmap package (version 1.0.12; http://cran.r-project.org/web/packages/pheatmap/index.html), using the default parameters, on R (version 3.3.0; http://cran.r-project.org/). The similarity between tumor samples was measured by Pearson correlation coefficients using R to determine if tumors within a subtype were more similar to each other than those from other subtypes. In addition, principal component analysis was used to investigate the methylation patterns of different tumor subgroups. Tumor subtype-specific lncRNAs were identified by ANOVA using R, as previously described (29). The RCircos package (version 1.2.1; http://cran.r-project.org/web/packages/RCircos/index.html) was used with R to display the distribution of methylation levels and the locations of tumor subtype-specific lncRNAs in the genome.

RNA-protein interaction data were downloaded from the RAID database (version 2.0) (30), which integrates experimental and computational prediction interactions from the literature and other database resources under one common framework. lncRNA-protein interaction associations were calculated by mapping subtype-specific lncRNAs to mRNAs.

A random forest algorithm (31) was performed using the DNA methylation value of the lncRNAs identified in the constructed network. The lncRNAs with a feature importance value >0, calculated using the random forest algorithm, were considered potential tumor subtype-specific signatures (32). All statistical analyses were performed using R (version 3.3.0; R project).

Tumor subtype-specific lncRNAs were annotated using their promoter regions by the Genomic Regions Enrichment of Annotations Tool (GREAT; version 3.3.0) using the default parameters (33); in addition, the lncRNA promoter regions were used as ‘background regions’. Specifically, the genomic coordinates of the lncRNAs used in GREAT were first converted into GRCh37 coordinates using the UCSC liftover tool, as aforementioned. Kyoto Encyclopedia of Genes and Genomes (KEGG; version 73.0) (34) enrichment analysis was performed on the PCGs using the Enrichr online tool (version 2.0), using the default parameters (35), and the significantly enriched KEGG pathways were identified (P<0.01).

To characterize lncRNA methylation patterns, a computational pipeline was used to annotate HM450 probes to lncRNA promoters. In total, 485,506 HM450 probes were successfully converted into GRCh38 coordinates, and 433,532 probes were obtained after filtering for SNPs and copy number-associated probes. The present analysis resulted in a set of 29,402 probes annotated in 6,967 lncRNA promoter regions (Fig. 1A). In total, 12,668 and 11,711 methylation probes were found in the promoters of long intergenic noncoding RNAs (lincRNAs) and antisense lncRNAs, respectively. Additionally, lncRNAs exhibited methylation states similar to those of PCGs around TSSs; both had low methylation values in proximity to the TSS and had higher methylation values as the distance from the TSS increased (Fig. 1B). However, the overall methylation level in lncRNAs was higher compared with PCGs (Fig. 1B). Cluster analysis suggested that tumor samples from patients with the same molecular subtype had similar methylation patterns in lncRNAs and that the methylation levels of various lncRNAs in the same class were not similar (Fig. 1C). Additionally, the majority of lncRNAs identified in the present study had a methylation value >0.4 (Fig. 1C). Certain lncRNAs, including those classified as ‘processed_transcript’ and ‘bidirectional_promoter_lncRNA’ lncRNAs, exhibited bimodally distributed methylation levels in tumor samples, whereas other lncRNAs, including ‘antisense’ and ‘sense_intronic’ lncRNAs, were hypermethylated in the majority of the patients (Fig. 1D).

The patterns of DNA methylation in the lncRNA genes were investigated in different tumor histopathological subtypes. In total, there were 57 patients with spinal cord tumors, 29 (50.9%) patients were in the SP-MPE group, whereas 21 (36.8%) and 7 (12.3%) were in the SP-EPN and SP-SE group, respectively. DNA methylation correlation values were calculated between each pair of samples, and tumor samples were clustered using their methylation correlation values. lncRNA methylation levels were found to be more similar within groups than between groups (Fig. 2A). In particular, patients in the SP-SE group exhibited the most consistent methylation status, with a mean methylation correlation value of 0.96. The mean methylation correlation values in the SP-EPN and SP-MPE were 0.93 and 0.92, respectively. This effect may be due to the little tumor heterogeneity among SE samples. According to the World Health Organization classification (36,37), SP-SE is considered a grade I tumor and its prognosis is more favorable compared with the majority of ependymal tumors (36). In addition, principal component analysis was performed to analyze the tumor samples based on the DNA methylation of the lncRNAs. The first three principal components had the highest proportion of variance values (Fig. 2B). These three components were used to characterize the DNA methylation features in these tumor subtypes, and the results suggested that patients with the same tumor subtype were likely to have similar lncRNA methylation levels (Fig. 2C). The present results suggested that the DNA methylation in lncRNA genes was correlated with histological characteristics of spinal cord ependymal tumors.

Tumor subtype-specific lncRNAs were identified based on their DNA methylation values using ANOVA. In total, 1,046 subtype-specific lncRNAs were identified [false discovery rate (FDR) <0.05; Figs. 3 and S1]. These lncRNAs tended to have higher methylation values (Fig. 3A), and the majority had methylation values ~0.9. Furthermore, the patients in the SP-EPN group had the highest methylation levels, whereas the patients in the SP-SE had the lowest methylation levels (Fig. 3B and C). GREAT was used to investigate the potential mechanisms of lncRNAs in spinal cord ependymal tumors. Among the significant Gene Ontology (38) terms (FDR q-value <0.01), three of them were associated with spinal cord development (Fig. 3D). ‘Ventral spinal cord development’ and ‘spinal cord development’ processes can affect cell differentiation of spinal cord cells, and a previous study demonstrated that cell differentiation is associated with tumorigenesis (39). Therefore, the distinct DNA methylation levels of various lncRNAs may affect the expression levels of lncRNAs in each tumor subtype, thus affecting genes associated with the nervous system and spinal cord development. The present results indicated the role of lncRNAs in human spinal cord development and suggested that they could be lncRNA signatures in spinal cord ependymal tumors.

Biological networks represent the interactions between molecules in vivo, and can be used to identify regulatory pathways and processes (40,41). lncRNAs may contribute to cancer by interacting with proteins (42), and an increasing number of studies have investigated lncRNA-protein interactions (43). lncRNA-protein interaction data in humans have been downloaded from the RAID database (30). In total, 32 subtype-specific lncRNAs exhibited associations with proteins (Fig. 4A). The network degrees of lncRNAs and proteins were separately analyzed (Fig. 4B and C). lncRNAs exhibited more associations compared with proteins. The present results indicated the important role of lncRNAs in interacting with proteins and suggested that these 32 lncRNAs may affect the disease status by affecting protein function in spinal cord ependymal tumors. Furthermore, in order to identify the lncRNA signatures in ependymal tumors, a random forest classification method was used. In total, 30 lncRNAs with importance values >0, ranging between 0.48 and 2.98, were identified (Fig. 4D). Using the methylation values of these 30 lncRNAs, the three tumor subtypes could be reliably distinguished with sensitivity and specificity values of 1 (Fig. 4E). These lncRNAs may be potential biomarkers associated with tumor subtypes, and may facilitate the diagnosis and treatment of patients with spinal cord ependymal tumors.

Long intergenic non-coding RNA 52 (LINC00052) is a lincRNA, and exhibited the highest importance value in the classification of the three tumor subtypes (Fig. 4D and F). LINC00052 has been previously shown to promote breast cancer and hepatocarcinoma development (44,45). HOX transcript antisense RNA (HOTAIR) is an oncogenic lncRNA in multiple types of cancer, including breast, gastric, colorectal and cervical cancer (46). Additionally, its DNA methylation status can serve as a biomarker in primary ovarian cancer (47). HOTAIR has been investigated as a prognostic factor in mesenchymal glioma, another type of tumor affecting the nervous system (48). In the network constructed in the present study, HOTAIR was identified to exhibit the features of a hub. HOTAIR interacted with 51 proteins and displayed tumor subtype-specific DNA methylation characteristics (Fig. 4A and G). The lncRNAs listed in Fig. 4D may be potential novel biomarkers for treating and diagnosing spinal cord ependymal tumors. Analyzing the DNA methylation values of these lncRNAs in patients with ependymal tumors may facilitate the classification of the tumors and the development of personalized therapies.

KEGG pathway enrichment analysis was performed in order to analyze the PCGs interacting with the identified lncRNA signatures. In total, 128 PCGs were analyzed. A total of 39 significantly enriched pathways (adjusted P<0.01; Fig. 5A) were identified. Out of 39 pathways, 17 were associated with cancer, including ‘transcriptional misregulation in cancer’, ‘p53 signalling pathway’, ‘colorectal cancer’ and ‘bladder cancer’. Furthermore, the pathway ‘proteoglycans in cancer’ was identified, which is involved in central nervous system-associated functions and diseases (49). Chondroitin sulphate proteoglycans (CSPGs) are enriched in the nervous system and contribute to neural cell migration and axon extension (50). A previous study demonstrated that CSPGs are upregulated after spinal cord injury (51). Additionally, the tumor necrosis factor (TNF) signaling pathway was identified, which is involved in various diseases (Fig. 5C). Interestingly, most of the lncRNAs identified to exhibit subtype-specific expression were hypomethylated compared with the total amount of lncRNAs identified (Figs. 1C and 5B).