Pervasive Hitchhiking at Coding and Regulatory Sites in Humans

1The neutral theory of molecular evolution postulates that adaptive substitutions occur so rarely that they can be safely ignored in most studies in population genetics or molecular evolution. This view has dominated the field of molecular evolution for the past 40 years. However, the past 4–6 years have seen a strong challenge to this view. This challenge comes not only from numerous studies detailing specific cases of molecular adaptation in a number of organisms (for example, see –) but also, and most compellingly, from a number of studies that indicate that adaptation might be common on the genomic scale –.

2High rates of adaptation on the genomic scale have been inferred from the excess of substitutions in functional regions relative to neutral expectations. The neutral expectations are derived from the polymorphism data at functional and putatively neutral sites and the divergence at the neutral sites using the reasoning of the McDonald-Kreitman (MK) test. The excess in the number of substitutions at functional sites over this expectation can be used to estimate the number of adaptive substitutions,. McDonald-Kreitman approaches can be modified to account for the presence of deleterious polymorphisms in the sample and the effects of demographic processes on polymorphism,,. The approach can also be extended to estimate rates of adaptation in regulatory regions,.

3McDonald-Kreitman analysis indicates that adaptive evolution in functional regions might be common in a range of organisms. In Drosophila, it has been estimated that from 30 to 60% of amino acid substitutions and ∼20% of substitutions in non-coding regions are adaptive,,, –. The rate appears similarly high in E. coli (>56% of amino acid substitutions are adaptive) but not in Arabidopsis (0–5% of amino acid substitutions are adaptive) and yeast.

4In humans, McDonald-Kreitman-based estimates have varied from zero to ∼35% of all amino acid substitutions being adaptive, –. A recent estimate by Boyko et al used information from the allele spectra of nonsynonymous and synonymous SNPs in human genes and the divergence with chimpanzee orthologs to estimate that ∼10% of amino acid substitutions between humans and chimpanzees have been fixed by positive selection. Thus, some of these studies suggest that adaptation might be fairly common in humans, although probably substantially less common than in Drosophila or E. coli.

5McDonald-Kreitman approaches are very powerful at detecting positive selection, however, they can be misleading for a variety of reasons,,. For example, if the strength of purifying selection over the evolutionary period separating two species has been different than it is in the present, McDonald-Kreitman-based approaches can either over- or underestimate the rate of adaptive evolution. As these estimates do not provide consistent answers about the prevalence of adaptation in humans and because they can be misleading under plausible demographic scenarios, reaching more reliable conclusions about the importance of adaptations in humans requires the investigation of other signatures of positive selection.

6An adaptive substitution reduces the level of polymorphism at neutral sites in its vicinity in a phenomenon known as a selective sweep. The width of the region in which the polymorphism is reduced is inversely proportional to the local recombination rate and directly proportional to the selection coefficient associated with the adaptive substitution –. The reduction of polymorphism is transient and the levels of polymorphism are expected to recover within roughly Ne generations. In addition to the reduction of the level of polymorphism, recurrent selective sweeps may also generate other signatures such as (i) an overabundance of low-frequency alleles,, (ii) a greater proportion of high-frequency derived alleles,, (iii) unusual haplotype structures,.

7A number of these expectations have been used to define signatures of positive selection for genome-wide scans for recent adaptation in humans: i.e., the detection of candidate regions that are likely to be experiencing a selective sweep at present or that have experienced one recently. For example, Nielsen et al. and Kelley et al used the deviation of the allele frequency spectrum from its background characteristics to detect candidate regions that may have experienced a sweep; several other methods have used summaries of haplotype structure and their deviation from the background to detect candidate regions that are undergoing a selective sweep,.

8Genomic scans for positive selection are primarily used to choose candidate regions for future investigation, but their application to the quantification of positive selection or even the establishment of its prevalence is problematic. To quantify the extent of positive selection based on the deviations of these signatures from the background requires a prior expectation about the likelihood of observing them under neutrality. These expectations, however, may be sensitive to the effects of non-equilibrium demography, –. As a result, it is difficult to generate robust a priori expectations for these statistics under neutrality. Therefore, scans for positive selection do not, by themselves, provide reliable quantification of the extent of positive selection in humans or establish that positive selection is prevalent in humans.

