Probabilities of Fitness Consequences for Point Mutations Across the Human Genome

General Features of the Prediction Problem

Information about genetic variation can be used to estimate probabilities of fitness consequences for moderately large groups of genomic positions but not for individual loci, owing to the sparsity of informative sites along the genome. This property of “groupwise” but not “individual” predictivity is common to many statistical problems, but it is complicated in our case by two additional features. First, an appropriate scheme for grouping or stratification is not clear a priori here, because genomic correlates of fitness consequences are incompletely understood. Second, the outcomes of interest in our problem—the fitness consequences of point mutations—are not directly evident from the data. To highlight these challenges, consider the simpler problem of estimating the expected risk of an automobile accident. This problem must also be addressed at the level of groups (either explicitly through stratification of drivers, or implicitly through regression), but in this case, the relevant features—such as the age, sex, and number of traffic violations of the driver—are mostly plain to the analyst. In addition, the outcomes of interest—the occurrences and costs of accidents—are directly observed. In our problem, the genomic “risk factors” for fitness-influencing mutations, particularly in unannotated noncoding regions of the genome, are much less clear. Furthermore, once a grouping is determined, it is still not possible to read off the associated fitness consequences of mutations; instead they must be inferred from patterns of genetic variation using an evolutionary model.

Calculation of FitCons Scores

We have addressed these challenges using the following strategy. Beginning with genome-wide functional genomic data sets obtained from each cell type (Fig. 1A), we first cluster genomic positions by their joint functional genomic “fingerprints” (Fig. 1B). We focus on three highly informative and largely orthogonal functional genomic data types—DNase-seq data, RNA-seq data, and ChIP-seq data describing histone modifications—which describe DNA accessibility, transcription, and chromatin states, respectively. We divide genomic positions into three levels of DNase-seq signal, four levels of RNA-seq signal, and 26 distinct chromatin states based on the ChromHMM method^31,33. In addition, we distinguish between sites that fall outside (0) or within (1) annotated protein-coding sequences (CDSs). We then consider all possible combinations of these four types of assignments, obtaining 3×4×26×2 = 624 distinct functional genomic classes. We apply this clustering step separately to three karyotypically normal cell types: human umbilical vein epithelial cells (HUVEC), H1 human embryonic stem cells (H1 hESC), and lymphoblastoid cells (GM12878), resulting in 443–447 usable classes of sites, with median numbers of 165 to 224 thousand sites per class (see Supplementary Table 1 and Methods for details).

Next, we use INSIGHT to estimate the probabilities of mutational fitness consequences within each of these classes based on patterns of polymorphism and divergence (Fig. 1C). This step yields an estimate of the fraction of sites under selection (ρ) for each of the analyzed classes, which serves as the fitCons score for that class. Finally, we assign to each nucleotide position in the genome the score estimated for the corresponding functional genomic class (Fig. 1D). Each genomic position is thus assigned a value between 0 and 1, representing the probability that the nucleotide at that position influences fitness, as estimated from patterns of variation at all genomic sites displaying the same functional genomic fingerprint. A vital property of these fitCons scores is that they integrate information from both evolutionary data and cell-type-specific functional genomic data.

Genomic Distribution of FitCons Scores

To obtain a general overview of the genomic distribution of fitCons scores, we first considered the composition and coverage of nucleotide sites of various annotation types as a variable threshold S was applied to the fitCons score, focusing on HUVEC (see Discussion for other cell types). When S is zero, all sites are considered and the composition of annotations reflects the overall genomic distribution (Fig. 2A). As S increases, however, the sites in known functional classes become strongly enriched relative to the intergenic and intronic sites. Regions such as 5’ and 3’ UTRs, promoters, and introns are most enriched at intermediate scores, reflecting moderate levels of natural selection in these regions, while CDSs dominate at the highest scores. The coverage properties (Fig. 2B) are best for CDSs, 3’ UTRs, and 5’ UTRs (in that order), but they are also considerably elevated above the intergenic background for promoters, transcription factor binding sites (TFBS), long intergenic noncoding RNAs (lincRNAs), and small noncoding RNAs (snRNAs). Notably, the enrichment for functionally annotated genomic regions at high scores occurs despite no use of genomic annotations in the scoring scheme (except for CDS annotations). Instead, these elevated scores reflect differences in patterns of polymorphism and divergence that arise naturally from the fitness consequences of mutations in these regions and become evident after clustering based on functional genomic data. The fitCons scores for each cell type are displayed genome-wide as tracks in our public mirror of the UCSC Genome Browser (Fig. 3; Supplementary Fig. 1).

