#
Estimating heterozygosity from a low-coverage genome sequence, leveraging data from other individuals sequenced at the same sites

An Investigation Submitted to Genetics

###### Abstract

High-throughput shotgun sequence data makes it possible in principle to accurately estimate population genetic parameters without confounding by SNP ascertainment bias. One such statistic of interest is the proportion of heterozygous sites within an individual’s genome, which is informative about inbreeding and effective population size. However, in many cases, the available sequence data of an individual is limited to low coverage, preventing the confident calling of genotypes necessary to directly count the proportion of heterozygous sites. Here, we present a method for estimating an individual’s genome-wide rate of heterozygosity from low-coverage sequence data, without an intermediate step calling genotypes. Our method jointly learns the shared allele distribution between the individual and a panel of other individuals, together with the sequencing error distributions and the reference bias. We show our method works well, first by its performance on simulated sequence data, and secondly on real sequence data where we obtain estimates using low coverage data consistent with those from higher coverage. We apply our method to obtain estimates of the rate of heterozygosity for 11 humans from diverse world-wide populations, and through this analysis reveal the complex dependency of local sequencing coverage on the true underlying heterozygosity, which complicates the estimation of heterozygosity from sequence data. We show filters can correct for the confounding by sequencing depth. We find in practice that ratios of heterozygosity are more interpretable than absolute estimates, and show that we obtain excellent conformity of ratios of heterozygosity with previous estimates from higher coverage data.

## 1 Introduction

Heterozygosity, or the fraction of nucleotides within an individual that differ between the chromosomes they inherit from their parents, is a crucial number for understanding human variation. Estimating this simple statistic from any type of sequence data is confounded by sequencing errors, mapping errors, and imperfect power for detecting polymorphisms. Obtaining an unbiased estimate is especially difficult for ancient genomes where the sequences have a higher error rate, or in cases of low-coverage sequence data where there is low power to detect heterozygous sites, or for hybrid capture where there may be additional biases due to the oligonucleotides used for fishing out sequences of interest.

Several methods for estimating individual heterozygosity have been proposed [\citeauthoryearJohnson and SlatkinJohnson and Slatkin2006, \citeauthoryearHellmann, Mang, Gu, Li, Francisco, Clark, and NielsenHellmann et al.2008, \citeauthoryearJiang, Tavaré, and MarjoramJiang et al.2009, \citeauthoryearLynchLynch2008, \citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010]. For an overview of these methods see [\citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010]. Haubold et al. 2010 [\citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010] describe mlRho, an implementation of a method that jointly infers , the scaled mutation rate, and , the scaled recombination rate for a shotgun sequenced genome. However, they examined performance of their method at 10X coverage and a small sequence error rate of , which is about four times lower than encountered currently in real data [\citeauthoryearShendure and JiShendure and Ji2008]. We developed a method that estimates the heterozygosity for an individual of interest by leveraging the genome-wide joint information across sequence reads from a panel of individuals. The advantage of leveraging the panel of individuals in our method is that it enables learning of the empirical distribution of alleles at heterozygous and homozygous positions, a distribution that encapsulates sequencing errors and the non-Bernoulli sampling of each allele at a heterozygous SNP. This allows one to disentangle the rate of heterozygosity from sequencing errors and other biases, and does not require explicit modeling of these platform–, batch–, and genome–specific (frequently unknown and unestimatable) error processes. As a result of including the allele or genotype information at other individuals, our method gains robustness to any unknown error sources that may also be present within the data.

