This is a question from /u/beneficii9 on reddit. The original post can be found here.

Through the Personal Genome Project, I have had my whole genome sequenced by Veritas, and have it in the form of a single VCF file for the whole genome and one BAS file for each chromosome. The reference genome associated with the VCF file is hg19. It has been helpful in health data; for example, I discovered I'm homozygous for the non-functional variant CYP-2D6 gene (rs3892097), which can render several common medications useless, and helps explain why some medicines didn't really work for me. My doctor has found this information very helpful.

Unfortunately, I can't find any way of looking at admixture or ancestry. I've tried setting everything up using a combination of VCFTools, Plink1.9, and ADMIXTURE, but I can't get it to work. I think for ADMIXTURE you have to have a bunch of genomes sorted by geographical origin to compare your genome against, but I'm not sure how to do that, and what's online isn't very clear to me. So scratch that one off.

I've tried converting the file to 23andme format (and at this /u/psychosomaticism has been very helpful). I did that (though it seems there were problems because of the way the VCF file was set up). But the websites that take the data want you to point them to your 23andme account, and that doesn't really work if you only have the file. 23andme doesn't provide for people who had their whole genomes sequenced. They want you to give them a saliva sample like everyone else.

So, what can I do?


A modified implementation of Vivek's answer.

peddy is a Python package that samples an input .vcf at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.

I downloaded the .vcf and .vcf.tbi for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom .ped file, NA12878.ped, with the contents:

NA12878 HG001 0 0 2 0

At the command line:

$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped

The output files all have the prefix myvcf., here's myvcf.pca_check.png myvcf.pca_check.png

| improve this answer | |
  • $\begingroup$ I think it might be worth adding: peddy requires coloredlogs and cyvcf2 the later of which will required you to have gcc&g++ and the headers (e.g. -dev versions in ubuntu) of libbz2, liblzma and curl. This knowledge would have saved a precious 15 precious moments of my life :) $\endgroup$ – George Feb 9 at 0:32
  • $\begingroup$ Also, now that it's installed I seem to get the error: AssertionError: error loading tabix index for b'HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz' (Also based on my understanding shouldn't the 2 be 0 or left to the user to set based on his gender ?) $\endgroup$ – George Feb 9 at 0:41

A suggestion from the UK Biobank QC methods white paper:

  1. Create a set of SNPs common to your VCF and the 1000 genomes phase 3 call set.
  2. Perform PCA of the 1000 genomes samples using eigenstrat smartpca. You might have to convert to binary plink format.
  3. Project your genotypes onto this pre-computed PCA space and visualize using ggplot to see which cluster you fall in.
| improve this answer | |

Ancestry testing is a tricky subject. I spent a good chunk of my PhD project on questions similar to this, and haven't really found a good answer for how to detect unmodelled ancestry.

The general idea of how I tried to approach ancestry determination was to create a model set containing well-defined groups of individuals with a known, specific ancestral background. A query individual (or individuals) were then added to this group, and an ancestry-estimating program was run (specifically Structure) to work out what proportion of the unknown individual was attributable to each of the known groups. I would expect that most genetic ancestry tests follow a similar approach to this, although perhaps with a little bit less care about the definition of the model population groups.

Problems arise when individuals are incorrectly assigned to a particular group, when groups are present in the model set that represent a larger proportion of individuals than most other groups, when groups are closely-related to other groups, and when an ancestral history is present in a test individual that doesn't match any of the model groups. And all that assumes that the marker set used for ancestry determination is perfect: no bias towards any particular group, and no systematic genotyping error.

This doesn't mean that ancestry testing won't work, but it's a good idea to take the results with a large grain of salt. There's a good example of a media personality in New Zealand who was told that she had a very high probability of being 100% Māori, despite having a good knowledge of her own family history that indicated one European ancestor a few generations back on both sides of her family.

| improve this answer | |

Adding another answer that, while much more involved, offers flexibility, walks through minimal steps, and provides visualizations.
I created a Dash app and repository that explores this type of analysis using published ancestry informative SNPs and 1000 Genomes Project Data.

  1. Identify published ancestry informative SNPs (AISNPs).
  2. Download 1000 Genomes Project genotype bcf data.
  3. Limit genotype data from (2) to AISNPs loci from (1).
  4. One-hot encode the genotypes.
  5. Perform dimensionality reduction (PCA, t-SNE, or UMAP)
  6. Plot components.

The original post seemed to want to project their genome against a reference sample. The plot_walkthrough.ipynb in the tgviz repo shows how this can be achieved.
Output from the notebook on my own data:

enter image description here

| improve this answer | |

Related: ethnicity check either from bam or vcf files

The missing piece of your puzzle is an ancestry-labeled reference dataset. 1000 Genomes phase 3 usually works well enough here, and it now has variant callsets for both GRCh37 and GRCh38.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.