The fitCons scores generally depend in expected ways on the marginal signals of functional genomic covariates, but they are also capable of capturing complex, non-additive relationships among covariates. For example, the scores outside of CDSs increase with marginal DNase-seq (Fig. 4A) and RNA-seq (Fig. 4B) signals, as expected; yet a closer examination reveals that the scores actually decrease with DNase-seq intensity in the presence of a high RNA-seq intensity, apparently due to an implicit partitioning of 5’ and 3’ UTRs (Fig. 4C). This example demonstrates that our exhaustive partitioning scheme allows the method to capture unanticipated relationships between functional genomic covariates and natural selection.

Predictive Power for Cis-Regulatory Loci

We evaluated the predictive power of the fitCons scores for known cell-type-specific regulatory elements in comparison with three widely used phylogenetic conservation scoring methods, the phastCons¹², phyloP¹⁵ and Genomic Evolutionary Rate Profiling (GERP)¹³ programs. In addition, we considered a new program, called Combined Annotation Dependent Depletion (CADD)³⁵, that estimates relative levels of pathogenicity of potential human variants using a support vector machine (SVM), many different genomic annotations, and simulations of nucleotide divergence rates. Where appropriate, we also considered RegulomeDB, a scoring system for the regulatory potential of variant sites based on combined experimental and computational data³⁶, and EnhancerFinder, a kernel-based predictor for developmental enhancers based on multiple data types³⁷. We evaluated the performance of these methods in predicting three types of functional elements that have putative roles in transcriptional regulation based on different data sets: (1) binding sites for various transcription factors supported by chromatin immunoprecipitation and sequencing (ChIP-seq) data from the ENCODE project^3,28; (2) high-resolution expression quantitative trait loci (eQTL) identified in a recent large-scale study⁶; and (3) enhancers based on characteristic chromatin marks³⁸ (see Methods for details).

To place the different predictors on equal footing, we plotted the basewise coverage of each type of regulatory element as a function of the total coverage of the noncoding genome, varying score thresholds to include 0–20% of noncoding sites (Fig. 5). This strategy allows us to measure the extent to which the elements of interest display signals that rise above the background of the noncoding genome, in a uniform manner across scoring methods. By this test, the fitCons scores showed dramatically better sensitivity for noncoding elements than almost all of the other methods considered. For example, at a total noncoding coverage of 2.5%, fitCons scores achieve nearly 70% coverage of TFBSs, whereas the other methods all have less than 20% coverage. Similarly, the coverage of enhancers was about 40% at 2.5% noncoding coverage, while most other scoring methods showed almost no signal above background. Only EnhancerFinder, which is specifically designed for this task, showed comparable prediction performance on enhancers. We also performed a more traditional evaluation of the trade-off between sensitivity and specificity using receiver operating characteristic (ROC) curves and found that fitCons scores were considerably better predictors of regulatory function than all other methods considered (see Methods and Supplementary Fig. 2).

The tests above were based on regulatory elements that are putatively active in the cell type for which the scores were produced, in order to highlight the benefits of using cell-type-specific functional data. To evaluate how well these advantages extended across cell types, we created an “integrated” fitCons score by combining information from three cell types (see Methods), and evaluated the performance of this score in predicting regulatory elements pooled from multiple cell types. We found that, in this less favorable setting, the fitCons scores still had better predictive performance for cis-regulatory elements than any of the other scoring methods (Supplementary Fig. 3).

