Joint Estimation of Contamination,Error and Demography for Nuclear DNA from Ancient Humans

When sequencing an ancient DNA sample from a hominin fossil, DNA from present-day humans involved in excavation and extraction will be sequenced along with the endogenous material. This type of contamination is problematic for downstream analyses as it will introduce a bias towards the population of the contaminating individual(s). Quantifying the extent of contamination is a crucial step as it allows researchers to account for possible biases that may arise in downstream genetic analyses. Here, we present an MCMC algorithm to co-estimate the contamination rate, sequencing error rate and demographic parameters—including drift times and admixture rates—for an ancient nuclear genome obtained from human remains, when the putative contaminating DNA comes from present-day humans. We assume we have a large panel representing the putative contaminant population (e.g. European, East Asian or African). The method is implemented in a C++ program called ‘Demographic Inference with Contamination and Error’ (DICE). We applied it to simulations and genome data from ancient Neanderthals and modern humans. With reasonable levels of genome sequence coverage (>3X), we find we can recover accurate estimates of all these parameters, even when the contamination rate is as high as 50%.     

The joint allele-frequency spectrum in closely related species

We develop the theory for computing the joint frequency spectra of alleles in two closely related species. We allow for arbitrary population growth in both species after they had a common ancestor. We focus on the case in which a single chromosome is sequenced from one of the species. We use classical diffusion theory to show that, if the ancestral species was at equilibrium under mutation and drift and a chromosome from one of the descendant species carries the derived allele, the frequency spectrum in the other species is uniform, independently of the demographic history of both species. We also predict the expected densities of segregating and fixed sites when the chromosome from the other species carries the ancestral allele. We compare the predictions of our model with the site-frequency spectra of SNPs in the four HapMap populations of humans when the nucleotide present in the Neanderthal DNA sequence is ancestral or derived, using the chimp genome as the outgroup.     

A method for estimating the intensity of overdominant selection from the distribution of allele frequencies

A method is proposed for estimating the intensity of overdominant selection scaled by the effective population size, S = 2Ns, from allele frequencies. The method is based on the assumption that, with strong overdominant selection, allele frequencies are nearly at their deterministic equilibrium values and that, to a first approximation, deviations depend only on S. Simulations verify that reasonably accurate estimates of S can be obtained for realistic sample sizes. The method is applied to data from several loci in the major histocompatibility complex (Mhc) in numerous human populations. For alleles distinguished by both serological typing and the sequence of the peptide-binding region, our estimates of S are comparable to those obtained by analysis of DNA sequences in showing that selection is strongest on HLA-B and weaker on HLA-A, HLA-DRB1, and HLA-DQA1. The intensity of selection on HLA-B varied considerably among populations. Two populations, Native American and Inuit, showed an excess rather than a deficiency in homozygosity. Comparable estimates of S were obtained for alleles at Mhc class II loci distinguished by serological reactions (serotyping) and by differences in the amino acid sequences of the peptide-binding region (molecular typing). A comparison of two types of data for DQA1 and DRB1 showed that serotyping led to generally lower estimates of S.     

Likelihood analysis of ongoing gene flow and historical association

Abstract.— We develop a Monte Carlo-based likelihood method for estimating migration rates and population divergence times from data at unlinked loci at which mutation rates are sufficiently low that, in the recent past, the effects of mutation can be ignored. The method is applicable to restriction fragment length polymorphisms (RFLPs) and single nucleotide polymorphisms (SNPs) sampled from a subdivided population. The method produces joint maximum-likelihood estimates of the migration rate and the time of population divergence, both scaled by population size, and provides a framework in which to test either for no ongoing gene flow or for population divergence in the distant past. We show the method performs well and provides reasonably accurate estimates of parameters even when the assumptions under which those estimates are obtained are not completely satisfied. Furthermore, we show that, provided that the number of polymorphic loci is sufficiently large, there is some power to distinguish between ongoing gene flow and historical association as causes of genetic similarity between pairs of populations.     

Incorporating genotypes of relatives into a test of linkage disequilibrium.

