1000 Genome Phase 3 Variant Set Mirrors

The 1000 Genomes consortium has released their final release of the Phase 3 variant set. This is a set of phased called variants in thousand genome individuals. Data is available in compressed .vcf.gz format along with the indexed tabix files. Looking for a way to download the data. See these FTP sites:

US Mirror:
ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/release/20130502/

UK Mirror:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/

The US 1000 Genomes mirror is substantially faster than the UK mirror when downloading files in the US.

How to Infer Ancestry from SNP Genotypes

Self-reported ancestry is poor metric to use when attempting to statistically adjust for the effects of ancestry since several individuals falsely report their ancestry or are simply unaware of their true ancestry.  Worse yet, sometimes you don't even have information collected on an individual's ancestry.  As many of you know, if you have SNP genotyping data you can rather precisely estimate the ancestry of an individual.  Classically, to do this one needed to combine genotypes from their study sample with genotypes from a reference panel (eg: HapMap or 1000 Genomes), find the intersection of SNPs in each dataset, and then run a clustering program to see which samples clustered with the reference ancestral populations.  Not a ton of work, but a minor annoyance at best.  Luckily, a relatively new program was just released that, in essence, does a lot of this ground work for you.  It is called SNPWEIGHTS and can be downloaded here.  Essentially, the program takes SNP genotypes as input, finds the intersection of the sample genotypes with reference genotypes, weights them based on pre-configured parameters to construct the first couple of principle components (aka. eigenvectors) and then calculates an individual's percentage ancestry for each ancestral population.  Here is how to run the program.

First, make sure you have Python installed on your system and that your genotyping data is in EIGENSTRAT format.  A brief tutorial to convert to EIGENSTRAT format using the convertf tool is here.

Next, download the SNPWEIGHTS package here and a reference panel. I usually use the European, West African and East Asian ancestral populations, but there are other options on the SNPWEIGHTS webpage as well.

Then, create a parameter file with directories of input files, your input population (designated "AA", "CO", and "EA"), and a output file.  An example of one is below:


Finally, run the program using the command inferancestry.py --par par.SNPWEIGHTS.  For the program to run correctly, make sure the inferancestry.info and snpwt.co files are in the same directory as your inferancestry.py file.

For more details, see the SNPWEIGHTS paper or the README file included in the SNPWEIGHTS zip folder.  For code on generating eigenvector plots with overlaid ancestry percentages and triangle plots according to percent ancestry, see this post.

How to get Phased Imputed Genotypes for Haplotype Analysis

I wanted to compare haplotypes from a set of cases genotyped on one platform to a set of reference controls genotyped on another platform for a small chromosomal region.  I looked into a few methods to do this and found that phasing the genotypes using SHAPEIT and then imputing off the 1000 Genomes reference panel using IMPUTE2 was the optimal approach to have a set of overlapping genotypes for a haplotype analysis.  Phasing from the .ped and .map files was easy using SHAPEIT.  The command I used with my .ped and .map files was:


Other file types can also be used as direct input to SHAPEIT (ex: .bed/.bim/.fam and .gen/.sample).  Genetic recombination maps can be downloaded here for build36 and build37.  Also the backslashes (\) are not necessary.  They just help organize the code by telling the computer to keep reading the input on the next line.

Once you have your .haps and .sample file generated, you are then ready to impute using IMPUTE2.  I used the command below:


Again the backslashes are just for a cleaner visualization of the code.  It is important when using the 1000 Genomes as your reference panel to have the SNP coordinates in hg19 coordinates.  LiftOver can help you convert from one build to another.  If you need the 1000 Genomes reference panel, it can be downloaded here.  The genetic map file is from the above SHAPEIT download.  Finally, if you want phased imputed results, the -phase command is essential to include.

You should now have your phased imputed genotypes to do whatever you wish with.  I chose to convert the phased imputed genotypes into a .fasta file with entries for each haplotype (see Python script).  Then I created parsimony trees with MEGA and visualized them in HapView.

Only Download Regional BAM File from Sequence Read Archive

Sequence Read Archive (SRA) stores raw next-generation sequencing data from a wide variety of sequencing platforms and projects.  It is the NIH's primary archive for high-throughput sequencing data containing well over 741,690,391,503,250 bases of publicly available sequence reads (that's a really big number).  Data is stored in one basic archived format, the SRA format, and can be downloaded for public use.  A toolkit is provided that supports conversion from .sra files to several popular data formats (see post for installation).  SRA recommends downloading .sra files using Aspera.

Not having sufficient space resources to download the entire .src files I would like regional .bam files for, I found two ways of downloading much smaller regional .bam files.

Web Based Interface
A somewhat cumbersome but practical solution is through one of SRA's web interfaces.  SRA offers a Run Browser which can be used to visualize data in the .src files.  After searching for the accession number and then choosing the reference chromosome and range, the run browser has the option of creating a .bam (as well as .fasta, .fastq, .sam, and .pileup) file for that region.  Just click to output the format to file and the download will begin in your web browser.

