Generating a consensus sequence with samtools/bcftools/vcftools

Say you have mapped reads to a reference genome in .bam format. If you would like to get a consensus sequence of the alignment, run the following (you will probably need to install/locate samtools, bcftools and vcftools):

samtools mpileup -AuDS -Q 10 -q 10 -d 10000 -f mito_eulemur.fasta AKF23_sorted.bam | bcftools view -cg - | /usr/share/samtools/vcfutils.pl vcf2fq | awk 'BEGIN{P=1}{if(P==1||P==2){gsub(/^[@]/,">");print}; if(P==4)P=0; P++}' | sed -n '/+/q;p' > AKF23MitoCons.fasta
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s