9To evaluate whether selective sweeps are common in the human genome, we require signatures that are unlikely to be generated by demography alone. The effects of recurrent selective sweeps (RSS) should be stronger in the regions of lower recombination and in regions of more frequent and selectively consequential adaptation. In Drosophila, for example, the level of neutral polymorphism is positively correlated with the recombination rate – and negatively with the rate and number of nonsynonymous substitutions in a region,,. These correlations, which are expected under models of RSS but should not be generated by demography alone, support the notion of high rates of adaptation in these taxa.

10Despite several compelling examples of adaptations, clear genome-wide signatures of RSS have been difficult to detect in humans. A relationship of diversity and recombination has been reported, but was attributed primarily to an association between recombination and mutation processes rather than to the effects of selected at linked sites –, with the possible exception of telomeric and centromeric regions. In turn, the relationships between levels of polymorphism and functional divergence have not yet been examined.

11If the recent MK estimates of the rate of adaptive evolution are correct and approximately 10% of amino acid substitutions are adaptive, we should expect to see a substantial number of recent selective sweeps in the polymorphism data. Indeed, ∼7×104 amino acid differences between human and chimpanzee proteins have accumulated over the past ∼14 million years. If 10% of these have been adaptive, then we can estimate that ∼7×103 adaptive amino acid substitutions have taken place over ∼14 million years. Assuming a constant rate of adaptation, this translates into ∼100 adaptive amino acid substitutions that occurred during the past Ne generations (Ne = ∼2×105 years). Moreover, if regulatory adaptations are common as well, then hundreds of recent selective sweeps should be detectable in the human polymorphism data.

12With these considerations in mind, we analyze genomic patterns of nucleotide polymorphism, recombination, functional density and functional divergence in humans using two independent, genome-wide SNP datasets. Consistent with the expectations of positive selection, we detect a positive correlation between levels of neutral polymorphism and recombination rate and a negative correlation between levels of nucleotide polymorphism and both functional density and functional divergence. These correlations remain intact after controlling for a number of possible covariates. The evidence is consistent with positive selection in both regulatory and protein-coding regions. We consider alternative explanations for these findings and argue that, in addition to recurrent selective sweep, only background selection (BS) (loss of neutral variants due to hitchhiking with linked deleterious mutations) can possibly generate most of these patterns. Hitchhiking of neutral polymorphisms with linked selected variants—either due to recurrent positive selection or background selection or possibly both—appears to be a substantial force determining levels of neutral polymorphism in the human genome.

13To study the effects of RSS, we separate the genomic sequences into two mutually exclusive sets of sequences: “functional” (genic and regulatory) and “nonfunctional”. Both sets of sequences are taken only from the internal parts of autosomes; specifically, we remove all sequences located within 10 Mbp of a telomere or a centromere. We further remove all sequences that cannot be aligned with the chimpanzee genome. The functional set is composed of several types of sequences (see Material and Methods). First, it contains all the genic regions, specifically those that (i) encode exons or are located within 1 kb of any predicted exon and (ii) are located within 5 kb from the starting and ending position of transcripts of protein-coding genes. Because many functional, noncoding sequences are located far from genes in the human genome –, we also take all the sequences that can be aligned between primates and zebrafish; sequences that can be aligned over such large evolutionary distances are very unlikely to be unconstrained (see Materials and Methods). The nonfunctional set contains all other sequences except for the repetitive sequences that are filtered out using RepeatMasker. We remove repetitive regions because both alignment and SNP discovery are more problematic in such regions. Hereafter, we will refer to the sequences in the primarily nonfunctional set (totaling ∼1,080 Mbp) as “neutral” sequences for brevity.