Genetic data from autosomal loci in diploids generally consist of genotype data for which no phase information is available, making it difficult to implement a test of linkage disequilibrium. In this paper, we describe a test of linkage disequilibrium based on an empirical null distribution of the likelihood of a sample. Information on the genotypes of related individuals is explicitly used to help reconstruct the gametic phase of the independent individuals. Simulation studies show that the present approach improves on estimates of linkage disequilibrium gathered from samples of completely independent individuals but only if some offspring are sampled together with their parents. The failure to incorporate some parents sharply decreases the sensitivity and accuracy of the test. Simulations also show that for multiallelic data (more than two alleles) our testing procedure is not as powerful as an exact test based on known haplotype frequencies, owing to the interaction between departure from Hardy-Weinberg equilibrium and linkage disequilibrium.     

The effective founder effect in a spatially expanding population

The gradual loss of diversity and the establishment of clines in allele frequencies associated with range expansions are patterns observed in many species, including humans. These patterns can result from a series of founder events occurring as populations colonize previously unoccupied areas. We develop a model of an expanding population and, using a branching process approximation, show that spatial gradients reflect different amounts of genetic drift experienced by different subpopulations. We then use this model to measure the net average strength of the founder effect, and we demonstrate that the predictions from the branching process model fit simulation results well. We further show that estimates of the effective founder size are robust to potential confounding factors such as migration between subpopulations. We apply our method to data from Arabidopsis thaliana. We find that the average founder effect is approximately three times larger in the Americas than in Europe, possibly indicating that a more recent, rapid expansion occurred.     

The geographic spread of the CCR5 Delta32 HIV-resistance allele

The Delta32 mutation at the CCR5 locus is a well-studied example of natural selection acting in humans. The mutation is found principally in Europe and western Asia, with higher frequencies generally in the north. Homozygous carriers of the Delta32 mutation are resistant to HIV-1 infection because the mutation prevents functional expression of the CCR5 chemokine receptor normally used by HIV-1 to enter CD4+ T cells. HIV has emerged only recently, but population genetic data strongly suggest Delta32 has been under intense selection for much of its evolutionary history. To understand how selection and dispersal have interacted during the history of the Delta32 allele, we implemented a spatially explicit model of the spread of Delta32. The model includes the effects of sampling, which we show can give rise to local peaks in observed allele frequencies. In addition, we show that with modest gradients in selection intensity, the origin of the Delta32 allele may be relatively far from the current areas of highest allele frequency. The geographic distribution of the Delta32 allele is consistent with previous reports of a strong selective advantage (>10%) for Delta32 carriers and of dispersal over relatively long distances (>100 km/generation). When selection is assumed to be uniform across Europe and western Asia, we find support for a northern European origin and long-range dispersal consistent with the Viking-mediated dispersal of Delta32 proposed by G. Lucotte and G. Mercier. However, when we allow for gradients in selection intensity, we estimate the origin to be outside of northern Europe and selection intensities to be strongest in the northwest. Our results describe the evolutionary history of the Delta32 allele and establish a general methodology for studying the geographic distribution of selected alleles.     

A population-genetic test of founder effects and implications for Ashkenazi Jewish diseases

A founder effect can account for the presence of an allele at an unusually high frequency in an isolated population if the allele is selectively neutral and if all copies are identical by descent with a copy that either was carried by a founder individual or arose by mutation later. Here, a statistical test of both aspects of the founder-effect hypothesis is developed. The test is performed by a modified version of a program that implements the Slatkin-Bertorelle test of neutrality. The test is applied to several disease-associated alleles found predominantly in Ashkenazi Jews. Despite considerable uncertainty about the demographic history of Ashkenazi Jews and their ancestors, available genetic data are consistent with a founder effect resulting from a severe bottleneck in population size between a.d. 1100 and a.d. 1400 and an earlier bottleneck in a.d. 75, at the beginning of the Jewish Diaspora. The relatively high frequency of alleles causing four different lysosomal storage disorders, including Tay-Sachs disease and Gaucher disease, can be accounted for if the disease-associated alleles are recessive in their effects on reproductive fitness.     

General Epistatic Models of the Risk of Complex Diseases

