Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Estimating genetic kin relationships in prehistoric populations

  • Jose Manuel Monroy Kuhn,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

    Current address: University of Freiburg, Institute of Biology I (Zoology), Hauptstrasse 1, D-79104 Freiburg, Germany

    Affiliation Uppsala University, Evolutionary Biology Centre, Department of Organismal Biology, Norbyvägen 18C, SE-752 36 Uppsala, Sweden

  • Mattias Jakobsson ,

    Roles Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing – review & editing

    mattias.jakobsson@ebc.uu.se (MJ); torsten.guenther@ebc.uu.se (TG)

    Affiliations Uppsala University, Evolutionary Biology Centre, Department of Organismal Biology, Norbyvägen 18C, SE-752 36 Uppsala, Sweden, Uppsala University, SciLifeLab, Norbyvägen 18C, SE-752 36 Uppsala, Sweden

  • Torsten Günther

    Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

    mattias.jakobsson@ebc.uu.se (MJ); torsten.guenther@ebc.uu.se (TG)

    Affiliation Uppsala University, Evolutionary Biology Centre, Department of Organismal Biology, Norbyvägen 18C, SE-752 36 Uppsala, Sweden

Abstract

Archaeogenomic research has proven to be a valuable tool to trace migrations of historic and prehistoric individuals and groups, whereas relationships within a group or burial site have not been investigated to a large extent. Knowing the genetic kinship of historic and prehistoric individuals would give important insights into social structures of ancient and historic cultures. Most archaeogenetic research concerning kinship has been restricted to uniparental markers, while studies using genome-wide information were mainly focused on comparisons between populations. Applications which infer the degree of relationship based on modern-day DNA information typically require diploid genotype data. Low concentration of endogenous DNA, fragmentation and other post-mortem damage to ancient DNA (aDNA) makes the application of such tools unfeasible for most archaeological samples. To infer family relationships for degraded samples, we developed the software READ (Relationship Estimation from Ancient DNA). We show that our heuristic approach can successfully infer up to second degree relationships with as little as 0.1x shotgun coverage per genome for pairs of individuals. We uncover previously unknown relationships among prehistoric individuals by applying READ to published aDNA data from several human remains excavated from different cultural contexts. In particular, we find a group of five closely related males from the same Corded Ware culture site in modern-day Germany, suggesting patrilocality, which highlights the possibility to uncover social structures of ancient populations by applying READ to genome-wide aDNA data. READ is publicly available from https://bitbucket.org/tguenther/read.

Introduction

An individual’s genome is a mosaic of different segments inherited from our various direct ancestors. These segments, shared between individuals, can be referred to as identical by descent (IBD). Knowledge about IBD segments has been used for haplotype phasing [1, 2], heritability estimation [3, 4], population history [5], inference of natural selection [6] and to estimate the degree of biological relationship among individuals [7]. A number of methods have been developed to estimate the degree of biological relationship by inferring IBD from SNP genotype or whole genome sequencing data. The methods for estimating relationship levels implemented in PLINK [8], SNPduo [9], ERSA [10, 11], KING [12], REAP [13] and GRAB [14] greatly benefit from genome wide diploid data, information about phase, recombination maps and population allele frequency, and are sometimes able to successfully infer relationships up to 11th degree [11].

Knowing whether a pair of individuals is directly related or not, and estimating the degree of relationship is of interest in various fields: Genome-wide association studies and population genetic analyses often try to exclude related individuals since they do not represent statistically independent samples; in forensics, archaeology and genealogy, individuals and their relatives can be identified based on DNA extracted from human remains [15, 16]; Breeders and conservation biologists are interested in the relatedness of mating individuals [17, 18]. Current methods present significant limitations for the analysis of degraded samples as they rely on diploid genotype calls, low proportions of missing data and sometimes even phase information. Especially in the fields of forensics and archaeology, the amount of endogenous DNA available for analysis is limited due to postmortem degradation [1921]. In archaeology, the analysis of IBD has the potential to provide an independent means to test kinship behavior and social organization [22, 23], but current methods would be restricted to exceptionally well-preserved samples. In forensic science and practice, the dominant approach has been to type several short tandem repeat (STR) markers, which in most cases provide sufficient information for relatedness assessment, but the STRs might be hard to type in degraded samples [24]. In addition to nuclear STRs, mitochondrial and Y-chromosome haplogroups have been widely used to infer family relationships (e.g. [15, 16, 25, 26]), although they can only exclude certain direct relationships since most mitochondrial and Y-chromosome haplogroups are relatively common among unrelated individuals. These uniparental markers can be typed from degraded samples, and can be used to exclude maternal or paternal relationships, but not to infer the actual degree of relationship. Genome-wide data, however, can be obtained from degraded samples at a higher success rate than STRs and it can be used to confidently identify individuals [27].

Single Nucleotide Polymorphism (SNP) data can be obtained from genotyping experiments (e.g. SNP arrays or RAD sequencing), targeted capture [28], and whole-genome shotgun sequencing (e.g. [29, 30]). The field of ancient DNA has developed rapidly over the last few years and allowed pivotal studies of the population history of Europe [23, 2842] and the peopling of the Americas [40, 43, 44]. However, both whole-genome shotgun sequencing (e.g. [30, 33, 34]) and genome-wide SNP capture (e.g. [28, 35]) usually achieve coverages <1x per informative site for most individuals which makes diploid genotype calls at all sites virtually impossible. Methods to infer relationships, however, rely on such ideal data to identify IBD blocks which is a major limitation for applying these methods to ancient DNA data.

However, even low coverage data contain information about the degree of relationship. To utilize this information, we developed READ (Relationship Estimation from Ancient DNA), a heuristic method to infer family relationships up to second degree from samples with extremely low coverage. The method is tested on publicly available data with known relationship, which we sub-sample to resemble the properties of degraded samples. We also apply our method to a number of ancient samples from the literature and confidently classify individual pairs as being related.

Results

Method outline

