SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 06-20-2010, 05:32 AM   #1
agc
Member
 
Location: Jerusalem

Join Date: May 2010
Posts: 26
Default 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?
agc is offline   Reply With Quote
Old 06-20-2010, 12:34 PM   #2
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by agc View Post
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.
nilshomer is offline   Reply With Quote
Old 06-21-2010, 06:33 AM   #3
GussowCarmelab
Junior Member
 
Location: Israel

Join Date: Apr 2010
Posts: 1
Default

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?
GussowCarmelab is offline   Reply With Quote
Old 06-21-2010, 10:59 AM   #4
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by GussowCarmelab View Post
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 View Post
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.
nilshomer is offline   Reply With Quote
Old 06-22-2010, 05:02 AM   #5
agc
Member
 
Location: Jerusalem

Join Date: May 2010
Posts: 26
Default

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 05:14 AM.
agc is offline   Reply With Quote
Old 06-22-2010, 03:37 PM   #6
bioinfosm
Senior Member
 
Location: USA

Join Date: Jan 2008
Posts: 482
Default

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
bioinfosm is offline   Reply With Quote
Old 06-27-2010, 07:55 AM   #7
agc
Member
 
Location: Jerusalem

Join Date: May 2010
Posts: 26
Default

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 is offline   Reply With Quote
Old 07-01-2010, 07:14 AM   #8
agc
Member
 
Location: Jerusalem

Join Date: May 2010
Posts: 26
Default

Quote:
Originally Posted by agc View Post
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.
agc is offline   Reply With Quote
Old 07-01-2010, 08:48 AM   #9
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by agc View Post
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.
nilshomer is offline   Reply With Quote
Reply

Tags
indels, pileup, samtools, tablet

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:50 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO