Using samtools/tablet to find indels

[agc]

As I understand, the following code should extract indels from a bam file and output them in pileup format:

samtools pileup -if <ref.fa> <file.sorted.bam>

This command outputs nothing, however, when viewing the BAM file in tablet I can clearly see sections of the reference with no reads mapped to them.

Am I running samtools correctly? Is there a different way I should be looking for indels?

[nilshomer]

As I understand, the following code should extract indels from a bam file and output them in pileup format:

samtools pileup -if <ref.fa> <file.sorted.bam>

This command outputs nothing, however, when viewing the BAM file in tablet I can clearly see sections of the reference with no reads mapped to them.

Am I running samtools correctly? Is there a different way I should be looking for indels?

It will call indels that are observed in the reads themselves. It wont do any type of re-assembly based on missing coverage etc. Hope that clears it up.

[GussowCarmelab]

I see. Thanks for the quick reply. So, in other words, hypothetically that pileup command would show me:

Read         : ACACGGGACAC
Reference  : ACACACAC

as an insertion of 3 Gs? Does that mean that indels can only be found via gapped aligners (IE bwa can find them, but not bowtie)?

Also, why wouldn't it output parts of the reference that no reads were aligned to? Aren't those likely to be deletions?

[nilshomer]

I see. Thanks for the quick reply. So, in other words, hypothetically that pileup command would show me:

Read         : ACACGGGACAC
Reference  : ACACACAC

as an insertion of 3 Gs? Does that mean that indels can only be found via gapped aligners (IE bwa can find them, but not bowtie)?

Yes, gapped aligners will find gaps, ungapped aligners will not.

Also, why wouldn't it output parts of the reference that no reads were aligned to? Aren't those likely to be deletions?

Missing coverage could be due to sampling, but are likely deletions (or SVs). Deletions can also be observed with paired-end or mate-pair reads, where the distance between the two ends are farther away than expected. There are tools (like Breakway) that try to find SVs from paired-end (mate-pair) data. SAMtools does not try to find such variants.

[agc]

Our data is not mate-pair (50bp fragment library), but I would like to find fairly large indels (>200bp). Is this possible? I've read up on gapped aligners (BFAST and BWA), but they seem to be limited to smaller indels.

Is there an aligner / software that can do this? Perhaps taking missing coverage into account?

Thanks.

[bioinfosm]

For such large indels, with short reads, one needs paired or mated data!

Using lack of coverage is an interesting concept, and IMO if there is a reference sample, you may be able to call indels based on comparison with coverage in reference sample. But simply taking missing coverage in a sample to be a deletion is very likely a false positive!

[agc]

Utilizing a reference sample, is there a program (or samtools command) that can output the coordinates and sequences on the reference sequence where there is no coverage by the reads? We can roughly estimate the size of the indel that we are looking for, and hopefully we can weed out false positives that way.

[agc]

Utilizing a reference sample, is there a program (or samtools command) that can output the coordinates and sequences on the reference sequence where there is no coverage by the reads? We can roughly estimate the size of the indel that we are looking for, and hopefully we can weed out false positives that way.

I apologize for bumping this, but if there is no such program, we will attempt to write our own, but if there is already a software option for this, we'd rather utilize it then put in the time and effort of rewriting something that already exists.

Thanks.

[nilshomer]

I apologize for bumping this, but if there is no such program, we will attempt to write our own, but if there is already a software option for this, we'd rather utilize it then put in the time and effort of rewriting something that already exists.

Thanks.

A quick perl script using the output of your pileup should be able to do this.