We use an Expectation-Maximization (EM) algorithm to estimate the most likely distribution of counts across the unknown underlying genotypic states, from which we obtain an estimate of the proportion of loci that are heterozygous in the target individual. An advantage of this method is that it returns an unbiased and accurate estimate of heterozygosity even when the individual has low sequence coverage. Our method learns the distribution of alleles directly from the sequence read data, and does not require modeling demographic relationships among the individuals nor genotype calls from the sequence reads. We validate our EM method on 1 Gb of simulated sequence data of 5X, 10X, and 20X coverage, and find that our method performs well at estimating the true heterozygosity even when the sequence error rate is extreme and mean coverage is low. As an empirical validation of the ability of our method to perform well on low-coverage datasets, we test our method on real high-coverage ( 30X) sequencing data, which we subsample to lower coverage, and verify that our estimates are consistent. In particular, we show that applying our method to a lower-coverage subsampling provides the same estimates of heterozygosity as those obtained on higher-coverage data, and are concordant with estimates of heterozygosity from other methods.

We apply our method to obtain estimates of heterozygosity for 11 individuals from many world-wide human populations, from [\citeauthoryearMeyer, Kircher, Gansauge, Li, Racimo, Mallick, Schraiber, Jay, Prüfer, de Filippo, et al.Meyer et al.2012]. Our finding underscores the need to compare ratios of heterozygosity across fixed genomic regions to infer the relative rates of diversity among individuals.

## 2 Materials and Methods

We apply our method to read data at sites with a target minimum coverage (for example, 5X coverage) for the sequenced diploid individual of interest, aligned to some reference genome of known sequence. We also use sequence read data from other individuals likewise aligned to the reference.

Let be the unknown diploid genotype of our target individual at some position in the genome, and be the aligned reference allele. Then the allele distribution in other individuals will depend on . Let be the vector of alleles by taking one randomly sampled read from each of the individuals. Let be the observed alleles from the reads for our individual. Both and are observed quantities for a given position in the genome for our individual, and we are interested in modeling the joint probability of , as the sum of the joint probabilities conditional on .

We assume conditional independence of on the true unobserved genotype . This assumption holds if the allele frequency spectrum of the panel of individuals depends only on the true underlying genotypic state of our individual, and not the allele counts we observe, and likewise the allele count distribution depends only on the true underlying genotypic state and not on the the alleles observed in the other individuals. From this conditional independence property, we then derive:

(1) | |||||

(2) | |||||

(3) | |||||

(4) |

which will later provide the leverage to infer the most likely values for the above probabilities, including , which gives us the genomic rate of heterozygosity.

For every site which has both sufficient coverage in our individual and for which we have complete information of the panel, we add this site to the corresponding bin of observed alleles and panel . This full matrix would be inconveniently large, so to simplify the data matrix of counts, we polarize our allele counts with respect to the reference, restricting to bi-allelic SNPs, which constitute the majority of sites. We denote the reference allele as , and allow only a single other variant per site, summarizing the observed alleles from the reads by the number of non-reference alleles. Thus, we denote genotypes as , which we will refer to as the homozygous ancestral, heterozygous, and homozygous derived states, respectively. If, for example, we consider only sites with a coverage of 4, then . We can also easily represent as a vector of 0’s and 1’s referring to reference or non-reference allele present in the randomly sampled reads, for example, where the length of is determined by the number of individuals we sample.

We create a count matrix of dimension corresponding to the number of observed sites with each particular combination of and . The counts of the number of loci where the individual is and the panel of individuals individuals comprise the alleles is represented by the corresponding row and column entry in the matrix .

From the matrix we estimate the true values of , and using the EM algorithm. Let be the observed counts of alleles in the matrix . Let be , the missing or unobserved counts of the alleles with the true parameter state . Then the likelihood of the data is:

(5) |

If the hidden variable, , corresponding to the true underlying genotypic state were observed, the log-likelihood would be:

(6) |

But this would require fitting different parameters per observed data point (i.e., count entry of ). This would require fitting three times as many parameters as there are data points. However, by relying on our conditional independence from equation (4) above we can reduce the number of parameters to be fitted from the data.

By EM theory, the function is given by:

where is the expected value of , which in our case derives from the multinomial distribution, under the posterior distribution calculated with the old parameters .

