I have probably a simple problem but which I cannot solve as of yet and puzzles me quite a bit. Using samtools mpileup and bcftools I am unable to call a small 6 base INSERT compared to my reference. Single SNPs are called ok. I have looked at many unanswered questions and potential answers here at BioStar and this question/answer gets closest but none actually has a working answer!
My play aorund data:
- Illumina PE 250 dataset of bacterial genome (2Mb at coverage 5k+)
- Assembled and reconstructed a genome from it
- If I check for SNPs by mapping the raw reads to the fresh assembly none are found (underlining the correct assembly)
- I created an artificial 2 position SNP and somewhere else I deleted 6 bases in the assembly
- I mapped the original reads which should detect the SNPs and the INSERT. Hoewever I only see the SNPs in the vcf and not the INSERT
My mapping is done with bowtie2 --local-sensitive settings to the ArtificialREFERENCE-6del sequence
Sam file is converted to bam and sorted using samtools
Then I call variants:
samtools mpileup -uf ArtificialREFERENCE-6del.fasta SORTED_BAMfile.bam -o SAMPLE_file.bcf bcftools call -vm -o SAMPLE_file.vcf SAMPLE_file.bcf
I tried different bcftools settings like -c and -p 0.5 ... no luck...
Any help would be appreciated! It is probably something very obvious but I fail to see it from the documentation.
NB: I tried samtools/bcftools version 1.3.1 but also an older 1.1
NB2: I tried the additional options and FLAGS as proposed in the above linked post without luck.
NB3: The below is an IGV snippet of the involved BAM file clearly showing there IS a 6 base insert.