SeqShop: Estimates of Genetic Ancestry Practical, June 2014
Introduction
See the tutorial slides for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.
The main purpose of this page is to provide step-by-step command lines for estimating ancestry of 6 targeted sequencing samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP).
For reference of more details about the options and usage of the software, please read the manual.
LASER workflow
Getting started
Create a working directory:
mkdir ancestry cd ancestry
Download and decompress software package:
wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz tar xzvf LASER-2.01.tar.gz
Set up to access data:
source /home/chaolong/LASER-Tutorial/setup.txt
What is in the setup.txt file:
export GC=/home/mktrost/seqshop/gotcloud export BAM=/home/chaolong/LASER-Tutorial/BAM export REF=/home/chaolong/LASER-Tutorial/reference export HGDP=/home/chaolong/LASER-Tutorial/HGDP
Preparing input files for LASER
Step 0: vcf --> geno
This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh".
Step 1: bam --> pileup
This step uses samtools to generate pileup files from bam files. Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps. It takes about 2 mins for each pileup job.
$GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101035.recal.bam > 121101035.recal.pileup & # $GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101043.recal.bam > 121101043.recal.pileup & # $GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101050.recal.bam > 121101050.recal.pileup & # $GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101052.recal.bam > 121101052.recal.pileup & # $GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101415.recal.bam > 121101415.recal.pileup & # $GC/bin/samtools mpileup -q 30 -Q 20 -f ./reference/hs37d5.fa.rz -l ./HGDP/HGDP_938.bed ./BAM/121101861.recal.bam > 121101861.recal.pileup &
Step 2: pileup --> seq
This step will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run. We used the pre-generated pileup files in the $BAM folder.
python ./LASER-2.01/pileup2seq/pileup2seq.py \ -m $HGDP/HGDP_938.site \ -b $BAM/AMD_roi_1-based.bed \ -i $BAM/AMD_hapmap_trios_id.txt \ -o hapmap_trios \ $BAM/121101035.recal.pileup \ $BAM/121101043.recal.pileup \ $BAM/121101050.recal.pileup \ $BAM/121101052.recal.pileup \ $BAM/121101415.recal.pileup \ $BAM/121101861.recal.pileup &
In the above command, -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID). -b and -i are optional.