#!/usr/bin/perl  -w

use strict;

use Benchmark::Timer;

use Bio::DB::Sam;



## The BAM in question
my $bam_file
  = '../Results/bowtie_combined.sorted.bam';

## Get in in BioPerl format
my $sam = Bio::DB::Sam->
  new( -bam  => $bam_file );

## Check that it's alive
print "there are ",
  $sam->n_targets, " target sequences (contigs) in '$bam_file'\n";



## Get a list of all the target sequences
my @seqs = $sam->seq_ids;



## Prove that we can count the number of reads using the longest
## target sequence...

my $seq = $seqs[0];
my $len = $sam->length($seq);

print
  join("\t", $seq, $len,
       scalar( $sam->segment( $seq, 1, $len )->features() )
      ), "\n";



## Do some benchmarking of the top 1000 contigs (most of the big
## ones)...

my $timer = Benchmark::Timer->new(skip => 1);

for my $i (0..1000){
  $seq = $seqs[$i];
  $len = $sam->length($seq);

  $timer->start('count reads');
  my @features = $sam->segment( $seq, 1, $len )->features();
  $timer->stop;
}

print $timer->report;