The input for READ are a set of TPED/TFAM files [8] containing pseudo-haploid genotypes for a population. The biallelic SNP sites included in that file would usually be from some externally ascertained SNP panel (e.g. Human Origins array or 1000 genomes). The data is assumed to be pseudo-haploid (i.e. one randomly sampled allele per individual and SNP site) as the low coverage in aDNA studies normally does not allow to call heterozygous genotypes. This procedure of randomly sampling one sequencing read per SNP site is widely used in aDNA studies of low coverage data, see e.g. [23, 2831, 3336, 39, 42, 45]. We divide the genome into non-overlapping windows of 1 Mbps each and for each pair of individuals calculate the proportion of non-matching alleles inside each window (P0). Before classifying the degree of relationship of a pair of individuals, we need to normalize P0 using the expected value for a randomly chosen pair of unrelated individuals from the same population in order to make the classification independent of within population diversity, SNP ascertainment and marker density. In most applications, that expected value is difficult to infer which is why several proxies can be used: a pair of unrelated individuals from the same population (similar to [46, 47]), a pair of individuals from a different population with similar expected diversity, or the median of all average pairwise P0 across all individuals which should correspond to a pair of unrelated individuals if the sample size is sufficient. The latter setting is the default option for READ and we are using it in all major simulations as well as the empirical data analysis of this study. Depending on the normalized proportion of shared alleles, each pair of individuals is classified as unrelated, second-degree (i.e. nephew/niece-uncle/aunt, grandparent-grandchild or half-siblings), first-degree (parent-offspring or siblings) or identical individuals/identical twins (Fig 1). As a method with the goal to classify pairs of individuals, READ always outputs the best fitting degree of relationship. This decision is based on the point estimate of the average P0 and we observe throughout our simulations that basing classifications on the point estimates has a low number of false positives. The user is provided with a graphical summary of the classification results which also includes the uncertainties of the different estimates. To additionally express the certainty of each categorization, the distance to the classification cutoffs are expressed as multiples of the standard error of the mean (Z).

thumbnail
Fig 1. Outline of the general READ workflow to estimate the degree of relationship between two individuals.

https://doi.org/10.1371/journal.pone.0195491.g001

Simulations based on modern data with known relationship

READ’s performance was tested on 1,326 individuals of 15 different populations from the phase 3 data of the 1000 genomes project [48]. A total of 86,336 pairwise comparisons were tested. The rates of false positives (unrelated individuals classified as related) and false negatives (related individuals classified as unrelated or as wrong degree) are highly dependent on the amount of data available for pairwise comparison. We measure the amount of available data as the number of SNP sites (out of a maximum of 1,156,468) with allelic information for both individuals. READ showed an overall good performance with false positive rates below three percent for as little as 1,000 overlapping SNP loci (Fig 2A). False negative rates are ≤ 10 percent across all tested pairs and are highest with the lowest amount of overlapping SNPs between the two individuals. Separating the error rates between first and second degree relatives shows that the false negative rate is much higher for pairs know to be second degree relatives who tend to be classified as unrelated (Fig 2C) while first degree relatives tend to be classified as second degree relatives if the classification is incorrect (Fig 2B). False positive rates are low for both degrees of relationship (Fig 2B and 2C). The rate of false negatives increases up to 7.5% for first degree relationships and 38% for second degree relationships when the number of SNPs is low (Fig 2B and 2C).

thumbnail
Fig 2. Simulation based estimates of false positive and false negative rates for different numbers of SNPs.

The analysis is based on pairs of individuals with known degree of relationship in the 1000 genomes data. (A) All degrees of relationship, (B) only first degree relatives and (C) only second degree relatives. For pairs known to be related but not classified correctly (“False negative”) we distinguish between pairs classified as unrelated and classified as related but to a wrong degree. Error rates were estimated for 1,000, 2,500, 5,000, 10,000 and 50,000 overlapping SNPs between the pair of individuals.

https://doi.org/10.1371/journal.pone.0195491.g002

Further complications in the analysis of empirical aDNA data are sequencing and mapping errors, contamination and post-mortem damage. Ultimately, these issues will increase the proportion of wrongly called alleles at the analyzed SNP sites, which means that READ would analyze a false allele instead of one of the two alleles actually carried by the individual. To see the effect of such allelic errors, we repeated the simulations with certain error rates meaning that alleles were randomly changed with a probability corresponding to a defined site specific allelic error rate. The results of this simulation are shown in Fig 3. Essentially, wrongly called alleles lead to an overestimation of genetic distance between individuals. As a consequence, pairs of individuals tend to get classified into more distant categories which can be seen by an increase in the false negative rate for higher rates of allelic error. False positive rates are not affected by wrongly called alleles but false negative rates increase substantially with more errors. In order to qualitatively investigate a situation where the normalization value is based on a data set with a different error rate than the data used for classification, we performed an additional set of simulations: The two populations IBS and YRI (the populations with the highest number of reported relatives) were split in two halves—one half with a simulated allelic error of 5%, the other with a simulated allelic error of 10%. For each half, a separate normalization value was estimated (the median across all pairs) which was then used for the normalization step when classifying related pairs in the other half of the data. A normalization value based on a lower allelic error rate resulted in an elevated false negative rate while a normalization value based on higher allelic errors caused an inflated false positive rate (S4 Table). These results highlight how important it is to keep the effects of contamination, post-mortem damage and other error sources low in empirical studies. Careful data curation as well as filtering should be able to minimize the rate of allelic errors, making error rates <5% realistic for most applications. We describe some important steps for preparing input data in the Discussion.

thumbnail
Fig 3. Effect of allelic errors on READ’s performance.

The simulations are identical to those conducted for Fig 2 but including a certain proportion of wrongly called alleles. The rates of false positives and false negatives were calculated accordingly. Error rates were estimated for 1,000, 2,500, 5,000, 10,000 and 50,000 overlapping SNPs between the pair of individuals.

https://doi.org/10.1371/journal.pone.0195491.g003

To illustrate how much sequencing would be needed to achieve the required SNP numbers, we estimate the expected number of SNPs covered by at least one read in both individuals depending on the sequencing coverage for each sample (Fig 4). This example assumes that the total number of SNP sites in the data set is 1,156,468 (as in the simulations above and similar to the 1.2 million SNP sites in the empirical aDNA data set studied below [35]) and that the read depth at each SNP locus follows a Poisson distribution with mean corresponding to the genome-wide average sequencing coverage [49].

thumbnail
Fig 4. Number of SNP sites covered in both individuals dependent on the sequencing coverage for each individual.

This figure shows expected number of SNP sites with overlapping data for two individuals for different combinations of sequencing depths. The contour lines mark different numbers of SNPs including those used in the simulations (see Fig 2). The total number of SNPs is set to 1,156,468, identical to what has been used in the simulations and similar to the 1.2 million SNPs used in the empirical data set [35]. The calculation assumes a Poisson distribution of sequencing coverage across the genome [49] and that coverage at each SNP site and individual is independent.

https://doi.org/10.1371/journal.pone.0195491.g004

Relationships among prehistoric Eurasians