SRA Toolkit
This is a very useful set of programs to access data in .sra files and can be used to access data remotely (without having to download .sra files locally!) through the NCBI site.  It is easy to install and provides clear documentation.  An essential included module for extracting sequence read data is sam-dump.  A genomic region can be selected and then output to .bam format using Samtools.  Here is an example script:

I used these approaches to get .bam files I could use to investigate depth of coverage by Complete Genomics in a few regions of interest.  These methods allowed me to get complete (rather than evidenceOnly or evidenceSupport) .bam files for a few genomic regions for a subset of individuals that had been whole-genome sequenced by Complete Genomics in the 1000 Genomes Project.  To get their accession numbers needed for Run Browser and SRA Toolkit, extract the SRR number out of the index file.  Hope this is helpful for those wanting to use SRA, but not having loads of space to store the .sra files locally.

Download One Thousand Genome Data for Haploview

Haploview has a built-in portal to download HapMap data, but Haploview development hasn't kept pace with developing a way to automatically download 1000G SNP data.  Searching for a way to visualize the higher density SNP coverage of the 1000G project, I found it was not all too difficult to do this manually.  It involves a couple of extra steps.

First, determine the genomic coordinates of region you are interested in.  This needs to be in hg19 coordinates.  If you have hg18 coordinates, liftOver is a useful tool to convert coordinates from one human genome build to another (liftOver format is chr:start-end, for example: chr8:1000-50000).

Next, go to this 1000G website, and plug in your genomic coordinates of interest.  Here the coordinates should not include chr in the chromosome name (for example: 8:1000-50000).  Then on the next page select ancestral populations you are interested in (you can select multiple populations by holding down Ctrl).  Give the website a few seconds to generate the files.  Eventually a link to a marker information file (.info) and linkage pedigree file (.ped) will appear.  Right click on each of these files and save them to your computer.

Now, fire up Haploview and select Open new data.  Go to the Linkage Format tab and browse for your .ped file in the Data File field and your .info file in the Locus Information File field (the .info file field is usually automatically generated after selecting your .ped file if your .ped and .info files have the same prefix).  Haploview will load the files and you should be ready to visualize the LD structure.  Enjoy using 1000G data in Haploview!

Download Region of Indexed VCF File

Tabix is a handy tool to download a region of interest from a .vcf file on an ftp site.  Below is an example of how to extract a region on chromosome 1 from a Complete Genomics .vcf file hosted on the 1000 Genomes ftp site.  The .vcf.gz file must have a matching .vcf.gz.tbi index file for tabix to work and the chromosome region of interest should not include "chr" with the chromosome number.

Thousand Genomes Complete Genomics Information

Recently I have been using the Complete Genomics high coverage sequencing data from the 1000 Genomes project.  Its a bit difficult to find information on the populations, samples, and available sequencing data since they are all stored in different places on their ftp server.  I decided to make a post that tried to combine all the useful information into one spot.  Here are links to files that may be of interest.

population file: gives a key of abbreviations used for the 1000 Genomes populations

pedigree file: provides relationship information on 1000 Genomes individuals that are related as well as gender and 1000 Genomes population

sample file: detailed spreadsheet that offers information on sample id, accession number, population, family, gender, relationship, sequencing center, and coverage for each 1000G sample.

Samtools Download BAM Region Only

Often times publicly available sequencing data can serve as a useful reference for a sequencing project. The 1000 Genomes project is a great source, especially with their newly released high coverage Complete Genomics data.  Here is an example UNIX script that shows how BAM files with genomic regions of interest can be created from a whole-genome BAM file that is hosted on an FTP server, without having to download the entire BAM file first.  The bai_file_list.txt file is a file that contains unique identifiers for each BAM file extracted from a previously selected list of BAI files of interest.  Here I am just extracting the BRCA1 and BRCA2 regions of the genome.  The extracted reads are sorted and then saved as a BAM file and an associated BAI index file is also created.  The final step removes excess BAI files that are downloaded and used by Samtools to extract the region of interest from the BAM files on the FTP server.

Thoudand Genomes Complete Genomics Index

As mentioned previously, 1000 Genomes has now made Complete Genomics whole genome sequencing publicly available for download.  They give an index file that mentions some of the individual high coverage .bam files that are available for download, but it seems to be missing a lot of the newer data released this April.  I was trying to find a way to efficiently search through the 1000 Genomes ftp site to get a better index of the available Complete Genomics data.  I am primarily interested in CEU samples and wanted the high coverage evidence support files.  Here is the code I used to search through the ftp. I chose just to search for and download .bai files from the site since they are quick to download and would create a useful index for downloading the bam files.

1000 Genomes Complete Genomics Data

1000 Genomes now has Complete Genomics high coverage sequencing data available for 57 samples online at the ftp sites: ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data and ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data.  The first ftp is best for US access and the second is best for access from Europe.  Here's a link to the 1000 Genomes page for more details: http://www.1000genomes.org/announcements/complete-genomics-data-available-2012-12-19.