Bio/Assembly/IO/sam.pm

   1 #
   2 # BioPerl module for Bio::Assembly::IO::sam
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Mark A. Jensen <maj -at- fortinbras -dot- us>
   7 #
   8 # Copyright Mark A. Jensen
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Assembly::IO::sam - An IO module for assemblies in Sam format *BETA*
  17
  18 =head1 SYNOPSIS
  19
  20  $aio = Bio::Assembly::IO( -file => "mysam.bam",
  21                            -refdb => "myrefseqs.fas");
  22  $assy = $aio->next_assembly;
  23
  24 =head1 DESCRIPTION
  25
  26 This is a (currently) read-only IO module designed to convert
  27 Sequence/Alignment Map (SAM; L<http://samtools.sourceforge.net/>)
  28 formatted alignments to L<Bio::Assembly::Scaffold> representations,
  29 containing .L<Bio::Assembly::Contig> and L<Bio::Assembly::Singlet>
  30 objects. It uses lstein's L<Bio::DB::Sam> to parse binary formatted SAM
  31 (.bam) files guided by a reference sequence fasta database.
  32
  33 B<NB>: C<Bio::DB::Sam> is not a BioPerl module; it can be obtained via
  34 CPAN. It in turn requires the C<libbam> library; source can be
  35 downloaded at L<http://samtools.sourceforge.net/>.
  36
  37 =head1 DETAILS
  38
  39 =over
  40
  41 =item * Required files
  42
  43 A binary SAM (C<.bam>) alignment and a reference sequence database in
  44 FASTA format are required. Various required indexes (C<.fai>, C<.bai>)
  45 will be created as necessary (via L<Bio::DB::Sam>).
  46
  47 =item * Compressed files
  48
  49 ...can be specified directly , if L<IO::Uncompress::Gunzip> is
  50 installed. Get it from your local CPAN mirror.
  51
  52 =item * BAM vs. SAM
  53
  54 The input alignment should be in (possibly gzipped) binary SAM
  55 (C<.bam>) format. If it isn't, you will get a message explaining how
  56 to convert it, viz.:
  57
  58  $ samtools view -Sb mysam.sam > mysam.bam
  59
  60 The bam file must also be sorted on coordinates: do
  61
  62  $ samtools sort mysam.unsorted.bam > mysam.bam
  63
  64 =item * Contigs
  65
  66 Contigs are calculated by this module, using the 'coverage' feature of
  67 the L<Bio::DB::Sam> object. A contig represents a contiguous portion
  68 of a reference sequence having non-zero coverage at each base.
  69
  70 The bwa assembler (L<http://bio-bwa.sourceforge.net/>) can assign read
  71 sequences to multiple reference sequence locations. The present
  72 implementation currently assigns such reads only to the first contig
  73 in which they appear.
  74
  75 =item * Consensus sequences
  76
  77 Consensus sequence and quality objects are calculated by this module,
  78 using the C<pileup> callback feature of C<Bio::DB::Sam>. The consensus
  79 is (currently) simply the residue at a position that has the maximum
  80 sum of quality values. The consensus quality is the integer portion of
  81 the simple average of quality values for the consensus residue.
  82
  83 =item * SeqFeatures
  84
  85 Read sequences stored in contigs are accompanied by the following
  86 features:
  87
  88  contig : name of associated contig
  89  cigar  : CIGAR string for this read
  90
  91 If the read is paired with a successfully mapped mate, these features
  92 will also be available:
  93
  94  mate_start  : coordinate of to which the mate was aligned
  95  mate_len    : length of mate read
  96  mate_strand : strand of mate (-1 or 1)
  97  insert_size : size of insert spanned by the mate pair
  98
  99 These features are obtained as follows:
 100
 101  @ids = $contig->get_seq_ids;
 102  $an_id = $id[0]; # or whatever
 103  $seq = $contig->get_seq_by_name($an_id);
 104  # Bio::LocatableSeq's aren't SeqFeature containers, so...
 105  $feat = $contig->get_seq_feat_by_tag(
 106             $seq, "_aligned_coord:".$s->id
 107          );
 108  ($cigar) = $feat->get_tag_values('cigar');
 109  # etc.
 110
 111 =back
 112
 113 =head1 TODO
 114
 115 =over
 116
 117 =item * Supporting both text SAM (TAM) and binary SAM (BAM)
 118
 119 =back
 120
 121 =head1 FEEDBACK
 122
 123 =head2 Mailing Lists
 124
 125 User feedback is an integral part of the evolution of this and other
 126 Bioperl modules. Send your comments and suggestions preferably to
 127 the Bioperl mailing list.  Your participation is much appreciated.
 128
 129   bioperl-l@bioperl.org                  - General discussion
 130 http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 131
 132 =head2 Support
 133
 134 Please direct usage questions or support issues to the mailing list:
 135
 136 L<bioperl-l@bioperl.org>
 137
 138 rather than to the module maintainer directly. Many experienced and
 139 reponsive experts will be able look at the problem and quickly
 140 address it. Please include a thorough description of the problem
 141 with code and data examples if at all possible.
 142
 143 =head2 Reporting Bugs
 144
 145 Report bugs to the Bioperl bug tracking system to help us keep track
 146 of the bugs and their resolution. Bug reports can be submitted via
 147 the web:
 148
 149   https://github.com/bioperl/bioperl-live/issues
 150
 151 =head1 AUTHOR - Mark A. Jensen
 152
 153 Email maj -at- fortinbras -dot- us
 154
 155 =head1 APPENDIX
 156
 157 The rest of the documentation details each of the object methods.
 158 Internal methods are usually preceded with a _
 159
 160 =cut
 161
 162 # Let the code begin...
 163
 164
 165 package Bio::Assembly::IO::sam;
 166 use strict;
 167 use warnings;
 168
 169 # Object preamble - inherits from Bio::Root::Root
 170
 171 use Bio::Seq::Quality;
 172 use Bio::Seq::PrimaryQual;
 173 use Bio::LocatableSeq;
 174 use Bio::Assembly::IO;
 175 use Bio::Assembly::Scaffold;
 176 use Bio::Assembly::Contig;
 177 use Bio::Assembly::Singlet;
 178 use Bio::SeqIO;
 179 use File::Spec;
 180 use File::Basename;
 181 use File::Temp qw(tempfile);
 182 use Carp;
 183 use Bio::Root::Root;
 184 use base qw(Bio::Root::Root Bio::Assembly::IO Bio::Root::IO);
 185
 186 our $HAVE_IO_UNCOMPRESS;
 187 BEGIN {
 188 # check requirements
 189     unless ( eval "require Bio::DB::Sam; 1" ) {
 190         Bio::Root::Root->throw(__PACKAGE__.' requires installation of samtools (libbam) and Bio::DB::Sam (available on CPAN; not part of BioPerl)');
 191     }
 192     unless ( eval "require IO::Uncompress::Gunzip; \$HAVE_IO_UNCOMPRESS = 1") {
 193         Bio::Root::Root->warn("IO::Uncompress::Gunzip is not available; you'll have to do your decompression by hand.");
 194     }
 195 }
 196
 197 my $progname = 'sam';
 198
 199 sub new {
 200     my $class = shift;
 201     my @args = @_;
 202     my $self = $class->SUPER::new(@args);
 203     my ($file, $refdb, $format) = $self->_rearrange([qw(FILE REFDB FORMAT)], @args);
 204     $self->file($file);
 205     $refdb && $self->refdb($refdb);
 206     $self->_init_sam() or croak( "Sam initialization failed" );
 207     $self->{_assigned} = {};
 208     return $self;
 209 }
 210
 211
 212 =head1 Bio::Assembly::IO compliance
 213
 214 =head2 next_assembly()
 215
 216     Title   : next_assembly
 217     Usage   : my $scaffold = $asmio->next_assembly();
 218     Function: return the next assembly in the sam-formatted stream
 219     Returns : Bio::Assembly::Scaffold object
 220     Args    : none
 221
 222 =cut
 223
 224 sub next_assembly {
 225     my $self = shift;
 226     my @contig_set;
 227     # get Bio::DB::Sam object
 228     # could add a refdb or fasfile attribute to contain the reference db name
 229     # iterate through seq_ids...
 230     my @refseq_ids = $self->sam->seq_ids;
 231     my $assembly = Bio::Assembly::Scaffold->new( -progname => $progname );
 232
 233     foreach my $id (@refseq_ids) {
 234         #### this is choice 1: all refseqs into one assy...###
 235         $self->_current_refseq_id( $id );
 236         # Load contigs and singlets in the scaffold
 237         while ( my $obj = $self->next_contig()) {
 238             # Add contig /singlet to assembly
 239             if ($obj->isa('Bio::Assembly::Singlet')) { # a singlet
 240                 $assembly->add_singlet($obj);
 241             } else { # a contig
 242                 $assembly->add_contig($obj);
 243             }
 244         }
 245     }
 246     return $assembly;
 247 }
 248
 249 =head2 next_contig()
 250
 251     Title   : next_contig
 252     Usage   : my $contig = $asmio->next_contig();
 253     Function: return the next contig or singlet from the sam stream
 254     Returns : Bio::Assembly::Contig or Bio::Assembly::Singlet
 255     Args    : none
 256
 257 =cut
 258
 259 sub next_contig {
 260     my $self = shift;
 261     if (!defined $self->_current_contig_seg_idx) {
 262         $self->_current_contig_seg_idx(0);
 263     }
 264     else {
 265         $self->_current_contig_seg_idx( 1+$self->_current_contig_seg_idx );
 266     }
 267     unless ($self->_current_refseq_id) {
 268         croak("No current refseq id set");
 269     }
 270     my $contig_segs = $self->_segset($self->_current_refseq_id);
 271     unless ($contig_segs && @$contig_segs) {
 272         croak("No contig segset for current id '".$self->_current_refseq_id."'")
 273     }
 274     # each segment in @$contig_segs represents a contig within the
 275     # current refseq
 276     my $contig_seg = $$contig_segs[$self->_current_contig_seg_idx];
 277     return if (!defined $contig_seg); # iterator finished
 278     # each 'feature' in $contig_seg represents a read;
 279     # $seqio lets us iterate efficiently over the reads:
 280     my $seqio = $contig_seg->features(-iterator => 1);
 281
 282
 283     # Contig and read related
 284     my $contigobj = $self->_store_contig($contig_seg);
 285     my $numseq = 0;
 286
 287     while ( my $read = $seqio->next_seq ) {
 288         if ($self->{_assigned}->{$read->name}) {
 289             next;
 290         }
 291         $self->{_assigned}->{$read->name}=1;
 292         $self->_store_read($read, $contigobj);
 293         $numseq++;
 294     }
 295     if ($numseq == 1) { # ooh! a singlet!
 296         $contigobj = $self->_store_singlet($contigobj);
 297     }
 298     return $contigobj;
 299 }
 300
 301 =head2 write_assembly()
 302
 303  Title   : write_assembly
 304  Usage   :
 305  Function: not implemented (module currrently read-only)
 306  Returns :
 307  Args    :
 308
 309 =cut
 310
 311 sub write_assembly {
 312     my $self = shift;
 313     $self->throw( "This module is currently read-only" );
 314 }
 315
 316 =head1 Internal
 317
 318 =head2  _store_contig()
 319
 320     Title   : _store_contig
 321     Usage   : my $contigobj = $self->_store_contig(\%contiginfo);
 322     Function: create and load a contig object
 323     Returns : Bio::Assembly::Contig object
 324     Args    : Bio::DB::Sam::Segment object
 325
 326 =cut
 327
 328 sub _store_contig {
 329     my ($self, $contig_seg) = @_;
 330
 331     # Create a contig
 332     my $contigobj = Bio::Assembly::Contig->new(
 333         -id     => 'sam_assy['.$self->_basename.':'.$self->_current_refseq_id.']/'.$contig_seg->start.'-'.$contig_seg->end,
 334         -source => $progname,
 335         -strand => 1
 336     );
 337     my $consobj = $self->_calc_consensus($contig_seg);
 338     my $consensus_seq = Bio::LocatableSeq->new(
 339         -id    => $contigobj->id,
 340         -seq   => $consobj->seq,
 341         -start => 1,
 342     );
 343     $contigobj->set_consensus_sequence($consensus_seq);
 344     my $consensus_qual = Bio::Seq::PrimaryQual->new(
 345         -id    => $contigobj->id,
 346         -qual  => $consobj->qual,
 347         -start => 1,
 348     );
 349     $contigobj->set_consensus_quality($consensus_qual);
 350
 351     # Add other misc contig information as subsequence feature
 352     #my @other = grep !/asmbl_id|end|qualobj|start/, keys %$contiginfo;
 353     #my %other;
 354     #@other{@other} = @$contiginfo{@other};
 355     #my $contigtags = Bio::SeqFeature::Generic->new(
 356     #    -primary     => '_main_contig_feature',
 357     #    -source      => $$contiginfo{'asmbl_id'},
 358     #    -start       => 1,
 359     #    -end         => $contig_seg->length,
 360     #    -strand      => 1,
 361     #    # dumping ground:
 362     #    -tag         => \%other
 363     #);
 364     #$contigobj->add_features([ $contigtags ], 1);
 365
 366     return $contigobj;
 367 }
 368
 369 =head2 _store_read()
 370
 371     Title   : _store_read
 372     Usage   : my $readobj = $self->_store_read($readobj, $contigobj);
 373     Function: store information of a read belonging to a contig in a contig object
 374     Returns : Bio::LocatableSeq
 375     Args    : Bio::DB::Bam::AlignWrapper, Bio::Assembly::Contig
 376
 377 =cut
 378
 379 sub _store_read {
 380     my $self = shift;
 381     my ($read, $contigobj) = @_;
 382     my $readseq = Bio::LocatableSeq->new(
 383         -display_id => $read->name,
 384         -primary_id => $read->name,
 385         -seq        => $read->dna,
 386         -start      => 1,
 387         -strand     => $read->strand,
 388         -alphabet   => 'dna'
 389     );
 390     my $qual = Bio::Seq::PrimaryQual->new(
 391         -id         => $read->name."_qual",
 392         -qual       => [$read->qscore]
 393     );
 394
 395     # add pair information
 396     my @pair_info;
 397     if ($read->proper_pair) { # mate also aligned
 398         @pair_info = (
 399             mate_start => $read->mate_start,
 400             mate_len   => $read->mate_len,
 401             mate_strand => $read->mstrand,
 402             insert_size => $read->isize
 403         );
 404     }
 405
 406     my $alncoord = Bio::SeqFeature::Generic->new(
 407         -primary    => $read->name,
 408         -start      => $read->start,
 409         -end        => $read->end,
 410         -strand     => $read->strand,
 411         -qual       => join(' ',$read->qscore),
 412         -tag        => { 'contig' => $contigobj->id,
 413                          'cigar'  => $read->cigar_str,
 414                          @pair_info }
 415     );
 416
 417     $contigobj->set_seq_coord($alncoord, $readseq);
 418     $contigobj->set_seq_qual( $readseq, $qual );
 419
 420     #add other misc read info as subsequence feature:
 421     #my @other = grep !/aln_(?:end|start)|seq(?:str)?|strand/, keys %$readinfo;
 422     #my %other;
 423     #@other{@other} = @$readinfo{@other};
 424     #my $readtags = Bio::SeqFeature::Generic->new(
 425     #    -primary     => '_main_read_feature',
 426     #    -source      => $readobj->id,
 427     #    -start       => $$readinfo{'aln_start'},
 428     #    -end         => $$readinfo{'aln_end'},
 429     #    -strand      => $$readinfo{'strand'},
 430     #    # dumping ground:
 431     #    -tag         => \%other
 432     #);
 433     #$contigobj->get_features_collection->add_features([$readtags]);
 434     #$contigobj->get_features_collection->add_SeqFeature($alncoord, $readtags);
 435
 436     return $readseq;
 437 }
 438
 439 =head2  _store_singlet()
 440
 441     Title   : _store_singlet
 442     Usage   : my $singletobj = $self->_store_singlet($contigobj);
 443     Function: convert a contig object containing a single read into
 444               a singlet object
 445     Returns : Bio::Assembly::Singlet
 446     Args    : Bio::Assembly::Contig (previously loaded with only one seq)
 447
 448 =cut
 449
 450 sub _store_singlet {
 451     my $self = shift;
 452     my ($contigobj) = @_;
 453
 454     my ($readseq) = $contigobj->each_seq;
 455
 456     my $singletobj = Bio::Assembly::Singlet->new( -id => $contigobj->id,
 457                                                   -seqref => $readseq );
 458
 459 # may want to attach this someday
 460 #    my $qual = $contigobj->get_qual_by_name($readseq->id);
 461
 462     return $singletobj;
 463 }
 464
 465 =head1 REALLY Internal
 466
 467 =head2 _init_sam()
 468
 469  Title   : _init_sam
 470  Usage   : $self->_init_sam($fasfile)
 471  Function: obtain a Bio::DB::Sam parsing of the associated sam file
 472  Returns : true on success
 473  Args    : [optional] name of the fasta reference db (scalar string)
 474  Note    : The associated file can be plain text (.sam) or binary (.bam);
 475            If the fasta file is not specified, and no file is contained in
 476            the refdb() attribute, a .fas file with the same
 477            basename as the sam file will be searched for.
 478
 479 =cut
 480
 481 sub _init_sam {
 482     my $self = shift;
 483     my $fasfile = shift;
 484     my $file = $self->file;
 485     my $sam;
 486     $fasfile ||= $self->refdb;
 487     $file =~ s/^[<>+]*//; # byebye parasitic mode chars
 488     my ($fname, $dir, $suf) = fileparse($file, ".sam", ".bam");
 489     $self->_set_from_args({ '_basename' => $fname },
 490                           -methods => [qw( _basename)],
 491                           -create => 1);
 492     if (!defined $fasfile) {
 493         for (qw( fas fa fasta fas.gz fa.gz fasta.gz )) {
 494             $fasfile = File::Spec->catdir($dir, $self->_basename.$_);
 495             last if -e $fasfile;
 496             undef $fasfile;
 497         }
 498     }
 499     unless (-e $fasfile) {
 500         croak( "Can't find associated reference fasta db" );
 501     }
 502     !$self->refdb && $self->refdb($fasfile);
 503     # compression
 504     if ($fasfile =~ /\.gz[^.]*$/) {
 505         unless ($HAVE_IO_UNCOMPRESS) {
 506             croak( "IO::Uncompress::Gunzip not available; can't decompress on the fly");
 507         }
 508         my ($tfh, $tf) = tempfile( UNLINK => 1);
 509         my $z = IO::Uncompress::Gunzip->new($fasfile) or croak("Can't expand: $@");
 510         while (<$z>) { print $tfh $_ }
 511         close $tfh;
 512         $fasfile = $tf;
 513     }
 514     if ($file =~ /\.gz[^.]*$/) {
 515         unless ($HAVE_IO_UNCOMPRESS) {
 516             croak( "IO::Uncompress::Gunzip not available; can't decompress on the fly");
 517         }
 518         my ($tfh, $tf) = tempfile( UNLINK => 1);
 519         my $z = IO::Uncompress::Gunzip->new($file) or croak("Can't expand: $@");
 520         while (<$z>) {
 521             print $tfh $_;
 522         }
 523         close $tfh;
 524         $file = $tf;
 525     }
 526     # sam conversion : just barf for now
 527     if (-T $file) {
 528         my $bam = $file;
 529         $bam =~ s/\.sam/\.bam/;
 530         croak( "'$file' looks like a text file.\n\tTo convert to the required .bam (binary SAM) format, run\n\t\$ samtools view -Sb $file > $bam\n");
 531     }
 532
 533     $sam = Bio::DB::Sam->new( -bam => $file,
 534                               -fasta => $fasfile,
 535                               -autoindex  => 1,
 536                               -expand_flags => 1);
 537     unless (defined $sam) {
 538         croak( "Couldn't create the Bio::DB::Sam object" );
 539     }
 540     $self->{sam} = $sam;
 541     # now produce the contig segments for each seq_id...
 542     for ($sam->seq_ids) {
 543         my $seg = $sam->segment(-seq_id=>$_, -start=>1, -end=>$sam->length($_));
 544         ${$self->{_segset}}{$_} = [$self->_get_contig_segs_from_coverage($seg)];
 545     }
 546
 547     return 1;
 548
 549 }
 550
 551 =head2 _get_contig_segs_from_coverage()
 552
 553  Title   : _get_contig_segs_from_coverage
 554  Usage   :
 555  Function: calculates separate contigs using coverage info
 556            in the segment
 557  Returns : array of Bio::DB::Sam::Segment objects, representing
 558            each contig
 559  Args    : Bio::DB::Sam::Segment object
 560
 561 =cut
 562
 563 sub _get_contig_segs_from_coverage {
 564     my $self = shift;
 565     my $segment = shift;
 566     unless ($self->sam) {
 567         croak("Sam object not yet initialized (call _init_sam)");
 568     }
 569     unless ( ref($segment) =~ /Bio::DB::Sam::Segment/ ) {
 570         croak("Requires Bio::DB::Sam::Segment object at arg 1");
 571     }
 572     my ($cov, @covdata, @rngs, @segs);
 573     ($cov) = $segment->features('coverage');
 574     unless ($cov) {
 575         croak("No coverage data!");
 576     }
 577     @covdata = $cov->coverage;
 578
 579     # calculate contigs:
 580     my $in_contig;
 581     my ($lf_end,$rt_end);
 582     for (0..$#covdata) {
 583         if ($covdata[$_]) {
 584             if ($in_contig) {
 585                 $rt_end = $_+1;
 586                 next;
 587             }
 588             else {
 589                 $in_contig = 1;
 590                 # push previous range
 591                 if (defined $lf_end && defined $rt_end) {
 592                     push @rngs, [$lf_end, $rt_end];
 593                 }
 594                 $lf_end = $_+1;
 595             }
 596         }
 597         else {
 598             $in_contig = 0;
 599         }
 600     }
 601     # last one
 602     push @rngs, [$lf_end, $rt_end] if (defined $lf_end and
 603                                        defined $rt_end and
 604                                        $lf_end <= $rt_end);
 605     unless (@rngs) {
 606         carp ("No coverage for this segment!");
 607         return;
 608     }
 609     for (@rngs) {
 610         push @segs, $self->sam->segment(-seq_id=>$segment->seq_id,
 611                                         -start=>$$_[0],
 612                                         -end=>$$_[1]);
 613     }
 614     return @segs;
 615 }
 616
 617 =head2 _calc_consensus_quality()
 618
 619  Title   : _calc_consensus_quality
 620  Usage   : @qual = $aio->_calc_consensus_quality( $contig_seg );
 621  Function: calculate an average or other data-reduced quality
 622            over all sites represented by the features contained
 623            in a Bio::DB::Sam::Segment
 624  Returns :
 625  Args    : a Bio::DB::Sam::Segment object
 626
 627 =cut
 628
 629 sub _calc_consensus_quality {
 630     # just an average over sites for now...
 631     my $self = shift;
 632     my $seg = shift;
 633
 634     my @quals;
 635     my $region = $seg->seq_id.':'.$seg->start.'..'.$seg->end;
 636     my $qual_averager = sub {
 637         my ($seqid, $pos, $p) = @_;
 638         return unless ($seg->start <= $pos and $pos <= $seg->end);
 639         my $acc = 0;
 640         my $n = 0;
 641         for my $pileup (@$p) {
 642             my $b = $pileup->alignment;
 643             $acc += $b->qscore->[$pileup->qpos];
 644             $n++;
 645         }
 646         push @quals, int($acc/$n);
 647     };
 648     $self->sam->pileup($region, $qual_averager);
 649     return @quals;
 650 }
 651
 652 =head2 _calc_consensus()
 653
 654  Title   : _calc_consensus
 655  Usage   : @qual = $aio->_calc_consensus( $contig_seg );
 656  Function: calculate a simple quality-weighted consensus sequence
 657            for the segment
 658  Returns : a SeqWithQuality object
 659  Args    : a Bio::DB::Sam::Segment object
 660
 661 =cut
 662
 663 sub _calc_consensus {
 664     # just an average over sites for now...
 665     my $self = shift;
 666     my $seg = shift;
 667
 668     my @quals;
 669     my $conseq ='';
 670     my $region = $seg->seq_id.':'.$seg->start.'-'.$seg->end;
 671
 672     my $weighted_consensus = sub {
 673         my ($seqid, $pos, $p) = @_;
 674         return unless ($seg->start <= $pos and $pos <= $seg->end);
 675         my %wt_tbl;
 676         my %n;
 677         for my $pileup (@$p) {
 678             my $b = $pileup->alignment;
 679             my $res = substr($b->qseq,$pileup->qpos,1);
 680             $wt_tbl{$res} += $b->qscore->[$pileup->qpos] || 0;
 681             $n{$res} ||= 0;
 682             $n{$res}++;
 683         }
 684         # really simple
 685         my $c = (sort { $wt_tbl{$b}<=>$wt_tbl{$a} } keys %wt_tbl)[0];
 686         $conseq .= $c;
 687         push @quals, int($wt_tbl{$c}/$n{$c});
 688     };
 689
 690     $self->sam->pileup($region, $weighted_consensus);
 691     return Bio::Seq::Quality->new(
 692         -qual => join(' ', @quals ),
 693         -seq => $conseq,
 694         -id => $region
 695     );
 696 }
 697
 698 =head2 refdb()
 699
 700  Title   : refdb
 701  Usage   : $obj->refdb($newval)
 702  Function: the name of the reference db fasta file
 703  Example :
 704  Returns : value of refdb (a scalar)
 705  Args    : on set, new value (a scalar or undef, optional)
 706
 707 =cut
 708
 709 sub refdb {
 710     my $self = shift;
 711
 712     return $self->{'refdb'} = shift if @_;
 713     return $self->{'refdb'};
 714 }
 715
 716 =head2 _segset()
 717
 718  Title   : _segset
 719  Usage   : $segset_hashref = $self->_segset()
 720  Function: hash container for the Bio::DB::Sam::Segment objects that
 721            represent each set of contigs for each seq_id
 722            { $seq_id => [@contig_segments], ... }
 723  Example :
 724  Returns : value of _segset (a hashref) if no arg,
 725            or the arrayref of contig segments, if arg == a seq id
 726  Args    : [optional] seq id (scalar string)
 727  Note    : readonly; hash elt set in _init_sam()
 728
 729 =cut
 730
 731 sub _segset {
 732     my $self = shift;
 733     return $self->{'_segset'} unless @_;
 734     return ${$self->{'_segset'}}{shift()};
 735 }
 736
 737 =head2 _current_refseq_id()
 738
 739  Title   : _current_refseq_id
 740  Usage   : $obj->_current_refseq_id($newval)
 741  Function: the "current" reference sequence id
 742  Example :
 743  Returns : value of _current_refseq (a scalar)
 744  Args    : on set, new value (a scalar or undef, optional)
 745
 746 =cut
 747
 748 sub _current_refseq_id {
 749     my $self = shift;
 750
 751     return $self->{'_current_refseq_id'} = shift if @_;
 752     return $self->{'_current_refseq_id'};
 753 }
 754
 755 =head2 _current_contig_seg_idx()
 756
 757  Title   : current_contig_seg_idx
 758  Usage   : $obj->current_contig_seg_idx($newval)
 759  Function: the "current" segment index in the "current" refseq
 760  Example :
 761  Returns : value of current_contig_seg_idx (a scalar)
 762  Args    : on set, new value (a scalar or undef, optional)
 763
 764 =cut
 765
 766 sub _current_contig_seg_idx {
 767     my $self = shift;
 768
 769     return $self->{'_current_contig_seg_idx'} = shift if @_;
 770     return $self->{'_current_contig_seg_idx'};
 771 }
 772
 773
 774 =head2 sam()
 775
 776  Title   : sam
 777  Usage   : $obj->sam($newval)
 778  Function: store the associated Bio::DB::Sam object
 779  Example :
 780  Returns : value of sam (a Bio::DB::Sam object)
 781  Args    : on set, new value (a scalar or undef, optional)
 782
 783 =cut
 784
 785 sub sam {
 786     my $self = shift;
 787     return $self->{'sam'};
 788 }
 789
 790 sub DESTROY {
 791     my $self = shift;
 792     undef $self->{'sam'};
 793     delete $self->{_segset}->{$_} foreach (keys %{$self->{_segset}});
 794 }
 795
 796 1;