This article is mentioned in:

Introduction

Cancer is caused by the accumulation of genomic alterations and consequent disruption of biological processes (1). The rapid progression and wide application of sequencing technologies has enabled the identification of hundreds of thousands of somatic variants in cancer (2). A significant issue in cancer genomics is the distinction of driver mutations, critical to oncogenesis, from passenger ones, which have little role in cancer initiation and progression (3). The development of reliable and efficient approaches to functionally annotate variants has been a consistent research focus in cancer-associated studies, and multiple computational tools have been investigated and widely utilized for the prediction of pathogenic mutations in the coding portion of the human genome, including the ‘sorting tolerant from intolerant’ algorithm (4) and the ‘polymorphism phenotyping’ tool (5). As an increasing number of non-coding pathogenic variants have been detected and annotated, there exists a great demand for the development of computational tools to prioritize non-coding drivers in the cancer genome (6,7). However, there have been few studies conducted in this field.

The recent completion of high-throughput projects, including the Encyclopedia of DNA Elements (ENCODE) (8), 29 Mammals Project (9) and Health Roadmap Epigenomics Project (10), has made non-coding variants interpretable. In particular, the ENCODE project has provided researchers with a genome-wide map of histone modification, DNase I hypersensitive sites, formaldehyde-assisted isolation of regulatory elements, transcription factor binding sites, RNA-seq and replication timing data across a number of cell lines (8). An increasing number of studies have taken advantage of these annotations of human functional elements to investigate non-coding disease-implicated variants or drivers in cancer, including RegulomeDB (11), HaploReg (12) and Funseq (13); the scoring systems that these approaches rely on are primarily empirical scoring algorithms, which are not scientifically rigorous and stringent (14).

Previous studies have taken advantage of machine-learning algorithms to better predict and score the functionality of non-coding variants (15–17). Kircher et al (18) contrasted the annotations of fixed or nearly fixed derived alleles in humans with those of simulated variants, and developed Combined Annotation-Dependent Depletion (CADD). CADD evaluates deleteriousness, which can be measured systematically across the genome assembly. Implementation of CADD as a support vector machine has successfully differentiated 14.7 million high-frequency human-derived alleles from 14.7 million simulated variants (18). Fu et al (19) developed a computational framework, FunSeq2, which processed large-scale genomics (including 1000 Genomes and ENCODE data) and cancer resources, and combined a high-throughput variant prioritization pipeline to annotate and prioritize somatic alterations, particularly regulatory non-coding mutations.

LncRNAs are a class of mRNA-like transcripts ranging from 200 bp to 100 kbp. They were regarded as transcription noise in the human genome, due to their lack of capability of protein translation. Over the previous decade, an increasing amount of evidence has indicated that lncRNAs have a variety of roles in numerous physiological processes (19–25). Despite a lack of capability of encoding proteins, lncRNAs may function through regulating gene expression at various levels, including chromatin architecture, transcription, RNA splicing, and protein translation and turnover (26,27). As a consequence, deregulation of lncRNAs may have a significant role in carcinogenesis (28–31).

In the present study, data concerning conservation information, regulatory features, expression and replication timing was collected, primarily from the ENCODE project, to create lung cancer-specific annotation and construct a logistic regression model based on ClinVar and HGMD pathogenic variants with the aim of functionally scoring non-coding variants in the lung cancer genome. This scoring system was applied to prioritize potential cancer-associated lncRNA candidates.

Materials and methods

Cancer mutation and pathogenic variant data

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).

Genome-wide data resources

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).

Logistic regression model training and validation

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 × 106).

Prioritization of cancer-associated lncRNA candidates

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.

RNA-seq data processing and expression analyses of lncRNAs

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−4 and absolute fold change ≥2.

Statistical analyses

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.

Results

Distinction of disease-associated non-coding variants from neutral ones with the logistic regression model

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.

Figure 1. Fitting and validation of the logistic regression model. (A) Densities of ClinVar and Human Gene Mutation Database pathogenic variants for all 25 non-coding features (red line, average density in the human genome). (B) Regression estimates for all features used in the logistic regression model. (C) Receiver operating characteristic curve for the model. (D) Scaled scores for GWAS, neutral SNPs (1 million random neutral SNPs), non-recurrent and recurrent mutations of lung cancer. CR, conserved region; TFBS, transcription factor binding site; cTFBS, conserved TFBS; UTR, untranslated region; HE, highly expressed region; SNP, single nucleotide polymorphism; Sensitive, known binding sites or motifs of transcription factors with high ratio of rare SNPs (allele frequency <0.01); ncExon, non coding Exon; H3K4me1, H3K9ac, etc., histone modification data; ER, early replicated region; Dnase, Dnase I hypersensitive site; LE, low expressed region; ECS, evolutionarily conserved structure; LR, late replicated region; TPR, true positive rate; FPR, false positive rate; GWAS, genome-wide association study.

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.