14We use two SNP datasets: (i) ∼1.2 million Perlegen “A” SNPs discovered using Perlegen chip technology in a panel of 71 individuals of mixed ancestry and (ii) ∼2.0 million SNPs discovered in the diploid sequence of James Watson (see Materials and Methods). In the remainder of the paper, we show the results derived from the analysis of the Perlegen dataset. The results derived from the analysis of the Watson SNPs are shown in the Supplementary Materials. All of the conclusions in the paper are supported by the analysis of either dataset.

15We measure the level of neutral nucleotide variation in a genomic window using the number of SNPs within the neutral regions divided by the total number of neutral sites (θneu) in a window (see Materials and Methods). This measure is proportional to the conventional Watterson's θ. In the remainder of the paper, all measurements are carried out over 400 kb windows. We have also carried out all of the analyses with two other window sizes, 200 and 600 kb; none of the conclusions change depending on the window size (Table S1, S2 and S3, Figures S7, S8, S9, and S10).

16The level of neutral polymorphism (θneu) depends both on the average time to coalescence within a particular genomic region and on the local constraint and mutation rate. For the purposes of detecting signatures of RSS, variation in constraint and mutation rate generates noise. We assess variability in constraint and mutation rate by measuring divergence per neutral site (dneu) within the neutral regions between the human and chimpanzee genomes (see Materials and Methods). We detect a positive correlation between dneu and θneu (Table 1), confirming that, as expected, constraint and/or mutation rate vary across the human genome. We control for the variation in neutral mutation rate either by carrying out partial correlations with dneu or by using a normalized measure of neutral variation,, where #SNPneu stands for the number of SNPs found in the neutral regions and Dneu stands for the number of divergent sites within neutral regions between humans and chimpanzee genomes. Pneu and θneu also correlate significantly with repeat density (RD) and GC content (GC) (Table 1, Table S1). Finally, in the case of the Watson data, we further carry out controls for the depth of sequence coverage (Table S4).

17The overall effect of RSS on the regional levels of neutral polymorphism should depend on (i) the regional rate of recombination, (ii) the number of recent sweeps (the rate of RSS), and (iii) the strength of positive natural selection associated with a typical adaptive substitution (the strength of RSS). The levels of neutral polymorphism across the genome should correlate positively with the rate of recombination and negatively with the rate and the strength of RSS.

18We take estimates of recombination rate from Myers et al., who used a statistical approach to infer recombination rates from linkage disequilibrium data in humans; these rates have been shown to be highly reliable by comparison to pedigree data. The levels of neutral polymorphism measured by both θneu and Pneu increase with the recombination rate (Figures 1, S1, and S2). The correlation remains when we control for possible confounders such GC content (GC), repeat density (RD), and divergence at neutral sites (dneu) separately (Table 2S) or together (Pearson r (θneu, RR|GC, RD, dneu) = 0.254, Pearson r (Pneu, RR|GC, RD) = 0.209, P<0.001 in both cases).

19Under a model of RSS regions experiencing more frequent or stronger selective sweeps should show lower levels of neutral polymorphism. Because positive selection should be more prevalent in regions of greater functional density, RSS is expected to generate a negative correlation between the degree of functional density and the level of neutral polymorphism. We measure functional density in two complementary ways. First, in each 400 kb window, we count the number of protein-coding codons (FDn) as a proxy of protein-coding density. In addition, we count the number of nongenic sites that can be aligned between primates and zebrafish (FDx) as a proxy of the number of conserved noncoding sites (CNRs) (see Materials and Methods for details).

20Consistent with the predictions under RSS, there are strongly negative correlations between either measure of functional density (FDn,FDx) and measures of neutral variability (Figures 2, S3, S4, and Tables 2, S3). After controlling for GC content (GC), recombination rate (RR), repeat density (RD), and divergence at putatively neutral sites (dneu) (in the case of θneu) the correlations become substantially weaker but do remain statistically significant (Tables 2, 3S). The correlations between FDn and both θneu and Pneu remain significant after we control for FDx; and similarly, the correlations between FDx and both θneu and Pneu are still significant when we control for FDn (Tables 2, 3S).

21The number of differences between humans and chimpanzee genomes at functional regions is likely to be a more direct proxy of the rate of positive selection than the functional density. Consistent with the expectations of RSS, we detect lower levels of θneu (Pneu) in regions of higher Dn (the count of divergent amino acid coding sites) or Dx (the count of divergent sites within conserved noncoding regions) (Tables 3, S3, Figures 3, S5, S6). These correlations remain significant when we control for GC content (GC), recombination rate (RR), repeat density (RD), and functional density (FDn, FDx, or both) (Tables 3, S4). The correlations between either Dn or Dx and either of the two measures of neutral variation (θneu or Pneu) remain statistically significant when we control/correct for the other measure of functional divergence (i.e. control for Dn in the case of correlations of neutral diversity with Dx and, similarly, control for Dx in the case of correlations of neutral diversity with Dn) (Tables 3, S4).

22The genome-wide patterns of nucleotide polymorphism in the human genome contain much information about the historical patterns of mutation, recombination, natural selection and population histories of modern humans. Here we search for traces of recurrent positive selection in the patterns of diversity at (mostly) neutral sites across the human genome. A number of studies argued that positive selection is reasonably common in humans,,,, although substantially less common than in Drosophila,,, –, and E. coli. A recent study estimated that ∼10% of all amino acid substitutions between humans and chimpanzees have been driven by positive selection. If true, then signatures of hundreds of recent selective sweeps should still be detectable in the pattern of neutral variation in the human genome.

23Because recurrent adaptive substitutions leave local (on the order of 0.1 s/ρ) and transient (on the order of Ne generations) dips in neutral polymorphism, persistent adaptation should lead to lower levels of neutral polymorphism in regions of lower recombination and regions where selective sweeps are more frequent and/or stronger on average. Here we have confirmed these predictions by showing that levels of SNP density are lower in the regions of lower recombination and in the regions of higher functional density and functional divergence.

24In addition to RSS, a number of other evolutionary forces can generate heterogeneous patterns of polymorphism: (i) variation in mutation rates and selective constraint, (ii) demographic events such as population structure, bottlenecks, and fast recent population growth, and (iii) hitchhiking of neutral variants with recurrent deleterious mutations (background selection (BS)). In addition, uneven ascertainment of SNPs across the genome could generate spurious variability in SNP density. Below we discuss the evidence in relation to these alternative possibilities and argue that hitchhiking—due to selective sweeps or background selection—needs to be invoked to explain the detected patterns.

25All SNP datasets suffer from ascertainment biases during the SNP discovery phase that can systematically under- or overestimate numbers of SNPs in particular genomic regions or at particular types of sites. We address this concern by using two very different SNP datasets that are likely to have different ascertainment biases: (i) the high quality (type A) SNPs from the Perlegen dataset and (ii) SNPs discovered in the sequenced diploid genome of James Watson. The type A SNPs were discovered using Perlegen oligo hybridization chip technology in a panel of 71 individuals of mixed ancestry. This set is biased against SNPs located in repetitive regions, given that it is difficult to design uniquely hybridizing oligonucleotides in such repetitive regions. The diploid genome of James Watson was sequenced using the 454 technology and does not suffer from the same technological problems as the Perlegen oligonucleotide chip hybridization technology.

26We obtain very similar results using both datasets, which argues that it is unlikely that specific ascertainment biases are responsible for the observed patterns. In addition, we also used the density of the repeats, GC content and functional density as variables in our statistical analyses and showed that all of the signatures of genetic hitchhiking in our data are robust to statistical controls for these variables. The depth of coverage in the Watson sequencing data also does not noticeably affect any of the detected correlations (Table S4).

27The demographic history of human populations in general, and specifically of the populations that have been used for SNP discovery and SNP typing in the Perlegen data, is very complex. Bottlenecks, quick population growth and complex patterns of admixture (for example in the African–American population) are expected to perturb levels of neutral polymorphism across the genome. Collectively, we will denote these forces as “demography”.

28The effects of recent demography undoubtedly generate much variation in neutral polymorphism; however, the correlations that we observe are likely to be weakened and unlikely to be generated by the demographic processes alone. For instance, the lower levels of neutral polymorphism in the regions that have large numbers of the protein-coding (Dn) and functional noncoding (Dx) differences are hard to explain by demography; demographic events cannot easily affect the longer-term rates of functional divergence that have been accumulating for ∼10–14 million years between chimpanzees and humans. On the other hand, it is clear that demography needs to be taken into account in order to use the detected signatures to evaluate the strength of hitchhiking in the human genome.

29Some of the variation in levels of polymorphism in the sequences that we use to measure levels of neutral polymorphism could be due to the variability in the rates of mutation and levels of selective constraint. We measure levels of neutral variation in the sequences that are less likely to be under selective constraint: they are noncoding, located far from exons, and cannot be aligned with distantly related species such as zebrafish. Nevertheless, some residual variation in constraint is likely to remain. Indeed, the positive correlation between our measures of the levels of neutral polymorphism (θneu) and divergence (dneu) (Table 1 and 1S) suggests that mutation rates and/or levels of constraint vary systematically in these regions. It is therefore important to control for the variation in the levels of selective constraint and mutation rate; we do so by using the levels of divergence (dneu) as a variable in partial correlation analyses or by using the measure Pneu.The levels of neutral variation correlate strongly with recombination rate, functional density and functional divergence after controlling for neutral divergence suggesting that these correlations are not due to the variation in mutation rate or constraint

30Partial correlations may not remove all of the effects of the variation in mutation rate and constraint, however. The variation in selective constraint among neutral regions should have a stronger effect on the levels of neutral divergence (dneu) than on the levels of neutral polymorphism (θneu) because deleterious mutations have a greater chance of segregating in the population than to become fixed. This implies that if the negative correlation between θneu and levels of functional density were entirely due to the variation in selective constraint (specifically higher remaining constraint in regions of higher functional density), then controlling for divergence (dneu) should make the partial correlation between neutral polymorphism (θneu) and functional density positive. Yet we see the opposite: the correlations between Pneu and functional density and the partial correlation between θneu and functional density with respect to dneu both remain strongly negative. This suggests that the variation in selective constraint is unlikely to generate the correlations between levels of neutral variation and recombination rate, functional density and functional divergence that we see in this study.

31On the other hand, variability in mutation rates might contribute to some of the observed patterns. Specifically, the positive correlation between neutral diversity and rates of recombination could be due to the mutagenic effects of recombination. Because rates of recombination at local scales (although not necessarily at the 200–600 kb scales relevant to this study) evolve fast,, –, mutagenic effects of recombination should have more pronounced effects on the levels of polymorphism than on the levels of divergence. If so, controlling for neutral divergence (dneu) may not entirely account for the higher mutation rates produced by recent recombination.

