[Bioperl-l] not able to use Bio::Root::IO method

Chris Fields cjfields at illinois.edu
Thu Jan 14 07:23:41 UTC 2010


You can remove separate 'use' directives if they are declared with 'use base' (they will be imported then).  Also, Bio::Root::IO inherits Bio::Root::Root, and Bio::Assembly::IO should inherit from Bio::Root::IO, so the only base module you should need is Bio::Assembly::IO.  It's possible having all three is confusing the interpreter.

chris

On Jan 13, 2010, at 11:35 PM, Dan Kortschak wrote:

> Thanks Mark, I'm not sure about that since @ISA still includes
> Bio::Root:IO when it's at the call, but it might be.
> 
> cheers
> Dan
> 
> Here is the entirety of the code (it reasonably short):
> 
> package Bio::Assembly::IO::bowtie;
> use strict;
> use warnings;
> 
> # Object preamble - inherits from Bio::Root::Root
> 
> use Bio::SeqIO;
> use Bio::Tools::Run::Samtools;
> use Bio::Assembly::IO;
> use Carp;
> use Bio::Root::Root;
> use Bio::Root::IO;
> use base qw( Bio::Root::Root Bio::Root::IO Bio::Assembly::IO );
> 
> our $HD = "\@HD\tVN:1.0\tSO:unsorted\n";
> our $PG = "\@PG\tID=Bowtie\n";
> 
> our $HAVE_IO_UNCOMPRESS;
> BEGIN {
> # check requirements
>    unless ( eval "require Bio::Tools::Run::Bowtie;") {
> 	Bio::Root::Root->throw("Bio::Tools::Run::Bowtie is not available - cannot extract refdb from index.");
>    }
>    unless ( eval "require IO::Uncompress::Gunzip; \$HAVE_IO_UNCOMPRESS = 1") {
> 	Bio::Root::Root->warn("IO::Uncompress::Gunzip is not available; you'll have to do your decompression by hand.");
>    }
> }
> 
> sub new {
> 	my $class = shift;
> 	my @args = @_;
> 	my $self = $class->SUPER::new(@args);
> 	my ($file, $index, $no_head, $no_sq) = $self->_rearrange([qw(FILE INDEX NO_HEAD NO_SQ)], @args);
> 	$file =~ s/^<//;
> 	$self->{'_no_head'} = $no_head;
> 	$self->{'_no_sq'} = $no_sq;
> 	# get the sequence so samtools can work with it
> 	my $inspector = Bio::Tools::Run::Bowtie->new( -command => 'inspect' );
> 	my $refdb = $inspector->run($index);
> 	my $bam_file = $self->_make_bam($self->_bowtie_to_sam($file, $refdb));
> 	my $sam = Bio::Assembly::IO->new( -file => "<$bam_file", -refdb => $refdb , -format => 'sam' );
> 	return $sam;
> }
> 
> sub _bowtie_to_sam {
> 	my ($self, $file, $refdb) = @_;
> 
> 	$self->throw("'$file' does not exist or is not readable.")
> 		unless ( -e $file && -r $file );
> 	my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$file);
> 	$self->throw("'$file' is not a bowtie formatted file.") unless $guesser->guess =~ m/^bowtie$/;
> 
> 	my %SQ;
> 	my $mapq = 255;
> 	my $in_pair;
> 	my @mate_line;
> 	my $mlen;
> 
> 	if ($file =~ m/\.gz[^.]*$/) {
> 		unless ($HAVE_IO_UNCOMPRESS) {
> 			croak( "IO::Uncompress::Gunzip not available, can't expand '$_'" );
> 		}
> 		my ($tfh, $tf) = $self->io->tempfile;
> 		my $z = IO::Uncompress::Gunzip->new($_);
> 		while (<$z>) { print $tfh $_ }
> 		close $tfh;
> 		$file = $tf;
> 	}
> 
>        open(my $fh, $file) or
> 		$self->throw("Can not open '$file' for reading: $!");
> 
> 	# create temp file for working
> 	my ($sam_tmp_h, $sam_tmp_f) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => '.sam' );
> 	
> 	while ($fh) {
> 		chomp;
> 		my ($qname,$strand,$rname,$pos,$seq,$qual,$m,$details)=split("\t",$_);
> 		$SQ{$rname} = 1;
> 		
> 		my $paired_f =  ($qname =~ m#/[12]#) ? 0x03 : 0;
> 		my $strand_f = ($strand eq '-') ? 0x10 : 0;
> 		my $op_strand_f = ($strand eq '+' && $paired_f) ? 0x20 : 0;
> 		my $first_f =  ($qname =~ m#/1#) ? 0x40 : 0;
> 		my $second_f = ($qname =~ m#/2#) ? 0x80 : 0;
> 		my $flag = $paired_f | $strand_f | $op_strand_f | $first_f | $second_f;
> 
> 		$pos++;
> 		my $len = length $seq;
> 		die unless $len == length $qual;
> 		my $cigar = $len.'M';
> 		my @detail = split(',',$details);
> 		my $dist = 'NM:i:'.scalar @detail;
> 		
> 		my @mismatch;
> 		my $last_pos = 0;
> 		for (@detail) {
> 			m/(\d+):(\w)>\w/;
> 			my $err = ($1-$last_pos);
> 			$last_pos = $1+1;
> 			push @mismatch,($err,$2);
> 		}
> 		push @mismatch, $len-$last_pos;
> 		@mismatch = reverse @mismatch if $strand eq '-';
> 		my $mismatch = join('',('MD:Z:', at mismatch));
> 
> 		if ($paired_f) {
> 			my $mrnm = '=';
> 			if ($in_pair) {
> 				my $mpos = $mate_line[3];
> 				$mate_line[7] = $pos;
> 				my $isize = $mpos-$pos-$len;
> 				$mate_line[8] = -$isize;
> 				print $sam_tmp_h join("\t", at mate_line),"\n";
> 				print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
> 				$in_pair = 0;
> 			} else {
> 				$mlen = $len;
> 				@mate_line = ($qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, undef, undef, $seq, $qual, $mismatch, $dist);
> 				$in_pair = 1;
> 			}
> 		} else {
> 			my $mrnm = '*';
> 			my $mpos = 0;
> 			my $isize = 0;
> 			print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
> 		}
> 	}
> 
> 	close($fh);
> 	$sam_tmp_h->close;
> 	
> 	return $sam_tmp_f if $self->{'_no_head'};
> 
> 	my ($samh, $samf) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => '.sam' );
> 
> 	# print header
> 	print $samh $HD;
> 	
> 	# print sequence dictionary
> 	unless ($self->{'_no_sq'}) {
> 		my $db  = Bio::SeqIO->new( -file => $refdb, -format => 'fasta' );
> 		while ( my $seq = $db->next_seq() ) {
> 			$SQ{$seq->id} = $seq->length if $SQ{$seq->id};
> 		}
> 	
> 		map { print $samh join("\t", ('@SQ', "SN:$_", "LN:$SQ{$_}")), "\n" } keys %SQ;
> 	}
> 	
> 	# print program
> 	print $samh $PG;
> 	
> 	open($sam_tmp_h, $sam_tmp_f) or
> 		$self->throw("Can not open '$sam_tmp_f' for reading: $!");
> 
> 	print $samh $_ while ($sam_tmp_h);
> 	
> 	close($sam_tmp_h);
> 	$samh->close;
> 	
> 	return $samf;
> }
> 
> sub _make_bam {
> 	my ($self, $file) = @_;
> 	
> 	$self->throw("'$file' does not exist or is not readable")
> 		unless ( -e $file && -r $file );
> 
> 	# make a sorted bam file from a sam file input
> 	my ($bamh, $bamf) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => '.bam' );
> 	my ($srth, $srtf) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => '.srt' );
> 	$_->close for ($bamh, $srth);
> 	
> 	my $samt = Bio::Tools::Run::Samtools->new( -command => 'view',
> 						   -sam_input => 1,
> 						   -bam_output => 1 );
> 
> 	$samt->run( -bam => $file, -out => $bamf );
> 
> 	$samt = Bio::Tools::Run::Samtools->new( -command => 'sort' );
> 
> 	$samt->run( -bam => $bamf, -pfx => $srtf);
> 
> 	return $srtf.'.bam'
> }
> 
> 1;
> 
> 
> On Thu, 2010-01-14 at 00:11 -0500, Mark A. Jensen wrote:
>> Hey Dan-- what does your constructor look like? I wonder if
>> something's getting 
>> lost in new() and _initialize() chaining spaghetti- MAJ
>> 
> 
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l





More information about the Bioperl-l mailing list