To investigate READ’s performance on empirical aDNA data, we analyzed a large published genotype data set of 230 ancient Eurasians from the Mesolithic, Neolithic and Bronze Age periods [35]. In accordance with the original publications [28, 30, 35], READ inferred RISE507 and RISE508 to be the same individual and all nine known relationships were correctly identified as first degree relatives (Table 1). In addition to those, READ identified one additional pair of first degree relatives as well as six new second degree relationships. All relatives are from the same location and their radiocarbon dates (if available) are overlapping.

thumbnail
Table 1. Pairs of relatives among the 230 individuals in the aDNA dataset as inferred by READ.

https://doi.org/10.1371/journal.pone.0195491.t001

READ identified an unknown pair of first degree relationship between two Srubnaya individuals (I0360 and I0354). Notably, Mathieson et al (2015) [35] have excluded I0354 since she was an outlier compared to other Srubnaya individuals. The classification of I0360 and I0354 as first degree relatives is probably genuine considering that READ has very low false positive rates. Fig 5 shows the results for all Srubnaya individuals and these two individuals clearly fall into the group of first degree relatives even when considering uncertainties in the P0 for this pair and for the normalization. If this prediction was a false positive, it would be very likely that they are at least second degree relatives as the fraction of unrelated individuals wrongly classified as first degrees is extremely low (Fig 2B). Furthermore, a highly distinct genetic background of one of the individuals should rather cause false negatives and not false positives, which increases the likelihood that the two individuals are in fact related. I0354 could have been a recent migrant to the region who produced offspring (I0360) with a local male, which would explain both the relationship between I0354 and I0360 and the genomic dissimilarity between I0354 and other Srubnaya individuals.

thumbnail
Fig 5. READ estimates for the Srubnaya sample.

(A) Sorted non-normalized average P0 values for all pairwise comparisons between Srubnaya individuals. Error bars show two standard errors of the mean. The solid horizontal line indicates the median value used for normalization. Dashed lines show the cutoffs used to classify the related individuals. The gray areas around dashed lines indicate 95% confidence intervals for the cutoffs accounting for the uncertainty in estimating the average P0 in the pair used to set the baseline for unrelated individuals. (B) A histogram of the non-normalized average P0 values, vertical dashed and solid lines indicate the same as in (A). A similar plot is produced as output when running READ.

https://doi.org/10.1371/journal.pone.0195491.g005

Particularly interesting is a group of five related males from the Corded Ware site in Esperstedt, Germany (Table 1, Fig 6). Mathieson et al (2015) [35] described two first degree relationships between I1540 and I1541 as well as between I1541 and I1538. Notably, READ missed the second degree relationship between I1540 and I1538, which is likely to be a false negative as the false negative rate for second degree relatives is known to be substantial with low amounts of data (Fig 2C) and the value for that pair (0.91) is only 0.097 standard errors above the threshold for second degree relatives (0.90625). Identical radiocarbon dates do not help to indicate a chronological order, but based on their Y-chromosomes (all likely R1a, S1 Table), one can suggest that they all represent a paternal line of ancestry. I1540 is classified as R1a1, but the Y-chromosomal marker this call is based on (L120) is missing in individuals I1538 and I1541, so they could all carry the same haplotype. In addition to these three individuals, I1534 is a second degree relative of I1538 and I1541, who was carrier of R(xR1b) but a more detailed classification was not possible due to the low coverage. I0104, who is a second degree relative to I1541, might also carry the same Y-chromosome as I1534, I1538, I1540 and I1541, but that cannot be determined due to low coverage in those individuals. Generally, the data would be consistent with all five individuals carrying the same Y-haplotype as there are no contradicting calls for R1a defining markers (S1 Table), which would suggest paternal relationship among them. In total, 13 Corded Ware individuals from Esperstedt were analyzed, nine of them were males. It is notable that all five related Esperstedt individuals discussed here were males and only one pair of related Corded Ware individuals from Esperstedt involved a female (I1539 and I1532; Table 1).

thumbnail
Fig 6. Kin-relationship among males at the Corded Ware site in Esperstedt, Germany.

The five individuals, their inferred degree of relationship and their uniparental haplogroups. The dashed line between I1540 and I1538 shows a second degree relationship missed by READ.

https://doi.org/10.1371/journal.pone.0195491.g006

Normalization in the aDNA data set

READ uses the average P0 from an unrelated pair of individuals to normalize the distribution for all test individuals. For our empirical data analysis, we assumed the median of all average P0 across pairs of individuals within a test population to represent unrelated individuals, as high values may be caused by recent migrants and low values by related individuals. Fig 7 shows the distributions of all average P0 before normalization highlighting that the populations exhibit different degrees of background diversity. It is also apparent how the pairs of related individuals (see Table 1) are outliers with lower pairwise differences (see also Fig 5). Most groups from similar geographic and cultural groups show similar medians. These include Neolithic groups (except Iberia_EN) and Yamnaya, and—to some degree—Late Neolithic and Bronze Age central Europeans. The latter set of populations could almost belong to two subgroups which cluster by data type (shotgun versus capture) instead of archaeological culture (Unetice, Corded Ware and Bell Beaker). This difference was not observed in Yamnaya for which both data types exist as well. The discrepancy highlights a potential risk of batch effects which has its consequences for the application of READ. Overestimating the distance between unrelated individuals could overestimate relationships in the test group and consequently cause false positives while underestimating the distance between unrelated individuals would have the opposite effect. The extent of the misclassification would be proportional to the ratio between true and used normalization value. For example, if the true value was 0.22 (e.g. Motala_HG, Fig 7) but 0.25 was used (e.g. Hungary_EN), an unrelated pair of individuals could be classified as second degree relatives (0.22/0.25 = 0.88 < 0.90625). Using the shotgun Bell Beaker median (0.245) to normalize the captured Bell Beaker data does not cause any changes in the classifications, whereas using the capture Bell Beaker median (0.257) for the shotgun data would classify RISE563 and RISE564 as second degree relatives. These two individuals might actually be related, but the value used for normalization would be higher than any pairwise comparison within the shotgun sequenced Bell Beakers. This violates the assumption that the normalization value represents the expectation for a pair of unrelated individuals so this result should be considered a false positive due to a batch effect.

thumbnail
Fig 7. Population distributions of average P0 before normalization.

The boxplots show all non-normalized average P0 scores (one per pair of individuals) per culture. CAP and SG indicate whether the individuals were subject to SNP capture or shotgun sequencing, respectively. A broader chronological/geographical context is shown on the left.

https://doi.org/10.1371/journal.pone.0195491.g007

Discussion

Applying READ to aDNA data

Several methods to estimate the degree of relationship between pairs of individuals have been developed. For genome-wide diploid data with low error rates, they successfully infer relationships up to 11th degree [11]. Since such data cannot be obtained from degraded samples, a loss in precision was expected. Estimation of second degree relationships (i.e. niece/nephew-aunt/uncle, grandparent-grandchild, half-siblings) is sufficient to identify individuals belonging to a core family which were buried together. We can show that obtaining as little as 2,500 overlapping common SNPs is enough to classify up to second degree relationships from effectively haploid data. The biggest limitations when using such low numbers of SNPs is the high rate of false negatives for second degree relatives. READ can be considered as a conservative tool that avoids false positives by having a relatively high false negative rate which can be decreased substantially with more data. Missing some second degree relationships seems preferable over wrongly inferring relationships for unrelated individuals. A consequent advantage of our method is that it is very unlikely that first degree relatives are classified as unrelated but some second degree relatives might be wrongly classified as unrelated. Shared uniparental haplotypes or a test result close to the threshold (e.g. less than one standard error difference) could raise such suspicions and might motivate additional sequencing of the samples in question. The amount of overlapping SNPs depends on the genome coverage of both individuals (Fig 4; e.g. two 0.1x individuals will have approximately the same amount of overlapping data as a 0.05x and a 0.2x individual or a 0.01x individual and a 1x individual). The number of SNPs required for a positive classification as first degree can be obtained by shotgun sequencing all individuals to an average genome coverage of 0.1x (Fig 4), which is in reach for most archaeological samples displaying some authentic DNA. More data would be beneficial to avoid false negatives in the case of second degree relatives. Recently developed methods for modern DNA, which use genotype-likelihoods to handle the uncertainty of low to medium coverage data require 1-3x genome coverage to estimate third degree relationships [5053]. A recent study successfully studied social organization in ancient DNA data for samples with ≥ 1x genome coverage [23]. Such approaches are promising for well-preserved samples but these coverages might not be within reach for most aDNA studies. Other methods specifically designed for ancient DNA data either require large reference data sets [47, 54] or are not directly designed to identify relatives and estimate their degrees [55]. A recent development [56] jointly estimates contamination, sequencing errors and relatedness coefficients for aDNA data, but it requires larger sample sizes than READ (N > 32 to accurately classify first degree relatives, N > 48 to classify second degree relatives [56]).

READ does not explicitly model aDNA damage and it only considers one allele at heterozygous sites. This implies that a careful curation of the data is required to avoid errors due to low coverage, short sequence fragments, deamination damage, sequencing errors and potential contamination. We recommend a number of well established filtering steps when working with low coverage aDNA data [2835, 39]. To avoid batch effects, all samples should be processed as similar as possible—at least the bioinformatic pipeline should be identical for all samples. Only fragments of 35 bp or longer should be mapped to the human genome as shorter fragments might represent spuriously mapping microbial contamination [57, 58]. All downstream analysis should be restricted to reads and bases with mapping and base qualities of 30 or higher to reduce the potential effects of mismapping and sequencing errors [58, 59]. To further reduce the effect of sequencing errors, most aDNA studies only consider biallelic SNPs known to be polymorphic in other populations, and call pseudo-haploid genotypes by randomly sampling one read covering that position. Deamination damage can be avoided during the data generation by enzymatic repair of damages [60], or later by computational rescaling of base qualities before SNP calling [61], or by excluding all transition SNPs as only those are affected by deamination damage. For humans, millions of polymorphic transversion sites are known across the genome [48, 62] still leaving substantial amounts of data for analyzing such data sets. Furthermore, a range of methods exist to estimate human contamination of a particular sample [6367] and the analysis could be restricted to fragments displaying characteristic damage to filter contamination [36, 68]. Finally, each study could simulate data exactly resembling the empirical data analyzed (fragment sizes, damages, contamination) to evaluate how these factors would affect the downstream analysis [58].

An important part of the READ pipeline is the normalization step. This step makes the classification independent of within population diversity, SNP ascertainment and marker density. This property, however, requires at least one additional and unrelated individual from the same population and ideally the same data type to avoid batch effects. The assignment of all individuals to a population can be checked with established methods as principal component analysis (PCA) or outgroup f3 statistics [44]. Alternatively, a pair of individuals from a different population with similar expected diversity could be used for normalization. Fig 7 shows that most (but not all) groups from similar cultural and geographical backgrounds have relatively similar normalization scores, but caution should be taken as strong misspecification of the normalization value can cause false negatives or false positives (see Results section). In practice, the relationships are not known a priori. For our data analysis, we assumed that the median across all pairs of individuals from a population of more than four samples represents a proxy of an unrelated pair (as the number of pairs is ; e.g. 10 pairs for a sample size of 5), which we also set as the default mode for READ. The implementation of READ also offers to use the maximum pairwise average P0 which should only be used in cases like supposed parent-child-trios (two first degree relationships, one unrelated), where the maximum value would represent the comparison between supposed mother and supposed father—the only unrelated pair in the sample. Other methods normalize by obtaining allele frequency data for a whole population [50, 54], which seems less feasible than obtaining just one unrelated individual (or a pair of unrelated individuals from a surrogate population). Furthermore, prehistoric populations are quite differentiated from modern groups [33, 39, 41] so using modern populations as references for the allele frequencies might introduce biases [23]. A certain limitation for all kinship estimation methods is if the sampled population itself cannot be considered homogeneous, for example due to varying degrees of admixture. Only quite recent developments in inferring relationships can efficiently deal with such cases for modern data [69].

Kinship in prehistoric populations

We successfully applied READ to data obtained from ancient individuals. READ confidently found all known relationships in the dataset. Furthermore, it identified a number of previously unknown relationships, mainly of second degree. The combination of genomic data, uniparental markers and radiocarbon dating allowed us to infer how two individuals were related to each other. Additional information such as osteological data on the age of the samples or stratigraphic information as burial location or depth could further help to assess the direction of a kinship. Of particular interest was a group of five males from Esperstedt in Germany who were associated with the Corded Ware culture—a culture that arose after large scale migrations of males [70] from the east [28, 30]. Around 50 Corded Ware burials, six of them stone cists, were excavated near Esperstedt in the context of road constructions in 2005 [28, 71]. Characteristic Corded Ware pottery was found in the graves and all male individuals had been buried on their right hand site [71]. Interestingly, the central individual of the group of related individuals (I1541, Fig 6) was buried in a stone cist approximately 700 meters from the graves of the other four individuals which were all close to each other [71]. The close relationship of this group of only male individuals from the same location suggest patrilocality and female exogamy, a pattern which has also been found from Strontium isotopes at another Corded Ware site just 30 kilometers from Esperstedt [15] and suggested for the Corded Ware culture in general [72]. This represents just one example of how the genetic analysis of relationships can be used to uncover and understand social structures in ancient populations. More data from additional sites, cultures and species other than humans will offer various opportunities for the analysis of relationships based on genome-wide data.

Materials and methods

Approach to detect related individuals

Our approach is based on the methodology used by GRAB [14] which was designed for unphased and diploid genotype or sequencing data. This approach divides the genome into non-overlapping windows of 1 Mbps each and compares for a pair of individuals the alleles inside each window. Each SNP is classified into three different categories: IBS2 when the two alleles are shared, IBS1 when only one allele is shared and IBS0 when no allele is shared. The program calculates the fractions for each category (P2, P1 and P0) per window and, based on certain thresholds, uses them for relationship estimation. GRAB can estimate relationships from 1st to 5th degree, but it has not been tested with data from different SNP panels or populations [14].

We assume that our input data stems from whole genome shotgun sequencing of an ancient sample resulting in low coverage sequencing data. In such situations, a common approach in many ancient DNA studies is to randomly sample one sequencing read per individual and SNP site and then use the allele carried on that read as pseudo-haploid information. Such approaches are obviously restricted to a set of biallelic SNPs ascertained in an external dataset. Consequently, we only expect to observe one allele per individual and SNP site which is either shared or not shared between the two individuals. READ does not model aDNA damage, so it is expected that the input is carefully filtered, e.g. by restricting to sites known to be polymorphic, by excluding transition sites or by rescaling base qualities before SNP calling [61]. Analogous to GRAB [14], we partition the genome in non-overlapping windows of 1 Mbps and calculate the proportions of haploid mismatches and matches, P0 and P1, for each window. Since P0 + P1 = 1, we can use P0 as a single test statistic. The average P0 is calculated from the genome-wide distribution. To reduce the effect of SNP ascertainment, population diversity and potential batch effects, each individual pair’s average P0 scores are then normalized by dividing all values by the average non-normalized P0 score from an unrelated pair of individuals from the same population ascertained in the same way as for the tested pairs. Such a normalization step is not implemented in GRAB [14] suggesting that GRAB might be sensitive to ascertainment bias and general population diversity. The normalization sets the expected score for an unrelated pair to 1 and we can define classification cutoffs which are independent of the diversity within the particular data set. We define three thresholds to identify pairwise relatedness as unrelated, second-degree (i.e. nephew/niece-uncle/aunt, grandparent-grandchild or half-siblings), first-degree (parent-offspring or siblings) and identical individuals/identical twins. The general work flow and the decision tree used to classify relationships is shown in Fig 1. There are four possible outcomes when running READ: unrelated (normalized P0≥0.90625), second degree (0.90625≥normalized P0≥0.8125), first degree (0.8125≥normalized P0≥0.625) and identical twins/identical individuals (normalized P0<0.625) (Fig 1). The cutoffs were chosen to lie halfway between the probabilities of one randomly chosen allele for an individual not being IBD to a randomly chosen allele from another individual considering their degree of relationship: 1/2 = 0.5 identical twins/identical individuals, 3/4 = 0.75 for first degree relatives, 7/8 = 0.875 for second degree relatives and 15/16 = 0.9375 for third degree relatives. We do not aim to classify higher degrees than second degree and, therefore, consider all relationships of third degree or higher as ‘unrelated’. This is a decision to keep the approach conservative and to allow for some variation within the group of unrelated individuals. Furthermore, the 1000 genomes data contains very few third degree relatives making it difficult to estimate error rates for this group. READ is implemented to classify pairs of individuals in certain categories, so it will always output the best fitting degree of relationship based on the point estimate of the average P0. As an assessment of confidence of that classification, we estimate the standard error of the mean of the distribution of normalized P0 scores ( where is the empirical standard deviation of P0 across all windows and n is the total number of windows) and calculate the distance to the cutoffs in multiples of the standard error (similar to a Z score also known as ‘standard score’). Furthermore, the user is provided with a graphical output (see Fig 5) showing the average P0 for each pair, their 95% confidence interval, and the cutoffs for classification together with their 95% confidence interval.

Relationship Estimation from Ancient DNA (READ) was implemented in Python 2.7 [73] and GNU R [74]. The input format is TPED/TFAM [8] and READ is publicly available from https://bitbucket.org/tguenther/read and as S1 File.

Modern data with reported degrees of relationships

Autosomal Illumina Omni2.5M chip genotype calls from 1,326 individuals from 15 different populations were obtained from the 1000 genomes project (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/hd_genotype_chip/) [48]. We used vcftools version 0.1.11 [75] to extract autosomal biallelic SNPs with a minor allele frequency of at least 10% (1,156,468 SNPs in total—similar to the aDNA data set used for the empirical data analysis [35]; see below) and to convert the data to TPED/TFAM files. The data set contains pairs of individuals that were reported as related, 851 of them as first degree relationships and 74 as second degree. We randomly sub-sampled 1,000, 2,500, 5,000, 10,000, and 50,000 SNPs and also randomly picked one allele per site in order to mimic extremely low coverage sequencing of ancient samples. In an additional simulation, we introduced different allelic error rates to the data to assess the possible effects of sequencing and mapping errors, contamination and post-mortem damage. Allelic errors were introduced by randomly changing alleles to the alternative based on a per site error rate, the per site error rates are aimed to reflect different error rates in different parts of the genome. Per site error rates were drawn from a Gaussian distribution with a mean corresponding to the average allelic error rate (0.05, 0.1, 0.15 or 0.2) and a standard deviation of 0.01.

READ was then applied to these data sets and the median of all average P0s per population was used to normalize scores assuming that this would represent an unrelated pair. Additionally, the data was tested employing a 10-fold cross-validation procedure, which allowed to infer the expected value of P0 for a pair of unrelated individuals from a different subset of the data than what was used to test the relationships avoiding potential circularity. The average value of P0 obtained for each pair was then used for classification.

To evaluate READ’s performance, we calculate false positive and false negative rates. Unrelated individuals classified as related were considered as false positives, related individuals classified as unrelated or as related but not at the proper degree were considered false negatives. READ’s performance was similar for both normalization approaches (median and cross-validation), so we present the results of using the median in the main text and the cross-validation approach results in Supplementary figures (S1 and S2 Figs). The cross validation approach would require large sample sizes per population which are not reached in most ancient DNA studies (see the empirical data set below for an example).

Ancient data

