Health Articles
remainder were retained for use in estimation.
Figure 2.8: MCMC output from one dataset. Nine independent chains are run from
overdispersed starting points and are shown here in individual colors. Displayed are the joint posterior
probability ("postProb"), µ ("X2Ns mean"), and σ ("X2Ns sd"), over the 50,000 iterations for the
case of UC regions using EA+AA pooled samples. Shown in the right hand panel are the marginal
distributions of each of the parameters. Convergence of our chains can be seen here by the relatively
consistent bands traced by the model parameters among chains.
Maximum Likelihood Estimates [AK]
The Hierarchical Bayesian MCMC model described in Section 2.12 treats the
mean selection coefficient µ as a random variable and allows calculation of its expected
last iteration in chain
last iteration in chain
last iteration in chain
Figure 2.9: Convergence plots for MCMC simulations. Shown here is the reduction in
the [12] scale reduction factor ("shrink factor") ˆ
R over the course of our MCMC iterations for the case
of UC regions using EA+AA pooled samples. Values close to 1.0 indicate convergence of the MCMC to
the posterior distribution of our parameters. Plot labels as in Figure 2.8.
value and variance. A less computationally intense and perhaps more conventional
approach can be used to derive the Maximum Likelihood (ML) point estimate of the
selection coefficient α where the likelihood of the data is, as before, the product of the
individual SNP probabilities given in Equation 2.3 but the SNPs are assumed to be
independent and identically distributed with a common selection coefficient α. That is,
we find the value of α that maximizes the likelihood
Following the methods of [13] we have numerically maximized Equation 2.7
to determine these ML estimates for our various genomic regions and population sub-
samples. The results comparing the ML estimates to the MAP estimates from the full
Bayesian model are shown in Table 2.5. It can be seen that the agreement between the
two methods is quite good so that our conclusion about the strength of selection on the
ultraconserved elements is not dependent on our Bayesian methods.
Table 2.5: Estimates of the selection coefficient using Bayesian and MaximumLikelihood methods MAP: Maximum A Posteriori estimate for mean value of selection coefficient
based on hierarchical Bayesian MCMC. ML: Maximum Likelihood point estimate for the selection
coefficient. Other symbols as in Table 2.3.
Effect of Linkage Disequilibrium
To investigate the possibility that linkage disequilibrium (LD) had an effect on
the derived allele frequency spectra, and hence my estimates of the selection coefficients
for the categories of sites examined, I performed calculations on subsets of the segre-
gating sites found in the pooled EA and AA samples analyzed above. For the Seattle
SNPs data, I chose 50 random subsets of the nonsynonymous and intron segregating
sites such that each subset had at most one segregating site per gene. For the UC and
UF sites I chose 50 random subsets such that each subset had at most one segregating
site per UX element. Thus the segregating sites in any of these subsets should have no
LD. For the UC sites I also chose 50 random subsets containing the same number (105)
of segregating sites as in the one-per-UX subsets of UC sites, but this time without
regard to the UX element of a site.
For each subset I calculated the ML estimate of the selection coefficient as
described in Section 2.14. The distribution across the subsets for each class is shown in
Figure 2.10. The median values of the subsets (listed in the caption of Figure 2.10) are
very close to the ML values for the EA+AA pooled samples (shown in Figure 2.10 and
listed in Table 2.5) with the largest effect being about a 15% reduction in the estimated
strong negative selection coefficient for the UC category. But this category has the
largest variance and the interquartile range of its subsets is completely consistent with
the its distribution in Figure 2.1. I am thus confident that my essential conclusion about
the selection coefficients in the different categories is not an artifact of LD.
ML Selection Coefficient
Figure 2.10: Selection Coefficient Estimates without Linkage Disequilibrium.
Boxplots of the distributions of Maximum Likelihood (ML) estimates of the selection coefficient, in units
of 2Nes, derived from 50 subsamples (for each category) of the segregating sites in the pooled EA+AA
samples are shown. The NonSyn and Intron subsamples include only one SNP per gene in the Seattle
SNPs data set. The UF and UC-noLD subsamples include at most one SNP per UX. The UC-rand
subsamples had the same number of SNPs as UC-noLD, but were chosen without replacement from UC
SNPs without regard to UX. The ML estimates for the full set of EA+AA samples in each category are
shown as thin red lines. The median values (-5.2, -4.4, -1.5, -1.3, -0.16) for the subsets are consistent
with the distribution of values in Figure 2.3, indicating that the latter are not strongly dependent on
linkage disequilibrium in the respective regions of segregating sites.
Ascertainment Effects [AK]
Ultraconserved elements of the human genome were identified through the
criterion of 100% sequence homology among species.
Since polymorphisms within
species are apt to be confused with fixed differences between species when single se-
quence comparisons are used, one would expect ultraconserved elements to contain
fewer polymorphisms than non-ascertained genomic regions even under the neutral
model. Furthermore, one would expect those polymorphisms that are present in ul-
traconserved elements of the genome to be segregating at lower frequencies than those
in non-ultraconserved regions, since intermediate and higher frequency polymorphisms
are more likely to be mistaken for fixed differences between species. Since our estimation
of selection coefficient uses information from the frequency spectrum of segregating poly-
morphism, these biases might mimic the effects of selection on ultraconserved regions
of the genome.
To quantify the strength of this ascertainment bias in our data, we performed
extensive coalescent simulation of the ultraconserved element identification process, to
examine the extent to which selection of ultraconserved elements alone can skew the
site frequency spectrum under neutral models. We used Hudson's "ms" program [32] to
perform the simulations under the standard neutral model (SNM). Our procedure con-
sisted of simulating 100,000 unlinked genomic regions from 96 individuals (representing
our human samples) and 1 outgroup lineage (representing a chimp sequence), with a
species divergence time of 5 units of 4Ne generations (about 5 million years assuming
Ne = 10,000 and a 25 year generation time), and a fixed value of θ (4Neu = 0.5).
We then split these 100,000 unlinked regions into those that were "ultraconserved" vs.
"non-ultraconserved" by comparing our outgroup sequence to one randomly selected
individual from our ingroup and finding those regions that had zero fixed differences in
the two-sequence comparison.
Distributions of Tajima's D [75] for our two classes of region are shown in the
left panel of Figure 2.11. The bias imposed by selection of the ultraconserved genomic
regions on Tajima's D results in a shift in the frequency spectrum towards rare alleles for
ultraconserved elements but the bias is quite small. Under the standard neutral model
we do not see values of Tajima's D that are as extreme as we observed in our data
(Table 2.3). Other demographic models incorporating population expansions yielded
similar results (data not shown).
Additionally we estimated the strength of selection α using ML estimates as
in Section 2.14, for the separately pooled sets of simulated "ultraconserved" and "non-
ultraconserved" polymorphisms. These simulations were repeated 10,000 times for the
standard neutral model (SNM) yielding the distributions shown in the right panel of
(Figure 2.11).
A fairly large difference between simulated ultraconserved and non-
ultraconserved sequences is observed. Ascertainment of ultras results in a negative shift
in estimates of the strength of selection, even under neutrality. In order to correct for
this ascertainment bias, we applied the procedure described in Section 2.17.
ML Selection Coefficient
Figure 2.11: Effect of Ascertainment Bias of Ultraconserved Elements onTajima's D and Estimates of Selection Coefficient.
Coalescent simulations were
performed using the standard neutral model (SNM) as described in the text, which were then split into
"ultraconserved" (ultras) and "non-ultraconserved" (non-ultra) subsets based on the absence/presence
of fixed differences in a two-sequence comparison between the outgroup sequence and a randomly chosen
ingroup sequence. Left panel shows boxplots of Tajima's D statistic, indicating that scertainment of the
ultraconserved elements has only a weak effect of biasing Tajima's D downward. Right panel shows box-
plots of Maximum Likelihood estimates of the selection coefficient, in units of 2Nes. If uncorrected, the
ascertainment bias for ultras results in an apparent negative selection coefficient for neutrally evolving
Correcting for Ascertainment Bias [AK]
As noted in Section 2.16, identification of the ultraconserved elements on the
basis of no divergence between species will impact the expected site frequency spec-
trum at those loci, even under the neutral model. To obtain accurate estimates of the
strength of selection at UC regions we modified our likelihood function (Equation 2.3
and Equation 2.4) to directly correct for this ascertainment bias. Following the general
framework of Nielsen et al. [51], we assume that our ascertainment was performed in
a subsample, size d, of our final sample, size n, where the derived allele frequency at
segregating site k was ik, such that the probability of ascertainment (probability of
choosing d ancestral alleles) is
In our case d = 1, corresponding to sampling the ancestral allele at our single genomic
reference sequence. With this probability in hand, using Bayes rule we can write down
the ascertainment corrected version of Equation 2.3 (suppressing the dependence on
sample size, but emphasizing that the bias is dependent on the site k). The probabil-
ity of observing a SNP segregating with derived allele frequency ik, conditional upon
ascertainment and the strength of selection α is
P r(ik α)P r(Asck ik, α)
P r(ik Asck, α) =
where the denominator is simply the normalization factor
P r(j α)P r(Asck j, α).
Our corrected likelihood function can then be expressed as
P (Data Asc, α) =
This ascertainment corrected likelihood function was tested using the coales-
cent simulation scheme described in Section 2.16. As can be seen in Figure 2.4, this
procedure returns the expected value of our estimates of the strength of selection back to
zero. Subsequently, all estimations done for UC regions, both Bayesian and Maximum
Likelihood, use this ascertainment corrected formulation.
A comparison of the ML estimates for the ultraconserved regions with and
without correction in various samples is given in Table 2.6. The correction reduces the
estimate for the UC regions by about 10-20%.
Table 2.6: Ultraconserved Element Ascertainment Bias Correction in Esti-mates of the selection coefficient. Maximum Likelihood point estimates for the selection
coefficient in the ultraconserved regions with (UC-corr) and without (UC-uncor) correction for the as-
certainment bias. The correction reduces the estimate by about 10-20%. Other symbols as in Table 2.3.
Evolution of Human Accelerated
Regions is GC-Biased
As noted in the Introduction, I enriched the genomic DNA from 11 individu-
als among the Yoruba from Ibidan, Nigeria (YRI samples) for 40kb neighborhoods of
the top 49 Human Accelerated Regions [59] (HARs) (hereinafter "harseq1-49"), and 13
similar control regions not containing HAR elements (hereinafter "ctlreg50-62") using
a microarray hybridization technique. Using ABI SOLiD [6] high throughput sequenc-
ing technology I obtained sufficient coverage in the target regions of the resulting short
(35bp and 50bp) sequencing reads to restrict my analysis to genotype calls at positions
where I had at least 35-fold coverage for an individual, to ensure accurate genotyping.
I determined the frequency of the derived (non-ancestral) allele in my set of 22 chro-
mosomal samples for all the resulting segregating sites (i.e. sites where not all samples
have the same allele).
I analyzed these segregating site frequencies in several ways. First, I compared
the frequency spectra of two subsets of my segregating sites: those comprising W2S
mutations (where the ancestral allele is either A or T, and the derived allele is either
G or C) versus the S2W mutations. Next, I compared the ratios of these different
classes of mutations to the same ratios for fixed substitutions that have arisen on either
the human or chimp lineages since our common ancestor ∼5 million years ago. For
both of these analyses I tested the strength of my results across the full scale of my
regions, and also performed the tests on the data pooled across my 49 genomic regions.
Finally, I examined the pattern of spatial variation of derived allele frequencies in search
of evidence for a selective sweep. Results for harseq regions and ctlreg regions were
compared to each other and to a set of 62 genic regions resequenced in the same 11 YRI
individuals by the Seattle SNPs project [4].
My key conclusions reinforce the notion that a bias in fixation of GC basepairs
is at play in the HAR neighborhoods. The use of the present tense ("is at play") in the
previous sentence is intentional, for I have two significant results. On the one hand, I
confirmed the historical GC-bias for fixation, finding in 10 of the 49 HAR neighborhoods
that fixed differences are more likely than polymorphic sites to show this bias. But in
my more novel result I also found that there is a similar GC-bias towards higher derived
allele frequencies at polymorphic sites alone. The high derived allele frequencies at these
GC-biased SNPs make them likely candidates for eventual fixation. I found this ongoing
bias in 11 cases, only 1 of which also shows historical GC-bias. Thus, almost half of HAR
regions (20) are currently undergoing, or have historically undergone, biased fixation of
GC alleles. This pattern strongly suggests that selectively neutral forces, such as biased
gene conversion, have helped to shape many of the genomic regions containing HARs
and could have influenced the fixation of a large number of GC alleles in the HAR
I also used my resequencing data to search for signals of a recent selective
sweep, which could have indicated that one or more of the HAR substitutions had been
subject to recent positive selection in the human population. In that analysis, I found
that only weak signals of that kind remain, possibly due to millions of years of decay.
It may be noteworthy that 3 of the 5 regions with those signals are also significant in
one of the two GC-bias tests.
The following sections provide more details of these results and describe the
methods used to derive them. In this work I collaborated with Andrew Kern who
implemented the coalescent simulations to determine p-values for the SweepFinder test
described below. Another collaborator in this work was Katherine Pollard who provided
the substitution rate matrix used in simulations to validate the p-values of the MK test
described below.
Results of HAR Neighborhood Analysis
Weak-to-Strong Mutations are Being Swept to Higher Derived
Allele Frequencies
When the HARs were first described, a strong W2S substitution bias was noted
in the human-specific substitutions in these elements [59, 60]. This bias was extremely
pronounced in HARs 1,2,3, and 5, but also noticeable as a general trend in the entire
set 1-49. This evidence suggested that biased gene conversion (BGC) could have had
an historical role in the evolution of the HARs. Here I analyze my list of segregating
sites in the 40kb HAR neighborhoods to determine if such bias is still ongoing in the
human population. In each region, I separately computed the derived allele frequency
spectra for the W2S mutations and the S2W mutations. I then tested for an offset in
the spectra between the two categories with a two-sided Mann Whitney U (MWU) test.
This test has been shown to have good power to detect fixation bias [3].
I found a significant (p ≤ 0.05) difference in 11 out of 49 harseq regions (Table
3.1). In all 11 significant harseq regions, the offset was for the W2S mutations to be
segregating at higher derived allele frequencies than S2W. This implies that regardless
of the rate of introduction of W2S or S2W mutations, it is the W2S mutations that are
more likely to reach high frequency and eventually fix in the human population. This
is certainly consistent with a mechanism of gene conversion that favors selection of G
or C alleles from a heterozygote or some other selective force generally favoring higher
GC content. The novel features of this result are that it indicates that this process is
ongoing and not confined to the core 100–400bp HAR elements.
This ongoing W2S fixation bias distinguishes the harseq regions from ctlreg
regions and Seattle SNPs regions. Significant MWU tests were observed at none of
ctrlreg50-62 (Table A.2) and five out of 62 Seattle SNPs regions (Table A.3), one of which
has higher derived allele frequencies in S2W compared to W2S mutations. Applying
the MWU test to simulations of a neutral coalescent model as described in Section 3.2.7
showed that the p-values from this test accurately reflect the fraction expected by chance
from a neutrally evolving locus.
The two harseq regions with the greatest offset were harseq21 and harseq34
(Figure 3.1). Note that across these two 40kb regions, the ratio of W2S to S2W segregat-
ing sites is not extreme; it is the W2S shift towards higher frequency in the population
that is significant. These ratios are consistent with the smoothed ratios reported in
the top several thousand conserved candidate HAR elements [59] even though my 40kb
regions do not comprise largely conserved regions.
It is natural to theorize that BGC will be a stronger driving force in areas
of high recombination and studies have shown there to be a good correlation with the
male recombination rate in particular [17]. I examined the recombination rates in the
enclosing 1Mb windows as determined by the deCODE project [41]. Harseq21 is an
outlier in that it is contained in a genomic region of extremely high recombination rate
(male 4.29 cM/Mb, sex-averaged 3.43 cm/Mb, in contrast to genome-wide averages of
0.93 cM/Mb and 1.29 cm/Mb respectively). But this is not true of harseq34 (0 male and
1.42 cm/Mb sex-averaged). For the remainder of the regions with a significant p-value
harseq21 segregating sites spectra. AT to GC vs. GC to AT mutations
fraction of sites
derived allele frequency (count of 22)
harseq34 segregating sites spectra. AT to GC vs. GC to AT mutations
fraction of sites
derived allele frequency (count of 22)
Figure 3.1: Frequency offset of weak-to-strong vs. strong-to-weak mutations.
MWU test is strongly significant for harseq21 (upper panel) and harseq34 (lower panel). This reflects
the fact that the derived allele frequency spectrum for weak-to-strong mutations (dark bars) is offset
towards higher frequencies compared to strong-to-weak mutations (light bars). N: count of segregating
sites of the indicated category in the region.
on the MWU test, the rates vary. An additional (not unrelated) factor that has been
strongly correlated with biased substitutions is chromosomal position near telomeres
[17]. With the exception of harseq1 this is not the case for the regions with significantly
shifted W2S spectra (Table 3.1).
I also performed the MWU test for the shift in W2S sites toward higher fre-
quencies after pooling all the segregating sites in my 49 harseq regions, and found a
p-value < 10−10. By contrast, the test was not significant (p = 0.26) for the pooled
data in my 13 control regions, thereby controlling for possible systematic bias in my
sequencing and genotyping techniques. I thus conclude that in many of the neighbor-
hoods of the top 49 HARs there is an ongoing force driving W2S mutations to higher
frequency in the human population.
Weak-to-Strong Mutations are More Likely to Have Become
Fixed Differences
Since the HARs were essentially defined based on fixed differences between
the human and chimp reference genomes in otherwise strongly conserved elements [59],
I compared such human/chimp fixed differences to the segregating sites in my human
samples. My set of fixed differences was based on high quality base calls from the
reciprocal best alignments [38] of human and chimp genomes. I further restricted this
list to the locations within my regions for which I had the above-mentioned 35-fold
coverage for an individual in my sample. Finally, I removed from this initial fixed
difference list those positions that I found to be segregating in my samples, or that
appeared in the dbSNP129 database [71] (see Section 3.2.6). The latter two filters
removed 6.7% of the fixed differences at high coverage positions. Since I do not have
information on sites that are segregating in the chimp population, I could not remove
those, but would expect the number to be similarly small.
I separated the mutations in my segregating site set into the categories: W2S,
S2W, and neither. I similarly divided the set of fixed differences, regardless of whether
the substitution occurred on the chimp or human lineage. As in reference [17], I per-
formed a variant of the McDonald-Kreitman (MK) test to compare the W2S:S2W ratios
in the sets of mutations and the sets of substitutions as described in Section 3.2.6.
I found a significant (p ≤ 0.05) difference between the substitution patterns
of segregating versus fixed sites in 11 of the 49 harseq regions (Table 3.2). In all but
one (harseq39) of these 11, the fixed substitutions had relatively more W2S mutations
(compared to S2W) from the ancestral form than did the segregating sites. Four of the
11 fell in the 40kb neighborhoods of the top 11 HARs, and indeed the strongest result
(p = 0.00015) was for harseq1. This is not surprising, since the HAR1 element has 18
fixed differences, all W2S [60]. By contrast, none of the ctlreg regions and nine out of
62 Seattle SNPs regions had a significant p-value on this test (Tables A.2 and A.3).
All of the nine significant Seattle SNPs regions had a higher W2S:S2W ratio in fixed
differences than in segregating sites. Thus, the harseq regions have a much stronger
signal for historical fixation bias than my control regions and a somewhat stronger
signal than the genic Seattle SNPs regions. Applying the MK test to simulations of a
neutrally evolving primate phylogeny (see Section 3.2.7) showed that the p-values from
this test accurately reflect the fraction expected by chance from a neutrally evolving
I also recapitulate and expand the finding [59] that this bias towards W2S fixa-
tion is associated with telomeres, since 7 of the 11 significant regions under the MK test
were found either in the karyotype band containing a telomere or the one immediately
adjacent. (This count includes the two cases of harseq19 and harseq36 from chromosome
2 that fall adjacent to the ancestral telomeric fusion event at 2q14.1.) Although many
of the 11 regions (including ones near telomeres) have elevated recombination rates, it
is also worth noting that 5 of the 11 are in regions with much lower than average male
recombination rates (including the two near the chromosome 2 fusion site, and harseq11
One noteworthy negative result from the MK test, (and the MWU test as
well) is harseq2. It had been noted [59] that the core HAR2 element showed a strong
bias towards W2S fixations, and that this extended to a region of ∼1kb. Here I find
no significant signal for either the MK or MWU test in my 40kb neighborhood of
that element (Table A.2). Another noteworthy negative result on these tests is the
harseq6 region, which has an extremely elevated rate of mutation as estimated either
by nucleotide diversity [76] or by the number of segregating sites [80] (Table A.1), but
apparently no strong bias towards weak-to-strong fixation (Table A.2).
Although the MK test, like the MWU test, is consistent with the mechanism
of BGC favoring fixation of G or C alleles, in fact only one of my regions (harseq1) had
significant results in both tests. Since the original identification of these elements as
HARs depended on substitutions fixed on the human lineage, it is not surprising under
the assumption of BGC to find that a bias towards W2S fixation is associated with
the HAR neighborhoods. The complementary evidence from the MWU test (a total of
20 of my 49 test regions are W2S-significant on one test or the other) indicates that
the ongoing bias in favor of W2S mutations has also probably led to human specific
substitutions in otherwise conserved elements.
The Scale of Fixation Bias towards Weak-to-Strong Mutations
Because BGC is posited to operate on a scale much smaller than the 40kb of my
target sequencing regions, perhaps operating at localized recombination hotspots, and
because the (100–400bp) HAR elements at the core of my target regions were suspected
of arising in part due to BGC [59], I wanted to test whether the MWU and MK signals
depended on these core elements. I therefore performed the same tests after masking
out the central 500bp, 1kb, 5kb or 10kb of each region.
I found that signals of both ongoing and historical fixation bias are fairly robust
to removing sequences including and flanking the core HAR element. For the MWU
test, all but two (harseq1, harseq18) of the 11 regions that were significant at the 5%
level for the MWU test were still significant at that level with the central 5kb omitted.
Seven of the 11 remained significant even with 10kb omitted (Table 3.1). For the MK
test, the results were slightly more sensitive to masking. Considering the 11 regions with
a significant result on that test (including the one favoring S2W fixation), 9(7) of these
were still significant with the central 1kb(5kb) masked out (Table 3.2). I conclude that
the evolutionary forces behind these results is not confined to the small HAR elements
themselves, but rather that any bias in the substitutions found in the HARs is likely a
byproduct of the forces acting at a larger scale.
To test whether there might be other localized elements within the 40kb regions
driving these results I performed the tests under a set of fifteen overlapping 5kb masks
(centered at a 2.5kb spacing along each region). Among the 11 MWU significant regions,
5 were still significant at 10% under this regime, while 3 completely lost significance (p
> 20%) for at least one such mask (Table 3.1). Of the 11 MK significant regions, 5
were still significant at 10% (including harseq1) under this regime, while 1 completely
lost significance (p > 20%) under at least one mask (Table 3.2). It is worth noting that
the harseq1 result reflects the fact that of the 105 segregating sites I found in that 40kb
neighborhood, 71 were S2W and only 19 W2S (Table A.2).
I conclude from this set of tests that the evolutionary forces behind W2S
fixation bias are not necessarily highly local. If fixation bias relies on recombination
hotspots and BGC, I have to posit that such hotspots extend over a long range of bases,
or are somehow temporally and spatially variable.
No Significant Evidence for Selective Sweeps
For each of the studied regions, I used the SFS to calculate two population
genetic statistics that can sometimes indicate positive selection: Tajima's D [75], which
is based on the folded SFS, and Fay and Wu's H [23]. Neither of these statistics exceeded
the value of ±2 in any region (Table A.1). I next compared the distributions of these
two statistics in the 49 harseq regions and the 13 ctlreg regions, to those for the same
population (YRI) in 104 genic regions resequenced by Seattle SNPs (see Section 3.2.7). I
found the ctlreg regions to be indistinguishable from the Seattle SNPs for these statistics,
while the harseq regions were only mildly more negative for Tajima's D (Wilcoxon rank
sum p-value = 0.08) and not significantly different for H (Figure 3.2). These results
strongly suggest that the site frequency spectrum in harseq regions is indistinguishable
from that found in either my control regions (ctlreg), or in the Seattle SNPs data set.
Thus there is no reason to believe that harseqs represent some kind of genomic outlier
with respect to recent selective events.
To test for evidence of a selective sweep, I analyzed the spatial variation of de-
rived allele frequencies at the segregating sites from my 22 chromosomal samples in each
of the target 40kb regions using the SweepFinder program. This software determines a
composite likelihood ratio (CLR) statistic comparing the hypothesis of a complete selec-
tive sweep at the location to the null hypothesis of no sweep using Test 2 from Nielsen
et al. [52]. I tested along a grid of 1000 points in each target region. This test has been
shown to be robust to demographic deviations from the standard neutral model in its
ability to use an arbitrary background site frequency spectrum [82]. I tested with two
such backgrounds: the first from the pooled set of all the data in the 49 harseq regions
plus 13 ctlreg regions, the second from the same population (YRI) as my samples but
with frequencies taken from the Seattle SNPs resequencing data [4] for a large set (104)
of genic regions. It should be noted that using the SFS from my data as the background
harseq1to49 vs p2seasnpgenes Wilcoxon 2−sided p−value = 0.08156
ctlreg50to62 vs p2seasnpgenes Wilcoxon 2−sided p−value = 0.9931
harseq1to49 vs p2seasnpgenes Wilcoxon 2−sided p−value = 0.1303
ctlreg50to62 vs p2seasnpgenes Wilcoxon 2−sided p−value = 0.6303
Figure 3.2: Comparison of Population Genetic Statistics in harseq regions andSeattle SNPs genes. Boxplots show the distribution of Tajima's D, Fay and Wu's H in my test
regions (harseq1-49), my control regions (ctlreg50-62) and YRI samples sequenced in 104 Seattle SNPs
to define the neutral model should be particularly conservative in that I am testing any
given region for deviations from that neutral model. To determine the significance of the
maximum CLR values, we performed coalescent simulations of each target region and
ran the SweepFinder program on each simulated set of segregating sites as described in
Section 3.2.5. I report as a p-value the fraction of simulations of each target that had
a CLR greater than or equal to the actual maximum CLR for that target (Table A.2)
The harseq regions with the most significant five SweepFinder p-values are
listed in Table 3.3. These are nearly all at the 95% confidence level for either of the two
background distributions used, but I note that none are individually significant after a
conservative Bonferroni correction, given the 49 harseq regions that were tested. Since
these may nevertheless harbor mutations that were selected for in the human lineage,
here I briefly note some of their characteristics that can be seen in tracks from the UC
Santa Cruz Genome Browser (Figure 3.3).
Table 3.3: SweepFinder Hits Target regions with most significant p-values for SweepFinder.
Har,Sea refer to null model background frequency spectrum derived from all target regions or Seat-
tleSNPs respectively. RecA and RecM are the sex-averaged or male only recombination rates at the
target region. karyo is the chromosomal karyotype band where the target region is found.
Unlike the other four SweepFinder hits, which all contain introns or exons of
coding genes, harseq25 is in a gene "desert". The nearest known gene, approximately
1Mb away on chromosome 4, is ODZ3, which is a transmembrane signaling protein most
23 q24 q25 4q26 27
A harseq25
SweepFinder CLR on harseq resequencing. bgrnd sfs = seasnp. gridsize = 1000
Segregating sites from harseq resequencing. Reddish for higher derived frequencies
Placental Mammal Basewise Conservation by PhyloP
Human Accelerated Regions
B harseq9
SweepFinder CLR on harseq resequencing. bgrnd sfs = seasnp. gridsize = 1000
Segregating sites from harseq resequencing. Reddish for higher derived frequencies
UCSC Genes Based on RefSeq, UniProt, GenBank, CCDS and Comparative Genomics
Placental Mammal Basewise Conservation by PhyloP
Human Accelerated Regions
C harseq11
SweepFinder CLR on harseq resequencing. bgrnd sfs = seasnp. gridsize = 1000
Segregating sites from harseq resequencing. Reddish for higher derived frequencies
UCSC Genes Based on RefSeq, UniProt, GenBank, CCDS and Comparative Genomics
Placental Mammal Basewise Conservation by PhyloP
Human Accelerated Regions
D harseq16
SweepFinder CLR on harseq resequencing. bgrnd sfs = seasnp. gridsize = 1000
Segregating sites from harseq resequencing. Reddish for higher derived frequencies
UCSC Genes Based on RefSeq, UniProt, GenBank, CCDS and Comparative Genomics
Placental Mammal Basewise Conservation by PhyloP
Human Accelerated Regions
E harseq24
SweepFinder CLR on harseq resequencing. bgrnd sfs = seasnp. gridsize = 1000
Segregating sites from harseq resequencing. Reddish for higher derived frequencies
UCSC Genes Based on RefSeq, UniProt, GenBank, CCDS and Comparative Genomics
Placental Mammal Basewise Conservation by PhyloP
Human Accelerated Regions
Figure 3.3: SweepFinder Composite Likelihood Ratio (CLR) in genomic con-text. For the indicated harseq regions, the SweepFinder results are shown along with the segregating
sites from my sequencing – more reddish for those with higher derived allele frequencies. Other tracks
show introns and exons of known genes and evolutionary conservation. The CLR peaks are not strongly
associated with the core HAR elements.
highly expressed in brain. Note that harseq25 also has a significant result on the MWU
test discussed above.
Evidence for a sweep in the harseq9 region is intriguing because it encompasses
the 42-codon long, second exon of the PTPRT gene, a phosphatase with possible roles
in the central nervous system. However, the human amino acid sequence of this exon
matches the other primates chimp, gorilla, and orangutan, except where chimp has an
obviously non-ancestral Thr>Ala substitution. The human sequence does have a single
G>A substitution near the 3' splice site just upstream of this exon, but it falls in a
position between the polypyrimidine (Py) tract and the AG acceptor site, for which the
consensus sequence across many splice sites is evenly divided among the 4 nucleotides.
The location of harseq11 on chromosome X places its evidence for a sweep in
the first intron of the 2.4Mb long dystrophin gene DMD. The evidence for a sweep in
harseq16 is offset to one end of its region, about 20kb from an apparent pseudogene
comprising a single coding exon with a 270-codon open reading frame (ORF) that is
probably derived from the Poly-A binding protein PABPC1.
The harseq24 region encompasses the second through fourth exons of the
SKAP2 gene with the strongest evidence for a sweep about 10kb from the closest exon,
but closer to a LINE transposable element that is present also in chimp, orangutan, and
rhesus macaque (but lost in gorilla).
Although I found only weak evidence for selective sweeps, and none that seems
to point directly to a mutation in a core HAR element based on the position of the
SweepFinder peak CLR values, it is important to note that while having the advantage
of robustness to demography and recombination rate, my tests would not likely have
power to detect sweeps that occurred beyond the last ∼200,000 years [82].
Methods for HAR Neighborhood Analysis
Sample Data Genomic DNA Selection
Genomic DNA for my samples was obtained from the NHGRI Sample Repos-
itory for Human Genetic Research distributed by the Coriell Institute for Medical Re-
search [Camden NJ] [1]. All of the 11 samples were chosen from the Yoruba from Ibidan
Nigeria (YRI) HapMap population. In particular, the samples were chosen as a subset
of the Seattle SNPs P2 panel [4]. The Seattle SNPs PGA-VDR (and Coriell repository
numbers were): DY01(NA18502), DY03(NA19223), DY04(NA19201), DY17(NA19143),
DY18(NA18517), DY19(NA18856), DY20(NA19239), DY21(NA18871), DY22(NA19209),
DY23(NA19152), DY24(NA19210). Sample preparation as described below was similar
for all samples, except that prior to processing, the DY01,DY03,DY04 samples were
subject to whole genome amplification (WGA) using the Repli-G Kit (Qiagen N.V, The
Netherlands) according to manufacturer's specifications. It was found that WGA caused
a loss of coverage in some isolated target regions, most notably in the 28kb section of
the harseq1 region with coordinates hg18:chr20:61,183,966–61,212,244, except for the
core HAR1 element at chr20:61,203,919-61,204,081.
Target Regions and Nimblegen Enrichment Array Design
The primary target regions consisted of 20kb extensions in both directions
from the 49 most statistically significant Human Accelerated Regions (HARs) identi-
fied as having a 5% false discovery rate [60]. Additional 40kb control regions were
chosen in neighborhoods of 13 of the set of 34,498 vertebrate conserved elements that
had extremely low likelihood ratio test (LRT) scores used to define the HARs. The
enrichment arrays were obtained from Nimblegen Systems (Madison WI) who designed
the probes on their 385K array based on my specifications of the coordinates of the 62
target genomic regions (Table A.1). Probes were chosen to tile the target regions from
both DNA strands. The design process avoided probes in highly repetitive sequences as
described previously [53, 5]. The fraction of bases in the target regions covered by the
probes ranged from 98% (harseq12) to 65% (harseq4) with a median of 88%. Details of
probed bases are available upon request.
At the time that I designed the enrichment array the main alternative tech-
nology for genomic DNA enrichment was to use oligonucleotide probes in solution. A
solution based method presumably would have better kinetic properties than the array
based method. However, the main drawback is cost. Synthesizing individual oligonu-
cleotides at a cost even as low as $0.10 per base would be prohibitively expensive to
cover target regions encompassing 2M bases. A commercially available alternative from
Agilent Technologies was also available based on their array technology in which they
separated the probes from the array after fabrication. This approach was successfully
used for enrichment of sets of exons in the human genome [61] but the cost was still
relatively high, since this technology provided 55,000 probes for about $10,000. Com-
pared to the 385,000 probes on the Nimblegen array (at about $1,000 per array) this
would have been much more expensive to cover my target regions, even considering the
fact that one order of Agilent probes would have sufficed to enrich DNA from all of my
SOLiD Barcoded Library Preparation
The library preparation of the samples for SOLiD (Applied Biosystems, Fos-
ter City, CA) sequencing generally followed the manufacturer's protocols for barcoded
SOLiD System 2.0 Fragment Library Preparation (samples DY01,DY03,DY04) and
SOLiD System 3.0 Barcoded Fragment Library Preparation (remaining samples), with
changes as necessary for enrichment on the Nimblegen arrays as noted in the following:
Samples were sheared to approximately 100bp using the Covaris S2 System Program B
(Covaris Inc, Woburn, MA). End repair was performed with End-It DNA End-Repair
Kit (Epicentre Biotechnologies, Madison, WI) per manufacturer's instructions. Single
stranded oligos for the P1 and P2-barcoded SOLiD adaptors were ordered from In-
vitrogen Corporation and annealed per the SOLiD protocol to form double stranded
adaptors, which were ligated to the end-repaired DNA fragments using the Quick Lig-
ation Kit (New England BioLabs, Ipswich, MA) per manufacturer's directions, leaving
a nick at the 3' end of each genomic DNA strand where the 5' end of the adaptor
was not phosphorylated. Samples DY01,DY03,DY04 were size selected to 150∼250bp
from a 6% polyacrylamide gel and purified via ethanol precipitation. This was followed
by nick translation and ligation mediated PCR (LMPCR) amplification (6 cycles) in a
combined reaction per the SOLiD protocols using Takara ExTaq polymerase (Takara
Bio, Madison, WI). After dividing into 10 aliquots, an additional 10 cycles of ExTaq
LMPCR amplification were performed on these samples in preparation for array hy-
bridization. Samples DY17-DY20 were first nick translated without amplification using
Pfu polymerase (Stratagene, La Jolla, CA). Samples DY21-DY24 were nick translated
and ExTaq LMPCR amplified (6 cycles). Quantitation using a DNA 2100 BioAnalyzer
(Agilent Technologies, Waldbronn, Germany) showed that the PCR associated with nick
translation prior to size selection mainly served to amplify the nick translated adaptors.
After the nick translation and any initial LMPCR, all of samples DY17-DY24 were size
selected on an E-Gel SizeSelect 2% agarose gel (Invitrogen) per manufacturer's instruc-
tions. These size selected samples were then ExTaq LMPCR amplified (6, 9 or 10 cycles)
in preparation for array hybridization.
Array hybridization to the Nimblegen arrays was performed on the Nimble-
gen Hybridization System 4 station under Mix Mode "B" for 64 to 70 hours using the
Nimblegen Sequence Capture Kit per manufacturer's instructions. Prior to hybridiza-
tion, samples DY21-24 were pooled to a total of 5.7µg. All the other samples were
hybridized individually, in amounts ranging from 1.9µg to 8.0µg per sample. For com-
petitive hybridization to probes on the array that might nonselectively bind repetitive
DNA, Human Cot-1 DNA (Invitrogen) that had been Covaris sheared to approximately
100bp was added to the hybridization mix in a 5:1 ratio by weight. Additionally, to
block the adaptor ends of the denatured, single stranded DNA fragments from binding
to each other, a 10:1 molar excess of adaptor oligos (P2 with an unused barcode sequence
and P1) was also added. At the completion of hybridization, the slide was washed with
the Nimblegen Sequence Capture Wash and Elution kit per manufacturer's directions,
and the enriched DNA was eluted with 350µL of 75◦C purified water using an affixed
SA200 SecureSeal Hybridization Chamber (Grace Bio-Labs, Bend OR). A secondary
elution with an additional 350µL was also taken. Quantitative real-time PCR (qPCR)
was performed on the eluted material to determine the rough fraction of target DNA
and to compare the primary and secondary elutions. For this purpose qPCR amplicons
within the target were compared to amplicons not in the target regions. It was gener-
ally found that the primary elution captured more than 95% of the target DNA (data
not shown). By normalizing to pre-enrichment material, and taking into account the
fact that the 2.1Mbp target region comprised approximately 0.1% of the entire human
genome, it was estimated that more than 35% of the eluted material fell in the target
regions (data not shown).
The eluted material was ExTaq LMPCR amplified (samples DY01 for 19 cy-
cles, samples DY03-04 for 15 cycles, pooled samples DY21-24 for 10 cycles, samples
DY17-20 for 12 cycles,) in preparation for the emulsion PCR step of SOLiD sequenc-
ing that was performed in the UC Santa Cruz Genome Sequencing Center. Samples
DY01,DY03,DY04 were processed with the SOLiD Version 2 system, producing 35 bases
of sequence information for each read. The remaining samples were processed with the
SOLiD Version 3 system, producing 50 bases of sequencing information for each read.
Sequence Mapping, Filtering, and Genotyping
The 35mer (DY01,DY03,DY04) or 50mer (remaining samples) sequencing reads
were mapped to the whole human genome using the bwa program [43] which generates
mappings and associated quality scores in the sam format [44] that can be processed
with the samtools suite. The bwa program is aware of the colorspace nature of the
SOLiD sequencing reads, and uses a dynamic programming algorithm to infer the best
nucleotide sequence for the read [43]. All reads were also mapped to the DNA of the
Epstein-Barr Virus, which was used to transform the Coriell cell lines from which the
supplied genomic DNA was extracted. For the female samples, the Y chromosome was
excluded from the mapping. For the male samples, the pseudo-autosomal region of the
Y chromosome was excluded from the mapping. The set of mappings for each sample
was then filtered to the regions covered by the probes on the Nimblegen enrichment
array described above. To eliminate spurious pileups caused by overamplification of
particular molecules in the library preparation process, the mapped reads were further
filtered to select at most 4 reads from each strand at a given genomic starting posi-
tion. Where there were more than 4, the 4 with the highest total read quality (not
the best mapping quality, which would bias against reads containing non-reference al-
leles) from the SOLiD instrument were selected. Between 40% and 60% of the reads
for a given sample were successfully mapped, and of those reads, between 33% and
48% mapped to bases covered by the probes on the enrichment array (Table 3.4). For
samples DY01,DY03,DY04 about 50% of the latter reads were lost in the "maximum 4
per strand" pileup elimination step, while only 11% to 23% were lost in this step in the
remaining samples (Table 3.4 and Figure 3.4). This was likely due to the difference in
LMPCR cycles used for the different samples as noted above.
To determine the coverage at each position in the target region and the con-
sensus genotypes for each sample, the command "samtools pileup -v" was used with
default parameters for its consensus calling model. Possible confounding of the geno-
types due to contamination by paralogous sequences was avoided in two ways. First,
as noted, only genotypes at positions delimited by the Nimblegen probes were used in
the analysis and these probes were designed to avoid repetitive sequences. Second, the
bwa mapping algorithm assigns low mapping quality to reads that are not genome-wide
unique, and the samtools consensus caller requires high mapping quality. To filter SNPs
from among the not homozygous reference genotypes, "samtools.pl varFilter" was run
with default parameters, except that the maximum read depth was set to 425, because
with up to 4 reads of length 50 on each strand, it was possible to get coverage of 400.
This command filters out potential SNPs when more than 2 fall within a 10bp window,
on the grounds that there might be an insertion/deletion event rather than separate
SNPs, and also filters out reads with RMS mapping quality value less than 25. A sim-
ilar quality filter was applied to the genotype calls that were homozygous reference.
Because of stochastic variation in the composition of reads from the two chromosomes
of each diploid individual, low coverage might cause an erroneous homozygous call in
a true heterozygote. Therefore a further filter restricted the subsequent analysis to the
SNP or homozygous reference calls made for a sample only at positions for which the
coverage was 35 or greater. For each target region, the count of the union of such posi-
tions across all samples is listed in Table A.1. As shown in Figure 3.5, the vast majority
of the segregating sites that remain after the application of my 35X coverage filter are
in Hardy-Weinberg Equilibrium, with only 0.3% having a p-value less than 0.05.
Table 3.4: HAR Region Mapping Statistics M: million reads. Total: the total number of
reads obtained for a given sample. Mapped: the subset of Total that mapped successfully to the human
genome. mapFrac is Mapped/Total. Probed: the subset of Mapped that mapped to the regions defined
by the probes on the Nimblegen enrichment array. probeFrac is Probed/Mapped. 4Best: the subset of
Probed that remain after filtering to allow at most 4 reads from each strand at each starting position in
the genome. 4bestFrac is 4Best/Probed. Sample DY02 was not included in the analysis because 97%
of its on-target reads were removed in the "4best" pileup elimination step.
Derived Allele Frequencies and SweepFinder
For purposes of all subsequent analysis, an ancestral allele at each position
in the target regions was determined from the Enredo-Pecan-Ortheus (EPO) pipeline
[54, 55] as published on the 1000 Genomes website [2]. This pipeline determines the
common ancestor of human and chimp at a locus by considering alignments of the
human, chimp, orangutan, and rhesus macaque genomes.
Figure 3.4: HAR Region Mapping Statistics by Sample. Graphical representation
of the data in Table 3.4. Although the starting number of reads ("Total") varied widely across the
samples, roughly the same number of reads ("4Best") were used for genotyping each sample. As noted
in Table 3.4 sample DY02 was excluded from the analysis.
Histogram of HWE p−values (two.sided)
N−total= 5835 N−pval<0.01= 1 N−pval<0.05= 17 N−pval<0.10= 56
umber of segsitesn
Figure 3.5: Histogram of Hardy-Weinberg Equilibrium Test. Histogram of the p-
values obtained in a 2-sided test for Hardy-Weinberg equilibrium across all of the SNPs used in the
analysis of my test (harseq) and control (ctlreg) regions. The vast majority conform to the HWE
From the sets of filtered genotype calls in the 11 diploid samples as described
above all the segregating sites were selected. A set of filters was applied to this list to
produce the final set of segregating site derived (i.e. non-ancestral) allele frequencies
(DAFs) for all downstream analyses. To avoid skewing the DAFs towards higher fre-
quencies, segregating sites with less than 8 chromosomal samples were eliminated. Also
eliminated were any positions with more than 2 alleles among the reference, ancestral,
or sample alleles, or where the ancestral allele was not determined by the EPO pipeline.
Lowercase values of the EPO ancestral allele, which result from various cases without
complete evidence in all species were not eliminated.
The SweepFinder program [52] was applied to the allele frequencies for the
final list of segregating sites to determine the composite likelihood ratio (CLR) of a
selective sweep at each one of a grid of 1000 positions across each target region. The
model used in this program requires a background derived allele frequency spectrum.
Two such backgrounds were used. First, all of the DAFs from all filtered segregating
sites in my sample were aggregated and used as input to the command "SweepFinder -
f", which accounts for missing data using a Broyden-Fletcher-Goldfarb-Shanno (BFGS)
algorithm. I refer to this as the "harseq" background. A second, presumably more
neutral background was obtained from the African-Derived (AD) YRI subset of 24
individuals in the Seattle SNPs P2 panel. The DAFs for all segregating sites in all 104
genes resequenced for this panel by Seattle SNPs [4] were included. As for the harseq
background, the "seasnp" background was obtained with the command "SweepFinder
-f" applied to these DAFs. The resulting "seasnp" frequency spectrum ran from 1 to 47
and was reduced to a spectrum running from 1 to 21, as needed by SweepFinder with my
data, by hypergeometric weighting the relevant components of the input Xj(j = 1.47)
allele frequencies at each target Yi(i = 1.21) allele frequency.
Dividing by the sum of the Yi in Eqn 3.1 produces a valid frequency spectrum
that sums to 1. A similar hypergeometric weighting was also required to reduce the
spectrum to a range from 1 to 19 for the harseq1,2 and ctlreg60 regions. In the latter
regions missing data reduced the maximum number of samples at the segregating sites
to 20. The reasoning behind the hypergeometric weighting can be understood from
an example. Consider a site for which, among 48 samples, there are 6 derived and
42 ancestral alleles.
At this site a subsample of 22 from the 48 samples has some
probability P3,6 of having 3 derived and 19 ancestral alleles that can be calculated as
follows. Think of the site as an urn with 6 "D" balls and 42 "A" balls. There are 48
possible subsamples. For each of these subsamples there are 6 ways to choose 3 "D"
balls and 42 ways to choose 19 "A" balls from the urn. So the probability P
hypergeometric factor 6 42/ 48. Now if there are M
6 such sites, the expected number
of such sites for which the subsample has 3 derived alleles is given by N3,6 = M6P3,6.
For other sites with at least 3 (of 48) derived alleles and at least 19 (of 48) ancestral
alleles similar calculations can be made and then summed to get the total expected
number of sites with 3 (of 22) derived alleles N3 = P MjP3,j. Dividing both sides by
the total number of sites yields equations in the form of Equation 3.1. The additional
normalization step is required because it is possible for many of the subsampled sites
to yield either 0 or 22 derived alleles, but we are only interested in sites which are
segregating in the subsample.
MWU and MK tests on Harseq and Ctlreg Regions
To determine if mutations from an ancestral weak (A or T) basepair to a
strong (G or C) basepair (W2S mutations) are more likely to spread in the popula-
tion represented by my samples, I compared W2S segregating sites to S2W segregating
sites for each target region and for the aggregate set of segregating sites. I performed
a Mann-Whitney U (MWU) test for a difference between the W2S and S2W derived
allele frequency spectra. The test was performed in the R language with the command
"wilcox.test (paired=FALSE, alternative=two.sided)". The resulting "location" param-
eter was normalized to 22 samples and is positive if W2S mutations are segregating at
higher frequencies than S2W. The resulting p-values are in Table A.2.
To determine if relatively more W2S mutations fixed along the human or chimp
lineages than are segregating in the human population represented by my samples, I
first determined the high mapping quality chimp reference bases that differ from the
human reference using reciprocal best alignments of the chimp and human genomes
[38]. This set was then restricted to the positions in my target regions for which I had
a genotype call for at least one sample with read depth of coverage of 35 or greater as
discussed above. From this set of fixed differences I removed any for which the EPO
ancestral allele was not determined as discussed above, or for which I had a segregating
site, or which appeared in dbSNP release 129 [71]. The remaining fixed differences as
well as the segregating sites were divided into W2S or S2W (or other). A McDonald
Kreitman-like (MK) test on the resulting 2x2 contingency table was performed in the
R language with the command "fisher.test (alternative=two.sided)". The resulting p-
values are in Table A.2. For the significant cases it was easy to determine from the data
in the contingency table if the fixed differences favored S2W mutations relatively more
than the segregating sites (column "S2W" in Table 3.2).
For the target regions with significant p-values on either the MWU or MK
tests, I tested whether the significance was due to a restricted locus within the region,
by removing all segregating sites and fixed differences under a mask of a given size at a
given position within the region and rerunning the test with the remaining data (Table
3.1 and Table 3.2).
MWU and MK Tests on Seattle SNPs Regions and Simulations
I downloaded the data for the 104 genes resequenced with the Seattle SNPs P2
panel. The genotypes for my 11 YRI samples (which are a subset of the P2 panel) and
the coordinates for the genotyped segregating sites were obtained from the global "pret-
tybase" file mapped to UCSC hg18 coordinates. The total set of positions genotyped
was obtained from the individual gene "genbank" files by excluding the features defined
as "Region not scanned for variation" and aligning the remaining regions to the hg18
coordinates of the full extent of the genic region sequenced specified in the associated
"ucscDataFile". Given these coordinates and the genotypes at the segregating sites,
the same techniques as described above for my resequencing data was applied to derive
ancestral alleles and human/chimp fixed differences, and to perform the MWU and MK
My data was derived from (probed) regions of a relatively tight size distribu-
tion (Table A.1). By contrast, the sizes of Seattle SNPs variation-mapped genic regions
varied widely. Some were rather small and contained few segregating sites. Therefore, I
included only genic regions with a minimum of 10kb variation-mapped and a minimum
of 40 segregating sites. Additionally, a small number of the genic regions were excluded
because of data missing from the "prettybase" file or because there were no associated
high quality reciprocal best human-chimp differences as described above, possibly be-
cause of paralogous genes in one or the other lineage. The remaining set of results for
the MWU and MK tests on 62 genic regions (Table A.3) were used for comparison to
my 49 harseq regions.
To determine if the p-values for the MWU and MK tests were accurate, I
also conducted the tests on sets of simulated data. For the MWU test, I performed
coalescent simulations using Hudson's ms program [32]. For each simulation I generated
22 samples at 85 segregating sites (the average number of W2S plus S2W segregating
sites in the 49 harseq target regions) and then randomly assigned the sites as either
W2S or S2W in Bernoulli trials using the W2S:S2W ratio from the 49 harseq regions
of 2057:2114. After calculating the MWU test p-value for each simulation, the fraction
of simulations with p-value less than a given value was computed, as well as the subset
of that fraction in which the W2S spectrum was offset towards higher derived allele
frequencies (Figure 3.6).
Fraction of Simulations by MWU−test p−value
all (slope = 0.998)
W2S (slope = 0.513)
Figure 3.6: Accuracy of MWU Test p-values. For a given p-value the fraction of neutral
model simulations having that p-value or better on the MWU test is shown. Red circles: all simulations
with a better p-value. Blue squares: the simulations with a better p-value for which the weak-to-strong
(W2S) spectrum is offset towards higher frequency. Regression lines are shown with the slope indicated
in the legend. For a given p-value the W2S-biased simulations comprise about half of all simulations.
For the MK test, for each simulation I separately derived a set of human-
chimp fixed differences and a set of segregating sites. The fixed differences were derived
using the phyloBoot program from the PHAST package [72]. I used a phylogenetic
model and substitution rate matrix derived from 4-fold degenerate amino-acid coding
synonymous sites across the genome as an unbiased neutral model. The equilibrium
GC-content of this model was adjusted to reflect the genome-wide average GC-content.
From the primate sequences so generated, I extracted positions containing a human-
chimp difference that could also be unambiguously assigned an ancestral allele based on
the macaque allele at that position. Each simulation used 335 such sites, (the average
number of fixed differences in the 49 harseq target regions) which were divided based on
whether they were W2S or S2W (or other). The segregating sites for each simulation
were derived from 101 (the average total number of segregating sites in the 49 harseq
target regions) Bernoulli trials, randomly dividing them as W2S or S2W (or other)
according the corresponding ratios in the fixed differences from all of the phyloBoot
simulations. After calculating the MK test p-value for each simulation, the fraction of
simulations with p-value less than a given value was computed, as well as the subset
of that fraction for which the ratio of W2S:S2W was higher for the simulated fixed
differences than for the simulated segregating sites (Figure 3.7).
Fraction of Simulations by MK−test p−value
all (slope = 0.865)
W2S (slope = 0.426)
Figure 3.7: Accuracy of MK Test p-values. For a given p-value the fraction of neutral
model simulations having that p-value or better on the MK test is shown. Red circles: all simulations
with a better p-value. Blue squares: the simulations with a better p-value for which the weak-to-strong
(W2S) fraction is relatively higher in fixed differences than segregating sites. Regression lines are shown
with the slope indicated in the legend. The p-values are conservative in that they overestimate the
fraction of simulations at a given significance level. For a given p-value the W2S-biased simulations
comprise about half of all simulations.
In the era of comparative genomics, strong signals of conservation across mul-
tiple species serve as signposts that can indicate regions where evolutionary forces may
be preserving functional elements that are subject to purifying selection. By contrast,
signals of positive selection pointing to adaptive changes in one lineage are harder to
find, often employing sets of polymorphic sequences from multiple individuals of the
same species.
In this work I have addressed both of these aspects of the evolutionary history
of the human genome. In my work on the Ultraconserved Elements (UCEs) I employed
resequencing data for a large number of samples that was obtained using classical meth-
ods based on PCR followed by capillary electrophoresis Sanger sequencing. From this
resequencing data I identified polymorphic positions and then derived and analyzed the
derived allele frequency (DAF) spectrum of these polymorphisms. The results of my
analysis clearly show that these elements, taken as a whole, are under strong selective
constraint. This constraint is in fact much stronger than that at play in protein coding
sequences. This result effectively laid to rest any speculation that these elements are
merely mutational cold spots, lacking any likely function. An important challenge that
was overcome in this work was the inherent ascertainment bias in the definition of these
elements wherein polymorphic positions at high derived allele frequency are likely to
prevent a region from being defined as ultraconserved. An especially striking aspect of
my results was that the bases immediately flanking the UCEs are apparently subject to
much weaker selective forces.
Although these conclusions about the UCEs are important, they are only ap-
plicable to the UCEs as an aggregate because the analysis of the DAF spectrum was
only possible using all the polymorphisms found across these elements. Most individual
UCEs have zero, one or two SNPs. Thus it was not possible to make statements about
the strength of selection of particular elements, nor to infer if any are in fact more likely
to be functional than others.
By contrast, my analyses of the HARs were applied separately to a 40kb neigh-
borhood of each region. I exploited the two recently developed techniques of genomic
enrichment and high throughput sequencing to characterize the polymorphism in a sin-
gle human population across these neighborhoods of the 49 HARs (harseq regions). I
investigated the harseq regions because the HARs were defined based on a presump-
tion that the human lineage specific fixed differences therein might have arisen due to
adaptive evolutionary forces. On the other hand, it has been emphasized by some that
the presumably evolutionarily neutral mechanism of biased gene conversion (BGC) can
influence the frequency spectra at polymorphic positions, or cause fixation of alleles in
a way that partially mimics the action of adaptive evolution. Indeed, fixation bias was
noted in connection with the limited set of human specific alleles for some of the HARs
when they were first described [59].
With the extensive novel polymorphism in my samples, I was able to carefully
characterize fixation bias — both historical and ongoing — in the harseq regions and
to conduct tests for recent selective sweeps across these regions. Like my resequencing
data for the UCEs, my harseq resequencing data eliminates issues of SNP ascertainment
bias that could have skewed previous investigations of polymorphism near HARs.
Consistent with published reports [59, 17, 26], I found evidence of historical
W2S fixation bias in harseq regions. Using a MK test, I compared the proportion of
W2S mutations among already fixed substitutions on the human or chimp lineage to
that among the still segregating sites in my samples. I found that 11 of the 49 regions
show statistically significant evidence of historical bias in allele fixation, with all but
one favoring W2S fixation. These results strengthen and expand previous findings by
identifying signals for W2S bias in much larger regions flanking the core HAR regions
in an ascertainment-free population sample.
But my work goes beyond previous approaches by also looking at ongoing
W2S fixation bias in the segregating site frequency spectrum. I performed a MWU
test using only sites that are still segregating in the human population, separating out
W2S from S2W mutations. This second test is designed to detect a phenomenon of bias
that is currently driving W2S mutations to higher frequency in the population than
S2W mutations. I found statistically significant evidence for this bias (and none in the
opposite direction) in the regions flanking 11 of 49 HARs.
For both tests, I showed that the core HAR element is generally not the main
source of the signal that was detected, since the signal usually remains strong even
when I mask out the central 1kb or even 5kb of the region. This is not consistent
with BGC due to a recombination hot spot that has remained in the same location for
millions of years, because the length scale of the effect of BGC is set by the length of
the heteroduplex tract formed during recombination that needs to be repaired, which
is thought to be ∼500bp (e.g. [11, 65, 25]). However, it is consistent with a model in
which the location of recombination hotspots drift fairly rapidly over evolutionary time
scales, but may be denser in some regions [64, 58, 48, 15, 57].
It is noteworthy that there was little overlap in the regions identified by these
two tests, one for older W2S fixations and the other for present day forces toward
fixation, with a total of 20 found in one or the other.
This is consistent with the
hypothesis that the regional focus of BGC, which may be recombination hot spots,
drifts significantly on a time scale of many hundreds of thousands or millions of years.
Another explanation for W2S fixation bias near HARs is selection for increased
GC-content or individual fitness-improving GC alleles. To investigate these hypotheses
and to attempt to disentangle the possible roles of BGC and positive selection in shaping
the HARs, I applied a recently developed powerful method for detecting selective sweeps.
Selection was previously investigated in much larger (∼500kb) regions using more sparse
polymorphic loci [82]. That study found 101 regions with strong evidence for a selective
sweep within 100kb of a known gene. Here, I found only weak evidence for selection
comprising 5 possible candidates for such sweeps among 49 target regions. As we are
dealing with a lineage-specific evolutionary period of about 5 million years, and these
tests can only see back a few hundred thousand years, it is quite possible that the original
signal for selective sweeps in these regions has already decayed beyond our ability to
recognize it in human population genetic data.
Consistent with the idea that HAR regions may have experienced positive se-
lection too long ago to be detected with population genetic methods, very few positively
selected regions in the human lineage have been identified to date, despite the existence
of numerous public databases. Selective sweeps that have been identified have typically
been the product of very recent events in human history, such as dairy farming affect-
ing the lactase gene [8] or climate differences influencing a salt sensitivity variant [78].
Such environmental or cultural changes result in differences in the genetic makeup of
disparate human populations, and such differences can be exploited to find evidence of
recent, possibly still ongoing, selective sweeps.
Despite the evidence that the unusually high level of recent substitution in the
more extreme HAR elements such as HAR1 is at least partly due to the non-selective
process of BGC, there is ample evidence that some of these genomic elements remain
functional, and thus that the effect of BGC was to mutationally stress but not destroy
these elements. HAR1 shows a very strong pattern of compensatory substitutions within
its RNA helix structures, indicating a selective force to maintain these helix structures.
The W2S substitutions all strengthen the RNA helices of HAR1, and in one case, a
substitution appears to extend one of them. Thus, although the nominally non-selective
GC-biased evolutionary force that I have documented may reduce the fitness advantage
of some HARs, it is also possible that it has enabled some of these regions to attain a
configuration subject to positive selection for higher fitness in humans.
In summary, I have investigated two sets of highly conserved genomic regions.
For both sets I have used polymorphism data from the present day human population
to elucidate evolutionary forces at work. I have shown that purifying selection is likely
preventing the UCEs from diverging in our species from the ancestral state, and that
some of the base substitutions now evident in the HARs are likely to be the result of a
process that favors fixation of mutations from weak to strong basepairs. As larger data
sets of human genetic variation become available it should be possible to expand these
results to the genome-wide scale.
Supplementary Tables
HAR Region Statistics. Characteristics of the target (harseq1-
49) and control (ctlreg50-62) regions in the HAR study. hg18pos: chromosome
coordinates in the hg18 human assembly. bpExp: the number of basepairs
used in the experimental analysis. This number starts with the count of bases
in the region that were included in probes in the Nimblegen array, but only
includes bases for which sequencing coverage of at least 35X was reached for
at least one sample, after removing pileups of more than 4 reads at the same
starting point. AncGC: the GC percentage of the ancestral bases in bpExp.
θπ: estimate of the neutral mutation parameter (= 4Neu x104) via nucleotide
diversity. θW : estimate of the neutral mutation parameter via the number of
segregating sites. tajD: measure of skew in the folded site frequency spectrum.
faywuH: measure of skew in the unfolded site frequency spectrum.
HAR Region Test Results
segsites (N): total number of
segregating sites from my samples. fixed diffs (N): number of positions for
which I had sufficient coverage to make a genotype call in at least one sample
that have a fixed difference between the human and chimp reference, exclud-
ing segregating sites and sites in dbSNP129.
w2s,s2w: subsets with weak-
to-strong or strong-to-weak mutations from the ancestral basepair. MWUp-
value,MWUoffset: p-value and estimated location offset (the latter normalized
to 22 samples and estimated under the assumption that w2s and s2w have
identical distributions) from Mann-Whitney U test comparing the frequency
spectra of segsites w2s vs s2w, with positive values indicating W2S mutations
shifted to higher derived allele frequencies. MKpval: p-value from McDonald-
Kreitman-like test comparing the counts in w2s,s2w categories of segregating
and fixed sites. SF: p-values for maximum CLR from SweepFinder run on the
region with a background SFS derived from segregating sites pooled from all
regions (Hpvalue) or from the same population segregating in Seattle SNPs
resequenced regions (Spvalue).
Seattle SNPs Gene Region MWU and MK Test Results
basepairs(Gene): total length of genic region. basepairs(Map): sum of base-
pairs in regions that were mapped for variation within genic region. AncGC:
the GC percentage of the ancestral bases in Map basepairs. segsites (N): to-
tal number of segregating sites in Map basepairs for my samples. fixed diffs
(N): number of positions in Map basepairs that have a fixed difference between
the human and chimp reference, excluding segregating sites and sites in db-
SNP129. w2s,s2w: subsets with weak-to-strong or strong-to-weak mutations
from the ancestral basepair. Other annotations as in Table A.2.
[1] Coriell cell repositories. http://ccr.coriell.org
[2] 1000 genomes: A deep catalog of human genetic variation. http://1000genomes.
org , August 2009.
[3] H. Akashi. Inferring the fitness effects of DNA mutations from polymorphism and
divergence data: statistical power to detect directional selection under stationarityand free recombination. Genetics, 151:221–238, Jan 1999.
[4] J.M. Akey, M.A. Eberle, M.J. Rieder, C.S. Carlson, M.D. Shriver, D.A. Nickerson,
and L. Kruglyak. Population history and natural selection shape patterns of geneticvariation in 132 genes. PLoS Biol., 2:e286, Oct 2004.
[5] T.J. Albert, M.N. Molla, D.M. Muzny, L. Nazareth, D. Wheeler, X. Song, T.A.
Richmond, C.M. Middle, M.J. Rodesch, C.J. Packard, G.M. Weinstock, and R.A.
Gibbs. Direct selection of human genomic loci by microarray hybridization. Nat.
Methods, 4:903–905, Nov 2007.
[6] Applied Biosystems. SOLiD System Brochure. System Brochure, Jun 2008.
[7] G. Bejerano, M. Pheasant, I. Makunin, S. Stephen, W.J. Kent, J.S. Mattick, and
D. Haussler. Ultraconserved elements in the human genome. Science, 304:1321–1325, May 2004.
[8] T. Bersaglieri, P. C. Sabeti, N. Patterson, T. Vanderploeg, S. F. Schaffner, J. A.
Drake, M. Rhodes, D. E. Reich, and J. N. Hirschhorn. Genetic signatures of strongrecent positive selection at the lactase gene. Am. J. Hum. Genet., 74:1111–1120,Jun 2004.
[9] C. P. Bird, B. E. Stranger, M. Liu, D. J. Thomas, C. E. Ingle, C. Beazley, W. Miller,
M. E. Hurles, and E. T. Dermitzakis. Fast-evolving noncoding sequences in thehuman genome. Genome Biol., 8:R118, 2007.
[10] M. Blanchette, W.J. Kent, C. Riemer, .L Elnitski, A.F.A. Smit, K.M. Roskin,
R. Baertsch, K. Rosenbloom, H. Clawson, E.D. Green, D. Haussler, and W. Miller.
Aligning multiple genomic sequences with the threaded blockset aligner. GenomeResearch, 14:708–715, 2004.
[11] H. L. Blanton, S. J. Radford, S. McMahan, H. M. Kearney, J. G. Ibrahim, and
J. Sekelsky. REC, Drosophila MCM8, drives formation of meiotic crossovers. PLoSGenet., 1:e40, Sep 2005.
[12] S.P. Brooks and A. Gelman. General methods for monitoring convergence of iter-
ative simulations. Journal of Computational and Graphical Statistics, 7:434–455,1997.
[13] C.D. Bustamante, J. Wakeley, S. Sawyer, and D.L. Hartl. Directional selection and
the site-frequency spectrum. Genetics, 159:1779–1788, Dec 2001.
[14] Francis S. Collins, Lisa D. Brooks, and Aravinda Chakravarti. A DNA polymor-
phism discovery resource for research on human genetic variation. Genome Re-search, 8:1229–1231, 1998.
[15] G. Coop, X. Wen, C. Ober, J. K. Pritchard, and M. Przeworski. High-resolution
mapping of crossovers reveals extensive variation in fine-scale recombination pat-terns among humans. Science, 319:1395–1398, Mar 2008.
[16] Jared A Drake, Christine Bird, James Nemesh1, Daryl J Thomas, Christopher
Newton-Cheh, Alexandre Reymond, Laurent Excoffier, Homa Attar, Stylianos EAntonarakis, Emmanouil T Dermitzakis, and Joel N Hirschhorn. Conserved non-coding sequences are selectively constrained and not mutation cold spots. NatureGenetics, 38:223–227, 2006.
[17] T.R. Dreszer, G.D. Wall, D. Haussler, and K.S. Pollard. Biased clustered substitu-
tions in the human genome: the footprints of male-driven biased gene conversion.
Genome Res., 17:1420–1430, Oct 2007.
[18] L. Duret and P.F. Arndt. The impact of recombination on nucleotide substitutions
in the human genome. PLoS Genet., 4:e1000071, May 2008.
[19] L. Duret and N. Galtier.
Comment on "Human-specific gain of function in a
developmental enhancer". Science, 323:714; author reply 714, Feb 2009.
[20] W. Enard, M. Przeworski, S. E. Fisher, C. S. Lai, V. Wiebe, T. Kitano, A. P.
Monaco, and S. Pbo. Molecular evolution of FOXP2, a gene involved in speechand language. Nature, 418:869–872, Aug 2002.
[21] B Ewing and P Green. Basecalling of automated sequencer traces using phred. II.
error probabilities. Genome Research, 8:186–194, 1998.
[22] B Ewing, L Hillier, M Wendl, and P Green. Basecalling of automated sequencer
traces using phred. I. accuracy assessment. Genome Research, 8:175–185, 1998.
[23] J. C. Fay and C. I. Wu. Hitchhiking under positive Darwinian selection. Genetics,
155:1405–1413, Jul 2000.
[24] R. A. Fisher. The Genetical Theory of Natural Selection. Clarendon Press, Oxford,
[25] L. Frisse, R. R. Hudson, A. Bartoszewicz, J. D. Wall, J. Donfack, and A. Di Rienzo.
Gene conversion and different population histories may explain the contrast be-tween polymorphism and linkage disequilibrium levels.
Am. J. Hum. Genet.,
69:831–843, Oct 2001.
[26] N. Galtier and L. Duret. Adaptation or biased gene conversion? Extending the
null hypothesis of molecular evolution. Trends Genet., 23:273–277, Jun 2007.
[27] N. Galtier, L. Duret, S. Glmin, and V. Ranwez. GC-biased gene conversion pro-
motes the fixation of deleterious amino acid changes in primates. Trends Genet.,25:1–5, Jan 2009.
[28] John H. Gillespie. Population Genetics: A Concise Guide. JHU Press, 2004.
[29] D Gordon, C Abajian, and Green P. Consed: a graphical tool for sequence finishing.
Genome Research, 8:195–202, 1998.
[30] E. Hodges, Z. Xuan, V. Balija, M. Kramer, M.N. Molla, S.W. Smith, C.M. Middle,
M.J. Rodesch, T.J. Albert, G.J. Hannon, and W.R. McCombie. Genome-wide insitu exon capture for selective resequencing. Nat. Genet., 39:1522–1527, Dec 2007.
[31] F. Hsu, W.J. Kent, H. Clawson, R.M. Kuhn, M. Diekhans, and D Haussler. The
UCSC known genes. Bioinformatics, 22:1036–46, 2006.
[32] Richard R. Hudson. Generating samples under a Wright-Fisher neutral model of
genetic variation. Bioinformatics, 18:337–338, 2002.
[33] Mark A. Jobling, Matthew Hurles, and Chris Tyler-Smith. Human Evolutionary
Genetics. Garland Publishing, 2004.
[34] N.L. Kaplan, R.R. Hudson, and C.H. Langley. The hitchhiking effect revisited.
Genetics, 123:887–899, Dec 1989.
[35] S. Katzman. Materials and Methods: Human genome ultraconserved elements are
ultraselected. Online Supporting Material, 2007. http://compbio.soe.ucsc.edu/papers/ultras meth v31.pdf
[36] S. Katzman, A.D. Kern, G. Bejerano, G. Fewell, L. Fulton, R.K. Wilson, S.R.
Salama, and D. Haussler. Human genome ultraconserved elements are ultraselected.
Science, 317:915, Aug 2007.
[37] Peter D. Keightley, Martin J. Lercher, and Adam Eyre-Walker.
widespread degradation of gene control regions in hominid genomes. PLoS Biology,3:e42, 2005.
[38] W. James Kent, Robert Baertsch, Angie Hinrichs, Webb Miller, and David Haus-
sler. Evolution's cauldron: duplication, deletion, and rearrangement in the mouseand human genomes. Proc Natl Acad Sci U S A, 100:11484–11489, 2003.
[39] A. D. Kern. Correcting the site frequency spectrum for divergence-based ascertain-
ment. PLoS ONE, 4:e5152, 2009.
[40] Y. Kim and W. Stephan. Detecting a local signature of genetic hitchhiking along
a recombining chromosome. Genetics, 160:765–777, Feb 2002.
[41] A. Kong, D. F. Gudbjartsson, J. Sainz, G. M. Jonsdottir, S. A. Gudjonsson,
B. Richardsson, S. Sigurdardottir, J. Barnard, B. Hallbeck, G. Masson, A. Shlien,S. T. Palsson, M. L. Frigge, T. E. Thorgeirsson, J. R. Gulcher, and K. Stefans-son. A high-resolution recombination map of the human genome. Nat. Genet.,31:241–247, Jul 2002.
[42] Gregory V Kryukov, Steffen Schmidt, and Shamil Sunyaev. Small fitness effect
of mutations in highly conserved non-coding regions. Human Molecular Genetics,14:2221–2229, 2005.
[43] H. Li and R. Durbin. Fast and accurate short read alignment with Burrows-Wheeler
transform. Bioinformatics, 25:1754–1760, Jul 2009.
[44] H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth,
G. Abecasis, and R. Durbin. The Sequence Alignment/Map format and SAMtools.
Bioinformatics, 25:2078–2079, Aug 2009.
[45] G. Marais. Biased gene conversion: implications for genome and sex evolution.
Trends Genet., 19:330–338, Jun 2003.
[46] J. Maynard Smith and J. Haigh. The hitch-hiking effect of a favourable gene.
Genet. Res., 23:23–35, Feb 1974.
[47] J. Meunier and L. Duret. Recombination drives the evolution of GC-content in the
human genome. Mol. Biol. Evol., 21:984–990, Jun 2004.
[48] S. Myers, C. C. Spencer, A. Auton, L. Bottolo, C. Freeman, P. Donnelly, and
G. McVean. The distribution and causes of meiotic recombination in the humangenome. Biochem. Soc. Trans., 34:526–530, Aug 2006.
[49] Masatoshi Nei and Takashi Gojobori. Simple methods for estimating the numbers
of synonymous and nonsynonymous nucleotide substitutions. Molecular Biologyand Evolution, 3:418–426, 1986.
[50] D.A. Nickerson, V.O. Tobe, and S.L Taylor. Polyphred: automating the detection
and genotyping of single nucleotide substitutions using fluorescence-based rese-quencing. Nucleic Acids Research, 25:2745–2751, 1997.
[51] R. Nielsen, M. J. Hubisz, and A. G. Clark. Reconstituting the frequency spectrum
of ascertained single-nucleotide polymorphism data. Genetics, 168:2373–2382, Dec2004.
[52] R. Nielsen, S. Williamson, Y. Kim, M.J. Hubisz, A.G. Clark, and C. Bustamante.
Genomic scans for selective sweeps using SNP data. Genome Res., 15:1566–1575,Nov 2005.
[53] D.T. Okou, K.M. Steinberg, C. Middle, D.J. Cutler, T.J. Albert, and M.E. Zwick.
Microarray-based genomic selection for high-throughput resequencing. Nat. Meth-ods, 4:907–909, Nov 2007.
[54] B. Paten, J. Herrero, K. Beal, S. Fitzgerald, and E. Birney. Enredo and Pecan:
genome-wide mammalian consistency-based multiple alignment with paralogs.
Genome Res., 18:1814–1828, Nov 2008.
[55] B. Paten, J. Herrero, S. Fitzgerald, K. Beal, P. Flicek, I. Holmes, and E. Birney.
Genome-wide nucleotide-level mammalian ancestor reconstruction. Genome Res.,18:1829–1843, Nov 2008.
[56] L. A. Pennacchio, N. Ahituv, A. M. Moses, S. Prabhakar, M. A. Nobrega,
M. Shoukry, S. Minovitsky, I. Dubchak, A. Holt, K. D. Lewis, I. Plajzer-Frick,J. Akiyama, S. De Val, V. Afzal, B. L. Black, O. Couronne, M. B. Eisen, A. Visel,and E. M. Rubin. In vivo enhancer analysis of human conserved non-coding se-quences. Nature, 444:499–502, Nov 2006.
[57] A. D. Peters. A combination of cis and trans control can solve the hotspot conver-
sion paradox. Genetics, 178:1579–1593, Mar 2008.
[58] M. Pineda-Krch and R. J. Redfield. Persistence and loss of meiotic recombination
hotspots. Genetics, 169:2319–2333, Apr 2005.
[59] K.S. Pollard, S.R. Salama, B. King, A.D. Kern, T. Dreszer, S. Katzman, A. Siepel,
J.S. Pedersen, G. Bejerano, R. Baertsch, K.R. Rosenbloom, J. Kent, and D. Haus-sler. Forces shaping the fastest evolving regions in the human genome. PLoSGenet., 2:e168, Oct 2006.
[60] K.S. Pollard, S.R. Salama, N. Lambert, M.A. Lambot, S. Coppens, J.S. Peder-
sen, S. Katzman, B. King, C. Onodera, A. Siepel, A.D. Kern, C. Dehay, H. Igel,M. Ares, P. Vanderhaeghen, and D. Haussler. An RNA gene expressed duringcortical development evolved rapidly in humans. Nature, 443:167–172, Sep 2006.
[61] G.J. Porreca, K. Zhang, J.B. Li, B. Xie, D. Austin, S.L. Vassallo, E.M. LeProust,
B.J. Peck, C.J. Emig, F. Dahl, Y. Gao, G.M. Church, and J. Shendure. Multiplexamplification of large sets of human exons. Nat. Methods, 4:931–936, Nov 2007.
[62] S. Prabhakar, J. P. Noonan, S. Pbo, and E. M. Rubin. Accelerated evolution of
conserved noncoding sequences in humans. Science, 314:786, Nov 2006.
[63] S. Prabhakar, A. Visel, J.A. Akiyama, M. Shoukry, K.D. Lewis, A. Holt, I. Plajzer-
Frick, H. Morrison, D.R. Fitzpatrick, V. Afzal, L.A. Pennacchio, E.M. Rubin, andJ.P. Noonan. Human-specific gain of function in a developmental enhancer. Science,321:1346–1350, Sep 2008.
[64] S. E. Ptak, A. D. Roeder, M. Stephens, Y. Gilad, S. Pbo, and M. Przeworski.
Absence of the TAP2 human recombination hotspot in chimpanzees. PLoS Biol.,2:e155, Jun 2004.
[65] S. J. Radford, S. McMahan, H. L. Blanton, and J. Sekelsky. Heteroduplex DNA
in meiotic recombination in Drosophila mei-9 mutants. Genetics, 176:63–72, May2007.
[66] M. V. Rockman, M. W. Hahn, N. Soranzo, F. Zimprich, D. B. Goldstein, and G. A.
Wray. Ancient and recent positive selection transformed opioid cis-regulation inhumans. PLoS Biol., 3:e387, Dec 2005.
[67] Steve Rozen and Helen J. Skaletsky. Primer3 on the www for general users and
for biologist programmers. In S Krawetz and S Misener, editors, BioinformaticsMethods and Protocols: Methods in Molecular Biology., pages 365–386. HumanaPress, Totowa, NJ, 2000.
[68] P.C. Sabeti, D.E. Reich, J.M. Higgins, H.Z. Levine, D.J. Richter, S.F. Schaffner,
S.B. Gabriel, J.V. Platko, N.J. Patterson, G.J. McDonald, H.C. Ackerman, S.J.
Campbell, D. Altshuler, R. Cooper, D. Kwiatkowski, R. Ward, and E.S. Lander.
Detecting recent positive selection in the human genome from haplotype structure.
Nature, 419:832–837, Oct 2002.
[69] P.C. Sabeti, P. Varilly, B. Fry, J. Lohmueller, E. Hostetter, C. Cotsapas, X. Xie,
E.H. Byrne, S.A. McCarroll, R. Gaudet, S.F. Schaffner, and E.S. Lander. Genome-wide detection and characterization of positive selection in human populations.
Nature, 449:913–918, Oct 2007.
[70] S. A. Sawyer and D. L. Hartl. Population genetics of polymorphism and divergence.
Genetics, 132:1161–1176, 1992.
[71] S.T. Sherry, M.H. Ward, M. Kholodov, J. Baker, L. Phan, E.M. Smigielski, and
K. Sirotkin. dbSNP: the NCBI database of genetic variation. Nucleic Acids Res.,29:308–311, Jan 2001.
[72] A. Siepel and D. Haussler. Phylogenetic estimation of context-dependent substitu-
tion rates by maximum likelihood. Mol. Biol. Evol., 21:468–488, Mar 2004.
[73] Elizabeth M. Smigielski, Karl Sirotkin, Minghong Ward, and Stephen T. Sherry.
dbsnp: a database of single nucleotide polymorphisms. Nucleic Acids Research,28:352–355, 2000.
[74] J.N. Strathern, B.K. Shafer, and C.B. McGill. DNA synthesis errors associated
with double-strand-break repair. Genetics, 140:965–972, Jul 1995.
[75] F. Tajima. Statistical method for testing the neutral mutation hypothesis by DNA
polymorphism. Genetics, 123:585–595, Nov 1989.
[76] Fumio Tajima. Evolutionary relationship of DNA sequences in finite populations.
Genetics, 105:437–460, 1983.
[77] K. Tang, K.R. Thornton, and M. Stoneking. A New Approach for Using Genome
Scans to Detect Recent Positive Selection in the Human Genome. PLoS Biol.,5:e171, Jun 2007.
[78] E. E. Thompson, H. Kuttab-Boulos, D. Witonsky, L. Yang, B. A. Roe, and
A. Di Rienzo. CYP3A variation and the evolution of salt-sensitivity variants. Am.
J. Hum. Genet., 75:1059–1069, Dec 2004.
[79] B.F. Voight, S. Kudaravalli, X. Wen, and J.K. Pritchard. A map of recent positive
selection in the human genome. PLoS Biol., 4:e72, Mar 2006.
[80] G. A. Watterson. On the number of segregating sites in genetical models without
recombination. Theoretical Population Biology, 7:256–276, 1975.
[81] S. H. Williamson, R. Hernandez, A. Fledel-Alon, L. Zhu, R. Nielsen, and C. D. Bus-
tamante. Simultaneous inference of selection and population growth from patternsof variation in the human genome. Proc. Natl. Acad. Sci. U.S.A., 102:7882–7887,May 2005.
[82] S. H. Williamson, M. J. Hubisz, A. G. Clark, B. A. Payseur, C. D. Bustamante,
and R. Nielsen. Localizing recent adaptive evolution in the human genome. PLoSGenet., 3:e90, Jun 2007.
[83] A. Woolfe, M. Goodson, D. K. Goode, P. Snell, G. K. McEwen, T. Vavouri, S. F.
Smith, P. North, H. Callaway, K. Kelly, K. Walter, I. Abnizova, W. Gilks, Y. J.
Edwards, J. E. Cooke, and G. Elgar. Highly conserved non-coding sequences areassociated with vertebrate development. PLoS Biol., 3:e7, Jan 2005.
[84] S Wright. Evolution and the Genetics of Populations, volume Vol. 2: The Theory
of Gene Frequencies. University of Chicago Press, Chicago, 1969.
Source: https://compbio.soe.ucsc.edu/theses/Sol-Katzman-PhD-2010.pdf
B. Raj Kumar et al. / JGTPS/ 5(3)-(2014) 1827- 1832 Journal of Global Trends in Pharmaceutical Sciences Journal home page: www.jgtps.com A NEW VALIDATED STABILITY-INDICATING RP-HPLC METHOD FOR THE DETERMINATION OF TELBIVUDINE B. Raj Kumar1* Dr. K. V. Subrahmanyam2 Telbivudine 5'–triphosphate inhibits HBV DNA polymerase (reverse
Evidence-Based RegenerativeInjection Therapy (Prolotherapy)in Sports Medicine K. Dean Reeves, MD; Bradley D. Fullerton, MD, FAAPMR;and Gaston Topol, MD use of prolotherapy, which is also called regenerative injection ther-apy (RIT), is expected to greatly accelerate in the next decade. This . The treatment of sports injuries to the point of restoration of full chapter will cover the pathology of injury; the current treatment meth-
