Skip to main content

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

Abstract

Background

We describe the distribution of indels in the 44 Encyclopedia of DNA Elements (ENCODE) regions (about 1% of the human genome) and evaluate the potential contributions of small insertion and deletion polymorphisms (indels) to human genetic variation. We relate indels to known genomic annotation features and measures of evolutionary constraint.

Results

Indel rates are observed to be reduced approximately 20-fold to 60-fold in exonic regions, 5-fold to 10-fold in sequence that exhibits high evolutionary constraint in mammals, and up to 2-fold in some classes of regulatory elements (for instance, formaldehyde assisted isolation of regulatory elements [FAIRE] and hypersensitive sites). In addition, some noncoding transcription and other chromatin mediated regulatory sites also have reduced indel rates. Overall indel rates for these data are estimated to be smaller than single nucleotide polymorphism (SNP) rates by a factor of approximately 2, with both rates measured as base pairs per 100 kilobases to facilitate comparison.

Conclusion

Indel rates exhibit a broadly similar distribution across genomic features compared with SNP density rates, with a reduction in rates in coding transcription and evolutionarily constrained sequence. However, unlike indels, SNP rates do not appear to be reduced in some noncoding functional sequences, such as pseudo-exons, and FAIRE and hypersensitive sites. We conclude that indel rates are greatly reduced in transcribed and evolutionarily constrained DNA, and discuss why indel (but not SNP) rates appear to be constrained at some regulatory sites.

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

Table 1 Indel density (for all 44 ENCODE regions)
Table 2 Indel density for annotation features (across all 44 ENCODE regions)
Table 3 Comparison of indel and SNP density by ENCODE experimental features

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 (sdnum/100 kb = 5.0 number of indels per 100 kb; sdbp/100 kb = 14.7 indel bp per 100 kb, where sdnum/100 kb and sdbp/100 kb refer to the standard deviation for number of indels and indel bp per 100 kb, respectively) to those in random regions (sdnum/100 kb = 4.0; sdbp/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 exhibit a significant reduction in indel rates, as assessed by identifying open reading frames (coding sequence [CDS] mean indel rate: 0.7 bp per 100 kb), transcription start sites (TSSs; 3.3 bp per 100 kb), rapid amplification of cDNA ends fragments (RACEfrags; 6.6 bp per 100 kb), and transcribed fragments (12.3 bp per 100 kb). Pseudo-exons (19.1 bp per 100 kb), 3' UTRs (23.6 bp per 100 kb), 5' UTRs (27.4 bp per 100 kb), and transcripts of unknown function (36.9 bp per 100 kb) all exhibit a reduction in indel rates compared with the overall mean for all ENCODE sequence, but these findings are not statistically significant.

Potential regulatory elements, assessed by measuring open chromatin sites, also reveal sequences with constrained indel rates (Table 2). Formaldehyde assisted isolation of regulatory elements (FAIRE) sites (23.8 bp per 100 kb) and DNAse hypersensitive sites (DHS; [NHGRI group] 19.7 bp per 100 kb and [Regulome group] 27.0 bp per 100 kb) both exhibit reduced indel rates. DHS are short regions of DNA that are relatively easily cleaved by DNAse I.

Acetylated histones are usually associated with transcriptionally active chromatin and deacetylated histones with inactive chromatin. Hence, histone modified regions often signify regulatory sites. Selected histone modifications and binding sites for RNA polymerase II and the general transcription factor TAF250 were assayed for the ENCODE regions (see ENCODE Project Consortium [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.

Table 4 Experimental feature definitions

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) also showed unexpectedly reduced indel rates. Cell cycle replicating segments (MidRepSeg: 43.2 bp per 100 kb) show no relationship with indel rates.

Figures 1 to 3 show the relationship between indel base pairs per 100 kb and measures of mammalian evolutionary constraint, human-primate evolutionary constraint, and AR rates, with each data point representing a summary score for each ENCODE region. The Pearson correlation coefficients relating to Figures 1 to 3 are statistically insignificant when all of the ENCODE region summary data points are considered. However, when outlying data points are identified and excluded using standard regression diagnostics, the correlations are of marginal statistical significance. Indel rates are (nonsignificantly) inversely correlated with mammalian MCS score (Figure 1; r = -0.25, P = 0.11 with outlier ENCODE region 10 excluded), and negatively associated with the primate genomic evolutionary rate profiling (GERP) score and GERP squared using multiple regression (Figure 2; multiple correlation coefficient: R = 0.32, P = 0.04). Indel rates are also observed to be marginally and negatively correlated with AR rates and AR squared (Figure 3; multiple correlation coefficient: R = -0.30, P = 0.06 with regions 8 and 15 identified as outliers).

Figure 1
figure 1

Indel rate versus MCS modest for human and 13 mammals. Indel rate and multi-species constrained sequences (MCS modest) 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.

Figure 2
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.

Figure 3
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 (bp per 100 kb) are strongly inversely correlated with MCS (Figure 4; r = -0.46, P < 0.002), but exhibit no relation with either human-primate or human-mammal GERP scores (plots not shown; GERP primate: r = 0.02, P = 0.91; GERP mammal: r = -0.03, P = 0.8). MCS and GERP constraint scores are positively correlated with one another in a curvilinear relationship (Figure 5; r = 0.42, P = 0.005), with the homeobox gene family HOXA cluster, ENCODE region 10, identified as a highly conserved outlier region on the MCS but not an outlier on either of the GERP scores.

Figure 4
figure 4

AR sequence rate versus MCS modest. Ancestral repeat (AR) sequence rate and multi-species conserved sequences (MCS modest) 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.

Figure 5
figure 5

MCS modest versus GERP human-primate score. Multi-species conserved sequences (MCS modest) 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.

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

Figure 6
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.

Figure 7
figure 7

Indel rates versus GC content. 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.

Table 3 compares the distribution of indel and validated SNP rates by experimental feature. In general, indel rates are lower than SNP rates, with a ratio of validated SNPs to indel event rates of 6.7 (102.4/15), or 2.4 (102.4/43.4) for validated SNPs:indel bp. The pattern of indel rates across genomic features 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).

For some features, however, the pattern of constraint for indel and SNP rates differ quite markedly (Table 3). Although indel rates are constrained in chromatin mediated transcription regulatory sites (FAIRE: 23.8 bp per 100 kb; DHS: 19.7 to 27.0 bp per 100 kb), SNP rates are not constrained for these features (FAIRE: 90 SNPs per 100 kb; DHS: 90 to 96 SNPs per 100 kb) as compared with the overall mean (102.4 SNPs per 100 kb).

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.

Table 5 Comparison of ENCODE and Bhangale et al. (ten ENCODE regions) indel data

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–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].

