← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to Bio/Assembly/IO/sam.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Assembly/IO/sam.pm
187
187
BEGIN {
188
188
# check requirements
189
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)");
 
190
        Bio::Root::Root->throw(__PACKAGE__.' requires installation of samtools (libbam) and Bio::DB::Sam (available on CPAN; not part of BioPerl)');
191
191
    }
192
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.");
 
193
        Bio::Root::Root->warn("IO::Uncompress::Gunzip is not available; you'll have to do your decompression by hand.");
194
194
    }
195
195
}
196
196
 
231
231
    my $assembly = Bio::Assembly::Scaffold->new( -progname => $progname );
232
232
 
233
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
 
        }
 
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
245
    }
246
246
    return $assembly;
247
247
}
259
259
sub next_contig {
260
260
    my $self = shift;
261
261
    if (!defined $self->_current_contig_seg_idx) {
262
 
        $self->_current_contig_seg_idx(0);
 
262
        $self->_current_contig_seg_idx(0);
263
263
    }
264
264
    else {
265
 
        $self->_current_contig_seg_idx( 1+$self->_current_contig_seg_idx );
 
265
        $self->_current_contig_seg_idx( 1+$self->_current_contig_seg_idx );
266
266
    }
267
267
    unless ($self->_current_refseq_id) {
268
 
        croak("No current refseq id set");
 
268
        croak("No current refseq id set");
269
269
    }
270
270
    my $contig_segs = $self->_segset($self->_current_refseq_id);
271
271
    unless ($contig_segs && @$contig_segs) {
272
 
        croak("No contig segset for current id '".$self->_current_refseq_id."'")
 
272
        croak("No contig segset for current id '".$self->_current_refseq_id."'")
273
273
    }
274
274
    # each segment in @$contig_segs represents a contig within the
275
275
    # current refseq
285
285
    my $numseq = 0;
286
286
 
287
287
    while ( my $read = $seqio->next_seq ) {
288
 
        if ($self->{_assigned}->{$read->name}) {
289
 
            1;
290
 
            next;
291
 
        }
292
 
        $self->{_assigned}->{$read->name}=1;
293
 
        $self->_store_read($read, $contigobj);
294
 
        $numseq++;
 
288
        if ($self->{_assigned}->{$read->name}) {
 
289
            next;
 
290
        }
 
291
        $self->{_assigned}->{$read->name}=1;
 
292
        $self->_store_read($read, $contigobj);
 
293
        $numseq++;
295
294
    }
296
295
    if ($numseq == 1) { # ooh! a singlet!
297
 
        $contigobj = $self->_store_singlet($contigobj);
 
296
        $contigobj = $self->_store_singlet($contigobj);
298
297
    }
299
298
    return $contigobj;
300
299
}
381
380
    my $self = shift;
382
381
    my ($read, $contigobj) = @_;
383
382
    my $readseq = Bio::LocatableSeq->new(
384
 
        -display_id => $read->name,
385
 
        -primary_id => $read->name,
386
 
        -seq        => $read->dna,
387
 
        -start      => 1,
388
 
        -strand     => $read->strand,
389
 
        -alphabet   => 'dna'
390
 
        );
 
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
    );
391
390
    my $qual = Bio::Seq::PrimaryQual->new(
392
 
        -id         => $read->name."_qual",
393
 
        -qual       => [$read->qscore]
394
 
        );
 
391
        -id         => $read->name."_qual",
 
392
        -qual       => [$read->qscore]
 
393
    );
395
394
 
396
395
    # add pair information
397
396
    my @pair_info;
398
397
    if ($read->proper_pair) { # mate also aligned
399
 
        @pair_info = (
400
 
            mate_start => $read->mate_start,
401
 
            mate_len   => $read->mate_len,
402
 
            mate_strand => $read->mstrand,
403
 
            insert_size => $read->isize
404
 
            );
 
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
        );
405
404
    }
406
 
            
 
405
 
407
406
    my $alncoord = Bio::SeqFeature::Generic->new(
408
 
        -primary    => $read->name,
409
 
        -start      => $read->start,
410
 
        -end        => $read->end,
411
 
        -strand     => $read->strand,
412
 
        -qual       => join(' ',$read->qscore),
413
 
        -tag        => { 'contig' => $contigobj->id,
414
 
                         'cigar'  => $read->cigar_str,
415
 
                         @pair_info }
416
 
        );
 
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
    );
417
416
 
418
417
    $contigobj->set_seq_coord($alncoord, $readseq);
419
418
    $contigobj->set_seq_qual( $readseq, $qual );
455
454
    my ($readseq) = $contigobj->each_seq;
456
455
 
457
456
    my $singletobj = Bio::Assembly::Singlet->new( -id => $contigobj->id,
458
 
                                                  -seqref => $readseq );
 
457
                                                  -seqref => $readseq );
459
458
 
460
459
# may want to attach this someday
461
460
#    my $qual = $contigobj->get_qual_by_name($readseq->id);    
488
487
    $file =~ s/^[<>+]*//; # byebye parasitic mode chars
489
488
    my ($fname, $dir, $suf) = fileparse($file, ".sam", ".bam");
490
489
    $self->_set_from_args({ '_basename' => $fname }, 
491
 
                          -methods => [qw( _basename)],
492
 
                          -create => 1);
 
490
                          -methods => [qw( _basename)],
 
491
                          -create => 1);
493
492
    if (!defined $fasfile) {
494
 
        for (qw( fas fa fasta fas.gz fa.gz fasta.gz )) {
495
 
            $fasfile = File::Spec->catdir($dir, $self->_basename.$_);
496
 
            last if -e $fasfile;
497
 
            undef $fasfile;
498
 
        }
 
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
        }
499
498
    }
500
499
    unless (-e $fasfile) {
501
 
        croak( "Can't find associated reference fasta db" );
 
500
        croak( "Can't find associated reference fasta db" );
502
501
    }
503
502
    !$self->refdb && $self->refdb($fasfile);
504
503
    # compression
505
504
    if ($fasfile =~ /\.gz[^.]*$/) {
506
 
        unless ($HAVE_IO_UNCOMPRESS) {
507
 
            croak( "IO::Uncompress::Gunzip not available; can't decompress on the fly");
508
 
        }
509
 
        my ($tfh, $tf) = tempfile( UNLINK => 1);
510
 
        my $z = IO::Uncompress::Gunzip->new($fasfile) or croak("Can't expand: $@");
511
 
        while (<$z>) { print $tfh $_ }
512
 
        close $tfh;
513
 
        $fasfile = $tf;
 
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;
514
513
    }
515
514
    if ($file =~ /\.gz[^.]*$/) {
516
 
        unless ($HAVE_IO_UNCOMPRESS) {
517
 
            croak( "IO::Uncompress::Gunzip not available; can't decompress on the fly");
518
 
        }
519
 
        my ($tfh, $tf) = tempfile( UNLINK => 1);
520
 
        my $z = IO::Uncompress::Gunzip->new($file) or croak("Can't expand: $@");
521
 
        while (<$z>) { 
522
 
            print $tfh $_;
523
 
            1;
524
 
        }
525
 
        close $tfh;
526
 
        $file = $tf;
 
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;
527
525
    }
528
526
    # sam conversion : just barf for now
529
527
    if (-T $file) {
530
 
        my $bam = $file;
531
 
        $bam =~ s/\.sam/\.bam/;
532
 
        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");
 
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");
533
531
    }
534
532
 
535
533
    $sam = Bio::DB::Sam->new( -bam => $file, 
536
 
                              -fasta => $fasfile,
537
 
                  -autoindex  => 1,
538
 
                              -expand_flags => 1);
 
534
                              -fasta => $fasfile,
 
535
                              -autoindex  => 1,
 
536
                              -expand_flags => 1);
539
537
    unless (defined $sam) {
540
 
        croak( "Couldn't create the Bio::DB::Sam object" );
 
538
        croak( "Couldn't create the Bio::DB::Sam object" );
541
539
    }
542
540
    $self->{sam} = $sam;
543
541
    # now produce the contig segments for each seq_id...
544
542
    for ($sam->seq_ids) {
545
 
        my $seg = $sam->segment(-seq_id=>$_, -start=>1, -end=>$sam->length($_));
546
 
        ${$self->{_segset}}{$_} = [$self->_get_contig_segs_from_coverage($seg)];
 
543
        my $seg = $sam->segment(-seq_id=>$_, -start=>1, -end=>$sam->length($_));
 
544
        ${$self->{_segset}}{$_} = [$self->_get_contig_segs_from_coverage($seg)];
547
545
    }
548
546
 
549
547
    return 1;
566
564
    my $self = shift;
567
565
    my $segment = shift;
568
566
    unless ($self->sam) {
569
 
        croak("Sam object not yet initialized (call _init_sam)");
 
567
        croak("Sam object not yet initialized (call _init_sam)");
570
568
    }
