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

[shuang]

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.

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

<result>
##fileformat=VCFv4.0
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT outdir/02A.bam
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.

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

<result>
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

[wwmm933]

same problem. No indel output!

[abi]

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]

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: