Using samtools/tablet to find indels

agc · 06-20-2010, 04:32 AM

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

Code:

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 · 06-20-2010, 11:34 AM

Quote:

Originally Posted by agc

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

Code:

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 · 06-21-2010, 05:33 AM

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

Code:

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 · 06-21-2010, 09:59 AM

Quote:

Originally Posted by GussowCarmelab

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

Code:

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.

Quote:

Originally Posted by GussowCarmelab

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 · 06-22-2010, 04:02 AM

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 · 06-22-2010, 02:37 PM

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 · 06-27-2010, 06:55 AM

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 · 07-01-2010, 06:14 AM

Quote:

Originally Posted by 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 · 07-01-2010, 07:48 AM

Quote:

Originally Posted by agc

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.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
samtools indels	Robby	Bioinformatics	3	11-08-2011 07:02 AM
find all snps/indels	prbndr	Bioinformatics	2	09-20-2011 10:43 AM
samtools (or other tools) to find larger indels	waalkes	Bioinformatics	0	06-15-2011 02:53 PM
indels from samtools	coonya	Bioinformatics	2	06-06-2011 04:12 PM
Tablet not showing pileup indels	agc	Bioinformatics	1	11-21-2010 07:35 AM

06-20-2010, 04:32 AM	#1
agc Member Location: Jerusalem Join Date: May 2010 Posts: 26	Using samtools/tablet to find indels As I understand, the following code should extract indels from a bam file and output them in pileup format: Code: 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?

06-21-2010, 05:33 AM	#3
GussowCarmelab Junior Member Location: Israel Join Date: Apr 2010 Posts: 1	I see. Thanks for the quick reply. So, in other words, hypothetically that pileup command would show me: Code: 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?

06-22-2010, 04:02 AM	#5
agc Member Location: Jerusalem Join Date: May 2010 Posts: 26	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. Last edited by agc; 06-22-2010 at 04:14 AM.

06-22-2010, 02:37 PM	#6
bioinfosm Senior Member Location: USA Join Date: Jan 2008 Posts: 482	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! __________________ -- bioinfosm

06-27-2010, 06:55 AM	#7
agc Member Location: Jerusalem Join Date: May 2010 Posts: 26	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.