571
569
    unless ( ref($segment) =~ /Bio::DB::Sam::Segment/ ) {
572
 
        croak("Requires Bio::DB::Sam::Segment object at arg 1");
 
570
        croak("Requires Bio::DB::Sam::Segment object at arg 1");
573
571
    }
574
572
    my ($cov, @covdata, @rngs, @segs);
575
573
    ($cov) = $segment->features('coverage');
576
574
    unless ($cov) {
577
 
        croak("No coverage data!");
 
575
        croak("No coverage data!");
578
576
    }
579
577
    @covdata = $cov->coverage;
580
578
    
582
580
    my $in_contig;
583
581
    my ($lf_end,$rt_end);
584
582
    for (0..$#covdata) {
585
 
        if ($covdata[$_]) {
586
 
            if ($in_contig) {
587
 
                $rt_end = $_+1;
588
 
                next;
589
 
            }
590
 
            else {
591
 
                $in_contig = 1;
592
 
                # push previous range
593
 
                if (defined $lf_end && defined $rt_end) {
594
 
                    push @rngs, [$lf_end, $rt_end];
595
 
                }
596
 
                $lf_end = $_+1;
597
 
            }
598
 
            
599
 
        }
600
 
        else {
601
 
            $in_contig = 0;
602
 
        }
 
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
        }
603
600
    }
604
601
    # last one
605
602
    push @rngs, [$lf_end, $rt_end] if (defined $lf_end and 
606
 
                                       defined $rt_end and
607
 
                                       $lf_end <= $rt_end);
 
603
                                       defined $rt_end and
 
604
                                       $lf_end <= $rt_end);
608
605
    unless (@rngs) {
609
 
        carp ("No coverage for this segment!");
610
 
        return;
 
606
        carp ("No coverage for this segment!");
 
607
        return;
611
608
    }
612
609
    for (@rngs) {
613
 
        push @segs, $self->sam->segment(-seq_id=>$segment->seq_id,
614
 
                                        -start=>$$_[0], 
615
 
                                        -end=>$$_[1]);
 
610
        push @segs, $self->sam->segment(-seq_id=>$segment->seq_id,
 
611
                                        -start=>$$_[0], 
 
612
                                        -end=>$$_[1]);
616
613
    }
617
614
    return @segs;
618
615
}
637
634
    my @quals;
638
635
    my $region = $seg->seq_id.':'.$seg->start.'..'.$seg->end;
639
636
    my $qual_averager = sub {
640
 
        my ($seqid, $pos, $p) = @_;
 
637
        my ($seqid, $pos, $p) = @_;
641
638
        return unless ($seg->start <= $pos and $pos <= $seg->end);
642
 
        my $acc = 0;
643
 
        my $n = 0;
644
 
        for my $pileup (@$p) {
645
 
            my $b = $pileup->alignment;
646
 
            $acc += $b->qscore->[$pileup->qpos];
647
 
            $n++;
648
 
        }
649
 
        push @quals, int($acc/$n);
 
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);
650
647
    };
651
648
    $self->sam->pileup($region, $qual_averager);
652
649
    return @quals;
673
670
    my $region = $seg->seq_id.':'.$seg->start.'-'.$seg->end;
674
671
 
675
672
    my $weighted_consensus = sub {
676
 
        my ($seqid, $pos, $p) = @_;
 
673
        my ($seqid, $pos, $p) = @_;
677
674
        return unless ($seg->start <= $pos and $pos <= $seg->end);
678
 
        my %wt_tbl;
679
 
        my %n;
680
 
        for my $pileup (@$p) {
681
 
            my $b = $pileup->alignment;
682
 
            my $res = substr($b->qseq,$pileup->qpos,1);
683
 
            $wt_tbl{$res} += $b->qscore->[$pileup->qpos] || 0;
684
 
            $n{$res} ||= 0;
685
 
            $n{$res}++;
686
 
        }
687
 
        # really simple
688
 
        my $c = (sort { $wt_tbl{$b}<=>$wt_tbl{$a} } keys %wt_tbl)[0];
689
 
        $conseq .= $c;
690
 
        push @quals, int($wt_tbl{$c}/$n{$c});
691
 
 
 
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});
692
688
    };
693
689
 
694
690
    $self->sam->pileup($region, $weighted_consensus);
695
691
    return Bio::Seq::Quality->new(
696
 
        -qual => join(' ', @quals ),
697
 
        -seq => $conseq,
698
 
        -id => $region
699
 
        );
 
692
        -qual => join(' ', @quals ),
 
693
        -seq => $conseq,
 
694
        -id => $region
 
695
    );
700
696
}
701
697
 
702
698
=head2 refdb()