Figures

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.

Data Availability: All relevant data are within the manuscript and Supporting Information files.

Funding: JMMK received scholarships from the Erasmus Mundus Master Programme in Evolutionary Biology (MEME) and the Consejo Nacional de Ciencia y Tecnología (CONACYT). MJ was supported by an European Research Council starting grant (no. 311413), a Swedish Research Council grant, Knut och Alice Wallenbergs Stiftelse, Knut och Alice Wallenbergs Stiftelse, and Vetenskapsrådet. TG was supported by a Wenner-Gren-Foundations fellowship. Computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under projects b2013203 and b2015094. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

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 [19–21]. 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, 28–42] 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, 28–31, 33–36, 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).

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).

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.

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.

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.

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].

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.

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.

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.

(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.

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).

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.

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.

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 [50–53]. 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 [28–35, 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 [63–67] 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.

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, 30–33, 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.

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

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.