Definition and characterization of high-scoring regions in the non-coding genome

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).

Figure 2. Characterization of high-scoring regions in lung cancer. (A) Fraction of high-scoring regions in various non-coding features. (B) Average score in protein-coding gene and lncRNA introns near 5′ splice site (left panel) and 3′ splice site (right panel). (C) Fraction of high-scoring regions in various gene classes. (D) Density plot of fraction of high-scoring regions in lncRNAs. lncRNA, long non-coding RNA; PC, protein-coding; 5′SS, 5′ splice site, 10 nucleotides from the 5′ end of introns of genes; 3′SS, 3′ splice site, 50 nucleotides from the 3′ end of introns of genes; UTR, untranslated region; LUCAT1, lung cancer associated transcript 1; MALAT1, metastasis associated lung adenocarcinoma transcript 1; GAS5, growth arrest-specific 5; HOTAIR, HOX transcript antisense RNA.

Prioritization of lung cancer-associated lncRNAs with the scoring system

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−4 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).

Figure 3. Characterization of functional lncRNA candidates in lung cancer. (A) Fraction of conserved regions in functional lncRNA candidates (candidate), control lncRNAs (control) and IR. (B) Average phastCons scores for functional lncRNA candidates (candidate), control lncRNAs (control) and IR. (C) Average densities of GWAS disease or trait-related SNPs for functional lncRNA candidates (candidate), control lncRNAs (control) and IR. (D) Average densities of somatic mutations for functional lncRNA candidates (candidate), control lncRNAs (control) and IR. (E) Relative expression (log scaled FPKM) for functional lncRNA candidates (candidate), control lncRNAs (control) and IR. lncRNA, long non-coding RNA; IR, intergenic regions; GWAS, genome-wide association studies; SNP, single nucleotide polymorphism; FPKM, fragments per kilobase.

Figure 4. Expression changes for differentially expressed lncRNA candidates between lung cancer and normal samples. lncRNA, long non-coding RNA.

Discussion

In the present study, a logistic regression model was presented and used to predict ‘high-impact’ somatic alterations, combining pathogenic variants from ClinVar and HGMD databases and lung-cancer specific features. There are two main advantages of the present scoring model: Firstly, the logistic regression model took into account all non-coding pathogenic variants from HGMD and ClinVar databases, which are two well-known databases of disease-associated variants worldwide, allowing for a complete assessment of the damaging impact of any non-coding variant in the human genome. Furthermore, a large number of features used in the annotation are lung-cancer specific, including histone modifications, TFBSs, replication timing and expression data, which facilitates the scoring of variants in a lung cancer-specific manner.

Non-coding features that most positively contributed to the model include conserved regions, early replicated regions, promoter, H3K36me3, H3K4me3, conserved TFBS, TFBS and sensitive regions. Among these features, H3K36me3 is associated with actively transcribed genes, and H3K4me3 is a hallmark of actively transcribed protein-coding promoters in eukaryotes (51). These findings support the fact that conserved and regulatory elements are critical to the formation and functionality of pathogenic variants in the non-coding genome (52). The area under the ROC curve was 0.89, which outperformed two well-known tools CADD and funSeq2 (14), however, more stringent comparison must be conducted to obtain a final conclusion. Furthermore, the present model successfully distinguished GWAS variants and recurrent cancer mutations from benign SNPs and non-recurrent mutations, demonstrating the reliability and efficient performance of the model.

Given that splicing sites and UTRs are more evolutionarily conserved across mammals (53), it was observed that these regions have a higher fraction of high-scoring regions and splicing sites have a higher score compared with intronic regions. With respect to the distribution of high-scoring regions in various gene classes, it was observed that known cancer genes and cancer-associated lncRNAs demonstrated increased enrichment of high-scoring regions compared with non-cancerous genes. Based on these findings, the present study combined the fraction of high-scoring regions and average score of each lncRNA to filter out a subset of functional lncRNA candidates, which contained a number of well-characterized cancer lncRNAs, for example, HOTAIR, the expression of which is elevated in lung cancer and correlated with metastasis and poor prognosis (54). MALAT1 has been implicated in tumorigenesis and progression in a variety of cancer types (55–57). A total of 104 functional lncRNA candidates were are differentially expressed in lung cancer and normal samples. This group of lncRNAs are important candidates for cancer researchers to conduct additional experimental validation and characterization in future studies.

In conclusion, the present scoring system provides an opportunity to identify cancer-driving mutations in the vast non-coding human genome, as well as prioritizes a number of lncRNA candidates for cancer research. This scoring system may assist with the identification of driver non-coding genes for improved clinical decision-making in the future.

Acknowledgements

The present research was made possible with financial support from the National Natural Sciences Foundation of China (Beijing, China; grant no., 81272142).

References