Technical Comments

Response to Comment on “Evidence of Abundant Purifying Selection in Humans for Recently Acquired Regulatory Functions”

Science  10 May 2013:
Vol. 340, Issue 6133, pp. 682
DOI: 10.1126/science.1233366

Abstract

Green and Ewing propose corrections to our methodology, which we incorporate and extend here. The improved methodology supports our initial conclusion of extensive lineage-specific constraint concentrated in ENCODE elements. We clarify that our estimate is dependent on the constrained and neutral references used, which can further increase the number of nucleotides involved, because a particularly stringent definition was initially used.

In our initial Report (1), we found reduced genetic diversity at noncoding genomic regions that have not been conserved across mammals but are biochemically active, suggesting that some fraction of these regions has experienced lineage-specific purifying selection, and proposed a method for estimating the proportion under human constraint (PUC). Green and Ewing (2) suggest that our PUC estimate is inflated by several technical artifacts and that a reduction in SNP density cannot be reliably distinguished from a reduction in mutation rate. In this Response, we incorporate and build upon their improved method by including additional quality filters, and we repeat our analysis, confirming the validity of each of their proposed corrections, but demonstrating continued support for our original conclusions. We address each of the corrections raised below.

First, Green and Ewing point out that masking CpG nucleotides formed by both reference and alternate alleles can lead to lower density measurements in GC-rich features and that including positions where the reference genome contains the derived allele leads to higher density measurements and higher derived allele frequency (DAF) measurements in less-constrained regions. We confirm that incorporating both corrections modestly decreases our density-based total estimate of recent constraint.

Second, they show that density-based estimates can be biased by variable mutation rate due to nucleotide composition and regional effects. However, we note that our analysis comparing heterozygosity levels at bound and unbound regulatory motifs would not be affected by such mutation rate differences, because the nucleotide composition is the same between bound and unbound instances. Moreover, the remaining results in our original paper are supported by DAF, which is not sensitive to mutation rate.

Third, they point out that the allele frequency estimates of the 1000 Genomes Project low-coverage pilot data are biased by sequencing depth, causing rarer alleles to be more efficiently called in higher-coverage regions, and thus lowering the mean DAF of higher-coverage regions. Because higher GC content is associated with both higher sequencing depth and with ENCODE-defined active regions, this could potentially lead to artificially lower DAF for ENCODE regions. Green and Ewing proposed a procedure for correcting this effect by calculating relative constraint within bins of equal sequencing depth, which led to a reduced, albeit still substantial, estimate of human constraint.

Because the rare part of the allele frequency spectrum (DAF < 2%) ought to show the strongest signal of selection, Green and Ewing also specifically focused on the fraction of variants that are rare. However, rare variants are also the most prone to genotype errors. To address these potential genotyping errors, we have extended Green and Ewing’s methodology by including a quality filter that requires at least 50 of the sampled individuals to have genotype calls from all three sequencing centers and exclude the top and bottom 5% of SNPs by total read depth to avoid technical artifacts. Finally, we use mammalian-conserved and unconserved regions as the reference points for DAF within each coverage bin, which are much more abundant annotations than nondegenerate conserved protein-coding nucleotides, to reduce sampling error.

We refer to the resulting value as coverage-corrected relative constraint (ccRC) to emphasize that it can take on both positive and negative values and that it represents an aggregate measure of overall constraint, rather than a partitioning between constrained and nonconstrained bases. The ccRC values of a test feature and two control features can then be used to produce a PUC value as described previously.

The resulting ccRC values (Table 1) confirm our original observation of extensive lineage-specific purifying selection concentrated at ENCODE elements, especially at regulatory motifs bound by their cognate proteins (3, 4) and at enhancers defined by histone modification patterns (5, 6). Moreover, the signal of selection is apparent both in the mean DAF and in the fraction of alleles with DAF less than 2% (Fig. 1), the additional criterion proposed by Green and Ewing, consistent with an increased accuracy at these rare sites.

Table 1

ccRC based on DAF at annotated genomic features. Values are in arbitrary units scaled between noncoding unconserved (zero) and noncoding conserved (unity).

View this table:
Fig. 1

Fraction of SNPs with DAF < 2% in annotated features, binned by depth of coverage. Blue, unconserved noncoding non-ENCODE; red, unconserved noncoding ENCODE; black, conserved noncoding.

Lastly, we stress that our original PUC estimates had made the assumption that constrained nucleotides have the same average selection coefficient as conserved coding positions, specifically in nondegenerate sites, which are the most constrained class of nucleotides we observe. These PUC estimates can change and become considerably higher if instead we model the constraint using different reference annotations.

For example, estimating PUC using the new ccRC values and our original report’s constrained and unconstrained references leads to an updated estimate of 51 Mb of the mappable non–mammalian-conserved genome being under human constraint (1.6% of the entire genome). When the constrained reference is the average of all mammalian-conserved regions, including both coding and noncoding elements, this estimate rises to 157 Mb (5.1% of the entire genome). We can also make a PUC estimate that quantifies the difference in constraint between ENCODE and non-ENCODE regions in the unconserved noncoding genome.

Using this approach, we can directly estimate the fraction of ENCODE nucleotides that are likely to be constrained at the level of mammalian-conserved noncoding elements, specifically testing noncoding ENCODE elements outside mammalian-conserved regions. This results in 129 Mb (11.5% of the unconserved noncoding ENCODE nucleotides) having the same level of constraint as mammalian-conserved noncoding nucleotides (our constrained reference), and the remaining 994 Mb (88.5%) having the same level of constraint as nonconserved noncoding regions outside ENCODE elements (the unconstrained reference).

It important to note that each of these estimates is a simplification of a much more complex picture, because each nucleotide has its own selection coefficient that we do not have power to estimate directly, and human constraint is due to a mixture of multiple classes of elements (not only two) that we do not seek to distinguish in our analysis.

In summary, we agree with the proposed corrections from Green and Ewing, which will also be important in interpreting related population genomics work estimating constraint at regulatory elements (79). However, these corrections should be coupled with the additional quality filters applied here when looking at low-coverage data. Our ccRC estimates, which incorporate both corrections, support our original conclusion of extensive lineage-specific constraint on regulatory elements and continue to be consistent with previous estimates by other groups (10, 11).

References and Notes

  1. Acknowledgments: We thank P. Green and B. Ewing for numerous suggestions, for sharing their manuscript prepublication, and for their helpful correspondence.
View Abstract

Subjects

Navigate This Article