The range of possible gene interactions in a multilocus model of a complex inherited disease is studied by exploring genotype-specific risks subject to the constraint that the allele frequencies and marginal risks are known. We quantify the effect of gene interactions by defining the interaction ratio, , where K_R is the recurrence risk to relatives with relationship R for the true model and is the recurrence risk to relatives for a multiplicative model with the same marginal risks. We use a Markov chain Monte Carlo (MCMC) procedure to sample from the space of possible models. We find that the average of C_R increases with the number of loci for both low frequency (p = 0.03) and higher frequency (p = 0.25) causative alleles. Furthermore, the probability that C_R > 1 is nearly 1. Similar results are obtained when more weight is given to risk models that are closer to the comparable multiplicative model. These results imply that, in general, gene interactions will result in greater heritability of a complex inherited disease than is expected on the basis of a multiplicative model of interactions and hence may provide a partial explanation for the problem of missing heritability of complex diseases.ALTHOUGH many genome-wide association studies (GWAS) have been performed and have found hundreds of SNPs associated with higher risk of complex inherited diseases, those SNPs so far account for only a small fraction of the inherited risk of those diseases (Altshuler et al. 2008). Several not mutually exclusive explanations have been proposed for the “missing heritability,” i.e., the heritability that is not yet accounted for by SNPs found in GWAS (Manolio et al. 2009): (i) common alleles of small effect that have not been found because GWAS done so far have been underpowered, (ii) low-frequency alleles of moderate effect that are difficult to find using HapMap SNPs, (iii) rare copy-number variants that are not in strong linkage disequilibrium (LD) with HapMap SNPs, (iv) inherited epigenetic factors that are not in strong LD with HapMap SNPs, and (v) interactions among causative alleles that conceal their true contribution to heritability. In this article we investigate the last possibility and determine the extent to which interactions may account for missing heritability.Our analysis is in the same spirit as that of Culverhouse et al. (2002). We assume that the risk of being affected by a complex disease is determined by an individual''s genotype at two or more loci and that the frequencies of causative alleles and the average risks for each one-locus genotype (the marginal risks) are known. Culverhouse et al. (2002) assumed the marginal risks were the same for all genotypes and all loci. In that case, causative alleles have odds ratios of 1; they contribute to risk only through their interactions. Culverhouse et al. found the risk function that maximized the heritability and showed that the maximum possible heritability attributable to interactions increased with the number of loci. They concluded that it is quite possible that interactions among loci that have no main effect could contribute substantially to the heritability of a complex disease and indeed could account for “virtually all the variation in affection status for diseases with any prevalence” (Culverhouse et al. 2002, p. 468).We generalize the analysis of Culverhouse et al. in three ways. First, we allow causative alleles to have odds ratios >1. Second, we explore the entire space of models instead of focusing only on the risk model that maximizes heritability. Third, we examine how the importance of gene interactions depends on the “distance” between a risk model and a comparable multiplicative model. We show that gene interactions can substantially increase the heritability of risk as measured by recurrence risk, K_R, and that the effect increases with the number of loci carrying causative alleles. Furthermore, we show that these results are true even if more weight is given to models that are closer to a comparable multiplicative model.Geometrically, the space of feasible genotype-specific risks subject to the aforementioned constraints (i.e., that the allele frequencies and marginal risks are known) corresponds to a high-dimensional convex polytope, and the computational problem of interest involves integrating a quadratic function over the polytope. The dimension of the polytope grows exponentially with the number of loci, and, therefore, analytic computation is intractable for more than two loci. Hence, we devise a Monte Carlo approach to tackle the problem. Note that, because of high dimensionality, rejection algorithms are not appropriate for this kind of problem. We instead employ a Markov chain Monte Carlo (MCMC) algorithm based on a random walk that always stays inside the polytope. We present empirical results for up to five loci and obtain a closed-form formula for the minimum of K_R over the polytope; the latter result applies to an arbitrary number of loci. Interestingly, the minimum of K_R decreases as the number L of loci increases, but the average of K_R over the polytope increases with L.     

The concordance of gene trees and species trees at two linked loci

The gene genealogies of two linked loci in three species are analyzed using a series of Markov chain models. We calculate the probability that the gene tree of one locus is concordant with the species tree, given that the gene tree of the other locus is concordant. We define a threshold value of the recombination rate, r*, to be the rate for which the difference between the conditional probability of concordance and its asymptotic value is reduced to 5% of the initial difference. We find that, although r* depends in a complicated way on the times of speciation and effective population sizes, it is always relatively small, <10/N4, where N4 is the effective size of the species represented by the internal branch of the species tree. Consequently, the concordance of gene trees of neutral loci with the species tree is expected to be on roughly the same length scale on the chromosome as the extent of significant linkage disequilibrium within species unless the effective size of contemporary populations is very different from the effective sizes of their ancestral populations. Both balancing selection and selective sweeps can result in much longer genomic regions having concordant gene trees.