Small indels are common and constitute approximately 15 insertions/deletions every 100 kb or, in terms of sequence length, 43 bp per 100 kb of the genome. The number of validated common SNPs is observed to be about seven times the number of small indels (indels per 100 kb) or twice the observed indel bp rate (bp per 100 kb). Indel rates are greatly reduced in regions associated with known functionality (largely coding DNA) and under evolutionary constraint. Compared with the overall mean, indel event rates are reduced by factors of about 20 for exon coding regions, about 5 for strict MCS sequence, and about 2 for measures of chromatin mediated regulatory sites. These observations are consistent with estimates from other studies [1, 2, 8]. The corresponding reduction in indel rates for these data compared with bulk DNA and when measured as indel bp per 100 kb rather than indel events, about 60 (CDS), about 10 (strict MCS), and about 2 (FAIRE and DHS).

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 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 human-primate 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 protein-coding 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 SSAHADIP 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 re-sequencing 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, MLAGAN, 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 leverage-residual diagnostics [36], and we assessed the sensitivity of results to outlier removal using Pearson product moment correlation coefficients (r) and adjusted linear R2 statistics (multiple correlation coefficient R).

Abbreviations

AR:

ancestral repeat

bp:

base pairs

CDS:

coding sequence

CI:

confidence interval

DHS:

DNAse hypersensitive sites

ENCODE:

Encyclopedia of DNA Elements

FAIRE:

Formaldehyde assisted isolation of regulatory elements

GERP:

genomic evolutionary rate profiling

HOXA:

HOXA cluster (homeobox gene family)

kb:

kilobases

MCS:

multi-species constrained sequence

RACEfrag:

rapid amplification of cDNA ends fragment

SNP:

single nucleotide polymorphism

TSS:

transcription start site

UTR:

untranslated region.

