|Home | About | Journals | Submit | Contact Us | Français|
Multiple computational tools have been widely applied to the detection of coding driver mutations in cancer; however, the prioritization of pathogenic non-coding variants remains a difficult and demanding task. The present study was performed to distinguish non-coding disease-causing mutations from neutral ones, and to prioritize potential cancer-associated long non-coding RNAs (lncRNAs) with a logistic regression model in lung cancer. A logistic regression model was constructed, combining 19,153 disease-associated ClinVar and Human Gene Mutation Database pathogenic variants as the response variable and non-coding features as the predictor variable. Validation of the model was conducted with genome-wide association study (GWAS) disease- or trait-associated single nucleotide polymorphisms (SNPs) and recurrent somatic mutations. High scoring regions were characterized with respect to their distribution in various features and gene classes; potential cancer-associated lncRNA candidates were prioritized, combining the fraction of high-scoring regions and average score predicted by the logistic regression model. H3K79me2 was the most negative factor that contributed to the model, while conserved regions were most positively informative to the model. The area under the receiver operating characteristic curve of the model was 0.89. The model assigned a significantly higher score to GWAS SNPs and recurrent somatic mutations compared with neutral SNPs (mean, 5.9012 vs. 5.5238; P<0.001, Mann-Whitney U test) and non-recurrent mutations (mean, 5.4677 vs. 5.2277, P<0.001, Mann-Whitney U test), respectively. It was observed that regions, including splicing sites and untranslated regions, and gene classes, including cancer genes and cancer-associated lncRNAs, had an increased enrichment of high-scoring regions. In total, 2,679 cancer-associated lncRNAs were determined and characterized. A total of 104 of these lncRNAs were differentially expressed between lung cancer and normal specimens. The logistic regression model is a useful and efficient scoring system to prioritize non-coding pathogenic variants and lncRNAs, and may provide the basis for detecting non-coding driver lncRNAs in lung cancer.
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.
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 × 106).
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−4 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−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).
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.
The present research was made possible with financial support from the National Natural Sciences Foundation of China (Beijing, China; grant no., 81272142).