The component maps of MDLP-World22 calculator

I would like to express my gratitude to the fellow members of ABF - Loxias and Wojewoda- for creating amazing maps of components and plots showing interrelationships between 22 components.

Since the readers of my blog might be interested in visual inspecting of components , i have decided to upload them on-line. 

The first batch of maps  includes "the component portraits", created by Loxias:

West-Asian

Atlantic-Mediterranean

East-Siberian

Amerind

Indian

Indo-Iranian

Indo-Tibetan

Paleo-Siberian

Pygmy

Sub-Saharan

Near-East

North-East-European

North-Siberian

North-Mesolithic-European

 The second batch includes the PCA scatter plot of components, created by Wojewoda. PCA scatter plots one PCA component vs. another PCA component of data obtained from ancestry coefficients  collected for each of MDLP World22 components: the spots are connected over time to produce a trajectory for each component.

 

Paint Me a Rainbow: Painting World 22 ancestral components

This update will be concerned with inter-related concepts of  "chromosome painting" and admixture. Modern genetics and personal genomics, particularly in the last 2-5 years, had paid very considerable attention to them, sometimes under the rubric of " determining ancestral origin of genomic segments".

In 2011 year, i performed a couple of experiments with LAMP software (cf. Experimental test III: (De)constructing ancestry with LAMP  and Experimental test: estimating ancestries at each locus in a population of admixed individuals (LAMP))).

Although the experiments met with moderate success, i wasn't satisfied with the results and decided to postpone the forthcoming experiments with chromosome painting to a future day. In so doing, i endorsed increasing  appreciation of the distinction between  population stratification's algorithms, implemented in LAMP and ADMIXTURE.

I have already discussed the differences between LAMP and Admixture, but in illustrating the idea of experiment, i need to turn to my previous explanation again: 

1) ADMIXTURE software  implements a model-based approach to estimate ancestry coefficients as the parameters of a statistical model. It is also important to add that the model-based approach in ADMIXTURE is based  on the global ancestry paradigm (i.e the goal of  ADMIXTURE/STRUCTURE analysis is to estimate the proportion of ancestry from each contributing population, considered as an average over the individual's entire genome).

2) LAMP software is built upon an efficient dynamic-programming algorithm WINPOP that infers locus-specific ancestries.Genome is partitioned into chromosome segments of definite ancestral origin (overlapping, contiguous windows of SNPs) and likelihood model optimized over each window. The goal then is to find the segment boundaries and assign each segment's origin.

I understand the problem in terms of the ancestry assignment. My experience shows that methods based on  the inference of locus-specific ancestries are usually very accurate for ancestral deconvolution of genotype data that has consistently been shown to do better than more popular statistical and PCA-based methods, while being able:
1) to handle more than two ancestral populations
2) to model the paths of recombination between ancestral segments.In June 2012, Jason Mezey Lab (Cornell University) released SupportMix - a machine learning algorithm for determining ancestral origin of genomic segments when analyzing individuals from a population with a recent or ancient history of admixture. As regards the accuracy of the software, the authors argued that SupportMix provides a robust tool for accurate and robust ancestral assignment by simultaneous analysis of a worldwide selection of ancestral populations. Such analyses will be critical for accurate assignment in the many world-wide admixed populations that are likely to have unexpected ancestry that reflects a richer history than known from anthropological or historical studies (from the provisional paper: Omberg et al.2012 "Inferring genome-wide patterns of admixture in Qataris using fifty-five ancestral populations"). The cited paper includes a number of other claims important for any analysis of genetic admixture: the accuracy of ancestry assignment was lower for more closely related populations but better than LAMP-ANC, a method that was been shown to consistently outperform other ancestry deconvolution methods.