To address possible deficiencies of these tests, we carried out two additional sets of validation experiments. First, we performed a second round of experiments on ChIP-seq-supported TFBSs that considered only the subset of nucleotide positions at which base preferences were especially strong, which should be enriched for bases having fitness consequences. The ROC curves based on this more stringent test were very similar to the original curves (Supplementary Fig. 4), demonstrating that the apparent performance of the fitCons scores was not artificially inflated by the coarse-grained nature of our scores and TFBSs. Second, we examined an alternative set of predicted enhancers for GM12878 cells based on characteristic patterns of divergent transcription initiation³⁹. Unlike the chromatin-based enhancer predictions described above, these predictions are based on data completely independent from those underlying the fitCons scores. Nevertheless, the fitCons scores still displayed excellent predictive power for this set, better than all other methods besides EnhancerFinder (Supplementary Fig. 5).

Proportion of the Human Genome Under Selection

The proportion of nucleotides in the human genome that directly influence fitness—sometimes called the “share under selection” (SUS)—has primarily been estimated using methods that consider divergence patterns among mammals, for which turnover of functional elements may be an important confounding factor^40–44. In addition to being useful as predictors of function, the fitCons scores could be useful in obtaining estimates of the SUS that are less sensitive to turnover because they measure natural selection over much shorter time scales.

An initial estimate of the SUS can be obtained by simply averaging the fitCons scores across all nucleotide positions in the genome. Because each score represents a probability that an individual nucleotide influences fitness, their average represents an expected fraction of nucleotides in the genome having fitness-influencing functions, or an expected SUS. This approach yields an estimate of 7.5% (±0.1%) for HUVEC, or 7.5–7.8% across cell types. These estimates are largely consistent with, but on the high end, of those based on cross-species divergence, which generally have fallen between 3 and 8%^{12,40,44–46}. Among sites under selection, we estimate that 9.0% are in CDS, 2.2% in 3’ UTRs, 35.2% in introns, 51.7% in intergenic regions, and <1% in each of several other noncoding annotation classes (Supplementary Table 2). Our estimates of the SUS are somewhat lower than previous estimates that have explicitly allowed for evolutionary turnover, most of which have been 2–3 times higher than the pan-mammalian estimates of ~5%^{26,44,46–48}. However, they are similar to a recent estimate of 7.1–9.2% based on improved alignments and a new model for turnover⁴⁹ (see Discussion).

Violations of modeling assumptions will tend to bias the fitCons scores upward, particularly for functional classes for which the true fraction is close to zero (see Supplementary Note). To address this problem, we performed a parallel calculation for “neutral” sites that intersect the large class of genomic positions having a “null” functional genomic fingerprint (i.e., no DNase-seq, RNA-seq, or histone modification signal). This calculation results in an estimate of 3.3%, which can be considered an upper bound on the contribution of error because these putatively neutral sites undoubtedly include some sites under selection. By subtracting this 3.3% from our naive estimate of 7.5%, we obtain an estimated lower bound for the SUS of 4.2%, with somewhat higher fractions of selected sites in CDSs and 3’ UTRs (Supplementary Table 2). (These estimates are for HUVEC, but the results for the other cell types were very similar.) Overall, our analysis of the SUS suggests that between 4.2% and 7.5% of nucleotides in the genome have direct fitness-influencing functions, and that the ratio of noncoding to coding functional sites is between 5.4 and 10.1.

Implications for Evolutionary Turnover of Functional Elements

