parsing bam

[taha.hojati]

Hi,

I recently concatenated many 200nt sequences to the end of each other; these sequences were from different parts of genome. Then I fed this to Geneious as a reference sequences and mapped sequencing reads (each read = 70nt) to it. I exported the results as a bam file. Now I need to write a Perl script that parses the bam file and finds out whether or not each of those initial sequences (the 200nt ones that were concatenated) is covered in sequencing reads. I have read the Bio:B::Sam tutorial on CPAN but I am still stuck. Here is all the code I have written so far:

use Bio::DB::Sam;
my $sam = Bio::DB::Sam -> new (-fasta => ... , -bam => ...) ;
my $segment;
($id)=$sam -> seq_ids;
while ($i < $sam->length($id)) { 
$segment =$ sam -> segment($id , $i ,$i +200) ;

# What should I write here?  

$i +=200; 
}

I would appreciate your help, since CPAN manual only confuses me.

[adaptivegenome]

What do you mean by covered? Do you mean the coverage across every position of a 200nt fragment is at least 1X?

If so its a pretty simple problem to solve. If you don't care about speed you could use mpileup to check if each of the 200 positions for a given range are covered by at least 1 70nt read.

If you like to program, then you could use any of the various samtools bindings to PERL, C, C++,etc. You would open the BAM and step through each position to see if there is at least 1 read returned. Maybe there is even a faster way to do it, but either of the above two suggestions would work.

[taha.hojati]

Thanks for the suggestion. I used the pileup method in Perl and it worked fine.