Functional constraint and small insertions and deletions in the ENCODE regions of the human genome

Indel rates were observed to be reduced approximately twenty-fold in exonic ENCODE regions, five-fold in sequence that exhibits high evolutionary constraint in mammals and up to two-fold in some classes of regulatory elements.


Background
Insertion-deletion polymorphisms (indels) have to date received less attention in the study of sequence variation than have single nucleotide polymorphisms (SNPs), despite their frequency (estimated at approximately 16% to 25% of all sequence polymorphism events) and their potential functional importance [1]. 5' Untranslated regions (UTRs) and gene coding regions have previously been observed to have lower indel rates compared with other regions, suggesting that the constraint may have arisen because of negative selection [2]. In general, indels that give rise to frame shifts in coding sequence are more disruptive than non frame-shifts and single point mutations, because of third base degeneracy [3]. As a result, coding sequence indels tend to have lengths that are multiples of three, whereas regulatory sequences tend to have more frequent indels that occur in distinct blocks [4]. The majority of indels are di-allelic and small, with allele length differences of relatively few (one to four) nucleotides [2,5,6]. Given their frequency, small indels could play an important role in contributing to phenotypic differences in humans, including susceptibility to diseases. It is therefore of interest to characterize indel distribution across the human genome, and to integrate indels into SNP marker maps in order to aid in the identification of natural genetic variation.
Recent theoretical work has considered the distribution of indels under neutrality and exploited the evolutionary imprint of sequence indels in order to pinpoint functional DNA regions that are subject to purifying selection [7]. Snir and Pachter [8] used Encyclopedia of DNA Elements (ENCODE) data and multiple primate sequences to study indel events between species. This work suggests that indel rates genome wide are not uniform and that indel events are not neutral; in particular, the work has identified indel hotspots in the human genome. A minority of insertions and deletions may also have plausibly played a major role in speciation events, including human-chimpanzee phenotypic differences [9,10]. An investigation of 2,000 human di-allelic indels found that the majority were monomorphic in chimpanzees and gorillas, indicating that most indels have arisen after the most recent common primate ancestor [6] and are lineage specific [5].
We used the small insertion and deletion ENCODE data [11] to address four questions. First, do the 14 manually selected regions have lower insertion and deletion rates compared with the 30 randomly selected regions? This might be expected to be the case if the selection process [12] for the manually selected ENCODE regions of interest were biased toward regions with greater density of genes or genes of evolutionary importance, with greater functional and evolutionary constraints. Second, do indel rates vary by genomic annotation feature (in turn reflecting varying levels of functional constraint)? Indels that arise in coding sequence are more likely to be deleterious and therefore subject to purifying selection. As a result, DNA sequences that encode proteins might be expected to have some of the lowest genomic indel rates, followed by a wide variety of functional features that are believed to regulate gene expression via an increasing number of previously unrecognized mechanisms [13][14][15][16][17].
Third, are indel rates negatively correlated with measures of evolutionary constraint? We expect indel rates to be negatively associated with evolutionary constraint scores (see Materials and methods, below) where DNA sequences are subject to purifying selection. To address this question, we also correlated indel rates with ancestral repeat (AR) sequence. AR sequences are mobile elements that inserted before the common ancestor of most mammals and have subsequently become inactive [18]. ARs are considered to be predominantly neutral sequences (not subject to purifying selection) and hence we would anticipate indels to accumulate in AR sequence regions with relatively little or no constraint. Based on the assumption that new indels have arisen in AR regions in the past at the same rate as elsewhere in the genome, observed indel rates might be expected to be positively correlated with AR sequence rates.
The fourth question we consider is how do ENCODE indel rates compare with SNP rates across genomic features and evolutionary constrained sequence?
Here we describe the distribution of small indels (ranging from 1 to 20 base pairs [bp]) in the manually and randomly selected ENCODE regions, their distribution in relation to genomic annotation features, and their relationship with measures of evolutionary constraint.

Results
All identified small indels (n = 4486) in the ENCODE regions were mapped onto physical coordinates for ENCODE functional features. The average indel length of identified small indels is 2.8 bp, ranging from 1 to 20 bp. The overall density is on average 15 indels per 100 kilobases (kb; 99% confidence interval [CI] 13.4 to 16.7) or, in terms of total indel length, 43.4 bp per 100 kb (99% CI 38.3 to 49.1). All results in Tables  1 to 3 are presented in two ways: as numbers of indel events (indels per 100 kb) and total indel length (indel bp per 100 kb). In the interests of brevity, indel rates are referred to in the text to as indel bp per 100 kb unless stated otherwise. This also facilitates comparison with SNP rates.
There are no substantial differences in indel or gene density between manually and randomly selected regions ( Table 1). The indel rates in manual regions are similarly variable (sd num/100 kb = 5.0 number of indels per 100 kb; sd bp/100 kb = 14.7 indel bp per 100 kb, where sd num/100 kb and sd bp/100 kb refer to the standard deviation for number of indels and indel bp per 100 kb, respectively) to those in random regions (sd num/ 100 kb = 4.0; sd bp/100 kb = 14.0), with no significant differences in the summary data (F [13,29] = 1.52, P = 0.34).
We observed a reduction in indel rates for coding sequence and annotation features that are believed to play a regulatory role in gene expression (Table 2). Compared with the overall mean (43.4 bp per 100 kb), ENCODE coding sequences all   [19] and Table 4 for details). These sites show modestly reduced indel rates (HisPolTAF: 32.4 bp per 100 kb), along with sites occupied by sequence specific binding proteins (all motifs: 35.8 bp per 100 kb), but neither finding is statistically significant.
Multi-species constrained sequence (MCS moderate; 11.2 bp per 100 kb) show greatly reduced indel rates (Table 2), similar to rates in coding regions. AR regions (26.5 bp per 100 kb)  Indel rate versus GERP score comparing human and primates Figure 2 Indel rate versus GERP score comparing human and primates. Indel rate is expressed as base pairs (bp) per 100 kilobases (kb). The solid line represents the fit from a cubic smoothing spline, whereas the dashed line is the fit from a robust linear regression. GERP, genomic evolutionary rate profiling. Indel rate versus all AR sequence rate Figure 3 Indel rate versus all AR sequence rate. Indel rate and ancestral repeat (AR) sequence rate are both expressed as base pairs (bp) per 100 kilobases (kb). The solid line represents the fit from a cubic smoothing spline, whereas the dashed line is the fit from a robust linear regression. Note that the same relationship is observed for indel rate versus long AR bp per 100 kb.  AR rates also exhibit a strong negative correlation with local GC content ( Figure 6; r = -0.55, P = 0.001). Indel rates show an overall positive correlation with GC content for the ENCODE regions (Figure 7), which illustrates that indel rates may be confounded by local GC content. In order to check the effect of GC content on indel rates, we recalculated the results presented in Table 2 including GC content as a confounder. For example, although indel events per 100 kb in AR sequence is observed to be about 7.9 (99% CI 6.7 to 9.2; see Table 2), the mean rates are about 4.7 (99% CI 3.5 to 6.4) and about 10.4 (99% CI 8.6 to 12.4) for AR sequence with GC content above 50% and GC content below 50%, respectively. However, the mean indel rates presented in Table 2 are not significantly altered when adjusted for local GC content at each annotational feature (data not presented).  AR sequence rate versus GC content Figure 6 AR sequence rate versus GC content. Ancestral repeat (AR) sequence rate is expressed as base pairs (bp) per 100 kilobases (kb). The reduced local GC content observed in AR sequence reflects the process of deamination of methylated CpG to TpG dinucleotides in vertebrate sequence over long evolutionary periods of time [3]. The solid line represents the fit from a cubic smoothing spline, whereas the dashed line is the fit from a robust linear regression. tures is broadly similar to SNP density. For example, as a percentage of their respective overall means, the indel rates for MCS evolutionary constraints of strict, moderate, and loose are 10%, 26% and 61%, compared with 29%, 43% and 55% for SNP rates. Similarly, the indel and SNP rates are reduced for many transcribed sequences (CDS, TSS, and RACEfrags).

AR sequence rate versus MCS modest
For some features, however, the pattern of constraint for indel and SNP rates differ quite markedly (  Table 5 compares indel rates by functional annotation for these data and the data presented by Bhangale and coworkers [20]. The overall indel rates are very similar for indel events (15 per 100 kb versus 13.8 per 100 kb for the data presented by Bhangale and coworkers [20]) and indel bp (43.4 bp per 100 kb versus 39.4 bp per 100 kb). The indel rates presented by Bhangale and coworkers [20] are also greatly reduced for coding DNA but not pseudo-exons or UTR sequence. Open chromatin indel rates are reduced in both datasets.

Discussion
This work represents the first systematic description of small insertion/deletion human polymorphism data in relation to functional and evolutionary annotation, which complements larger scale structural variation data across the genome [2,[21][22][23][24]. In order to understand the potential contribution made by indels to human genetic variation, we contrasted small indel rate variation by type of ENCODE region (manual or random selection), indel rates by functional annotation features, and indel rates by evolutionary constraint scores and neutral (AR) sequence; finally, we compared indel and SNP rates and their relative pattern of distribution across genomic features.
Overall, indel rates do not vary significantly between manual and randomly selected regions, suggesting that the ENCODE selection criteria for manual regions (the presence of well studied genes and availability of substantial comparative sequence) do not preclude similar genomic profiles for manual and random regions, with stratified randomly selected regions designed to be representative of a broad range of the genome [11]. Approximately 5% of the ENCODE sequence is estimated to be subject to moderate evolutionary constraint across mammalian species (Table 2), but only a minority of these constrained sequences are estimated to overlap with known protein coding exons and their associated UTRs (about 40%). The majority either overlap with known noncoding functional features (20%) or are suspected to be associated with previously unrecognized (40%) noncoding transcription [25].
As expected, coding (CDS, TSS, and RACEfrags) and constrained sequence (MCS) show the most constrained indel rates, followed by noncoding transcripts (transcriptionally active regions/transcribed fragments) and regulatory features (FAIRE sites, DHS, and HisPolTaf). To the extent that indels arise in functional sequence, in general indels appear to be subject to purifying selection, with indel rates negatively correlated with past evolutionary constraint across mammal Indel rates versus GC content and primate sequences (MCS human-mammal and GERP human-primate scores; Figures 1 and 2).
An apparent exception to the negative relationship between indel rates and constraint score is the HOXA cluster (ENCODE region 10), which runs counter to this trend. This region simultaneously exhibits the highest evolutionary constraint in the comparison of mammalian sequence (MCS) and the third highest indel rate for all the ENCODE regions (Figure 1). However, the HOXA cluster is in the centre of the region and is surrounded by gene deserts with limited evidence of evolutionary constraint. Hence, the explanation for this potentially counterintuitive observation is probably that the indel polymorphisms are largely confined to the gene deserts, whereas the constrained sequence is confined to the central portion of the HOXA cluster.
AR sequence rates are negatively correlated with mammalian sequence constraint (MCS; Figure 5), which is expected because AR sequence is neutral and not subject to natural selection. However, AR is not associated with GERP humanprimate and GERP human-mammal scores (data not shown), because AR sequence was defined and identified in relation to broad mammalian sequence comparisons and not specifically primate sequence.
Multi-species constrained scores for mammals (MCS modest) and GERP for human-primate comparisons are strongly negatively correlated ( Figure 6). The nonlinear relationship also reflects the fact that relatively recent (human-primate) sequence constraint comparisons fail to discriminate between the shared, more highly conserved sequences, which are only observed using broader phylogenetic comparisons.
Indel and SNP rates do not vary by the timing of DNA sequence replication during S-phase (the synthesis of DNA in preparation for mitosis) when classified as early, mid, and late S-phase replication timing [19].
Based upon two assumptions, we anticipated AR sequence rates to be positively correlated with indel rates across the ENCODE regions. What we in fact observe is a negative correlation between AR and indel rates ( Figure 3). This unexpected result initially suggests that one or both of the assumptions may be false. The first assumption is that AR sequence is effectively functionless (and therefore neutral sequence), and the second is that indel mutations arise at the same rate in AR sequence as elsewhere in the genome.
Although there is evidence that interspersed repeats in mammalian genomes may acquire functional roles as both proteincoding and transcriptional regulatory regions, only about 5% of the total amount of nonexonic constrained sequence (GERP) in the ENCODE regions is estimated to overlap with AR sequence [26]. This indicates that most AR sequence is still likely to be neutral and, for the most part, is unlikely to be subject to selection. By contrast, a lack of uniform indel mutation rates across the genome is more plausible [8]. Just as nucleotide point mutation rates [27] and segmental duplications [28] vary widely across the genome, it has been shown that the rate (and perhaps mechanism) of indel generation also varies widely across the genome [8].
Alternatively, the observed reduction in AR indel rates could in part arise from confounding caused by local GC content or experimental ascertainment bias. We observed AR rates to be negatively correlated with GC content (Figure 6; r = -0.55, P = 0.001) and ENCODE indel rates to be positively correlated with GC content (Figure 7). The overall indel event rate when adjusted for mean centered GC content remains unaltered, at 15 events per 100 kb (99% CI 13.4 to 16.7), whereas AR indel rates are about 4.7 events (99% CI 3.5 to 6.4) and about 10.4 events (99% CI 8.6 -12.4) for sequence with GC content above 50% and GC content below 50%, respectively. However, although indel rates are associated with local GC content, the latter only partly accounts for reduced AR indel rates, because the indel rate for AR with reduced GC content (10.4 per 100 kb for sequence with <50% GC content) is still lower than indel rates for bulk DNA (15 per 100 kb).
One final possibility for the observed indel rate reduction for these data in AR regions, could also be an artefact of the data generation process. Ascertainment bias could arise against AR sequence with common indels because of the nature of identifying ancestral repeats common to mammalian species. Indels arising in AR sequence would reduce the alignment score used to identify ancestral repeats, so that true AR containing indels would be less likely to be identified as AR. The problem could be exacerbated if indel rates are elevated in regions that have also experienced increased rates of indel changes throughout mammalian evolution (which is likely to be the case because lineage-specific rates of indel divergence between mammals are strongly correlated with genomic region). Both of these mechanisms would give rise to experimental artefact and apparent reduction in AR indel rates.
Some of the annotation features that show significantly reduced indel rates in this analysis also show reduced levels of nucleotide substitutions (for example, CDS, TSS, RACEfrags, and MCS), indicating that selective constraint is acting to both reduce SNP as well as indel density. However, other categories such as noncoding transcription sites (transcriptionally active regions/transcribed fragments, pseudo-exons) and chromatin regulatory elements, as assessed by DHS (activated cis-regulatory elements in mammalian genomes) and FAIRE sites, appear to show reduced indel rates, but not reduced SNP density. This observation is consistent with experimental data that show DNA regulation of nucleosome stability to be diffuse, cumulative across base pairs, and apparently on the scale of a single nucleosome, at about 200 bp [14]. In this context, indels may have important implications for understanding genome function and variation, because chromatin composition plays a central role in regulating all DNA templated processes, including transcription, recombination, repair, and replication.
There are two potential limitations of the present study. The first relates to the completeness and accuracy of the indel and genomic annotation data [19], ensuring which is a continuing exercise for coding and noncoding transcript features [29].
Although the complete accuracy of annotations is essential to the future success of genomic and complex trait research [30], in this study we have deliberately taken a conservative statistical approach to investigating the distribution of indels in relation to annotation features, in order to account for inherent uncertainty, both in terms of biology and experimental measurement error. We also used independent data from Bhangale and coworkers [20] to compare indel rates by functional annotation and evolutionary constraint. Table 5 shows similar overall rates and reduction in indel rates for coding sequence (CDS, TSS, and RACEfrags), but not for pseudo-exons or UTRs. The data from Bhangale and coworkers also show reduced rates for open chromatin features (FAIRE and hypersensitive sites).
The second potential limitation is that, for most of our analyses, we have used summary measures for each ENCODE region, and it is likely that some effects of interest in small sequences will therefore be overlooked. Nevertheless, relatively crude summary measures by region and annotation feature still reveal clear trends between indel rates and indirect (experimental and computational) measures of functional and evolutionary constraint. We assessed the robustness of our results to various potential biases by conducting several sensitivity analyses. For instance, some of the encode regions (ENm010, ENm013, ENm014, ENr112, ENr113, ENr123, ENr131, ENr213, ENr232, and ENr321) were genotyped more intensively than others, but we found no evidence that these regions yielded substantially different results in our analyses. Indels are predominantly (58%) 1 bp in length, and we repeated analyses with only those indels with lengths in excess of 1 bp, and found that the trends in our analysis do not substantially alter (data not shown). We also repeated the analyses for insertions and deletions separately and reached the same conclusions.

Conclusion
Small indels that arise in functional sequence are likely to be subject to negative selection, as shown by the reduced indel rates in transcribed DNA, evolutionarily constrained sequence, and -to a lesser extent -regulatory elements. Although reduced indel and SNP rates are both clearly related to coding sequence constraints, constrained indel rates in regulatory regions may reflect that indels are more likely than SNPs to moderate the structural function of regulatory elements. Indels may play a more important role than SNPs in contributing to natural genetic variation at regulatory sites, and hence they could be an important source of variation in gene expression levels.

Materials and methods
The ENCODE project aims to identify and catalog all functional elements, including coding sequences of genes and noncoding DNA, in the human genome. A pilot study phase considered 44 discrete regions that encompass 30 megabases, or about 1% of the human genome, with 14 of these regions (about 15 megabases) selected manually and the remainder randomly [11]. Small indels in the ENCODE regions were called from shotgun re-sequencing reads and traces of the SNP discovery efforts from both the SNP consortium and the HapMap (see the report from the ENCODE Project Consortium [19] for details of discovery and validation procedures). The shotgun technology used identified indels with a maximum length of 20 bp. Whole genome sequence data were generated totalling onefold coverage of the human genome from DNA derived from a pool of cell lines from eight unrelated adult African Americans (four male and four female) enrolled in Houston, Texas, USA [31]. The SSA-HADIP software package, a modification of SSAHASNP [32], was used to align these reads to build 35 of the human reference sequence, generating polymorphism calls, while keeping track of the total bases aligned for each read. In brief, the neighbourhood quality standard base alignment method was adapted to identify indels by requiring the inserted/deleted bases and the flanking five bases on either side of the indel to exceed a minimum Phred quality score of 22. If these minima were not met, then the indel was not reported.
For this study, indels and SNPs were called using the eight Baylor samples in order to facilitate comparison. Only validated SNPs (those with heterozygosity scores) were used.
As part of the HapMap project [33], ten ENCODE regions had in-depth SNP discovery by polymerase chain reaction resequencing on 48 individuals in four populations; this dataset represents the deepest multi-megabase resequencing data currently available and is about three times as dense as phases I and II of the HapMap project. Experimental data, sequence conservation, and feature definitions were obtained from the three experimental groups of the ENCODE Consortium and the multiple sequence analysis group [19]. All data used in our work are available at the ENCODE project at University of California, Santa Cruz [34] and are available from the corresponding author. The full set of ENCODE indels can be downloaded directly [35].
To evaluate the potential contributions of insertion and deletion events to functional variation, we calculated indel density as a percentage of nucleotides for each ENCODE region and classification feature. Genomic coordinates for features and ENCODE regions were used to estimate two summary measures of indel density: the number of indels per 100 kb of the region length or total feature, and the number of indel bp per 100 kb for the region length or total feature. The densities were analyzed using a negative binomial model with the number of indels or base pairs as the response, the lengths of sequence as an offset, and data aggregated to the region level [36]. This approach allowed us to calculate 99% confidence intervals for indel and SNP densities, compensated for potential over-dispersion, and provided a conservative framework for testing for differences between manually and randomly selected regions and genomic features. Comparisons across genomic features are also likely to be conservative, because confidence intervals were generated using aggregate summary measures across ENCODE regions, rather than raw data. The SNP densities correspond to a measure of heterozygosity (1 SNP per 100 kb, corresponds to a heterozygosity of 1 × 10 -5 ).
Comparative sequence analysis has become a key bioinformatics tool for identifying noncoding functional DNA [37]. We use two derived scores that attempt to measure the relative evolutionary constraint of DNA sequences: the lengths of MCSs determined from the multiple sequence alignments comparing human with 13 species of mammal [25], and rejected substitution or GERP scores from comparisons between humans and primates [26].
Evolutionarily constrained sequences were identified using three independent sequence conservation constraint programs (binCons, phastCons, and GERP) for three different multiple sequence alignments generated using TBA, MLA-GAN, and MAVID for 14 mammalian species [25]. Each alignment used human sequence as the reference. Three levels of MCS were defined: strict, in which sequences are constrained in all alignment/conservation combinations; moderate, in which sequences are constrained in at least two of three alignments, and from two of three conservation programs; and loose, in which sequences are constrained in at least one alignment/conservation combination. We found results did not alter qualitatively between use of the three scores and we present results using moderate MCS in this report. Note that we refer to 'constrained' rather than 'conserved' sequence because conservation per se does not imply function, whereas constraint does. GERP identifies regions at high resolution that exhibit nucleotide substitution deficits, and measures these deficits as 'rejected substitutions'. Rejected substitutions reflect the intensity of past purifying selection and are used to rank and characterize constrained elements. GERP scores are positive in constrained regions and negative in neutral DNA [26], and MCS scores are high in constrained regions and low in neutral DNA [25].
To illustrate potential relationships between indel rates and constraint scores, summary data for the 44 ENCODE regions were plotted using cubic smoothing splines and robust linear regression using an M estimator [36]. The latter approach is robust to potential outliers but conservative. Potential outliers were also identified using standard regression leverageresidual diagnostics [36], and we assessed the sensitivity of results to outlier removal using Pearson product moment correlation coefficients (r) and adjusted linear R 2 statistics (multiple correlation coefficient R).