32Mutagenic effects of recombination are expected to affect levels of polymorphism proportionately to the rate of recombination in the area, whereas hitchhiking (RSS or BS) is expected to affect levels polymorphism in regions of very low recombination much more substantially. We observe a mostly linear effect of recombination on divergence (dneu) suggestive of the mutagenic effect of recombination and further arguing that the regional recombination rates at the level of our analysis (200 to 600 kb) do not evolve as fast as the location of recombination hotspots. In contrast, the effect of recombination on the levels of polymorphism (θneu and Pneu) is curvilinear, with most of the effect limited to the regions of the lowest recombination rates (Figure 1 and S1). Indeed, when we split the data by the median value of recombination rate (RR = 1.040 cM/MBp), the correlation between the levels of neutral divergence (dneu) and recombination rate (RR) for the two halves of the data are of similar strength (r (dneu, RR|RR<1.040) = 0.197 and r (dneu, RR|RR>1.040) = 0.220). However, the correlations between recombination rate and levels of polymorphism (θneu or Pneu) are much stronger in the low recombination regions than in the high recombination regions ((r (θneu, RR|RR<1.040) = 0.249 versus r (θneu, RR|RR>1.040) = 0.045; r (Pneu, RR|RR<1.040) = 0.194 versus r (Pneu, RR|RR>1.040) = −0.0241). These considerations suggest that most of the positive correlation between recombination rates and levels of neutral polymorphism, and especially the reduction at lower recombination rates, is caused by some form of hitchhiking. These results are consistent with the findings of Hellman et al who detected lower levels of polymorphism in the areas of low recombination close to centromeres and telomeres. Note that in our study we explicitly excluded telomeric and centromeric regions (see Materials and Methods), making our findings complementary to those of Hellman et al.

33Background selection (BS) is the process of hitchhiking of neutral or weakly deleterious polymorphism with linked strongly deleterious polymorphisms –. BS should be more efficacious and lead to lower levels of neutral polymorphism in regions of lower recombination. It is thus quite possible that the positive correlation between neutral polymorphism and recombination rate is due in part to BS. In addition, BS should be stronger in the more constrained genomic regions because such regions should experience higher rates of deleterious mutation (e.g. ). Therefore BS is likely to contribute to the negative correlation between levels of neutral polymorphism and functional density as well. Because regions of higher functional density also exhibit higher rates of functional divergence (Tables 1 and S1), BS could contribute to the negative correlation between levels of neutral polymorphism and functional divergence as well.

34It is less clear whether BS could generate the negative correlation between the levels of neutral polymorphism and functional divergence after controlling for levels of functional density (Tables 3, S4, Figure S6). Two regions of equal functional density can differ in the strength of BS if they differ in the rate of deleterious mutations in the functional sequences. The higher level of deleterious mutations should lead to stronger BS and therefore lower levels of polymorphism in the linked neutral sequences. At the same time, the higher rate of deleterious mutations is likely to come at the expense of neutral mutations at functional sites and thus should lead to lower levels of protein and regulatory divergence. The reduction of neutral mutation rate in the regions of higher deleterious mutation should lead to a positive correlation between levels of neutral polymorphism and functional divergence after controlling for functional density—the opposite of what is seen. On the other hand, the increase in the rate of fixation of weakly deleterious mutations, also expected in the regions of stronger BS, counteracts the reduction of the rate of functional divergence due to the reduction of neutral mutation rate. The combined effect is difficult to estimate given that we do not have information about the distribution of the rates of mutations of different selective effects along the genome.

35There is another pattern we observed that is not naturally predicted by BS. The correlations between functional density (FDn or FDx) and neutral polymorphism weaken very substantially and in some cases become nonsignificant when we control for functional divergence at replacement (Dn) and conserved noncoding sites (Dx) (Table 2). Functional density is likely to be a better of proxy of regional constraint than functional divergence. If BS is indeed the dominant force in the generation of the observed patterns, we might have expected correlations between neutral polymorphisms with FDn and/or FDx to be the most robust.

36Without a better understanding of the distribution of selective effects and rates of new mutations, we cannot reject the possibility that BS contributes substantially to all of the detected patterns. It appears, however, that only specific distributions of selective effects of new mutations would generate all of the observed patterns. Whether such a distribution exists in principle and whether the distribution of selective effects of human mutation satisfies these requirements in fact remains to be determined.

37The arguments above suggest strongly that some form of hitchhiking, either BS or RSS, needs to be invoked to explain the results presented in this paper. These results also suggest that natural selection at both coding and regulatory sites affect linked neutral polymorphism. This is because the measures of the rate of functional evolution at coding and regulatory sites appear to influence levels of neutral polymorphism independently of each other. Specifically, divergence at coding sites and divergence at regulatory sites correlate negatively with the levels of neutral polymorphism after controlling for each other and for the variation in levels of functional divergence (Table 3, S3). To the extent that this is due to recurrent adaptation selection at both coding and regulatory sites, this would echo results of McDonald-Kreitman analyses of adaptation in Drosophila.

38Levels of neutral polymorphism correlate stronger with divergence at coding than at non-coding regions, possibly implying that either a higher proportion of nonsynonymous changes are adaptive compared to changes in regulatory regions or that the nonsynonymous adaptations have higher selective coefficients. It is also possible and even likely that Dx is a noisier measure than Dn due to greater difficulties in identification of regulatory regions and the noise in estimating Dx due to misalignments. This pattern may also be due to different rates or distributions of the selective effects of deleterious mutations located in coding and regulatory regions, leading to varying effects of BS on linked neutral polymorphism and functional divergence.

39These results can also be used to assess the importance of hitchhiking (either RSS or BS) in affecting patterns of neutral polymorphism. The levels of neutral polymorphism appear to be ∼50% lower in the regions of high Dn or Dx (Figures 3, S5) relative to the regions of zero functional divergence (Dn or Dx = 0). If we assume that this effect is entirely due to hitchhiking, then by using the observed correlation between θneu and Dn, we estimate that the levels of polymorphism genome-wide are reduced by 6% genome-wide (Materials and Methods). This reduction is much more pronounced in the more gene-rich regions. For instance, in the 50% of the most gene-rich regions (regions that have greater than the median density of codons (FDn)), the neutral polymorphism is reduced by 11%, while in the regions that contain 50% of the genes (regions that have greater than the mean density of codons (FDn)), the neutral polymorphism is reduced by 13%.

40It is clear that hitchhiking has left a significant imprint on the patterns and levels of neutral variability in the human genome and that the effects of natural selection at linked sites cannot be ignored in the analysis of polymorphism data in humans. The challenge for the future is to use these signatures to answer a number of outstanding questions. What are the selective effects and genomic distributions of adaptive and deleterious changes responsible for RSS and BS? What is the biological nature of these changes? What is the relative importance of RSS and BS? Can we estimate parameters of adaptive evolution in the presence of BS? The availability of whole genome sequences in a large number of humans may provide the necessary data to answer these questions. What is needed now are the models and tools to harness these data to provide a cogent picture of the effects of natural selection on human genome and human evolution.

41All analyses have been carried out using two SNP data sets—Perlegen data and Watson data. Perlegen data were downloaded from http://genome.perlegen.com. These data were annotated based on the NCBI build 35 of the human genome sequence. We updated all the genomic positions of the SNPs to match the latest NCBI build 36, according to the rs number of SNPs in the dbSNP build 127. During the processing, 1,361 SNPs were discarded because they could not be uniquely mapped to the human genome. Perlegen data contain three classes of SNPs: (A) array-based genomic resequencing, (B) reliable external SNP collections, and (C) unvalidated, lower confidence sources (see Supplementary text of ). We excluded class B and C SNPs and retained 1,235,057 class A SNPs located on autosomes for our analysis. The Watson data were downloaded from http://jimwatsonsequence.cshl.edu/. The genome of James Watson was sequenced at 6× coverage using 454 Life Sciences Technology and matched to the human genome project's published reference sequence. In the Watson DNA sequence, heterozygous sites, in which each site was sequenced multiple times and both forms of the base were found in the diploid genome, were ascertained as SNPs. Homozygous sites of Watson's DNA sequence that have been sequenced multiple times and that differ from the reference sequence of the human genome were also ascertained as SNPs. In total our Watson dataset consisted of 2,020,767 SNPs.

42Whole-genome alignments of human (H), chimpanzee (C), and zebrafish (Z) sequences were obtained from the Ensembl compara database through the Ensembl Application Program Interfaces (APIs). We defined the “neutral” genomic regions of the human genome if the regions were: (1) H-C aligned, (2) not H-C-Z aligned, (3) located at least 5 kb away from the starting and ending position of transcripts of protein-coding genes and at least 1 kb away from any exons, (4) located on autosomes at least 10 Mbp away from the boundaries of centromeres and the ends of telomeres, (5) not located in the simple repetitive regions of the human genome. The chromosomal coordinates of exons, transcripts and simple repeats were obtained from the finished and annotated human chromosome sequence from the Ensembl database (build 36).

43Neutral divergence was assessed from H-C alignments. The accuracy of estimation of neutral divergence may be influenced by the misaligned sequences. Indeed, we discovered some short (2 kb on average) neutral genomic regions having extremely high levels of divergent sites, which may result from misalignments (data not shown). To minimize the possible influence of misalignments, we only counted “isolated” substitutions that are flanked by two monomorphic positions on each side (i.e. no substitutions or SNPs were mapped to these sites). We denoted the number of isolated substitutions between human and chimpanzee sequences as Dneu, and the number of isolated substitutions per neutral site, dneu. To measure neutral polymorphism, we counted the number of SNPs in neutral regions and denoted the number of SNPs per site as θneu. Alternatively, we measured neutral polymorphism with. Data manipulation was done using Matlab functions based on PGEToolbox and MBEToolbox.

44We used four metrics as proxies of the rate of adaptive evolution for a given region in the human genome. Functional density was measured using FDn, the number of codons, and FDx, the number of aligned bases in the H-C-Z three-way alignments. Functional divergence was measured using Dn, the number of codons involved in nonsynonymous substitutions between H-C orthologous gene pairs, and Dx, the number of H-C substitutions in H-C-Z alignments that are located in noncoding human genomic regions. For each pair of genes, the amino-acid sequences were extracted and aligned using CLUSTALW with the default parameters. The corresponding nucleotide sequence alignments were derived by substituting the respective coding sequences from the protein sequences. The synonymous substitution rate (Ks) was then estimated by the maximum-likelihood method implemented in the CODEML program of PAML. Insertions and deletions within alignments were discarded. Poorly aligned orthologous pairs, as indicated by Ks>5, were excluded. The codons containing nonsynonymous substitutions were mapped back onto the human genome and positions were recorded. For simplicity we counted the numbers of codons causing amino-acid changes instead of the numbers of single nucleotide replacement substitutions. In calculation of Dx, we excluded “tri-allelic” sites where the bases of H, C and Z all differ from each other.

45We used 400 kb (as well as 200 and 600 kb) sliding window with a step of 100 kb to scan along the human genome. For each window, two measures of neutral polymorphism (θneu and Pneu) and four proxies of the rate adaptive evolution (FDn, FDx, Dn, and Dx) were estimated. To reduce noise arising from small sample size, we also discarded the windows with Dneu<500 and the ones with the total amount of “neutral” sequence less than 2 kb. 22,553 400 kb windows have been used for the correlation analysis. Spearman rank correlation or Kendall's correlation coefficients have been calculated in all cases. To visualize correlations between variables, we used scatter plots with regression lines superimposed. We also pooled the data points of neutral polymorphism by the values of the proxy of adaptation under consideration (e.g. Dn). To do this, we ranked all the data points of the neutral polymorphism by the values of the proxy and then pooled them into 100 bins such that each bin had equal size (i.e., 1%) of the data points. We then computed average values of the proxies of adaptation and the average value of neutral polymorphism for each bin, and superimposed them onto the scatter plots.

46To control for confounding variables, we calculated Spearman partial correlation coefficients between variables X and Y controlling for Z, using the function partialcorr in the Matlab statistic toolbox. Recombination rate estimated by using the coalescent method of were downloaded from http://hapmap.org/downloads/recombination/. The density of simple repeats was computed as the proportion of bases of simple repeats in the given region. Chromosomal coordinates of simple repeats in the human genome, identified by RepeatMasker, were obtained from the UCSC genome browser.

47We also calculated the partial correlation coefficients between variables X and Y by calculating the correlation between the two sets of residuals formed by two linear models X∼Z and Y∼Z (see also ) where Z stands for either one or a series of variables. The distribution of Dn, Dx, FDn, and FDx values is approximately exponential, which is a problem in a least squares linear model framework in controlling for a third variable, Z. The linear model used to regress out Z is sensitive to the highly non-normal distribution of variables, and the residuals will be highly non-normal, making the results difficult to interpret. Therefore, we quantile-normalized values, replacing the original estimates with their theoretical quantiles based on a normal distribution. Then, we fitted linear models, using as the response variable quantile-normalized Dn, Dx, FDn, or FDx, and using as the predictor variables various combinations of recombination rate, GC content, and the density of simple repeats.

48Estimation of the effect of hitchhiking on the level of neutral polymorphism was calculated using the regression between θneu on Dn, using the formulawhere q is the reduction of polymorphism due to hitchhiking, i is a window count for the subsets of windows used in the analysis (e.g. FDn>median (FDn)), b is the intercept of the regression of θneu on Dn.