As the doctor gone rogue

April 6, 2012

Convert plink format to vcf

Filed under: data management, plink, plinkseq, vcf, vcftools — hypotheses @ 12:20 pm

So, I recently got to work with vcf file format for the first time. It is been quite a mess to start with, but after all the helps I can gather around here (Thank you very much, YK & HS !!!) I finally came down to a solution.

You probably have heard about “plink”, “plink/seq”, and “vcftools”. Among the three programs, “plink” has been around the longest, and one of the reason may be because it is versatile, has good documentation, and quite easy to use. Although plink/seq was inspired by plink, its documentation is still during development (I tried to be optimistic here).

For this short tutorial, I am going to use PLINK/BED file as a medium for file conversion, since it is the fastest file to manipulate and the most efficient to use in terms of disc space and memory requirement (from my experience).

“So, my question is how do I convert my GWAS data to vcf format, with a specific reference allele that I want to use i.e. not the default ‘major allele’ as a reference as we typically do”

So, I now have PLINK/Binary format file call csOmni25 and a file containing reference allele csOmni25.refAllele

As a short note: csOmni25.refAllele looks like this


-------- csOmni25.refAllele ----------------

rs12345 A

rs12958 G

rs10596 C

rs18607 T

...

...

---------end csOmni25.refAllele--------------

  1. First, we will convert PLINK/Binary format file so that A1 [reference allele] correspond to the reference allele that we want
    NOTE: when you create the reference allele file, make sure that all reference alleles are in UPPER CASE.
  2. Second, we will import plink/bed to plink/seq and write out vcf format file

Here’s the script that I use. It’s that simple. And for your reference, it took me a week to figure this out, and tested it.


#!/bin/sh

##-- SCRIPT PARAMETER TO MODIFY--##

PLINKFILE=csOmni25

REF_ALLELE_FILE=csOmni25.refAllele

NEWPLINKFILE=csOmni25Ref

PLINKSEQ_PROJECT=csGWAS

## ------END SCRIPT PARAMETER------ ##

#1. convert plink/binary to have the specify reference allele

plink --noweb --bfile $PLINKFILE --reference-allele $REF_ALLELE_FILE --make-bed --out $NEWPLINKFILE

#2. create plink/seq project

pseq $PLINKSEQ_PROJECT new-project

#3. load plink file into plink/seq

pseq $PLINKSEQ_PROJECT load-plink --file $NEWPLINKFILE --id $NEWPLINKFILE

#4. write out vcf file, as of today 4/6/2012  using vcftools version 0.1.8, although the documentation says that you can write out a compressed vcf format using --format BGZF option, vcftools doesn't recognize what this option is. So, I invented my own solution

pseq $PLINKSEQ_PROJECT write-vcf | gzip > $NEWPLINKFILE.vcf.gz

At the end, this will create a compressed vcf file “csOmni25Ref.vcf.gz” with the specified reference alleles.

Advertisements

April 22, 2011

Split Whole-Genome PLINK Binary Files to Small Chunks

Filed under: bash, data management, genetics, plink, Uncategorized — hypotheses @ 1:03 pm

I prefer to store my GWAS data in PLINK binary format. The PLINK binary file format gives us many advantages. With tons of data, binary file format is a very efficient way to save disk space. To give you some number, you can store the data of 2000 people from Illumina Human 1M chip with ~1 million marker in 1 CD (~700Mb). I prefer to keep all the genetic markers in the same file. This is convenient when you want to extract a list of SNPs from several chromosome. You can just go to this one file and extract them out.
However, keeping everything in one single file is not always a good solution to all the problems. For most of the genome-wide association analysis, we can speed up the analysis using parallelization. The simplest way to do this is to split your data into several smaller chunks and submit each part for analysis on several computers (preferably computing cluster). Although this is so-called "poor man" parallelization, it gets the job done in no time. 
For this purpose, from a single PLINK binary file format, I wrote a script to split this file by chromosome (in by_chr directory), and for each chromosome into a smaller chunk of your choice (I used 5000 for the example below).

#!/bin/bash
## generate PLINK script to split whole genome PLINK BED TO smaller chunks
set -e
#set -o verbose
split_by_chromosome=0
prefix=$1
chunk=$2
postfix=$3
if [ ! `type -p plink` ];then
echo "Cannot find PLINK. Make sure PLINK is in your \$PATH";
exit 1
fi
if [ "${prefix}" == "" ]; then
echo "Usage: ./splitPlink.sh [prefix_plink_input] [chunk_size] [prefix_plink_output]"
echo "Example: $sh splitPlink.sh /research/imp_data/shared/washuDataReformatted/blackR21 10000 blackSR21"