References

  1. Mills RE, Luttig CT, Larkins CE, Beauchamp A, Tsui C, Pittard WS, Devine SE: An initial map of insertion and deletion (INDEL) variation in the human genome. Genome Res. 2006, 16: 1182-1190.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  2. Bhangale TR, Rieder MJ, Livingston RJ, Nickerson DA: Comprehensive identification and characterization of diallelic insertion-deletion polymorphisms in 330 human candidate genes. Hum Mol Genet. 2005, 14: 59-69.

    Article  PubMed  CAS  Google Scholar 

  3. Strachan T, Read AP: Human Molecular Genetics. 2006, Garland Science/Taylor & Francis Group, ,

    Google Scholar 

  4. Wasserman WW, Sandelin A: Applied bioinformatics for the identification of regulatory elements. Nat Rev Genet. 2004, 5: 276-287.

    Article  PubMed  CAS  Google Scholar 

  5. Initial sequence of the chimpanzee genome and comparison with the human genome. Nature. 2005, 437: 69-87.

  6. Weber JL, David D, Heil J, Fan Y, Zhao C, Marth G: Human diallelic insertion/deletion polymorphisms. Am J Hum Genet. 2002, 71: 854-862.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Lunter G, Ponting CP, Hein J: Genome-wide identification of human functional DNA using a neutral indel model. PLoS Comput Biol. 2006, 2: e5-

    Article  PubMed  PubMed Central  Google Scholar 

  8. Snir S, Pachter L: Phylogenetic Profiling of Insertions and Deletions in Vertebrate Genomes. Research in Computational Molecular Biology. 2006, Berlin, Heidelberg: Springer, 265-280.

    Chapter  Google Scholar 

  9. Kuroki Y, Toyoda A, Noguchi H, Taylor TD, Itoh T, Kim DS, Kim DW, Choi SH, Kim IC, Choi HH, et al: Comparative analysis of chimpanzee and human Y chromosomes unveils complex evolutionary pathway. Nat Genet. 2006, 38: 158-167.

    Article  PubMed  CAS  Google Scholar 

  10. Li WH, Saunders MA: News and views: the chimpanzee and us. Nature. 2005, 437: 50-51.

    Article  PubMed  CAS  Google Scholar 

  11. ENCODE Project Consortium: The ENCODE (ENCyclopedia Of DNA Elements) Project. Science. 2004, 306: 636-640.

    Article  Google Scholar 

  12. ENCODE project target selection process and target regions. [http://genome.gov/10506161]

  13. Eddy SR: Non-coding RNA genes and the modern RNA world. Nat Rev Genet. 2001, 2: 919-929.

    Article  PubMed  CAS  Google Scholar 

  14. Giresi PG, Gupta M, Lieb JD: Regulation of nucleosome stability as a mediator of chromatin function. Curr Opin Genet Dev. 2006, 16: 171-176.

    Article  PubMed  CAS  Google Scholar 

  15. Kellis M, Patterson N, Endrizzi M, Birren B, Lander ES: Sequencing and comparison of yeast species to identify genes and regulatory elements. Nature. 2003, 423: 241-254.

    Article  PubMed  CAS  Google Scholar 

  16. Maquat LE: Nonsense-mediated mRNA decay: splicing, translation and mRNP dynamics. Nat Rev Mol Cell Biol. 2004, 5: 89-99.

    Article  PubMed  CAS  Google Scholar 

  17. Pagani F, Baralle FE: Genomic variants in exons and introns: identifying the splicing spoilers. Nat Rev Genet. 2004, 5: 389-396.

    Article  PubMed  CAS  Google Scholar 

  18. Waterston RH, Lindblad-Toh K, Birney E, Rogers J, Abril JF, Agarwal P, Agarwala R, Ainscough R, Alexandersson M, An P, et al: Initial sequencing and comparative analysis of the mouse genome. Nature. 2002, 420: 520-562.

    Article  PubMed  CAS  Google Scholar 

  19. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature. 2007, 447: 799-816.

  20. Bhangale TR, Stephens M, Nickerson DA: Automating resequencing-based detection of insertion-deletion polymorphisms. Nat Genet. 2006, 38: 1457-1462.

    Article  PubMed  CAS  Google Scholar 

  21. Conrad DF, Andrews TD, Carter NP, Hurles ME, Pritchard JK: A high-resolution survey of deletion polymorphism in the human genome. Nat Genet. 2006, 38: 75-81.

    Article  PubMed  CAS  Google Scholar 

  22. Locke DP, Sharp AJ, McCarroll SA, McGrath SD, Newman TL, Cheng Z, Schwartz S, Albertson DG, Pinkel D, Altshuler DM, et al: Linkage disequilibrium and heritability of copy-number polymorphisms within duplicated regions of the human genome. Am J Hum Genet. 2006, 79: 275-290.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  23. Redon R, Ishikawa S, Fitch KR, Feuk L, Perry GH, Andrews TD, Fiegler H, Shapero MH, Carson AR, Chen W, et al: Global variation in copy number in the human genome. Nature. 2006, 444: 444-454.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  24. Wong KK, deLeeuw RJ, Dosanjh NS, Kimm LR, Cheng Z, Horsman DE, MacAulay C, Ng RT, Brown CJ, Eichler EE, et al: A comprehensive analysis of common copy-number variations in the human genome. Am J Hum Genet. 2007, 80: 91-104.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  25. Margulies EH, Cooper GM, Asimenos G, Thomas DJ, Dewey CN, Siepel A, Birney E, Keefe D, Schwartz AS, Hou M, et al: Analyses of deep mammalian sequence alignments and constraint predictions for 1% of the human genome. Genome Res. 2007, 17: 760-774.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  26. Cooper GM, Stone EA, Asimenos G, Green ED, Batzoglou S, Sidow A: Distribution and intensity of constraint in mammalian genomic sequence. Genome Res. 2005, 15: 901-913.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  27. Ellegren H, Smith NG, Webster MT: Mutation rate variation in the mammalian genome. Curr Opin Genet Dev. 2003, 13: 562-568.

    Article  PubMed  CAS  Google Scholar 

  28. Bailey JA, Eichler EE: Primate segmental duplications: crucibles of evolution, diversity and disease. Nat Rev Genet. 2006, 7: 552-564.

    Article  PubMed  CAS  Google Scholar 

  29. Guigo R, Flicek P, Abril JF, Reymond A, Lagarde J, Denoeud F, Antonarakis S, Ashburner M, Bajic VB, Birney E, et al: EGASP: the human ENCODE Genome Annotation Assessment Project. Genome Biol. 2006, 7 (Suppl 1): S2-S31.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Reese MG, Guigo R: EGASP: introduction. Genome Biol. 2006, S1-S3. Suppl 1

  31. A haplotype map of the human genome. Nature. 2005, 437: 1299-1320.

  32. Ning Z, Cox AJ, Mullikin JC: SSAHA: a fast search method for large DNA databases. Genome Res. 2001, 11: 1725-1729.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  33. The International HapMap Project. [http://www.hapmap.org/]

  34. The ENCODE project at UCSC. [http://genome.ucsc.edu/cgi-bin/hgTables]

  35. The ENCODE project at UCSC. [http://hgdownload.cse.ucsc.edu/goldenPath/hg17/encode/database/encodeIndels.txt.gz]

  36. Venables WN, Ripley BD: Modern Applied Statistics with S. 2002, New York: Springer

    Book  Google Scholar 

  37. Hardison RC: Conserved noncoding sequences are reliable guides to regulatory elements. Trends Genet. 2000, 16: 369-372.

    Article  PubMed  CAS  Google Scholar 

  38. Giresi PG, Kim J, McDaniell RM, Iyer VR, Lieb JD: FAIRE (Formaldehyde-Assisted Isolation of Regulatory Elements) isolates active regulatory elements from human chromatin. Genome Res. 2007, 17: 877-885.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

Download references

Acknowledgements

We wish to thank all those within the ENCODE Consortium who have assisted directly and indirectly with this work, and Tushar Bhangale for providing data. TC was supported by the Medical Research Council (UK). TA is supported by the Department of Health.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Taane G Clark.

Additional information

Authors' contributions

TC and TA analysed the indel data and wrote the manuscript. GC generated the evolutionary constraint scores and JM generated the shotgun indel data. EM, GC, JM and DB provided detailed advice and commented on the text. All authors read and approved the final version of the manuscript.

Taane G Clark, Toby Andrew contributed equally to this work.

Authors’ original submitted files for images

Rights and permissions

Reprints and permissions

About this article

Cite this article

Clark, T.G., Andrew, T., Cooper, G.M. et al. Functional constraint and small insertions and deletions in the ENCODE regions of the human genome. Genome Biol 8, R180 (2007). https://doi.org/10.1186/gb-2007-8-9-r180

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/gb-2007-8-9-r180

Keywords