This overly optimistic conclusion influenced my choice between LAMP-ANC and SupportMix in favor of latter. To be honest, i'm not the first genome blogger to use Supportmix - in July of 2012 Polako from Eurogenes carried out a loci-specific analysis of Finnish genomes using SupprotMix. I decided to repeat this experiment. There is, of course, a significant  difference between Polako's analaysis and my experiment. While Polako was using the modern populations as 'putative' donor populations, the final goal of my design was to imitate the results of the long-awaited update of 23andme's Ancestry Painting, which is is being updated to offer more detailed results based on approximately 20 world regions, drawn from both customer data and academic reference populations. In order to do so, i have used the dummy set of 22 simulated putative ancestral populations simulated from the allele frequencies of the World-22 calculator.

The experiment

SupportMix requires at least three input files. One file for each of the putative ancestral populations and one file containing the genetic information of the admixed individuals with additional requirement that he markers have to be phased. Each population should be represented by two files in Plink transposed format, a .tped file and a .tfam file.

The markers (80751 SNPs) were phased per each chromosome using default settings in BEAGLE software. While the original Plink format does not specify the order of the alleles in in the file, SupportMix works with phased data. Keeping that in mind, i have converted BEAGLE-phased dataset directly into Plink's tped format without pre-processing the dataset in Plink (hint: Kantele's beagle_to_tped script). Then i used UNIX text processing utilities to extract the genotypes from ancestral populations into the corresponding subsets ('references') and split the project dataset (93 individuals) into 9 subgroups.

Finally, i have interpolated the genetic map position of each SNP along chromosome using Rutgers genetic maps.
 SupportMix was run by specifying a configuration file with the default options: (window_size =400, generations_from_admixture_event=6).

Below are some screenshots of  SupportMiX output for the MDLP project participants:

Please note that each 'recipient' (i.e the project participant) is represented by two phased chromosomes, i.e V199_a and V199_b. The color legend of the components used in the analysis, has been attached to the right side of the plot.

The complete set (Chr.1-22 for all project participants) in tar.gz format (14.6 Mb) could be downloaded here.


SupportMix results: what to do next.

First of all, i encourage every participant of my project to compare their results to "chromosome painting" in MDLP World-22 calculator on John Olson's Gedmatch site. I have to contemplate the possibility that the painting on Gedmatch site might be different from that one produced by SupportMix. I would also suggest to compare SupportMix's paintings to other calculators' paintings and 23andme's Ancestry Finder, etc.

If you are familiar with basic techniques of image editing software, then it is a good idea to have your chromosomes cut&merged into the composite image (see example below):