To better understand the differences between fitCons scores and conventional divergence-based scores, we devised an alternative scoring system (denoted “fitConsD”) based on the same site clusters but an estimator of the fraction of nucleotides under selection that instead considers nucleotide divergence patterns across primates (see Methods). Thus, the fitCons and fitConsD scores both represent probabilities of fitness consequences per nucleotide but over two different evolutionary time scales. Overall, these two measures are remarkably well correlated, with R²=0.88 (Fig. 6A). Furthermore, a measure based on the difference between fitConsD and fitCons scores suggests relatively low amounts of turnover across annotation classes, accounting for no more than about 10% of all functional sites (Fig. 6B). These observations suggest that the main signal for selection has been maintained over long evolutionary time periods and that turnover has been modest during primate evolution, but that there are some classes of sites that show stronger recent than ancient natural selection.

Functional Genomic Data

RNA-seq and DNase-seq data for HUVEC, H1 hESC, and GM12878 were downloaded from the University of California, Santa Cruz (UCSC) Genome Browser. Chromatin states for the same three cell types were downloaded from the European Bioinformatics Institute ftp site (see Supplementary Table 3). For DNase-seq, we considered two replicate experiments from University of Washington (UW) data for each cell type. However, only one UW replicate was available for H1 hESC so additional DNase-seq data for this cell line was obtained from Duke University. For each replicate DNase-seq experiment, we downloaded broad and narrow peak calls. For RNA-seq, we selected a single replicate from the Caltech poly-A+ 75 bp paired-end read data, after examining several alternative data sets. For chromatin states, we used the 25-state ChromHMM segmentation generated in December, 2012³³.

Clustering Approach

We produced a separate partitioning for each cell type based on the functional genomic data. The broad and narrow DNase-seq peaks were used to partition sites in the genome into three mutually exclusive classes: sites that fall in a narrow peak in both replicate experiments (2); sites that fall in a broad peak in at least one replicate and do not fall in a narrow peak in both replicates (1); and sites that fall outside of all called peaks (0). This three-level scheme allows for both high sensitivity (class 1) and high specificity (class 2). For H1 hESC, only one set of broad peak calls was available to define class 1. For the RNA-seq data, we partitioned sites in the genome into four mutually exclusive classes (0–3) based on numbers of reads aligned at each position. Read depth thresholds were set separately for each cell type through a process that aims to minimize the conditional entropy of concentrations of predicted sites under selection (see Supplementary Note). Chromatin states were defined directly from the 25 states in ChromHMM, with a 26th state containing sites not assigned to any chromatin class. The Cartesian product of these partitions, together with the partition into coding and noncoding sequences (see below), resulted in 3×4×26×2 = 624 distinct functional classes.

Running INSIGHT

The INSIGHT method infers the fraction of nucleotide sites under selection (ρ) for a given collection of sites by comparing patterns of within-species polymorphism and between-species divergence within those sites and within putatively neutrally evolving sites nearby. A detailed description of the method, sequence data, and data-quality filters is given in reference [34]. For each non-empty fitCons site cluster, INSIGHT was used to estimate ρ, which was used as the fitCons score of all sites in that cluster. To reduce sensitivity to estimates with high uncertainty, we filtered out clusters for which the estimated standard error was larger than 40% of the estimated value of ρ. To increase computational efficiency, clusters larger than 20 Mb were partitioned into smaller subclasses, and estimates of ρ were computed as weighted averages (by numbers of informative sites) across subclasses.

Neutral Sites

The collection of sites predicted to be free from the influence of natural selection (“neutral” sites) was derived from a set identified previously^28,34,53. Briefly, this set is obtained by eliminating from the genome sites likely to be under direct natural selection, including (1) exons of annotated protein-coding genes and the 1000 bp flanking them on either side; (2) RNA genes from GENCODE v11 and 1000 bp flanks; and (3) conserved non-coding elements (identified by phastCons) and 100 bp flanks. This set was used both by INSIGHT and for our power analysis (see below).

GENCODE Annotations

Transcript annotations from GENCODE v15⁵⁴ were downloaded from the Sanger Institute FTP server and used to define eight site classes: coding sequences (CDS), 5’ untranslated regions (UTRs), 3’ UTRs, promoters, introns, long intergenic noncoding RNAs (lincRNA), short noncoding RNAs (sncRNA), and intergenic (sites not falling within any protein-coding transcription unit). Transcripts annotated with feature type=“CDS” and gene type=“protein coding” were used to define the CDS set for fitCons. For subsequent analysis, we used a slightly more conservative set, obtained by additionally requiring feature type=“gene”, gene status=“KNOWN”, transcript status=“KNOWN”, and the identification of both start and stop codons within the transcript. UTRs were defined from transcripts having feature type=“UTR” and gene type=“protein coding” and designated as 5’ or 3’. Introns were defined by positions that fall within a protein-coding transcript but outside of the CDS and UTR regions. Promoters were defined as the 1000 bp immediately upstream of the first (i.e., most upstream) transcription start site for each protein-coding gene. A similarly defined alternative set of 100 bp promoter regions was used in assessing differences between cell types (see Supplementary Fig. 9). LincRNAs were identified by transcripts with feature type=“exon” and gene type=“lincRNA”. Similarly, sncRNAs consisted of transcripts with feature type=“exon” and gene type ∈ {“miRNA”, “snRNA”, “snoRNA”}. Positions in the more inclusive CDS set were removed from all noncoding classes.

Cis-Regulatory Elements

Transcription factor binding sites (TFBSs) were drawn from a set for 78 transcription factors, based on ChIP-seq data from ENCODE²⁸ downloadable from our Genome Browser mirror. This set contains roughly 1.4 million binding sites of mean length 11 bp, each of which is associated with the cell types in which it was detected. For some tests, we considered only the subset of nucleotide positions inside these TFBSs that corresponded to motif positions with strong base preferences, defined as those positions at which the consensus allele appeared in at least 90% of all binding sites (according to the inferred motif model). For enhancers, we used the distal regulatory modules described in reference [³⁸]. We downloaded the file enets4.Distal cell line.txt from the Gerstein Lab and extracted from it a total of 19,005 enhancer-transcript associations, covering 5,834 unique autosomal loci with a mean length of 888 bp, along with the cell types associated with each predicted enhancer. Expression quantitative trait loci (eQTL) described in reference [6] were downloaded from the European Bioinformatics Institute’s E-GEUV-1 data set. We used the four files {EUR373, YRI89}.{exon, gene}.cis.FDR5.best.rs137.txt.gz to identify 6,760 distinct autosomal positions and the associated transcripts, removing all positions overlapping CDSs.

Identifying Active Elements per Cell Type

In several analyses, we considered the subset of elements in each annotation class for which we had evidence of activity in a given cell type. To identify the cell types in which TFBS and enhancers were active, we used the cell type designations provided in the corresponding annotation files (see above). For other classes of elements, we defined the active elements using a set of GENCODE transcripts and genes that showed significantly elevated levels of RNA transcription in the Caltech RNA-seq data. Those were transcripts (or genes) for which the 95% confidence interval of the normalized read count in a given cell type fell within the top one third of normalized read counts for transcripts (or genes) across all three cell types (threshold 1.477 for transcripts and 4.966 for genes). Active eQTLs were identified via associated active genes using the GENCODE gene identifier specified for each eQTL. Active promoters, UTRs, CDS, and introns, were identified via associated active transcripts. For the comparison between cell types (Supplementary Fig. 9) we also used collections of eQTLs and promoters found to be inactive in a given cell type. Those were defined in a similar way, by using transcripts and genes falling in the bottom third of the distribution of normalized read counts.

Comparison with Other Scores

