Sample collection

All AncestryDNA samples included in this study were collected from AncestryDNA customers, who have agreed to the informed consent for the Independent Review Board approved Ancestry Human Diversity Project (Quorum Review #26168/1)—an AncestryDNA sponsored research project. Samples obtained from external sample collections were included only in the ADMIXTURE reference panel and were used in accordance with any applicable restrictions. Following sample quality-control steps described in Supplementary Methods, we obtain a final panel of 774,516 genotyped samples consented to participate in research. Consult Supplementary Methods for more details on sample collection procedures.

Genealogical data

We compile statistics from pedigree data linked to genotyped individuals to better understand the historical and geographical context for the IBD clustering. After quality-control steps (Supplementary Methods), we obtain a final set of 432,611 genotyped samples linked to non-private pedigrees. We only include pedigree nodes corresponding to ancestors of the genotyped individuals. We use two types of information in our analyses of associated pedigree data: birth year and birth location (longitude and latitude). Among DNA samples linked to pedigrees, 322,683 (96% of reported birth locations) were born in the United States, 13,748 (4% of reported birth locations) were born outside the United States and 96,180 (22% of all DNA samples linked to pedigrees) have unreported birth locations. See Supplementary Methods for additional details.

Genotyping

Customer saliva sample accessioning, DNA extraction, and genotyping were completed by the Illumina FastTrack Microarray Services labs. Customer genotype data for this study was generated using the Illumina Human OmniExpress platform. This genotyping array assays 730,525 SNPs (709,358 SNPs are on autosomal chromosomes) and small indels across the genome. The SNPs on this array were carefully selected to capture the majority of common genetic variation in European and other worldwide populations. (Note that genotype data on sex-linked chromosomes are not used in this study). Genotypes were called by Illumina technicians using the GenomeStudio platform.

Genotype quality-control procedures

We perform extensive quality-control checks for each processed AncestryDNA customer sample. We have developed an almost entirely automatic workflow to signal problematic genotype samples. This process includes the following steps: (1) identifying and removing any duplicate samples (these samples are possibly identical twins); (2) identifying and removing samples with a per-sample call rate <98%, since low call rate can be indicative of either a poor-quality DNA sample or a technical failure in the genotyping array; (3) identifying discrepancies between gender inferred from the DNA sample and self-reported gender recorded by the customer during test activation; and (4) flagging an unusually high rate of genotype heterozygosity, which can indicate sample cross-contamination. We tolerate small numbers of samples failing checks 3 and 4. Additional tests are applied to each batch of 96 samples on a genotyping plate to identify and remove cases of large-scale laboratory mistakes such as incorrect plate rotation. If more than two samples on a single plate have mismatching genders or high heterozygosity, the entire plate of samples is withheld for manual examination. Samples that pass all quality-control tests proceed to the analysis pipeline; samples that fail one or more of the above tests must be collected again from customers, or manually cleared for analysis by lab technicians.

SNP quality-control procedures

As these data were generated over a period of 3 years, the DNA product has gone through three revisions, mainly due to beadpool depletion and updates to array geometry. To minimize the impact of different versions of the array on any inferences made from the genotype data, we have developed a versioning quality-control procedure. Briefly, we re-genotype a collection of 180 reference samples for each new array version, and assess genotyping performance in terms of SNP presence/absence, genotype call rate and genotype discordance. If a SNP fails at least one of these assessments for a particular version of the array—for example, the SNP is absent due to probe dropout in re-manufacturing, its call rate is <95%, or its genotype discordance is >1%—the genotypes for that SNP are treated as missing for all samples processed on this version of the array. In the worst case, the genotypes of ∼2.6% SNPs were treated as missing. Any genotypes treated as missing were later imputed (see below).

Genotyping phasing

Our method for phasing genotypes is similar to BEAGLE51, but with substantial improvements to achieve high phasing accuracy on a large scale. Our strategy is to first learn BEAGLE haplotype models from a reference panel of 217,722 genotype samples previously phased at 633,299 autosomal SNPs, which includes 558 phased samples from Phase 1 of the 1000 Genomes project34,51, then we use these models to quickly phase new samples. See Supplementary Methods for details on design of the phasing reference panel, modifications to BEAGLE to accommodate large data sets, and empirical evaluation of the phasing algorithm.

Detecting IBD

For computational scalability, we have developed our own distributed-processing implementation of the GERMLINE algorithm52 to detect long chromosome segments suggestive of inheritance from a recent common ancestor (‘IBD segments’). We apply this method to the phased genotypes to detect pairs with total shared IBD >5 cM. Although other methods have been developed for accurate detection of IBD53, most of these methods scale quadratically with the number of samples, hence are not suitable for our data. GERMLINE has been shown to be particularly inaccurate for IBD segments <4 cM28, but this is less of a concern here since these short IBD segments contribute little or no weight to the IBD network.

Our algorithm produces the same output as GERMLINE, but offers two computational advantages that allow us to efficiently handle hundreds of thousands of phased genotypes: first, it distributes the computation over a Hadoop compute cluster; second, it stores the phased genotypes in a database so that new samples can be efficiently compared with previously processed samples. The inputs to our GERMLINE implementation are the phased genotypes, and the genetic distances (in cM) between consecutive pairs of SNPs. (We use interpolated HapMap genetic distances54). The algorithm output reproduces the same result as GERMLINE 1.5.1 with the following command-line arguments: germline -bits 96 -err_hom 0 -err_het 0 -min_m 5. In an effort to reduce computational cost and the rate of false positives in identifying IBD segments >5 cM among all pairs of >700,000 individuals, we set the ‘bits’ parameter to a value slightly larger than some have recommended55, and we do not tolerate homozygous or heterozygous mismatches. The empirical distribution of IBD grows exponentially with decreasing IBD length, as expected (Supplementary Fig. 28).

Visualizing geographic patterns of IBD

To get an initial suggestion that IBD in our large, un-curated sample might be informative of demography within the United States, in a preliminary phase of our study we compiled IBD statistics aggregated by US state, in which samples are assigned to states based on self-reported birth location. More precisely, we tabulate the total amount of IBD shared by individuals born in the same state (‘within-state IBD’) and between states (‘cross-state IBD’) from 322,683 genotypes with self-reported US birth locations. We only count pairs if total IBD is >12 cM. To summarize the state-level IBD data in a single figure (Fig. 1), we use kernel PCA56, implemented in R package kernlab57. We project US states onto the first two PCs, in which the kernel matrix is defined by total cross-state IBD, normalized to eliminate the effect of variation in total within-state IBD on the projection. More precisely, we define an n × n kernel matrix with entries K(i, j)=t(i, j)/{d(i) d(j)}1/2, where n=51 is the number of states including Washington, DC, t(i, j) is the cross-state IBD for i≠j and within-state IBD for i=j, (Supplementary Data 1), and d(i)=t(i, 1)+t(i, 2)+…+t(i, n) is the total IBD between state i and all other states, including within-state IBD t(i, i).

Estimating admixture proportions from genotype data

Based on an individual’s genotypes at 112,909 SNPs, we estimate the proportion of their genome that is attributed to different ancestral populations. We have subdivided the global human population into K=26 regions (Supplementary Fig. 7; Supplementary Table 2). We use the program ADMIXTURE12, together with a curated reference panel of 3,000 putatively single-origin (‘labelled’) genotype samples, to jointly estimate admixture proportions in the unlabelled customer samples from their genotypes. This panel includes 855 samples from the Human Genome Diversity Project58,59, samples from a proprietary AncestryDNA reference collection, and putatively single-origin samples from consenting AncestryDNA customers.

To validate our admixture estimates against public collections where the genotype samples have been annotated by experts, we estimated admixture proportions in 1,043 Human Genome Diversity Project (HGDP) samples (of which 855 are also included in the reference panel) and 1,816 unrelated samples genotyped as part of the 1000 Genomes project33 using the Illumina OMNI 2.5M genotyping chip. (Note that none of these 1000 Genomes samples are included in the reference panel). These results are summarized in Supplementary Data 3. Overall, the admixture estimates align closely with the expert-provided population labels; that is, individuals are assigned higher proportions in the appropriate ancestral populations. For example, the Han HGDP samples are attributed 86–96% to the Asia East population. Many of these test samples are expected to be admixed (for example, ASW and MXL), and exhibit admixture in the expected proportions. Two of the defined ancestral populations, ‘Mali’ and ‘European Jewish’, could not be validated from these data because these ancestral populations contribute at most a small proportion to any of the test samples. We note that admixture proportions for populations that are less genetically differentiated, such as Great Britain and Europe West, are expected to exhibit less accurate admixture estimates using ADMIXTURE; these are typically combined in our calculations when reporting final admixture statistics for each cluster (Table 1; Supplementary Data 2). Finally, since few samples included in our IBD analyses have large contributions from individual West African ancestral populations, we collapse admixture proportions for all West African regions into a single statistic, for a total of 20 reported populations. Consult Supplementary Methods for more details on design of the reference panel, selection of SNPs for admixture calculations, and ADMIXTURE parameter settings.

Constructing the IBD network

To define the IBD network, we apply a weight function, w[e(i, j)] € [0,1], to each edge e(i, j). We define w[e(i, j)] as the proportion of total IBD lengths observed in simulated genotypes that are due to relationships separated by at most eight reproductive events, or meioses (corresponding to common ancestors at most four generations back), although it may reflect more distant relationships for subpopulations that conform less closely to our simulations. This empirical distribution is fit to the Beta cumulative density function, and this fitted distribution (with scale parameters α=2, β=200) defines the weights for all edges in the network (see Supplementary Fig. 2 for more details). We remove all edges corresponding to pairs with total IBD <12 cM since they signal the target familial relationships <6% of the time, and therefore contribute little weight to the network. See Supplementary Methods and Supplementary Figs 29,30 for a description of the simulations, and a detailed rationale for this choice of edge weight function.

Hierarchical clustering of IBD network

To identify network modules, we employ a simple and fast heuristic algorithm, the multi-level or Louvain method29, implemented in the igraph R package60, which heuristically maximizes the modularity by recursively merging subgraphs. (Note that the multi-level algorithm internally generates a hierarchy as it iteratively optimizes the modularity, but we do not use this internal hierarchy in our results). In an attempt to reduce clustering of ‘extended families,’ before running the clustering algorithm we remove all edges in the network corresponding to total shared IBD >72 cM. Since this represents only 0.2% of all edges, removing these edges has little effect on our ability to detect larger modules. In addition to the multi-level method, we tested two other methods, both implemented in igraph, that have low computational complexity—O(m), where m is the number of edges—and so could feasibly be applied to our network: the Infomap method61, and the label propagation method of Raghavan et al.62. Although all three provided similar higher-level clustering, only the multi-level method was able to identify substructure within the portion of the network that represents the vast majority of the sample—samples primarily of European or African descent (African Americans). Even though the multi-level method partitions this sub-network of 687,470 genotyped individuals into only two clusters, no other tested algorithms were able to identify non-trivial structure within this sub-network. This constitutes the main motivation for using the multi-level community detection algorithm. See Supplementary Methods for more details, including the procedure for recursively subdividing the IBD network using the multi-level algorithm.

Spectral analysis of IBD network

The spectral analysis is based on the Laplacian eigenmaps method30, which has close connections to spectral clustering24. Here we briefly describe computation of the spectral embedding; our procedure for identifying ‘stable subsets’ from the spectral embedding is described in the Supplementary Methods. The Laplacian eigenmaps method is derived from a spectral decomposition of the (normalized) Laplacian matrix, L=D−1/2 WD−1/2, where W is the n × n weighted adjacency matrix with entries W(i, j)=w[e(i, j)], and D is the n × n diagonal matrix in which diagonal entry D(i, i) is equal to the degree of node i, or the sum of the edge weights w[e(i, j)] for individual i. Here we define W(i, i)=1 for all i, so that there is always a nonzero probability of remaining at the same occupancy state in a random walk of the graph24. We define the spectral embedding as the first m eigenvectors of the normalized Laplacian. Here we limit each spectral embedding to the top m=40 eigenvectors, primarily for manageability of the analysis procedure. We cannot use the ‘eigengap’ heuristic26 to choose m because it only appropriate to use if the network contains well-pronounced modules24, which is not the case here.

Once we have completed the spectral analysis of the completely connected graph with 769,444 vertices, we compute a second spectral embedding from a subgraph with 586,147 vertices that is obtained by first removing the small sets of individuals and the clusters that project away from the origin in the initial spectral embedding. This step is taken because the spectral decomposition captures the most dominant modular structure in the network, and possibly obscures other, more subtly disconnected subsets. Briefly, to define stable subsets, we visualize the spectral embedding, labelled by the hierarchical clustering, and extract subsets with the same label that project away from the origin. Using this method (described in detail in the supplement and in Supplementary Fig. 31), we identify 18 stable subsets from the spectral embedding (Supplementary Data 2): 10 in the initial spectral embedding, and 8 more in the subgraph embedding. Projecting new samples not included in the IBD network onto the spectral embedding—in particular, the 1000 Genomes samples used as a validation—is described in the supplement.

Historical and geographic interpretation of clusters

Once we have completed hierarchical clustering and spectral analysis of the IBD network, we use the available annotations to investigate how the clusters relate to demography. To accomplish this, we identify features that distinguish members of the cluster, then we deduce a likely demographic scenario from these distinguishing features. For this analysis, we rely on two sets of features: (1) admixture proportions in 20 global populations estimated from the genotypes; and (2) ancestral birth dates and locations from pedigrees associated with some genotyped individuals. To simplify the presentation of admixture summary statistics, some population labels used to define the summary statistics are taken as combinations of ancestral populations; for example, we define ‘Europe West’ as Ireland, Great Britain, Scandinavia, and the region containing Germany and France.

We generate birth location maps by converting each birth location, within a specified range of generations, to the nearest coordinate on a two-dimensional grid, with grid points every 0.5° of latitude and longitude. Then, we count the number of birth locations at each grid point. The location of each grid point plotted on the map is the mean latitude and longitude over all the annotations assigned to that grid point. By scaling the area of each grid point by the number of birth annotations at that location, the maps yield population density estimates, and highlight large urban areas, at different time periods in the United States and Europe. All maps are produced in the same way, differing only in the granularity of the grid and the scale of the plotted points. The distribution of ancestral birth locations by generation recapitulate broad population trends in the United States and in Europe, such as increasing concentration in urban areas over time, increases in population density west of the Mississippi River reflecting westward expansion of European settlement, and migration trends from Europe to the United States (Supplementary Figs 16,17).

To discover geographic features characteristic of a given cluster, we compile statistics from genealogical data specific to each cluster. Specifically, we compute, for each grid point, the odds ratio (OR) for a given cluster—the odds that the grid point is associated with a cluster member over the odds that the grid point is associated with a non-member—then we visualize the distribution of map locations with the largest odds ratios. One rationale for using the OR statistic is that it is informative of cluster prediction accuracy; if we label all map locations with OR>x as ‘ground-truth cluster locations’, then cluster assignments will yield a higher rate of true positives (recovered ground-truth cluster locations) for larger x (assuming the map location frequency remains the same). To highlight the geographic concentration of individual clusters in Fig. 4, we plot only locations satisfying OR>x, with x chosen separately for each cluster. All plotted map locations require a minimum of 10 birth locations associated with cluster members. In some cases, the geographic concentration of birth locations becomes more apparent when the OR calculations are restricted to certain pedigree generations; for example, the birth locations of ancestors 0–5 generations ago associated with the Utah cluster are more highly concentrated in Utah, and other ancestral generations are more dispersed across the eastern United States (Supplementary Fig. 24). Although this strategy is useful for characterizing most clusters, in the Supplementary Discussion we point out some limitations in using the genealogical data to interpret the clusters. For example, our ability to accurately infer demographic trends depends on the composition of the AncestryDNA database, and the availability of genealogical records. Also, note that ORs are more relevant than evidence for enrichment (for example, P value); for example, consider that we would often expect strong evidence for enrichment in large US cities even when the OR is only slightly >1 (for example, association of Chicago with African American cluster), but these locations would provide relatively little information for interpreting the cluster given that many other groups have typically settled in large US cities.

Genetic differentiation between clusters

We calculate pairwise F ST to assess genetic differentiation in common variation between clusters. We include 611,560 SNPs on autosomal chromosomes in the F ST calculations, a subset of the 633,299 SNPs used in phasing that have genotype call rate >95% in the sample. We use the ratio-of-averages formula that also adjusts for differences in sample size63.

Data availability

The HGDP58,59 genotype samples included in the ancestry reference panel were obtained from the HGDP website ( http://www.hagsc.org/hgdp). 1000 Genomes Project Phase 3 (ref. 3333) genotype samples that were used for validation of admixture estimates and the spectral analysis were downloaded from the NCBI FTP site ( ftp://ftp-trace.ncbi.nih.gov/1000genomes). For the purpose of ensuring reproducibility, we will share the IBD network topology, edge weights, and cluster labels on request and subject to relevant data use policies. Although we cannot make the genealogical and genotype data widely available to the academic community in light of our commitment to our customers, we are interested to pursue research collaboration opportunities. Please contact C.A.B. ( cball@ancestry.com) for guidelines on submitting a research proposal.