Chromosome Painting (in 23andme's  style)

Chromosome Painting (an imitation of 23andme's Ancestry Finder)

 


 

Behind the Curtains: MDLP World 22 showcase

Preliminary remarks

As you all may know, the MDLP  blog hasn't been updated since February 2012.

Half of year ago i promised myself that i would stop writing new posts on the MDLP blog before i'll finally get my scientific report on blog written.  Since I had to prioritize the completion of a scientific paper over the routine of blog posting,I was unable to continue updating the blog on a regular basis due to a lack of time, and had to make a change in how I conducted my research. So i decided to abstain from posting on the MDLP blog for a couple of month, being focused on more important matters. Despite of all limitations, i kept  working secretly on the MDLP project, collecting necessary data and performing different 'genomic' experiments in order to achieve my final goal (publishing of paper).  The results of secret experiments with new genomic samples and tools eventually leaked to the curious public, spawning immense interest in my project. After releasing a new version of my own modification of DIYDodecad calculator on Gedmatch.com, i was literally flooded by emails from Gedmatch.com users asking me  questions they wanted me to answer.

I understood the strategical mistake of releasing poorly documented data/analysis on Internet  and felt obliged  to explain details. Obviously, i will start new series of  the blog posts by covering the project feature the people most interested in, i.e the MDLP World22 calculator.
 

The population dataset of MDLP World22 calculator.

The reference population dataset of the calculator was assembled in PLINK by intersecting and thinning the samples from different data sources: HapMap 3 (the filtered dataset CEU,YRI,JPT,CHB), 1000genomes, Rasmussen et al. (2010),  HGDP (Stanford) (all populations),  Metspalu et al. (2011),Yunusbayev et al. (2011), Chaubey et al. (2010) etc. Furthermore i handpicked random 10 individuals from each European country panel in POPRES dataset, or the maximum number of individuals available otherwise, to select the POPRES European individuals to be included in our study. Finally, in order to evaluate the correlation between the modern and the ancient genetic diversity, i have also included ancient DNA genomic samples of Ötzi,(Keller et al.(2012)) Swedish Neolithic samples Gök4, Ajv52, Ajv70, Ire8, Ste7 (Skoglund et al. (2012)) and 2 La Braña individuals from the Mesolithic sites of the Iberian Peninsula (Sánchez-Quinto et al.(2012)). Then i added 90 samples of individuals-participants of our MDLP project.  After merging the aforementioned datasets and thinning the SNP set with PLINK command to exclude SNPs with missing rates greater than 1% and minor alleles, i filtered out duplicates, the individuals with high pairwise IBD-sharing  (estimated in Plink as as the average fraction of alleles shared between two individuals over all loci) and the individuals with kinship coefficient suggesting relatedness (kinship coefficients were estimated in KING software). Also i had to filter out individuals with more more tham 3 standard deviations from the population averages. Since kinship coefficient is robustly estimated by HWE (Hary-Weinberg expectations) among SNPs with the same underlying allele frequencies, SNPs showing strong deviation (p < 5.5 x10⁻⁸) from Hardy-Weinberg expectations were removed from the merged and filtered dataset. After that I filtered to keep the list of common SNPs present in Illumina/Affymetrix chips and performed  linkage disequilibrium based pruning using a window size of 50, a step of 5 and r^2 threshold of 0.3.

This complex sequence of consequent operations with the initial reference and project datasets yielded a final dataset which included 80751 SNPs  in 2516 individuals from 225 populations.


ADMIXTURE analysis

 As always, the final dataset in PLINK linked format was further processed in ADMIXTURE software. Sketching the plan for the design of ADMIXTURE test,  i had to face the difficult problem: as it has been shown in (Patterson et al.2006) the number of markers needed to resolve populations in ADMIXTURE analysis is inversely proportional to the genetic distance (Fst ) betweeen the populations. According to ADMIXTURE best practice, it is believed that 10,000 markers  are suffice to perform GWAS correction for continentally separated populations (for example, African, Asian, and European populations FST > .05) while more like 100,000 markers are necessary when the populations are within a continent (Europe, for instance, FST < 0.01).

To increase the accuracy of ADMIXTURE results i decided to use a method proposed by  Dienekes' for converting allele frequencies into 'synthetic individuals'(see also Zack's example). The idea is fairly simple: run an unsupervised ADMIXTURE analysis once to generate allele frequencies for your K ancestral components; then generate zombie populations using these allele frequencies; whenever you want to estimate admixture proportions in new samples run supervised ADMIXTURE analysis using the zombie populations.  Like any genome blogger engaged in the task of evaluating admixtures in samples, i must grapple with obvious question of the reliability of this approach.  Although i am aware of  methodological controversies in using simulated individuals, i would rather concur with Dienekes who considered "synthetic individuals" the best abstract proxies for the ancient ancestral populations. But my purpose is served if i can use the approach used by Dienekes and Zack to obtain meaningful results. To begin with, i routinely ran unsupervised ADMIXTURE  K=22 analysis (assuming 22 ancestral populations) which yielded the admixture proportions of individuals from these K populations, as well as the allele frequencies for all SNPs for each of 22 ancestral populations (below are conventional names for each of inferred components in order of appearance):

Pygmy
West-Asian
North-European-Mesolithic
Tibetan
Mesomerican
Arctic-Amerind
South-America_Amerind
Indian
North-Siberean
Atlantic_Mediterranean_Neolithic
Samoedic
Proto-Indo-Iranian
East-Siberean
North-East-European
South-African
North-Amerind
Sub-Saharian
East-South-Asian
Near_East
Melanesian
Paleo-Siberean
Austronesian

Therefore i took the allele frequencies which were computed earlier in unsupervised Admixture K=22 for the merged dataset, pooled them into PLINK and generated 10 "synthetic individuals per ancestral component) using PLINK command --simulate.  When the simulation had been  finished, i visualized the distance between simulated individuals using multi-dimensional scaling:



 

As a next step,i included simulated individuals you have as part of a new reference population (including 220 simulated individuals in 22 simulated populations).Then, I ran ADMIXTURE anew, this time in “supervised” mode for K = 22 (with simulated individuals being 'reference' individuals). The Admixture K=22 converged in 31 iterations (37773.1 sec) with final loglikelihood:-188032005.430318 (below are Fst divergences between estimated 'ancestral' populations):

The Fst distance/divergence matrix was used for inferring a most probable NJ-based topology of component distance tree (outgroup: South-African):

 The individual 'supervised' ADMIXTURE  results (in Excel spreadsheet) for the project participants have been uploaded to GoogleDocs (please note that the average results for reference populations is also available on special request).

MDLP World22 DIYcalculator

The output files of Admixture K=22 supervised run (average values of admixture coefficients in reference populations and FsT values)  were used for designing a new version of  the MDLP DIYcalculator, which is better known by its codename "World22" (online version is available in AdMix-Utilities section of Gedmatch under MDLP project). MDLP DIYcalculator itself is based on the code of Dodecad DIY calculator (c)ourtesy of Dienekes Pontikos and was developed as part of the Dodecad Ancestry Project. In its Gedmatch implementation MDLP 'World22' DIYcalculator is paired by MDLP 'World22' Oracle, also based on Dienekes' and Zack's code (Harappa/DodecadOracle). The 'Oracle' is designed to find in a single population mode your closest (closest in terms of similarity) population from MDLP ''Word22' admixture results. In a mixed mode, Oracle considers all pairs of populations, and for each one of them calculates the minimum Fst-weighted distance to the sample in consideration, and the admixture proportions that produce it.

Please notice: 'ancestral' populations (i.e 'simulated populations' from the previous step - see above) are labeled in Oracle results as (anc), while the 'real world' modern and ancient populations are marked as "derived".

If you have troubles with understanding/interpreting the results of Oracle and DIYcalculcator, please consult the corresponding topics on Dodecad and HarappaWorld blogs. It is not of avail to repeat in this blog everything they wrote in their own blogs.

What the heck are MDLP Word-22 components?

 One of those questions that i usually keep getting in emails is what do the various reference populations and ancestral components for my World K=12 and World-22 analyses mean. I've already provided hints to the answer  earlier, but - as old Chinese proverb says - one picture is worth ten thousand words. That's why i decided to display the admixture coefficients spatially on the globe surface. Following Francois Olivier, who proposed to use the graphical library of the statistical software R to display  spatial interpolates of the admixture coefficients (Q matrix) in two dimensions (where spatial coordinates are recorded as longitude and latitude), i created  2 contour maps per component.

Pygmy (modal in Biaka and Mbuti population)

West-Asian (bimodal component with peaks in Caucasian populations and south-western part of Iran, equal to Dienekes' Caucasian/Gedrosia component) 

 North-European-Mesolithic (local component with peaks in European Mesolithic samples of La_Brana and  modern North-European Saami population).

 Tibetan (Indo-Burmese) component (Himalay, Tibet)

Mesomerican (major genetic component in Native Americans from Mesoamerica)

North-Amerind (the 'native' component in North American Natives)

South-Amerind (the 'native' component in South American Natives)

  Atlantic-Mediterranean-Neolithic (the main genetic component  in Western and South-Western Europe)

  The rest contour maps for all components could be downloaded here.