The estimates for that maximize , also derived from the MLE estimates for the multinomial distribution, are:

Further, by Bayes theorem this expands to,

By basic EM theory these re-estimated values of will generate a non-decreasing sequence of values for the log likelihood . Finally, we obtain the parameter of interest after convergence.

### 2.1 Implementation

In practice, without constraining the parameters we reach local but not consistently global likelihood maxima, which do not necessarily correspond to the genotypic state parameters we wish to obtain. To improve the ability of the EM to achieve global maxima, we fit Beta-Binomial distributions (effectively an over-dispersed Binomial distribution) to the probabilities of number of non-reference alleles for each possible genotypic state . This constraint, as well as the choice of reasonable starting parameters for the EM initialization, in practice improves our convergence to global maximum corresponding to the homozygous ancestral, heterozygous and homozygous derived genotypic states.

Like the Beta distribution, the MLE estimates for the Beta-Binomial distribution do not have a closed form, though they can be found using direct numerical optimization [such as a fixed-point iteration or a Newton-Raphson iteration]. However, instead, we estimate the two parameters () using MOM estimators for the Beta-Binomial, by setting:

In the case of under-dispersed data, it is possible to obtain MOM estimates that are invalid. Though unlikely to occur in the read data, for this contingency, we instead fit a Binomial distribution to the data.

There are several challenges in implementing the EM for our problem. The first is that the likelihood is poorly defined when any of the parameters we are interested in estimating approach 0, as then the likelihood also goes to zero. So to avoid this situation, we add a “prior” to the likelihood calculation which adds a small count value in the step calculating the parameters to avoid probabilities reaching 0. So instead, we calculate the posterior:

rather than the maximum likelihood estimation, hence, we obtain a MAP (maximum a posteriori) estimate, which is a Bayesian method that incorporates a prior over the distribution to be estimated (in this case, a small uniform prior). We choose a small prior (less than in total counting one site across all possible matrices) that does not impact our estimates. In general, our estimates are robust to choice of this prior, within a range examined of 1e-10 to 1e-50, and we will continue to refer to our method as an EM implementation though in fact we use a non MLE method. In effect, our equations for each step remain the same, except that in the M-step of the EM (where we estimate the parameters) we instead estimate the MAP using the prior. Specifically, we estimate:

In practice, we set to 1e–20 which does not alter estimates of the probabilities while preventing ill-defined likelihoods.

As with all likelihood calculations, our probabilities approach very small numbers. To avoid numerical error due to underflow of small likelihoods and parameter estimates, we implement the EM storing all probabilities and likelihoods in the log form.

Lastly, likelihood maximization occurs on an arbitrary base, so to avoid numerical issues due to any remaining underflow of the likelihood calculation, we compute a factor at the start of the EM. For each iteration, we compute the likelihood of the data minus this constant factor, which is a standard practice and does not affect the computation of the maximum. This is equivalent to calculating the log odds:

for some constants . In practice, we set to be the likelihood at initialization of the EM. We then iterate the EM until both the change in parameters and the change in the likelihood is smaller than our chosen threshold, which in practice we set as 1e–50.

It should be noted that any form for tallying read counts may be used, including the allele profile used in [\citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010]; our choice of was motivated by a choice of dimensionality that is a compromise between simplicity and capturing relevant information. Our method is highly generalizable to any choice of count data, and could be implemented assuming that a reasonable starting position for the EM could be proposed, such that the iterations are likely to converge to a global maximum corresponding to the genotypic states.

### 2.2 Proof of principle 1: Application to simulated data

We generated simulated sequence data and applied our EM method for estimating heterozygosity to assess the accuracy of our estimation procedure.

#### 2.2.1 Generating coalescent simulations

