A total of 1,623,250 somatic mutations detected by whole genome sequencing of 24 pairs of lung cancer and normal specimens were obtained from the supplementary data files of a previous study (32). Recurrent mutation represents two or more mutations that have the same mutation site across multiple samples (n=14,515 mutations). Non-recurrent mutation denotes mutations that only occur once in all patients. Germline polymorphism data comprising 38,248,779 single nucleotide polymorphisms (SNPs) was downloaded from the 1000 Genome project pilot 1 (www.1000genomes.org) (33). SNPs with derived allele frequencies >0.01 were considered to be neutral SNPs; rare SNPs denote those whose allele frequencies were <0.01. Disease-associated variants data from ClinVar (www.ncbi.nlm.nih.gov/clinvar) and the Human Gene Mutation Database (HGMD; www.hgmd.cf.ac.uk) are known (published) gene variants responsible for human inherited diseases (34,35). Trait or disease-associated SNPs were obtained from genome-wide association studies (GWAS; www.gwascentral.org) (36).

Human genome annotations were obtained from Gencode (www.gencodegenes.org/) (37), including protein coding genes, exons, introns, untranslated regions (UTRs) and non-coding exons (37). lncRNA annotation was primarily acquired from three different sources, Gencode (37), Human Body Map large intergenic non-coding RNAs and transcripts of uncertain coding potential generated from 4 billion RNA-seq reads across 24 tissues and cell types (38) and Refseq annotation (www.ncbi.nlm.nih.gov/refseq/) (39). In total, there were 39,952 lncRNA annotations collected from these three different databases. The 5′ splicing sites were 10 nucleotides from the 5′ end of introns of genes (40). The 3′ splicing sites were 50 nucleotides from the 3′ end of introns of genes (41). Evolutionarily conserved bases were identified using a recently published analysis of 46 mammalian genomes (42). A genome-wide phastCons score was obtained from Siepel et al's study (16) (hgdownload.cse.ucsc.edu/goldenPath/phastConsPaper/vertebrate-scores/). Sensitive regions from Khurana et al (13) consisted of binding sites or motifs of important transcription factors and contained an increased fraction of rare SNPs. Evolutionarily conserved structures were RNA secondary structures predicted using comparative structure prediction algorithms based on multiple genomes (42). Promoters, defined as regions 2.5 kb from transcription start sites (TSS), were generated from the Gerstein lab (http://funseq.gersteinlab.org/data) (13). RNA-seq data in bam format, transcription factor binding sites (TFBS), DNase I hypersensitive sites and histone modification data (H3K4me1, H3K9ac and others) of the A549 cell line were acquired from ENCODE (8). Conserved TFBS were transcription factor binding sites conserved in the human/mouse/rat alignment and obtained from University of California, Santa Cruz directly (41). The expression level was calculated by counting the number of reads per kilobase per million reads (RPKM) for each protein coding gene and lncRNA. Genes whose RPKM was >20 or <0.25 were defined as high and low expressed regions, respectively. A wavelet-smoothed, weighted average signal was used, and the high and low signal values corresponded with early and late replication during the S phase, respectively (genome.ucsc.edu/ENCODE, ‘Repli-seq track’) (8). Genome-wide replication timing was mapped to protein coding genes and lncRNAs. An early-to-late ratio was calculated as (G1b+S1)/(S4+G2) for each protein coding gene and lncRNA (43). When the ratio (G1b+S1)/(S4+G2) was >1, genes were considered to be early replicated, while late replicated genes had an early-to-late ratio <1.

Cancer lncRNAs containing 25 lncRNAs are a collection of mammalian long non-coding transcripts that have been experimentally demonstrated to be associated with a variety of cancer types. A list of cancer census genes was obtained from the current release of the catalogue of somatic mutations in cancer version 71 (COSMIC; cancer.sanger.ac.uk/cosmic) (44).