In addition to the modern data, published ancient data was obtained from the study of Mathieson et al. (2015) [35]. The data set consisted of 230 ancient Europeans from a number of publications [28, 3033, 76] as well as new individuals from various time periods during the last 8,500 years. The data set consisted of haploid data for up to 1,209,114 SNPs per individual. We extracted only autosomal data for all individuals and applied READ to each cultural or geographical group (as defined in the original data set of Mathieson et al (2015) [35]) with more than four individuals separately. Shotgun sequencing data was also analyzed separately from SNP capture data to avoid batch effects. The median of all average P0s per group was used for normalization assuming that this would represent an unrelated pair. Mathieson et al (2015) [35] report nine pairs of related individuals and they infer all of them to be first degree relatives without providing details on how they were classified. Y-chromosome haplotypes of the five individuals shown in Fig 6 were checked using samtools [77] (applying a minimum mapping and base quality of 30) and marker information for the haplotypes R1a and R1b from the International Society of Genetic Genealogy (http://www.isogg.org, accessed January 16, 2017). The results are shown in S1 Table.

Supporting information

S1 Fig. Simulation based estimates of false positive and false negative rates for different numbers of SNPs estimated using a cross validation scheme.

Compare Fig 2.

https://doi.org/10.1371/journal.pone.0195491.s001

(PDF)

S2 Fig. Effect of allelic errors on READ’s performance, simulations conducted employing a cross validation scheme.

Compare Fig 3.

https://doi.org/10.1371/journal.pone.0195491.s002

(PDF)

S1 Table. Y chromosome calls for haplogroup R defining markers in the five individuals shown in Fig 6.

https://doi.org/10.1371/journal.pone.0195491.s003

(XLS)

S2 Table. Pairs of first degree related individuals in the 1000 genomes data.

https://doi.org/10.1371/journal.pone.0195491.s004

(XLS)

S3 Table. Pairs of second degree related individuals in the 1000 genomes data.

https://doi.org/10.1371/journal.pone.0195491.s005

(XLS)

S4 Table. Classification performance for YRI and IBS when one half of the data had a different allelic error rate than the other half.

https://doi.org/10.1371/journal.pone.0195491.s006

(DOC)

S1 File. READ methodology implemented in R and Python scripts.

https://doi.org/10.1371/journal.pone.0195491.s007

(GZ)

Acknowledgments

We thank Federico Sanchez-Quinto, Jan Storå, Rita Peyroteo Stjerna, four anonymous reviewers, and those who commented on the preprint for constructive feedback on the manuscript as well as Gülşah Merve Dal Kılınç and Mehmet Somel for discussions on the approach. Computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under projects b2013203 and b2015094.

References

  1. 1. Kong A, Masson G, Frigge ML, Gylfason A, Zusmanovich P, Thorleifsson G, et al. Detection of sharing by descent, long-range phasing and haplotype imputation. Nat Genet. 2008;40(9):1068–1075. pmid:19165921
  2. 2. Palin K, Campbell H, Wright AF, Wilson JF, Durbin R. Identity-by-descent-based phasing and imputation in founder populations using graphical models. Genet Epidemiol. 2011;35(8):853–860. pmid:22006673
  3. 3. Zuk O, Hechter E, Sunyaev SR, Lander ES. The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci U S A. 2012;109(4):1193–1198. pmid:22223662
  4. 4. Browning SR, Browning BL. Identity by descent between distant relatives: detection and applications. Annu Rev Genet. 2012;46:617–633. pmid:22994355
  5. 5. Ralph P, Coop G. The geography of recent genetic ancestry across Europe. PLoS Biol. 2013;11(5):e1001555. pmid:23667324
  6. 6. Albrechtsen A, Moltke I, Nielsen R. Natural selection and the distribution of identity-by-descent in the human genome. Genetics. 2010;186(1):295–308. pmid:20592267
  7. 7. Weir BS, Anderson AD, Hepler AB. Genetic relatedness analysis: modern data and new challenges. Nat Rev Genet. 2006;7(10):771–780. pmid:16983373
  8. 8. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–575. pmid:17701901
  9. 9. Roberson EDO, Pevsner J. Visualization of shared genomic regions and meiotic recombination in high-density SNP data. PLoS One. 2009;4(8):e6711. pmid:19696932
  10. 10. Huff CD, Witherspoon DJ, Simonson TS, Xing J, Watkins WS, Zhang Y, et al. Maximum-likelihood estimation of recent shared ancestry (ERSA). Genome Research. 2011;21(5):768–774. pmid:21324875
  11. 11. Li H, Glusman G, Hu H, Shankaracharya , Caballero J, Hubley R, et al. Relationship estimation from whole-genome sequence data. PLoS Genet. 2014;10(1):e1004144. pmid:24497848
  12. 12. Manichaikul A, Mychaleckyj JC, Rich SS, Daly K, Sale M, Chen WM. Robust relationship inference in genome-wide association studies. Bioinformatics. 2010;26(22):2867–2873. pmid:20926424
  13. 13. Thornton T, Tang H, Hoffmann TJ, Ochs-Balcom HM, Caan BJ, Risch N. Estimating kinship in admixed populations. The American Journal of Human Genetics. 2012;91(1):122–138. pmid:22748210
  14. 14. Li H, Glusman G, Huff C, Caballero J, Roach JC. Accurate and robust prediction of genetic relationship from whole-genome sequences. PLoS One. 2014;9(2):e85437. pmid:24586241
  15. 15. Haak W, Brandt G, de Jong HN, Meyer C, Ganslmeier R, Heyd V, et al. Ancient DNA, Strontium isotopes, and osteological analyses shed light on social and kinship organization of the Later Stone Age. Proc Natl Acad Sci U S A. 2008;105(47):18226–18231. pmid:19015520
  16. 16. King TE, Fortes GG, Balaresque P, Thomas MG, Balding D, Delser PM, et al. Identification of the remains of King Richard III. Nature Communications. 2014;5.
  17. 17. Oliehoek PA, Windig JJ, van Arendonk JAM, Bijma P. Estimating relatedness between individuals in general populations with a focus on their use in conservation programs. Genetics. 2006;173(1):483–496. pmid:16510792
  18. 18. Habier D, Fernando R, Dekkers J. The impact of genetic relationship information on genome-assisted breeding values. Genetics. 2007;177(4):2389–2397. pmid:18073436
  19. 19. Pääbo S. Ancient DNA: extraction, characterization, molecular cloning, and enzymatic amplification. Proceedings of the National Academy of Sciences. 1989;86(6):1939–1943.
  20. 20. Briggs AW, Stenzel U, Johnson PL, Green RE, Kelso J, Prüfer K, et al. Patterns of damage in genomic DNA sequences from a Neandertal. Proceedings of the National Academy of Sciences. 2007;104(37):14616–14621.
  21. 21. Sawyer S, Krause J, Guschanski K, Savolainen V, Pääbo S. Temporal patterns of nucleotide misincorporations and DNA fragmentation in ancient DNA. PLoS One. 2012;7(3):e34131. pmid:22479540
  22. 22. Ensor BE. The archaeology of kinship: Advancing interpretation and contributions to theory. University of Arizona Press; 2013.
  23. 23. Sikora M, Seguin-Orlando A, Sousa VC, Albrechtsen A, Korneliussen T, Ko A, et al. Ancient genomes show social and reproductive behavior of early Upper Paleolithic foragers. Science. 2017; p. eaao1807.
  24. 24. Canturk KM, Emre R, Kınoglu K, Başpınar B, Sahin F, Ozen M. Current status of the use of single-nucleotide polymorphisms in forensic practices. Genet Test Mol Biomarkers. 2014;18(7):455–460. pmid:24754266
  25. 25. Deguilloux M, Pemonge M, Mendisco F, Thibon D, Cartron I, Castex D. Ancient DNA and kinship analysis of human remains deposited in Merovingian necropolis sarcophagi (Jau Dignac et Loirac, France, 7th–8th century AD). Journal of Archaeological Science. 2014;41:399–405.
  26. 26. Cui Y, Song L, Wei D, Pang Y, Wang N, Ning C, et al. Identification of kinship and occupant status in Mongolian noble burials of the Yuan Dynasty through a multidisciplinary approach. Philos Trans R Soc Lond B Biol Sci. 2015;370(1660):20130378. pmid:25487330
  27. 27. Hughes-Stamm SR, Ashton KJ, van Daal A. Assessment of DNA degradation and the genotyping success of highly degraded samples. International Journal of Legal Medicine. 2011;125(3):341–348. pmid:20419381
  28. 28. Haak W, Lazaridis I, Patterson N, Rohland N, Mallick S, Llamas B, et al. Massive migration from the steppe was a source for Indo-European languages in Europe. Nature. 2015;522(7555):207–211. pmid:25731166
  29. 29. Skoglund P, Malmström H, Raghavan M, Storå J, Hall P, Willerslev E, et al. Origins and genetic legacy of Neolithic farmers and hunter-gatherers in Europe. Science. 2012;336(6080):466–469. pmid:22539720
  30. 30. Allentoft ME, Sikora M, Sjögren KG, Rasmussen S, Rasmussen M, Stenderup J, et al. Population genomics of Bronze Age Eurasia. Nature. 2015;522(7555):167–172. pmid:26062507
  31. 31. Lazaridis I, Patterson N, Mittnik A, Renaud G, Mallick S, Kirsanow K, et al. Ancient human genomes suggest three ancestral populations for present-day Europeans. Nature. 2014;513(7518):409–413. pmid:25230663
  32. 32. Gamba C, Jones ER, Teasdale MD, McLaughlin RL, Gonzalez-Fortes G, Mattiangeli V, et al. Genome flux and stasis in a five millennium transect of European prehistory. Nat Commun. 2014;5:5257. pmid:25334030
  33. 33. Skoglund P, Malmström H, Omrak A, Raghavan M, Valdiosera C, Günther T, et al. Genomic diversity and admixture differs for Stone-Age Scandinavian foragers and farmers. Science. 2014;344(6185):747–750. pmid:24762536
  34. 34. Günther T, Valdiosera C, Malmström H, Ureña I, Rodriguez-Varela R, Sverrisdóttir ÓO, et al. Ancient genomes link early farmers from Atapuerca in Spain to modern-day Basques. Proc Natl Acad Sci U S A. 2015;112(38):11917–11922. pmid:26351665
  35. 35. Mathieson I, Lazaridis I, Rohland N, Mallick S, Patterson N, Roodenberg SA, et al. Genome-wide patterns of selection in 230 ancient Eurasians. Nature. 2015;528(7583):499–503. pmid:26595274
  36. 36. Fu Q, Hajdinjak M, Moldovan OT, Constantin S, Mallick S, Skoglund P, et al. An early modern human from Romania with a recent Neanderthal ancestor. Nature. 2015;524(7564):216–219. pmid:26098372
  37. 37. Cassidy LM, Martiniano R, Murphy EM, Teasdale MD, Mallory J, Hartwell B, et al. Neolithic and Bronze Age migration to Ireland and establishment of the insular Atlantic genome. Proceedings of the National Academy of Sciences. 2016;113(2):368–373.
  38. 38. Hofmanová Z, Kreutzer S, Hellenthal G, Sell C, Diekmann Y, Díez-del Molino D, et al. Early farmers from across Europe directly descended from Neolithic Aegeans. Proceedings of the National Academy of Sciences. 2016; p. 201523951.
  39. 39. Lazaridis I, Nadel D, Rollefson G, Merrett DC, Rohland N, Mallick S, et al. Genomic insights into the origin of farming in the ancient Near East. Nature. 2016.
  40. 40. Slatkin M, Racimo F. Ancient DNA and human history. Proceedings of the National Academy of Sciences. 2016;113(23):6380–6387.
  41. 41. Günther T, Jakobsson M. Genes mirror migrations and cultures in prehistoric Europe—a population genomic perspective. Current Opinion in Genetics & Development. 2016;41:115–123.
  42. 42. Saag L, Varul L, Scheib CL, Stenderup J, Allentoft ME, Saag L, et al. Extensive Farming in Estonia Started through a Sex-Biased Migration from the Steppe. Current Biology. 2017;27(14):2185–2193.e6. pmid:28712569
  43. 43. Rasmussen M, Anzick SL, Waters MR, Skoglund P, DeGiorgio M, Stafford TW Jr, et al. The genome of a Late Pleistocene human from a Clovis burial site in western Montana. Nature. 2014;506(7487):225–229. pmid:24522598
  44. 44. Raghavan M, Skoglund P, Graf KE, Metspalu M, Albrechtsen A, Moltke I, et al. Upper Palaeolithic Siberian genome reveals dual ancestry of Native Americans. Nature. 2014;505(7481):87–91. pmid:24256729
  45. 45. Yang MA, Gao X, Theunert C, Tong H, Aximu-Petri A, Nickel B, et al. 40,000-Year-Old Individual from Asia Provides Insight into Early Population Structure in Eurasia. Current Biology. 2017;27(20):3202–3208.e9. pmid:29033327
  46. 46. Kennett DJ, Plog S, George RJ, Culleton BJ, Watson AS, Skoglund P, et al. Archaeogenomic evidence reveals prehistoric matrilineal dynasty. Nature Communications. 2017;8. pmid:28221340
  47. 47. Martin MD, Jay F, Castellano S, Slatkin M. Determination of genetic relatedness from low-coverage human genome sequences using pedigree simulations. Molecular Ecology. 2017;.
  48. 48. Consortium GP, Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68–74.
  49. 49. Wheeler DA, Srinivasan M, Egholm M, Shen Y, Chen L, McGuire A, et al. The complete genome of an individual by massively parallel DNA sequencing. Nature. 2008;452(7189):872–876. pmid:18421352
  50. 50. Lipatov M, Sanjeev K, Patro R, Veeramah K. Maximum likelihood estimation of biological relatedness from low coverage sequencing data. bioRxiv. 2015; p. 023374.
  51. 51. Korneliussen TS, Moltke I. NgsRelate: a software tool for estimating pairwise relatedness from next-generation sequencing data. Bioinformatics. 2015;31(24):4009–4011. pmid:26323718
  52. 52. Vieira FG, Albrechtsen A, Nielsen R. Estimating IBD tracts from low coverage NGS data. Bioinformatics. 2016;32(14):2096–2102. pmid:27153648
  53. 53. Waples RK, Albrechtsen A, Moltke I. Allele frequency-free inference of close familial relationships from genotypes or low depth sequencing data. bioRxiv. 2018; p. 260497.
  54. 54. Fernandes D, Sirak K, Novak M, Finarelli JA, Byrne J, Connolly E, et al. The Identification of a 1916 Irish Rebel: New Approach for Estimating Relatedness From Low Coverage Homozygous Genomes. Scientific Reports. 2017;7.
  55. 55. Vohr SH, Najar CFBA, Shapiro B, Green RE. A method for positive forensic identification of samples from extremely low-coverage sequence data. BMC Genomics. 2015;16(1):1034. pmid:26643904
  56. 56. Theunert C, Racimo F, Slatkin M. Joint estimation of relatedness coefficients and allele frequencies from ancient samples. Genetics. 2017; p. genetics–117. pmid:28396504
  57. 57. Meyer M, Arsuaga JL, de Filippo C, Nagel S, Aximu-Petri A, Nickel B, et al. Nuclear DNA sequences from the Middle Pleistocene Sima de los Huesos hominins. Nature. 2016;531(7595):504–507. pmid:26976447
  58. 58. Renaud G, Hanghøj K, Willeslev E, Orlando L. gargammel: a sequence simulator for ancient DNA. Bioinformatics. 2016; p. btw670.
  59. 59. Schubert M, Ginolhac A, Lindgreen S, Thompson JF, Al-Rasheid KA, Willerslev E, et al. Improving ancient DNA read mapping against modern reference genomes. BMC Genomics. 2012;13(1):178. pmid:22574660
  60. 60. Briggs AW, Stenzel U, Meyer M, Krause J, Kircher M, Pääbo S. Removal of deaminated cytosines and detection of in vivo methylation in ancient DNA. Nucleic Acids Research. 2010;38(6):e87–e87. pmid:20028723
  61. 61. Jónsson H, Ginolhac A, Schubert M, Johnson PL, Orlando L. mapDamage2.0: fast approximate Bayesian estimates of ancient DNA damage parameters. Bioinformatics. 2013; p. btt193.
  62. 62. Mallick S, Li H, Lipson M, Mathieson I, Gymrek M, Racimo F, et al. The Simons genome diversity project: 300 genomes from 142 diverse populations. Nature. 2016;538(7624):201. pmid:27654912
  63. 63. Green RE, Malaspinas AS, Krause J, Briggs AW, Johnson PL, Uhler C, et al. A complete Neandertal mitochondrial genome sequence determined by high-throughput sequencing. Cell. 2008;134(3):416–426. pmid:18692465
  64. 64. Rasmussen M, Guo X, Wang Y, Lohmueller KE, Rasmussen S, Albrechtsen A, et al. An Aboriginal Australian genome reveals separate human dispersals into Asia. Science. 2011;334(6052):94–98. pmid:21940856
  65. 65. Jun G, Flickinger M, Hetrick KN, Romm JM, Doheny KF, Abecasis GR, et al. Detecting and estimating contamination of human DNA samples in sequencing and array-based genotype data. The American Journal of Human Genetics. 2012;91(5):839–848. pmid:23103226
  66. 66. Fu Q, Mittnik A, Johnson PL, Bos K, Lari M, Bollongino R, et al. A revised timescale for human evolution based on ancient mitochondrial genomes. Current Biology. 2013;23(7):553–559. pmid:23523248
  67. 67. Renaud G, Slon V, Duggan AT, Kelso J. Schmutzi: estimation of contamination and endogenous mitochondrial consensus calling for ancient DNA. Genome Biology. 2015;16(1):224. pmid:26458810
  68. 68. Skoglund P, Northoff BH, Shunkov MV, Derevianko AP, Pääbo S, Krause J, et al. Separating endogenous ancient DNA from modern day contamination in a Siberian Neandertal. Proceedings of the National Academy of Sciences. 2014;111(6):2229–2234.
  69. 69. Moltke I, Albrechtsen A. RelateAdmix: a software tool for estimating relatedness between admixed individuals. Bioinformatics. 2014;30(7):1027–1028. pmid:24215025
  70. 70. Goldberg A, Günther T, Rosenberg NA, Jakobsson M. Ancient X chromosomes reveal contrasting sex bias in Neolithic and Bronze Age Eurasian migrations. Proceedings of the National Academy of Sciences. 2017;114(10):2657–2662.
  71. 71. Leinthaler B, Bogen C, Döhle H. Von Muschelknöpfen und Hundezähnen—Schnurbandkeramische Bestattungen in Esperstedt. In: Meller H, Dresely V, editors. Archäologie auf der Überholspur. Ausgrabungen an der A38. Landesamt für Archäologie Sachsen-Anhalt; 2006. p. 59–82.
  72. 72. Sjögren KG, Price TD, Kristiansen K. Diet and Mobility in the Corded Ware of Central Europe. PLoS One. 2016;11(5):e0155083. pmid:27223117
  73. 73. Van Rossum G, et al. Python Programming Language. In: USENIX Annual Technical Conference. vol. 41; 2007.
  74. 74. R Core Team. R: A Language and Environment for Statistical Computing; 2016. Available from: https://www.R-project.org/.
  75. 75. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–2158. pmid:21653522
  76. 76. Keller A, Graefen A, Ball M, Matzas M, Boisguerin V, Maixner F, et al. New insights into the Tyrolean Iceman’s origin and phenotype as inferred by whole-genome sequencing. Nat Commun. 2012;3:698. pmid:22426219
  77. 77. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–2079. pmid:19505943