We generated 100 replicate datasets of sequence data using MaCS [\citeauthoryearChen, Marjoram, and WallChen et al.2009]. Each replicate dataset contains 10 chromosomes of length 100Mb for a total of 1Gb of sequence, for 1 chimpanzee chromosome, 7 African chromosomes, 5 European chromosomes, and 5 East Asian chromosomes, using demographic parameters fitted by [\citeauthoryearGutenkunst, Hernandez, Williamson, and BustamanteGutenkunst et al.2009]. We include a chimpanzee chromosome assuming a constant population size of 50,000 individuals and a split time from humans of 6Mya, using the same generation time as humans of 25 years per generation.

#### 2.2.2 Adding simulated error

We simulate sequence data from the true genotypes by adding errors to reads. First, for all variable loci in the target individual, we randomly choose which allele is on a read, then adds errors to the each read (with high error rate of 0.002) to generate the total number of derived reads for the individual at the locus out of the total sequencing depth. For each other sequenced chromosome, we add errors with a lower error rate of 0.0001 (since we assume the panel is composed of higher-quality genomes), then add a count for the final simulated locus in the appropriate hash bin. Lastly, for each invariant locus, we add errors to the target individuals reads, and at a lower rate, adds errors to the other sequenced chromosomes, and input these counts into the hash bin. With an error rate of 0.001 we add errors to the chimpanzee chromosome which inverts the ancestral and derived reads.

### 2.3 Proof of principle 2: Downsampling high-coverage genomes

To assess the efficacy of our method at lower coverages on real sequence data, we obtain estimates of heterozygosity for a San individual sequenced to higher coverage using Illumina’s Genome Analyzer IIx next-generation sequencing technology, which we then downsample to varying levels low coverage. We use this dataset of sequence reads to explore the ability of our method to perform on low-coverage sequence data, and the lower bound of coverage at which we are able to obtain accurate estimates of heterozygosity. We compare the performance of our method to the estimate of obtained from MlRho [\citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010].

### 2.4 Application to 11 world-wide human genomes

We align sequence data from 11 human genomes from world-wide populations and an archaic Denisovan genome to the chimpanzee reference genome to avoid introducing human-reference-population biases. For details on the populations, samples and the sequencing performed, see [\citeauthoryearMeyer, Kircher, Gansauge, Li, Racimo, Mallick, Schraiber, Jay, Prüfer, de Filippo, et al.Meyer et al.2012]. We generate a counts matrix for each of the genomes using a panel generated by a single read sampled from each of the other 11 genomes.

We include only sites where there is a chimpanzee reference allele, exclude sites where two or more non-reference bases are equally present or if there are more than 5 reads showing a third (non-variant and non-reference) allele. We also exclude CpG sites, as well as sites where any individual from the panel has no coverage or sites that have insufficient coverage for the target individual.

To demonstrate the relationship between sequencing coverage and the true rate of heterozygosity of different regions, we generate count data for each bin of 5X coverage ranging between 5X and 50X, where for each bin dataset, we only include sites where the coverage of the individual falls within the target range. We downsample coverage at each bin (where possible) to 5X, 10X, and 20X, and compare results stratified by downsampling, as well as by genomic coverage.

Lastly, we produce estimates for the 11 present-day genomes and the archaic Denisova genome on a fixed set of sites, and compare to previous estimates for these samples [\citeauthoryearMeyer, Kircher, Gansauge, Li, Racimo, Mallick, Schraiber, Jay, Prüfer, de Filippo, et al.Meyer et al.2012].

## 3 Results

### Simulation results

We obtain accurate estimates of heterozygosity across a variety of coverage levels (5X to 30X) (see Figure 1). We note a tiny bias of 0.3% (in relative terms) from the true rate for 5X coverage read data, but with higher coverage this bias goes to zero.

### Downsampling results

Figure 2 illustrates that our EM estimation method and MlRho give consistent estimates of heterozygosity for the HGDP San individual starting at about 10X coverage and higher. However, at lower coverage (about 4X) our method significantly outperforms MlRho, giving a slightly biased, but nearly convergent, estimate.

### Heterozygosity estimates for 11 present-day and Denisovan genomes

We present our initial estimates of heterozygosity, downsampled to three different depths, for each sequencing coverage bin (normalized by individual mean sequence coverage) in Figure 3A. Our estimates of heterozygosity are consistent, independent of count matrix (i.e., downsampling) size, as would be expected from our simulated downsampling results shown in Figure 2. However, we find a strong signal that the estimates of heterozygosity are correlated to sequencing coverage of the region. We note that this is not an artifact of the larger amount of data available at higher coverage, since each bin is calculated after being downsampled to the same depth. Instead, the U-shaped curves in Figure 3B indicate that the apparent next-generation sequencing coverage is dependent on properties of the underlying genomic sequence. In particular, we find that regions of lower coverage and higher coverage (relative to the mean sequencing depth) show higher heterozygosity.

We witnessed increased heterozygosity at regions of higher coverage, which we suspected was due to perceived genetic diversity due to cryptic segmental duplications. To explore this hypothesis, we restricted our analyses to regions of the genome which have been identified as unlikely to contain segmental duplications, available on the Eichler Laboratory website (http://eichlerlab.gs.washington.edu/database.html). We find that this filter strongly reduces the effect of regions with likely segmental duplications on our estimates of heterozygosity (Figure 3C), confirming that unidentified segmental duplications, which result in a net higher sequencing coverage of the region, result in an high estimate of heterozygosity for such regions. Roving these regions with known segmental duplications reduces this effect at regions with higher sequencing coverage. However, the increase in heterozygosity at higher coverage still is present even after this correction (see Figure 3C), suggesting that there may be further individual or population-specific duplications remaining.

Using only data that passed the segmental duplication filter, we obtain estimates for the sequenced genomes on the same set of regions, restricting to regions with sequencing coverage between 20X and 40X. Using the EM, we estimate the total genome-wide fraction of heterozygosity for each individual, and we also can extract estimates of the allelic distribution of heterozygous and homozygous sites (Figure 4). We present the absolute estimates we obtain in Table 1, as well as the ratio of heterozygosity in the Denisova genome relative to the other individuals. We find the highest estimates of heterozygosity for the San African individual, and next highest estimates of heterozygosity for other African individuals from Mandenka, Yoruba, Mbuti, and Dinka populations. The next highest levels of heterozygosity are in individuals from European populations (French, Sardinian), followed by East Asian populations (Dai, Han). We find the lowest estimates of heterozygosity in the individuals from Melanesia (Papuan) and from a Native American population (Karitiana).

Individual | Heterozygosity estimate (%) | Ratio |
---|---|---|

Denisova | 0.0165 | – |

San | 0.0721 | 23% |

Mandenka | 0.0686 | 24% |

Yoruba | 0.0649 | 25% |

Mbuti | 0.0657 | 25% |

Dinka | 0.0635 | 26% |

Sardinian | 0.0490 | 34% |

French | 0.0473 | 35% |

Dai | 0.0465 | 35% |

Han | 0.0454 | 36% |

Papuan | 0.0386 | 43% |

Karitiana | 0.0353 | 47% |

## 4 Discussion

We have shown that our heterozygosity estimation method both performs well in low coverage simulated sequence data and provides consistent estimates on real low-coverage data downsampled from higher coverage. In particular, our method outperforms other methods on data that has been sequenced at less than 10X coverage, and provides reasonable estimates for as low as 4X coverage.

Our estimates for 11 world-wide human genomes and the archaic Denisovan genome provide important insights into the distribution of heterozygosity across human populations. Furthermore, our results show that estimates of heterozygosity are strongly affected by genomic properties such as copy-number variability, and these properties affect sequencing coverage. Hence, we show that the heterozygosity is not independent of sequencing coverage even within one genome, and is elevated in both regions with low coverage (relative to the mean sequencing depth) as well as regions with high coverage. This is an unexpected result if one assumes a the “Lander-Waterman” Poisson distribution of read depth [\citeauthoryearLander and WatermanLander and Waterman1988, \citeauthoryearWeber and MyersWeber and Myers1997]. Furthermore, even after excluding regions with known copy-number variable regions, an increase in heterozygosity is still present at the more extreme levels of sequence coverage, suggesting that other correlations of sequence diversity with coverage, or possibly individual-specific segmental duplications, still remain. Implications from our result suggest that using the higher tail of sequencing coverage for population genetic inference may result in a biased set of genomic regions with selectively higher heterozygosity, possibly due to population and individual segmental duplications.

Our absolute estimates of heterozygosity are lower than those reported for these genomes in other papers using other methods [\citeauthoryearMeyer, Kircher, Gansauge, Li, Racimo, Mallick, Schraiber, Jay, Prüfer, de Filippo, et al.Meyer et al.2012]. The absolute estimates are notably lower, but consistent with the relative diversity previously documented in these populations, and with other patterns of genetic diversity such as decay of linkage disequilibirium [\citeauthoryearJakobsson, Scholz, Scheet, Gibbs, VanLiere, Fung, Szpiech, Degnan, Wang, Guerreiro, Bras, Schymick, Hernandez, Traynor, Simon-Sanchez, Matarin, Britton, van de Leemput, Rafferty, Bucan, Cann, Hardy, Rosenberg, and SingletonJakobsson et al.2008]. Because the regions of the genome that pass our filters (and in particular, the copy-number variable filter) are likely to be lower in complexity and substantially biased towards lower diversity due to alignment biases, the lower absolute values of heterozygosity are expected. However, our relative heterozygosity estimates are consistent with previously documented levels of genetic diversity, with African populations showing highest levels of diversity, and decreasing levels with distance away from Africa [\citeauthoryearLi, Absher, Tang, Southwick, Casto, Ramachandran, Cann, Barsh, Feldman, Cavalli-Sforza, and MyersLi et al.2008, \citeauthoryearJakobsson, Scholz, Scheet, Gibbs, VanLiere, Fung, Szpiech, Degnan, Wang, Guerreiro, Bras, Schymick, Hernandez, Traynor, Simon-Sanchez, Matarin, Britton, van de Leemput, Rafferty, Bucan, Cann, Hardy, Rosenberg, and SingletonJakobsson et al.2008]. Our estimates confirm previous findings that the archaic Denisovan genome shows substantially lower levels of heterozygosity than any of the other present-day populations, with only a fraction of the rate of heterozygosity.

More generally, we emphasize that absolute heterozygosity is not a well-defined quantity in the analysis of genomic data, as it strongly depends on the particular filters that are used to select the regions being analyzed, and may be an implausible concept in highly repetitive regions (such as centromeres and telomeres) and copy-number-variable regions. The absolute value of heterozygosity can vary based on the regions chosen to be examined, but the relative heterozygosity estimates or ratios among individuals (using the same regions and filters) are consistent. Hence, in practice, heterozygosity estimates are most meaningful when viewed as relative ratios among individuals for the same regions of the genome, and not as absolute values inherent to diploid genomes.

## 5 Acknowledgments

KB gratefully acknowledges that this investigation was support by the National Institutes of Health under Ruth L. Kirschstein National Research Service Award #5F32HG006411.

## References

- [\citeauthoryearChen, Marjoram, and WallChen et al.2009] Chen, G., P. Marjoram, and J. Wall, 2009 Fast and flexible simulation of DNA sequence data. Genome research 19(1): 136–142.
- [\citeauthoryearGutenkunst, Hernandez, Williamson, and BustamanteGutenkunst et al.2009] Gutenkunst, R., R. Hernandez, S. Williamson, and C. Bustamante, 2009 Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS genetics 5(10): e1000695.
- [\citeauthoryearHaubold, Pfaffelhuber, and LynchHaubold et al.2010] Haubold, B., P. Pfaffelhuber, and M. Lynch, 2010 mlRho–a program for estimating the population mutation and recombination rates from shotgun-sequenced diploid genomes. Molecular Ecology 19(s1): 277–284.
- [\citeauthoryearHellmann, Mang, Gu, Li, Francisco, Clark, and NielsenHellmann et al.2008] Hellmann, I., Y. Mang, Z. Gu, P. Li, M. Francisco, A. Clark, and R. Nielsen, 2008 Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals. Genome research 18(7): 1020–1029.
- [\citeauthoryearJakobsson, Scholz, Scheet, Gibbs, VanLiere, Fung, Szpiech, Degnan, Wang, Guerreiro, Bras, Schymick, Hernandez, Traynor, Simon-Sanchez, Matarin, Britton, van de Leemput, Rafferty, Bucan, Cann, Hardy, Rosenberg, and SingletonJakobsson et al.2008] Jakobsson, M., S. W. Scholz, P. Scheet, J. R. Gibbs, J. M. VanLiere, H. C. Fung, Z. A. Szpiech, J. H. Degnan, K. Wang, R. Guerreiro, J. M. Bras, J. C. Schymick, D. G. Hernandez, B. J. Traynor, J. Simon-Sanchez, M. Matarin, A. Britton, J. van de Leemput, I. Rafferty, M. Bucan, H. M. Cann, J. A. Hardy, N. A. Rosenberg, and A. B. Singleton, 2008 Genotype, haplotype and copy-number variation in worldwide human populations. Nature 451(7181): 998–1003.
- [\citeauthoryearJiang, Tavaré, and MarjoramJiang et al.2009] Jiang, R., S. Tavaré, and P. Marjoram, 2009 Population genetic inference from resequencing data. Genetics 181(1): 187–197.
- [\citeauthoryearJohnson and SlatkinJohnson and Slatkin2006] Johnson, P. and M. Slatkin, 2006 Inference of population genetic parameters in metagenomics: a clean look at messy data. Genome research 16(10): 1320–1327.
- [\citeauthoryearLander and WatermanLander and Waterman1988] Lander, E. and M. Waterman, 1988 Genomic mapping by fingerprinting random clones: a mathematical analysis. Genomics 2(3): 231–239.
- [\citeauthoryearLi, Absher, Tang, Southwick, Casto, Ramachandran, Cann, Barsh, Feldman, Cavalli-Sforza, and MyersLi et al.2008] Li, J. Z., D. M. Absher, H. Tang, A. M. Southwick, A. M. Casto, S. Ramachandran, H. M. Cann, G. S. Barsh, M. Feldman, L. L. Cavalli-Sforza, and R. M. Myers, 2008 Worldwide human relationships inferred from genome-wide patterns of variation. Science 319(5866): 1100–4.
- [\citeauthoryearLynchLynch2008] Lynch, M., 2008 Estimation of nucleotide diversity, disequilibrium coefficients, and mutation rates from high-coverage genome-sequencing projects. Molecular biology and evolution 25(11): 2409–2419.
- [\citeauthoryearMeyer, Kircher, Gansauge, Li, Racimo, Mallick, Schraiber, Jay, Prüfer, de Filippo, et al.Meyer et al.2012] Meyer, M., M. Kircher, M. Gansauge, H. Li, F. Racimo, S. Mallick, J. Schraiber, F. Jay, K. Prüfer, C. de Filippo, and others, 2012 A high-coverage genome sequence from an archaic denisovan individual. Science.
- [\citeauthoryearShendure and JiShendure and Ji2008] Shendure, J. and H. Ji, 2008 Next-generation DNA sequencing. Nature biotechnology 26(10): 1135–1145.
- [\citeauthoryearWeber and MyersWeber and Myers1997] Weber, J. and E. Myers, 1997 Human whole-genome shotgun sequencing. Genome Research 7(5): 401–409.