The disease-implicated set of variants was composed of 19,153 non-coding pathogenic variants from the ClinVar and HGMD databases. For the control sets, the present study used neutral variants whose minor allele frequency was ≥1% to reduce the possibility of including functional rare SNPs. A total of 15,789,242 potential control SNPs were included in the model. In the logistic regression model, a matrix of 425,565 rows was formed throughout the non-coding genome, and each row represented one unique combination of features. Disease-causing variants from HGMD and ClinVar databases and neutral SNPs were used as the binary response variables, and the 25 genomic features served as the predictor variables to predict the likelihood of a variant being disease-associated. The logistic regression model was constructed with the general linear model. The receiver operating characteristic (ROC) curve was generated with a R script (version 2.15.3; www.r-project.org). Scores were predicted with the model for GWAS, neutral SNPs, and non-recurrent and recurrent somatic mutations of lung cancer and subsequently scaled using the following formula: scaled score=log(predicted score × 10⁶).

Cancer-associated lncRNA candidates were determined with the following criteria. Firstly, the logistic regression model was used to score each nucleotide of the lncRNAs and the average score was calculated for each lncRNA. Secondly, 100 Mb non-coding regions whose scores were >8.4149 were defined as high scoring regions, and the fraction of high scoring regions for each lncRNA was calculated. Subsequently, the final subset of lncRNA candidates was determined by identifying the overlap between the top 10% of lncRNAs with the highest average score and the top 10% of lncRNAs with the highest fraction of high scoring regions.

A total of 161 RNA-seq data samples, including 76 normal lung samples and 85 cancerous samples, were obtained from the Ju et al (45) study at the European Bioinformatics Institute. Reads were mapped to the hg19 genome using the Star aligner (https://github.com/alexdobin/STAR/releases) (46). Read counts were calculated with bedtools version 2.22.1 (bedtools.readthedocs.org/en/latest/#) for each lncRNA (47). The expression level in FPKM was calculated with Cufflinks version 2.2.1 (cole-trapnell-lab.github.io/cufflinks/) (48) and log scaled for each lncRNA. DESeq2 Release version 3.0 (bioconductor.org/packages/release/bioc/html/DESeq2.html) (49) was used to identify differentially expressed transcripts between tumor and normal pairs, with a cutoff of false discovery rate (FDR) ≤10⁻⁴ and absolute fold change ≥2.

Data are presented as the mean ± standard deviation. Differences between different groups were drawn with the two-sided Mann-Whitney U test or Fisher's exact test in R (version 2.15.3; www.r-project.org). P<0.05 was considered to indicate a statistically significant difference.

Estimates of the densities of ClinVar and HGMD disease-causing variants revealed that the densities of disease-associated variants varied greatly across various non-coding features (Fig. 1A). Certain features, including conserved regions, conserved TFBS, UTRs, promoters and highly-expressed regions, demonstrated the highest enrichment of pathogenic variants; however, features including H3K9me3, late replicated regions, H3K27me3, evolutionarily conserved structures and H2az had low densities of disease-causing variants, suggesting that different non-coding features have importance to the functionality of non-coding variants. It was observed that conserved regions, early replicated regions, promoters, H3K36me3 and conserved TFBSs most positively contributed to the model, while H3K79me2, H3K4me2, H3K9me3, H3K9ac and low-expressed regions were the most negatively informative for the model (Fig. 1B). It was demonstrated that the area under the ROC curve was 0.89 for the logistic regression model (Fig. 1C), which indicated that the model was able to discriminate between disease-implicated and control variants with a high specificity and sensitivity.

To investigate whether the present model could be applied to prioritize candidate functional variants, the disease or trait-associated variants from GWAS were selected for an independent validation. It was observed that non-coding GWAS SNPs had a significantly higher average score compared with 1 million random, neutral SNP control variants (mean, 5.9012 vs. 5.5238; P<0.001, two-sided Mann-Whitney U test; Fig. 1D). Recurrence is considered to be a potential sign of positive selection among tumors and is more likely to be associated with driver events (50). Subsequently, the present study evaluated recurrent mutations that occurred at the exact same site across >2 samples, as well as non-recurrent mutations, identified by whole-genome sequencing of 24 lung cancer samples. It was identified that the same-site recurrent mutations (n=14,515 mutations) had significantly higher scores compared with the non-recurrent mutations (mean, 5.4677 vs. 5.2277; P<0.001, Mann-Whitney U test; Fig. 1D), which suggested that this approach may be useful for the identification of non-coding driver mutations in lung cancer.

The present study defined 100 Mb non-coding regions, which were scored >8.4149 as high-scoring regions, and analyzed fractions of high-scoring regions in a variety of feature types. The 5′ and 3′ splice sites and UTRs were among the features that contained the highest fraction of high-scoring regions; by contrast, intergenic regions, lncRNA introns and lncRNA demonstrated the lowest fraction of high-scoring regions (Fig. 2A). The present model assigned a higher average score to splicing sites compared with adjacent intronic regions in protein coding genes (mean, 9.4374 vs. 8.3959; P<0.001, Mann-Whitney U test; Fig. 2B) and lncRNAs (mean, 8.1802 vs. 7.8146; P<0.001, Mann-Whitney U test; Fig. 2B). Subsequently, the present study sought gene classes with various fractions of high-scoring regions and identified that known cancer genes from COSMIC had a significantly increased fraction of high-scoring regions compared with non-cancerous ones (mean, 0.0817 vs. 0.0596; P<0.001, Fisher's exact test; Fig. 2C). Cancer-associated lncRNAs that were collected from recent publications demonstrated a significantly increased fraction of high-scoring regions compared with non-cancerous ones (mean, 0.1112 vs. 0.0590; P<0.001, Fisher's exact test; Fig. 2C), for example, HOX transcript antisense RNA (HOTAIR), metastasis associated lung adenocarcinoma transcript 1 (MALAT1), growth arrest-specific 5 (GAS5) and lung cancer associated transcript 1 are among the top 10% of lncRNAs with respect to high-scoring coverage (Fig. 2D).

Regarding prioritization of lung cancer-implicated lncRNAs, the fraction of high-scoring regions and average score were calculated for each lncRNA. Subsequently, overlapping lncRNAs were determined between the top 10% of lncRNAs with the highest fraction of high scoring regions and the top 10% of lncRNAs with the highest average score. A total of 2,679 lncRNAs were filtered out as functional candidates, including some experimentally characterized cancer-associated lncRNAs, including MALAT1, HOTAIR and GAS5. In the present study it was demonstrated that this subset of lncRNA candidates had a significantly increased fraction of conserved regions (mean, 0.1741 vs. 0.0528; P<0.001, Mann-Whitney U test; Fig. 3A) and average phastCons score (mean, 0.2770 vs. 0.2602; P<0.001, Mann-Whitney U test; Fig. 3B) compared with control lncRNAs, indicating that they were more conserved relative to control lncRNAs. It was additionally observed that this subset of lncRNAs had an increased enrichment of disease or trait-associated GWAS SNPs (mean, 6.2106 vs. 4.0618 SNPs/Mb; P<0.001, Fisher's exact test; Fig. 3C) and a lower somatic mutation density compared with the control lncRNAs (mean, 329.8380 vs. 573.2742 mutations/Mb; P<0.001, Fisher's exact test; Fig. 3D). RNA-seq data of 76 normal lung samples and 85 cancer samples were obtained from Ju et al's (45) study, which is publicly available from the European Bioinformatics Institute. Read alignment was conducted with a Star aligner and coverage was calculated for each lncRNA with bedtools software. DESeq2 was used to investigate the differential expression of lncRNAs between lung cancer and normal samples. It was observed that the lncRNA candidates showed significantly increased expression compared with control lncRNAs in cancerous and normal samples (log scaled FPKM, 1.8924 vs. 1.1386; P<2.2e-16, Mann-Whitney U test; Fig. 3E). Differentially expressed lncRNAs were determined based on the criteria that lncRNAs have cutoff FDR <10⁻⁴ and absolute fold change >2. The number of differentially expressed lncRNAs was 2,208, and 104 of them were among the list of potentially cancer-associated lncRNAs (Fig. 4).