Can not find idels!!

[scami]

Hi guys
I know that there will be a trivial answer to my question but I really can not understand where the problem is. I am aligning a 20x coverage Illumina paired end reads by using soap aligner. The command I am using is the following:

./soap -D ./reference.fasta.index  -a ./s_2_1_sequence.txt -b ./s_2_2_sequence.txt -o s_2_paired_mapped -u s_2_unpaired_reads -2 s_2_unpaired_mapped -v 3 -m 100 -x 400 -r 1 -g 10

I then concatenate the paired_mapped and unpaired_mapped (the real name should be single_mapped) and convert it with the tool soap2sam.pl available at the soapdenovo website. I then use the sam file to gather information about indel by using samtools:

1) samtools import (to get a bam file)
2) samtools sort (to sort it)
3) samtools pileup -f reference.fasta -i file_name.bam

With the last command I can not get any indel. I am not an expert but I think that no indel at all is not possible. Have you got any idea on what is going on here?

thanks a lot for the help



UPDATE: by looking at the soap aligner first output I realised that indels are present. This make me thing that the problem is the perl script soap2sam. I will try to call snp and indel directly with soapsnp.

[nilshomer]

Try the "-c" option in the "samtools pileup" command to call the consensus sequence.

[scami]

Just added -c before the -i flag in the above samtools pileup command. Unfortunately nothing changed!

cheers

[lh3]

use another aligner

[zee]

I agree,

BWA, Novoalign and BFAST all do a good job with finding indels from short read alignments and all of these programs support native SAM output.

use another aligner