SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics
Similar Threads
Thread Thread Starter Forum Replies Last Post
No output from Samtools pileup pravee1216 Bioinformatics 1 12-05-2012 03:59 AM
samtools pileup BioSlayer Bioinformatics 9 01-16-2012 03:10 AM
samtools pileup raghu1990 Bioinformatics 0 06-28-2011 02:26 AM
samtools pileup wisosonic Bioinformatics 0 05-11-2011 04:57 AM
Samtools version 1.0.7 and 1.0.10 different pileup mrxcm3 Bioinformatics 1 11-19-2010 07:24 AM

Reply
 
Thread Tools
Old 05-17-2010, 08:33 AM   #1
Graham Etherington
Member
 
Location: Norwich, England.

Join Date: Apr 2010
Posts: 22
Default Samtools Pileup Parser

Hi,
I'm trying to parse the Samtools Pileup format in Perl.

I've created a paired-end assembly using BWA (aln, sampe) and then used Samtools to create a pileup (using import - sort - index - pileup).

What I want to do is to open up the pileup (which is obviously in the Samtools Pileup format) and iterate over each position of the pileup, examining differences in base-calls.

Does anyone know of anything available, or have and sample code they could post?
Many thanks.
Graham Etherington is offline   Reply With Quote
Old 05-17-2010, 08:57 AM   #2
Pepe
Member
 
Location: Germany

Join Date: Mar 2009
Posts: 28
Default

You cannot open a file in Pileup format with Bio:B::sam

You should convert your SAM file to BAM format (samtools view) and sort it (samtools sort), maybe index it too? I'm not sure (samtools index).

Then follow the instructions in the CPAN website, which are very, veryhelpful. Basically:

use Bio:B::Sam;
my $sam = Bio:B::Sam->new(-bam => "my_sorted_bam_file.bam",
-fasta => "my_ref.fasta");
my @targets = $sam->seq_ids;
foreach my $chr (@targets){
$sam->pileup($chr,$my_subroutine);
}

where $my_subroutine is a subroutine that you'll need to get to do what you want: snp calling, coverage calculation...
Pepe is offline   Reply With Quote
Old 05-18-2010, 12:01 AM   #3
Graham Etherington
Member
 
Location: Norwich, England.

Join Date: Apr 2010
Posts: 22
Default

Quote:
Originally Posted by Pepe View Post
You cannot open a file in Pileup format with Bio::DB::sam

You should convert your SAM file to BAM format (samtools view) and sort it (samtools sort), maybe index it too? I'm not sure (samtools index).

Then follow the instructions in the CPAN website, which are very, veryhelpful. Basically:

use Bio::DB::Sam;
my $sam = Bio::DB::Sam->new(-bam => "my_sorted_bam_file.bam",
-fasta => "my_ref.fasta");
my @targets = $sam->seq_ids;
foreach my $chr (@targets){
$sam->pileup($chr,$my_subroutine);
}

where $my_subroutine is a subroutine that you'll need to get to do what you want: snp calling, coverage calculation...
Hi,
Thanks for your answer. I'd realised that Bio::DB:Sam couldn't parse the pileup, so this is why I am asking if anyone has, or knows about, anything that could parse it, before I spent time writing a Samtools parser myself.
I have a 'Parser.pm' module that parses various pileup formats (maq, bowtie, novoalign, etc), but I haven't added one for Samtools format yet.
Cheers,
Graham
Graham Etherington is offline   Reply With Quote
Old 05-18-2010, 05:14 AM   #4
nekrut
Member
 
Location: Penn State

Join Date: Apr 2009
Posts: 22
Default

Galaxy's (http://usegalaxy.org) pileup parser might do the trick:
http://bit.ly/cXtqD9
nekrut is offline   Reply With Quote
Old 05-18-2010, 07:58 AM   #5
Graham Etherington
Member
 
Location: Norwich, England.

Join Date: Apr 2010
Posts: 22
Default

Quote:
Originally Posted by nekrut View Post
Galaxy's (http://usegalaxy.org) pileup parser might do the trick:
http://bit.ly/cXtqD9
Thanks. That should give me enough to get started with.
Best wishes,
Graham
Graham Etherington is offline   Reply With Quote
Old 08-24-2012, 07:15 AM   #6
cschu
Junior Member
 
Location: Earth

Join Date: Aug 2012
Posts: 1
Default

I noticed that the referenced Galaxy code does not check for the '>' and '<' characters. Can anybody please provide some help how to process those? I.e., I don't really understand what "reference skip" means and would think that such a symbol means that the read does not cover the the respective position. Any help would be appreciated, thanks!
cschu is offline   Reply With Quote
Reply

Tags
pileup, samtools

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 05:28 PM.


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