Base-wise scores from the Genomic Evolutionary Rate Profiling (GERP)¹³ method, the Combined Annotation Dependent Depletion (CADD) method, and phastCons¹² and phyloP¹⁵ methods were downloaded from the respective websites (see URLs; file hg19.GERP scores.tar.gz generated in August, 2010 for GERP, file whole genome SNVs.tsv.gz, downloaded in September, 2013 for CADD, and the 46 placental mammals conservation tracks for phastCons and phyloP). CADD scores are specified for each genomic position and each of the three possible variant bases in that position. We took the maximum of these three scores, which yielded the best performance for CADD in our comparisons. We also used RegulomeDB³⁶ (downloaded in January, 2013) to rank single nucleotide polymorphisms (SNPs), such as eQTLs, into one of 13 categories according to evidence from functional genomic data. Finally, we obtained EnhancerFinder scores³⁷ for 1500 bp windows tiled across the genome directly from the authors. We used the general, non-tissue-specific scores, and averaged them at positions contained in multiple overlapping windows.

Receiver Operating Characteristic (ROC) Curves

We used ROC curves to measure the ability of each scoring scheme to discriminate between functional and nonfunctional regulatory elements. For TFBSs and enhancers, we used the annotations described above as “true positives” and defined “true negatives” from our filtered, putatively neutral sites (see above). For eQTLs, our negative set consisted of all 9.8 million variants tested in reference [6], excluding indels, non-simple variants, and positions that showed possible associations at a threshold of nominal p < 0.05 (7.6 million SNPs remained). In all three cases, we additionally removed any sites in the positive set from the negative set. A point on a ROC plot indicates the fraction of annotated genomic positions with scores higher than a given score (true positive rate) versus the fraction of control genomic positions with scores higher than that score (false positive rate). Positions with no scores are ignored when computing fractional coverage.

Integrating fitCons Scores across Cell Types

We generated a series of fitCons scores that integrate functional genomic data across the three cell types by using the original 624 fingerprints (see above) and altering the rule by which sites are assigned to clusters to reflect information from multiple cell types. Our approach attempts to select a fingerprint for each site that was likely to be most informative about the site’s function, while avoiding a bias toward higher scores with increasing cell types. See Supplementary Note for details.

Share Under Selection

Assume a partitioning of the genome into K mutually exclusive and exhaustive clusters, C₁,…,C_K, and a corresponding set of fitCons scores, ρ(C₁),…,ρ(C_K). Note that the expected number of genomic positions under selection in cluster C_i is given by ρ(C_i)|C_i|, because ρ is an estimate of the fraction of sites under selection. For an arbitrary collection of sites, S, the expected number of sites in S that are under selection is given by sel(S) = ∑_iρ(C_i) |C_i∩S|, and the average fitCons score for S is given by ρ(S) = sel(S) / |S|. To avoid under-estimation of ρ(S), we do not filter out fitCons scores with high uncertainty in these calculations, as we do for other analyses (see above). In addition, to account for possible overestimation of ρ in very large clusters having low fractions of sites under selection, we ran INSIGHT on the intersection of our neutral sites (see above) and all non-coding sites in the ‘Quiescent’ chromatin state with no DNase-seq or RNA-seq signal. We then subtracted the estimated value of ρ, denoted ρ_neut, from the raw fitCons score to obtain a conservative lower bound, ρ(S)–ρ_neut, for the fraction of sites under selection in S.

FitConsD and Evolutionary Turnover

To make the comparison between fitCons and fitConsD as direct as possible, fitConsD scores were computed using the same pipeline we developed for fitCons (see Fig. 1), except that in step C we replaced the INSIGHT model with an evolutionary model that considers sequence divergence between the human, chimpanzee, orangutan, and rhesus macaque genomes. The fitConsD scores are based on an estimate s_i for the relative evolutionary rate of each cluster C_i compared with a neutral model globally estimated for the four-primate phylogeny. This relative rate is then compared with the relative rate $s_i^{\text{neut}}$ estimated for the putative neutral regions flanking sites in C_i and a divergence-based estimate of the fraction of sites under selection in cluster C_i is given by ${\rho }_{\text{div}}(C_i)=1-s_i/s_i^{\text{neut}}$ (see Supplementary Note).