Go Back   SEQanswers > Bioinformatics > Bioinformatics

Similar Threads
Thread Thread Starter Forum Replies Last Post
NO indel output from Samtools by using nearly default settings wwmm933 Bioinformatics 1 12-05-2012 04:18 AM
How to make sense of Tophat's output file 'junctions.bed' gsinghal RNA Sequencing 4 09-03-2012 07:49 AM
what does "hetero SNP/INDEL" look like in BWA-SW+Samtools shuang Bioinformatics 5 10-11-2011 09:44 AM
MAQ Simulator SNP/Indel Output peveralldubois Bioinformatics 1 02-12-2011 07:30 PM
samtools pileup output for SNP detection Mansequencer Bioinformatics 1 09-08-2010 07:36 PM

Thread Tools
Old 10-04-2011, 08:32 AM   #1
Senior Member
Location: IL

Join Date: Jul 2011
Posts: 100
Default How to make Samtools to output both SNP and INDEL variants?

How to make Smatools to output both SNP and INDEL variants?

I used bwa-sw to align my Sanger reads with the reference genome. Then I used Samtools mpileup to do variant calling (see commands below). However, it did NOT show INDELs variants.

$samtools mpileup -uf Sorbil.fasta 1A.bam > 1A.pileup
$bcftools view -bvcg 1A.pileup > 1A.bcf
$bcftools view 1A.bcf > 1A.vcf

chr_8 50759327 . T C 4.77 . DP=1;AF1=0.9999;CI95=0.5,1;DP4=0,0,0,1;MQ=60 PL:GT:GQ 33,3,0:1/1:41
chr_8 50759520 . T C 10.4 . DP=1;AF1=0.9999;CI95=0.5,1;DP4=0,0,0,1;MQ=60 PL:GT:GQ 40,3,0:1/1:41

When I tried "pileup", instead of "mpileup", it only showed INDELs, no SNPs.

$samtools pileup -vcf Sorbil.fasta 1A.bam > 1A.pileup

chr_8 50759483 a A 30 0 60 1 ,+1t M
chr_8 50759483 * +T/+T 40 0 60 1 +T * 1 0 0 0 0
shuang is offline   Reply With Quote
Old 12-12-2011, 12:37 PM   #2
Location: ct

Join Date: Jan 2010
Posts: 17

same problem. No indel output!
wwmm933 is offline   Reply With Quote
Old 05-20-2013, 06:36 PM   #3
Location: Blacksburg

Join Date: May 2013
Posts: 18

Here is what you should do to get SNPs and InDels

The BAM files for alignment map were sorted and indexed.

Then, here is what I did for pileup for processing for a small region of interest:

samtools mpileup -r 0:4,000,679-5,000,894 -uf e.coli.fa sorted.bam > sorted.region.mpileup

Then to obtain SNPs and InDels from mpileup file, here is what I did:

bcftools view -bvcg contigalignsorted.region.mpileup > sorted.region.bcf

bcftools view sorted.region.bcf > sorted.region.vcf

However, when i see my alignment of BAM file on IGV for this region, I notice that the final VCF file obtained through commands above, has a lot of SNPs and InDels missing! I can show a snapshot of a picture of SNP/InDel missed if required.

What went wrong in the calls ?
abi is offline   Reply With Quote
Old 05-21-2013, 12:02 PM   #4
Location: Blacksburg

Join Date: May 2013
Posts: 18

More importantly , it even shows up wrong SNPs which otherwise is not a SNP as per the bam file loaded on IGV.

0 4009334 . A C 125 . DP=5;VDB=3.277706e-02;RPB=-1.291774e+00;AF1=0.5;AC1=1;DP4=1,1,2,1;MQ=59;FQ=74;PV4=1,0,1,0.38 GT:PL:GQ 0/1:155,0,101:99
0 4009336 . T A 17.1 . DP=5;VDB=6.080000e-02;RPB=1.291774e+00;AF1=0.5;AC1=1;DP4=2,1,1,1;MQ=59;FQ=20.1;PV4=1,1,0.14,1 GT:PL:GQ 0/1:47,0,154:50

As above 4009334 is a SNP but not as per IGV picture below:
Attached Images
File Type: png igv_snapshot.png (22.4 KB, 4 views)
abi is offline   Reply With Quote

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

All times are GMT -8. The time now is 11:01 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO