Analysis of sequence conservation at nucleotide resolution

Abstract

One of the major goals of comparative genomics is to understand the evolutionary history of each nucleotide in the human genome sequence, and the degree to which it is under selective pressure. Ascertainment of selective constraint at nucleotide resolution is particularly important for predicting the functional significance of human genetic variation and for analyzing the sequence substructure of cis-regulatory sequences and other functional elements. Current methods for analysis of sequence conservation are focused on delineation of conserved regions comprising tens or even hundreds of consecutive nucleotides. We therefore developed a novel computational approach designed specifically for scoring evolutionary conservation at individual base-pair resolution. Our approach estimates the rate at which each nucleotide position is evolving, computes the probability of neutrality given this rate estimate, and summarizes the result in a Sequence CONservation Evaluation (SCONE) score. We computed SCONE scores in a continuous fashion across 1% of the human genome for which high-quality sequence information from up to 23 genomes are available. We show that SCONE scores are clearly correlated with the allele frequency of human polymorphisms in both coding and noncoding regions. We find that the majority of noncoding conserved nucleotides lie outside of longer conserved elements predicted by other conservation analyses, and are experiencing ongoing selection in modern humans as evident from the allele frequency spectrum of human polymorphism. We also applied SCONE to analyze the distribution of conserved nucleotides within functional regions. These regions are markedly enriched in individually conserved positions and short (<15 bp) conserved "chunks." Our results collectively suggest that the majority of functionally important noncoding conserved positions are highly fragmented and reside outside of canonically defined long conserved noncoding sequences. A small subset of these fragmented positions may be identified with high confidence.

Figure 1. Examples of SCONE p-Value Scores for Coding (A), Highly Conserved Noncoding (B), and Nonconserved Regions

Positions likely to be conserved (p < 0.05) are in light green; other positions are dark. Below each plot is the portion of the multiple sequence used to generate scores for each sequence region. Deviations from human sequence (green) are indicated in red. (A) A portion of an exon from the MET gene (chr7:115,933,744–115,933,793). The pattern of conserved positions is indicative of the triplet structure of the genetic code. (B) A highly conserved intronic sequence in the FOXP2 gene (chr7:113,646,877–113,646,926). (C) An intergenic region near the AXIN1 gene (chr16:343,046–343,095) showing little overall conservation, but containing a significant number of individually conserved positions.

Figure 2. Rare Derived Allele Frequency in Conserved versus Nonconserved Sites

Positions are partitioned according to (i) ENCODE MCS elements for all ENCODE positions, (ii) SCONE conservation score for all ENCODE positions, and (iii) SCONE conservation score for all ENCODE positions outside of MCS elements. p-Values are calculated using Fisher's exact test.

Figure 3. Conservation in Various Functional Classes

For each functional class, the fraction of positions with SCONE p-value ≤ 0.005 is shown, both including (dark) and excluding (light) positions falling within MCS elements. Ancestral repeats are included as a control.

Figure 4. Islands of Conservation in Functionally Annotated Regions

(A) Localization of short (5–12 bp) conserved islands in functionally annotated regions. Shown is the fraction of all islands that fall within a particular region. Nongenic transcribed regions were omitted to preserve scale, but contain 93% of short conserved islands. (B) Fold excess of short (5–12 bp) conserved islands in functionally annotated regions compared to ancestral repeat regions. Shown is the ratio of the density in each region (number of clusters divided by total number of positions in the region) to the density in ancestral repeat positions.