UPDATE: With fixed Reference 3.
Let's start with some admixture analysis with our new reference 3 dataset. We already did K=2 admixture run.
Here's the results spreadsheet for K=3.
You can click on the legend to the right of the bar chart to sort by different ancestral components.
Fst divergences between estimated populations for K=3 in the form of an MDS plot.
And the numbers:
C1 C2
C2 0.114
C3 0.157 0.197
UPDATE: With fixed Reference 3.
Let's start with some admixture analysis with our new reference 3 dataset.
Here's the results spreadsheet for K=2 (i.e. two ancestral populations).
You can click on the legend to the right of the bar chart to sort by different ancestral components.
Fst divergences between estimated populations for K=2:
C1
C2 0.168
With the increase in the number of groups to more than 160, the bar chart has gotten too busy and is hard to figure out. Any suggestions to improve its readability? I have two ideas: One is to have a dropdown menu to select regions and thus have separate charts for each region. The other, which I like better, is to show only the top 50 groups for the ancestral component you are sorting by. Any other ideas are welcome.
I have fixed the problem with Reference 3 but if you notice any strange results, do let me know.
While the Reference 3 admixture results were generally good (and I have some nice surprises on the way I hope), the Reich et al populations had some weird behavior. From one K value to the next, their admixture would swing wildly especially among the minor components.
For example, for Chenchu, the 2nd component after South Asian was Southwest Asian (42%) at K=6, European (45%) at K=7 and American (32%) at K=8. That just didn't make any sense. It was similar for other Reich et al populations, but all the other reference populations seemed pretty stable.
The issue was that when I was creating Reference 3, I had to juggle lists of SNPs to figure out a way to include Reich et al with a large (>100,000) number of SNPs in the dataset since Reich doesn't have as many SNPs in common with the other datasets plus 23andme (v2 and v3) and FTDNA. In that effort where I was doing lots of SNP set intersections and unions I messed up. I used 217,000 SNPs. While these SNPs were present in all the other datasets, Reich et al had only 102,000 SNPs common with that set. Ouch! This was a royal mess as the high missing rate of Reich et al caused weird instability in its admixture results even though the rest of the results were mostly stable.
Now, I have pared down Reference 3 to 118,000 SNPs. These have a low missing rate in all the datasets. So I don't expect the same problems.
I am redoing the admixture runs with this new data and will have some of the results up soon.
Behar as in the Behar et al paper/dataset and not the Indian state of Bihar. The Behar dataset contains 4 samples of Paniya, which apparently is a Dravidian language of some Scheduled Tribes in Kerala.
I had always been suspicious of those four samples since one of them had admixture proportions similar to other South Indians but the other three were like Southeast Asians.
When I got the Austroasiatic dataset, I found out that they had the four Paniyas from Behar et al in their data. However, only one of those four was the same as Behar. The other three were different. So I now had 7 Paniya samples.
Let's look at the K=12 admixture results for these Paniyas.
Behar's GSM536916 was the one which was the same as Austroasiatic's D36 and it has regular South Indian results. The other three Behar Paniyas are very Southeast Asian (yellow in the plot) while the three Paniyas from Austroasiatic data are similar to GSM536916/D36.
Since the Austroasiatic Paniya samples originated from Behar et al, I guess at some point before the Behar data being submitted to the GEO database the Paniyas got mislabeled.
I am now excluding the four Paniyas from Behar et al dataset and only using the Paniya samples from Austroasiatic dataset.
Even though the Pan-Asian dataset is not public, there was a request for my script to convert the data to Plink's PED format.
Here is how I convert the Pan-Asian data to Plink's transposed file format.
#!/usr/bin/perl -w $file="Genotypes_All.txt"; open(INFILE,"<",$file); open(TFAM,">","panasian.tfam"); open(TPED,">","panasian.tped"); $line = <INFILE>; chomp $line; @first = split('\t',$line); foreach my $sample (5..$#first) { print TFAM "0 $first[$sample] 0 0 0 -9\n"; } my $alleles; while(<INFILE>) { chomp; @lines = split('\t',$_); my ($major,$minor) = split('/',$lines[4]); print TPED "$lines[2] $lines[1] 0 $lines[3]"; foreach my $snp (5..$#lines) { if ($lines[$snp] == 0) { $alleles = "$major $major";} elsif ($lines[$snp] == 1) { $alleles = "$major $minor";} elsif ($lines[$snp] == 2) { $alleles = "$minor $minor";} else { $alleles = "0 0";} print TPED " $alleles"; } print TPED "\n"; } close(INFILE); close(TFAM); close(TPED); |
Again, no guarantees! It's Perl though, so it should be more stable across various operating systems.
I have withdrawn the Admixture results for Reference 3 for now while I figure out why a few of them were weird and unstable.Далматин
I will report back on what I find and will have fixed results soon.
Following mallu's request, here is the code I used to convert Xing et al data to Plink's PED format.
#!/bin/bash dos2unix *.csv head --lines=1 JHS_Genotype.csv > header.txt awk -F, '{for (i=2;i<=NF;i++) print "0",$i,"0","0","0", "0"}' header.txt > xing.tfam sed '1d' JHS_Genotype.csv > genotype.csv sort -t',' -k 1b,1 genotype.csv > genotype_sorted.csv sort -t',' -k 1b,1 JHS_SNP.csv > snp_sorted.csv join -t',' -j 1 snp_sorted.csv genotype_sorted.csv > xing_compound.csv awk -F, '{printf("%s %s 0 %s ",substr($2,4),$1,$3); for (i=6;i<=NF;i++) printf("%s %s ",substr($i,1,1),substr($i,2,1)); printf("\n")}' xing_compound.csv > xing.tped plink --tfile xing --out xing --make-bed --missing-genotype N --output-missing-genotype 0 |
I make no guarantees that it will work for you. I used it on my Ubuntu box, but I am sure it'll have trouble on Mac OS.
Having collected 12 datasets, I have gone through them and finally selected the samples and SNPs I want to include in my new dataset, which I'll call Reference 3.
It has 3,889 individuals and 217,957 SNPs. Since this is a South Asia focused blog, there are a total of 558 South Asians in this reference set (compared to 398 in my Reference I).
You can see the number of SNPs of various datasets which are common to 23andme version 2, 23andme version 3 and FTDNA Family Finder (Illumina chip).
The following datasets had more than 280,000 SNPs common with all three platforms and hence were included in Reference 3:
Reich et al had about 100,000 SNPs in common with 23andme (v2 & v3 intersection) and 137,000 with FTDNA, but there was not a great overlap. Only 59,000 Reich et al SNPs were present in all three platforms. Since I really wanted Reich et al data in Reference 3, I included it but the SNPs used for FTDNA comparisons won't be the same as for the 23andme comparisons.
Of the datasets I could not include, I am most disappointed about the Pan-Asian dataset since it has a good coverage of South and Southeast Asia. Unfortunately, it has only 19,000 SNPs in common with 23andme v2 and 23,000 with 23andme v3. I am going to have to do some analyses with the Pan-Asian data but it just can't be included in my Reference 3.
I am also interested in doing some analysis with the Henn et al African data with about 52,000 SNPs for personal reasons.
Xing et al has about 71,000 SNPs in common with 23andme v3, so some good work could be done with that, though I'll have to use only 23andme version 3 participants.
The information about the populations included in Reference 3 is in a spreadsheet as usual.
As part of my effort to create one big reference dataset for my use, I have been going over all the datasets I have and make sure there's no duplicates or relatives or any other strange things that could cause issues with my analysis.
So I went back to the Reich et al Indian dataset.
The dataset doesn't have any duplicate or likely relative samples itself. However, there are two Kharia samples that are the same as the Austroasiatic dataset. Since Austroasiatic dataset has more SNPs in common with 23andme, I removed these two samples from Reich et al.
The IBS/IBS analysis and the sample IDs are in a spreadsheet as usual.
I computed the IBS similarity matrix for the Harappa participants HRP0001 to HRP0080 over 500,000 SNPs. This is exactly the same thing as the genome-wide gene comparison at 23andme.
Then, I converted the similarity matrix to a dissimilarity/distance matrix with the standard formula:
dij = sqrt(2 - 2 * sij)
where sij is the similarity between individuals i and j and dij is the distance/dissimilarity between the two.
Using the dissimilarity matrix, I classified all the participants (excluding close relatives) using hierarchical clustering with complete linkage. You can see the dendrogram below.
Then I used the same dissimilarity matrix to calculate 6-dimensional MDS. You can see the MDS plots below. The numbers on the plots are your Harappa IDs.
MDS Dimensions 1 & 2:
MDS Dimensions 3 & 4:
MDS Dimensions 5 & 6:
As you can see I (HRP0001) and my sister (HRP0035) are far away in the first four dimensions.
I'll let you guys speculate on what each dimension represents.
Now why create an MDS this way instead of directly using Plink's MDS functionality? Well, I needed to check if I could do it using only the similarity matrix because that would be really useful for something else. Tune in on my other blog for more later this week.
Recent Comments