echo "Exit no input file"
exit 1
fi
if [ "${chunk}" == "" ]; then
echo "Usage: ./splitPlink.sh [prefix_plink_input] [chunk_size] [prefix_plink_output]"
echo "Example: $sh splitPlink.sh /research/imp_data/shared/washuDataReformatted/blackR21 10000 blackSR21"
echo "Exit no chunk size"
exit 1
fi
if [ "${postfix}" == "" ]; then
echo "Usage: ./splitPlink.sh [prefix_plink_input] [chunk_size] [prefix_plink_output]"
echo "Example: $sh splitPlink.sh /research/imp_data/shared/washuDataReformatted/blackR21 10000 blackSR21"
echo "Need prefix of output"
exit 1
fi
echo "File ${prefix} by ${chunk} SNPs chunk"
echo "Spliting PLINK ${prefix} into ${chunk} file size"
if [ ! -e ${prefix}.bim ];then
echo "File ${prefix}.bim not found"
exit 1
fi
chrs=`awk '{print $1}' ${prefix}.bim | sort -nur`
echo "Chromosome in Map File: ${chrs}" | tr "\n" " "
echo ""

## 0. Create script directory if not exist
if [ ! -d by_chr ]; then mkdir -p by_chr; fi
if [ ! -d script ]; then mkdir -p script; fi

# 1. Split by chromosome first
for chr in ${chrs}; do
cat < script/plinkSplit.chr${chr}.sh
#!/bin/sh
# extracting chromosome $chr from ${prefix}
plink --noweb --bfile ${prefix} --chr ${chr} --make-bed --out by_chr/${postfix}_chr${chr}
EOF
if [ $split_by_chromosome -eq 1 ]; then
echo "Processing map of chromosome $chr"
grep -w "^$chr" ${prefix}.bim > by_chr/${postfix}_chr${chr}.bim
fi
sh script/plinkSplit.chr${chr}.sh
done

## 2. Split chromosome to a chunk of ${chunk} SNPs
for chr in $chrs;
do
if [ ! -d c${chr} ]; then mkdir -p c${chr}; fi
mapfile=by_chr/${postfix}_chr${chr}.bim
## wait till finish extracting bim by chromosome to submit the following jobs
## default wait_duration = 15 minutes (change value below for longer wait time
wait_duration=5 ## Wait for 1 hour before stopping
count=1
while [ "$count" -le "${wait_duration}" ];
do
if [ ! -e by_chr/${postfix}_chr${chr}.bed ]; then
echo -n -e "\r                                                   \r$count. Waiting for file by_chr/${postfix}_chr${chr}.bed"
let "count += 1"
sleep 60
else
let "count = wait_duration + 1"
n_end=`wc -l by_chr/${postfix}_chr${chr}.bim | awk '{print $1}'`
echo "Chromosome $chr has $n_end SNPs"
snp_start=1
let "snp_end = snp_start + ${chunk} -1"
if [ $snp_end -gt ${n_end} ];then snp_end=${n_end};fi
pc=1
echo "Extracting chromosome $chr from $snp_start to ${snp_end}"
while [ "$snp_end" -le "${n_end}" ];
do
first_snp=`tail -n +${snp_start} ${mapfile} | head -n 1| awk '{print $2}'`
last_snp=`tail -n +${snp_end} ${mapfile} | head -n 1 | awk '{print $2}'`
echo "First SNP $first_snp Last SNP $last_snp"
cat <<EOF > script/plinkSplit.chr${chr}.pc${pc}.sh
#!/bin/sh
echo "Extracting chromosome $chr from $snp_start to ${snp_end}"
plink --noweb --bfile by_chr/${postfix}_chr${chr} \
--from $first_snp --to $last_snp --make-bed \
--out c${chr}/${postfix}_chr${chr}_p${pc}
exit 0
EOF
sh script/plinkSplit.chr${chr}.pc${pc}.sh
let "pc += 1"
let "snp_start += ${chunk}"
let "snp_end += ${chunk}-1"
if [ $snp_end -gt ${n_end} ];then let "snp_end= n_end + 1" ;fi
done ## while loop all snps by chromosome
fi ## if checking whether chromosome has been extracted or not
done ## while loop waiting for split by chromosome
## Terminate if PLINK has not finished extracting file
if [ ! -e by_chr/${postfix}_chr${chr}.bed ]; then
echo "Timeover: PLINK failed to \extract by_chr/${postfix}_chr${chr}.bed"
exit 1
fi
done ## for loop all chromosomes
exit 0

So, this script will first split your file into a plink binary format by chromosome (in by_chr). Then, it will split each chromosome into a smaller file. The files by_chr may not be very useful for parallelization, but believe me, there're times that you will get to use them.

Create a free website or blog at WordPress.com.