1
$\begingroup$

Posted a similar question on Biostars but got no response. Not sure if I'm allowed to link to it? Basically I want to pull genotype frequency data for a population group (such as CEU) instead of allele frequency data, via the Perl API for 1000 genomes. I have tabix and perl API installed. This is for 100,000+ SNPs so the solution should hopefully not involve manually downloading genotype and calculating the frequency manually in a for-loop. My understanding also is no solution exists in BioMart (according to the Biostars answer).

Following the instructions here, I can see how to get specific genotypes for all individuals as a list.

Eg: Input -> Some snp, CEU... Output -> G/G 0.87, G/A 0.13 (frequency of some snp for CEU population).

I need to do this for 100,000+ SNPs, so I imagine manually pulling all genotypes for each SNP and calculating the frequency manually in a loop would not be practical.

Improve this question
$\endgroup$
6
  • $\begingroup$ Please add a link to the data repository and show a dummy example of the data format. Questions about a specific programming language are considered off-topic. So you should edit the question to make it "Language-neutral", focussing only on the methodology rather than the actual script. $\endgroup$
    – WYSIWYG
    Aug 17, 2015 at 7:03
  • $\begingroup$ Hi, I think Perl is the only way to access the 1000 genomes data genotype frequency, according to this post: biostars.org/p/153923/#153929 Unfortunately I got no further answer there on how to actually do this. $\endgroup$
    – user4779
    Aug 17, 2015 at 8:41
  • $\begingroup$ Can you add a link to the actual repository? You may be able to simply download the entire data (from here perhaps) and parse for the desired data. Moreover, the API mentioned in your link would allow you to customize the data that you want to download. You may have to do it one by one using a loop. $\endgroup$
    – WYSIWYG
    Aug 17, 2015 at 8:46
  • $\begingroup$ I'm new to 1000 genomes so not sure how to access the data directly, I've just been following the perl API tutorial for homo sapiens as per the link. That link you provided seems helpful but it looks like I'll need to manually iterate through the hundreds of genomes and calculate the frequency in the loop as you suggest, rather than there being an actual functional call like there is for allele frequency? I think if that's the case I'll just stick with the HapMap genotype frequency data. $\endgroup$
    – user4779
    Aug 17, 2015 at 9:05
  • $\begingroup$ You don't have to manually do it. Get the VCF or GVF file. A small script will parse the desired co-ordinates in a few minutes (or less than that). $\endgroup$
    – WYSIWYG
    Aug 17, 2015 at 9:13

2 Answers 2

Reset to default
3
$\begingroup$

If you want population specific allele frequencies you have three options: * For a single variant you can look at the population genetics page for a variant in our browser. This gives you piecharts and a table for a single site. * For a genomic region you can use our allele frequency calculator tool which gives a set of allele frequencies for selected populations * If you would like sub population allele frequences for a whole file, you are best to use the vcftools command line tool.

This is done using a combination of two vcftools commands called vcf-subset and fill-an-ac

An example command set using files from our phase 1 release would look like

grep CEU integrated_call_samples.20101123.ALL.panel | cut -f1 > CEU.samples.list

vcf-subset -c CEU.samples.list ALL.chr13.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz | fill-an-ac | bgzip -c > CEU.chr13.phase1.vcf.gz

Once you have this file you can calculate your frequency by dividing AN (allele number) by AC (allele count)

Please note that some early VCF files from the main project used LD information and other variables to help estimate the allele frequency. This means in these files the AF does not always equal AC/AN. In the phase 1 and phase 3 releases, AC/AN should always match the allele frequency quoted.

Improve this answer
$\endgroup$
2
$\begingroup$

Another possibility is to use glactools. I will show this by using data from the 1000 genomes.

First, we download the chromosome names and length for the reference:

 wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai

We then require information on the panels from the VCF file:

  wget https://personal.broadinstitute.org/armartin/ginger/integrated_call_samples_v3.20130502.ALL.panel.txt

 grep -v ^sample integrated_call_samples_v3.20130502.ALL.panel.txt  | cut -f 1,3  > panel.txt

You can run the following to transform VCF into allele counts format (ACF):

 tabix -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr2.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz 2:136486829-136653337  |glactools vcfm2acf --onlyGT --fai human_g1k_v37.fasta.fai  -  |glactools meld -f panel.txt -  | glactools view -h -

This will print the allele counts:

#chr    coord   REF,ALT root    anc AFR AMR EAS EUR SAS
2   136486850   G,T 0,0:0   0,0:0   1322,0:0    694,0:0 1008,0:0    1006,0:0    974,4:0
2   136486967   C,T 0,0:0   0,0:0   1321,1:0    694,0:0 1007,1:0    1006,0:0    978,0:0
2   136487007   C,T 0,0:0   0,0:0   1322,0:0    694,0:0 1007,1:0    1006,0:0    978,0:0
2   136487181   C,T 0,0:0   0,0:0   1322,0:0    693,1:0 1008,0:0    1005,1:0    978,0:0
2   136487214   G,A 0,0:0   0,0:0   1321,1:0    694,0:0 1008,0:0    1006,0:0    978,0:0
2   136487246   G,A 0,0:0   0,0:0   1282,40:0   693,1:0 1008,0:0    1006,0:0    978,0:0
2   136487336   G,T 0,0:0   0,0:0   1322,0:0    694,0:0 1008,0:0    1006,0:0    977,1:0
2   136487417   G,A 0,0:0   0,0:0   1321,1:0    693,1:0 1008,0:0    1006,0:0    978,0:0
2   136487504   A,C 0,0:0   0,0:0   1316,6:0    694,0:0 1008,0:0    1006,0:0    978,0:0
Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .