File Coverage

lib/Fry/Lib/BioPerl.pm

Criterion	Covered	Total	%
statement	7	9	77.7
branch			n/a
condition			n/a
subroutine	3	3	100.0
pod			n/a
total	10	12	83.3

package Fry::Lib::BioPerl;

use strict;
#use Benchmark qw/:all/;
use Time::HiRes qw/tv_interval gettimeofday/;
use Bio::SeqIO;
use Bio::AlignIO;
use Bio::Tools::Run::Alignment::Clustalw;
use Bio::Tools::Run::Alignment::TCoffee;
use Bio::DB::GenBank;
#use Bio::Graphics;
#Bio::Seq::RichSeq;
#Bio::Graphics::Feature;

our $VERSION = 0.15;
our $Verbose = 1;
$Bio::Root::Root::DEBUG = 0;

#these are defaults
our $seqio_file = "$ENV{HOME}/.cpanplus/5.8.1/build/bioperl-1.4/t/data/test.genbank";
#our $seqio_file = "t/data/test.genbank";
our $seqio_format = "genbank";
our $alnio_format =  "clustalw";
our $alnio_file = "$ENV{HOME}/.cpanplus/5.8.1/build/bioperl-1.4/t/data/testaln.aln";
#our $alnio_file = "t/data/testaln.aln";
our $default_seqio = "seqio1";
our $default_alnio = "alnio1";
our %clustal_params = (qw/matrix BLOSUM/);
our %tcoffee_params = ();
our %db_genbank_params = ();
my $seqcount = 0;
my $alncount = 0;
#hackish variable
our %temp;

sub _default_data {
	return {
		cmds=>{
			readSeq=>{qw/a brs/,d=>'Read a sequence from a file',
				u=>'$format$file$obj'},
			readAln=>{qw/a bra/,d=>'Read an alignment from a file',
				u=>'$format$file$obj'},
			startClustalw=>{qw/a bsc/,d=>'Initialize clustalw',u=>'%options'},
			startTCoffee=>{qw/a bst/,d=>'Initialize tcoffee',u=>'%options'},
			startDB=>{qw/a bsd/,d=>'Initialize database',u=>'%options'},
			clustalAlign=>{qw/a bca/,d=>'Align with clustalw',u=>'@seq'},
			tcoffeeAlign=>{qw/a bta/,d=>'Align with tcoffee',u=>'@seq'},
			compareAlign=>{qw/a bcom/,d=>'Compare two alignment methods',
				u=>'@seq'},
			alnioFormat=>{qw/a baf/,d=>'Change alnioout format',u=>'$format'},
			seqioFormat=>{qw/a bsf/,d=>'Change seqioout format',u=>'$format'},
			listSeqs=>{qw/a bls/,d=>'Read all sequences from a file',
				u=>'$seqio'},
			#nextSeq=>{qw/a ,/,d=>'Get the next Seq'},
			writeAln=>{qw/a bwa/,d=>'Write alignments',u=>'@aln'},
			displayAln=>{qw/a bda/,d=>'Display alignment info',u=>'@aln'},
			getAlnSeqs=>{qw/a bas/,d=>'Save alignment seqs',u=>'$aln'},
			writeSeq=>{qw/a bws/,d=>'Write sequences',u=>'@seq'},
			translateSeq=>{qw/a bts/,d=>'Translate sequences',u=>'$seq'},
			revcomSeq=>{qw/a bts/,d=>'Reverse complement sequences',u=>'$seq'},
			getSeqById=>{qw/a bgsi/,d=>'Get sequence by id',u=>'$id'},
			#featureImage=>{qw/a bfi/,d=>''},
			libCmdAttr=>{qw/a lca/,d=>"Display a library's cmds' attributes",
				u=>'$attr$lib'},
		},
		objs=>{
			#seqio1=>{}
		}
	}
}
sub _initLib {
	my ($cls) = @_;
	#reads in sequence file
	$cls->newSeqIoInObj;
	#sequence output
	$cls->newSeqIoObj('seqioout',-fh=>\*STDOUT,-format=>$seqio_format);
	#alignment output
	$cls->newAlnIoObj('alnioout',-fh=>\*STDOUT,-format=>$alnio_format);
}

#new*Obj
	#these methods create objects and index them as Fry::Obj objects as well
	sub newSeqIoInObj {
		my ($cls,%arg) = @_;

		#defaults
		my $objId = delete $arg{obj} || $default_seqio;
		$arg{-file} = $arg{-file} || $seqio_file;
		$arg{-format} = $arg{-format} || $seqio_format;

		$cls->newSeqIoObj($objId,%arg);

		#index seqs of seqio
		my ($id,@ids);
		my @newSeqs = $cls->ripStream($cls->obj->get($objId,'obj'));
		for (@newSeqs) {
			$id = $cls->newSeqObj($_);
			push(@ids,$id);
		}
		$cls->obj->set($objId,seqs=>\@ids);
	}
	sub newSeqIoObj {
		my ($cls,$objId,%arg) = @_;
		#print Dumper \%arg;

		my $obj =  $cls->seqioNew(%arg) ;
		$cls->obj->setOrMake($objId=>$obj);
	}
	sub newSeqObj {
		my ($cls,$obj) = @_;
		$seqcount++;
		my $id = "seq$seqcount";
		$cls->obj->manyNew($id=>{obj=>$obj});
		$cls->view("Created sequence '$id'\n") if ($Verbose);
		return $id;
	}
	sub newAlnObj {
		my ($cls,$obj) = @_;
		$alncount++;
		my $id = "aln$alncount";
		$cls->obj->manyNew($id=>{obj=>$obj});
		$cls->view("Created alignment '$id'\n") if ($Verbose);
		return $id;
	}
	sub newAlnIoObj {
		my ($cls,$objId,%arg) = @_;
		my $obj =  $cls->alnioNew(%arg) ;
		$cls->obj->setOrMake($objId=>$obj);
	}
	sub newAlnIoInObj {
		my ($cls,%arg) = @_;
		my $objId = delete $arg{obj} || $default_alnio;
		$arg{-file} = $arg{-file} || $alnio_file;
		delete $arg{-format};
		#$arg{-format} = $arg{-format} || "clustalw";

		$cls->newAlnIoObj($objId,%arg);

		#index seqs of seqio
		my ($id,@ids);
		my @newAlns = $cls->ripStream($cls->obj->get($objId,'obj'));
		for (@newAlns) {
			$id = $cls->newAlnObj($_);
			push(@ids,$id);
		}
		$cls->obj->set($objId,alns=>\@ids);
	}
#wrappers around classes to easily substitute your own subclasses of any of these
	sub seqioNew { shift; return Bio::SeqIO->new(@_) }
	sub alnioNew { shift; return Bio::AlignIO->new(@_) }
	sub biodbNew { shift; return Bio::DB::GenBank->new(@_) }
	sub clustalwNew { shift; return Bio::Tools::Run::Alignment::Clustalw->new(@_) }
	sub tcoffeeNew { shift; return Bio::Tools::Run::Alignment::TCoffee->new(@_) }
	sub objObj { shift->obj->get($_[0],'obj') }
#commands
	#h: will probably break in other releases
	sub libCmdAttr {
		my ($cls,$attr) = @_;
		my $lib = "Fry::Lib::BioPerl";

		my @ids =  @{$cls->lib->get($lib,'cmds') };
		$cls->printGeneralAttr('cmd',$attr,@ids);
	}
#constructors
	sub readSeq {
		#d: initializes $seqio
		my ($cls,$format,$file,$objId) = @_;
		$cls->newSeqIoInObj(-file=>$file,-format=>$format,obj=>$objId);
	}
	sub readAln {
		my ($cls,$format,$file,$objId) = @_;
		$cls->newAlnIoInObj(-file=>$file,-format=>$format,obj=>$objId);
	}
	sub startClustalw {
		my ($cls,%arg) = @_;
		#$cls->verify_params(\%arg);
		my %params = (%clustal_params,%arg);
		my $obj = $cls->clustalwNew(%params);
		$cls->obj->manyNew('cw'=>{obj=>$obj});
	}
	sub startDB {
		my ($cls,%arg) = @_;
		my %params = (%db_genbank_params,%arg);
		my $obj = $cls->biodbNew(%params);
		$cls->obj->manyNew('db1'=>{obj=>$obj});
	}
	sub startTCoffee {
		my ($cls,%arg) = @_;
		my %params = (%tcoffee_params,%arg);
		my $obj = $cls->tcoffeeNew(%params);
		$cls->obj->manyNew('tc'=>{obj=>$obj});
	}
#Alignment
	##$tcoffee
	sub tcoffeeAlign {
		my ($cls,@seqs) = @_;
		@seqs = $cls->idToObj(@seqs);

		my $t0= [ gettimeofday()];
		my $alnobj = $cls->obj->get('tc','obj')->align(\@seqs);
		my $elapsed = tv_interval($t0);
		$cls->newAlnObj($alnobj);
		$temp{tc} = {time=>$elapsed,alnid=>"aln$alncount"};

		$cls->view("This alignment took $elapsed seconds\n") if $Verbose;
	}
	##$clustalw
	sub clustalAlign {
		my ($cls,@seqs) = @_;
		@seqs = $cls->idToObj(@seqs);

		my $t0= [ gettimeofday()];
		my $alnobj = $cls->obj->get('cw','obj')->align(\@seqs);
		my $elapsed = tv_interval($t0);
		$cls->newAlnObj($alnobj);
		$temp{cw} = {time=>$elapsed,alnid=>"aln$alncount"};

		$cls->view("This alignment took $elapsed seconds\n") if $Verbose;
	}
	##macros
	sub compareAlign {
		my ($cls,@seqs) = @_;
		@seqs = $cls->idToObj(@seqs);

		$cls->clustalAlign(@seqs);
		$cls->tcoffeeAlign(@seqs);

		$cls->view("Alignment from Clustalw\n");
		$cls->view("------------------------------------\n");
		$cls->displayAln($temp{cw}{alnid});
		$cls->view("Time: ",$temp{cw}{time},"\n");
		$cls->view("\nAlignment from TCoffee\n");
		$cls->view("------------------------------------\n");
		$cls->displayAln($temp{tc}{alnid});
		$cls->view("Time: ",$temp{tc}{time},"\n");
	}
#IO
	##$alignio
	sub alnioFormat {
		my ($cls,$format) = @_;
		$format ||= $alnio_format;
		$cls->newAlnIoObj('alnioout',-fh=>\*STDOUT,-format=>$format);
	}
	##$seqio
	sub seqioFormat {
		my ($cls,$format) = @_;
		$format ||= $seqio_format;
		$cls->newSeqIoObj('seqioout',-fh=>\*STDOUT,-format=>$format);
	}
#Sequences
	##$biodb
	sub getSeqById {
		my ($cls,$id) = @_;
		my $method = "get_Seq_by_id";
		my $seq = $cls->objObj('db1')->$method($id);
		$cls->newSeqObj($seq);
	}
	##$seq
	sub listSeqs {
		#d: readSeqs
		my ($cls,$seqioId) = @_;
		$seqioId ||= $default_seqio;
		my @seqs = @{$cls->obj->get($seqioId,'seqs') || []};
		my $output = "Sequence ids of $seqioId are:\n";
		$output .= $cls->seqIds(map {$cls->obj->get($_,'obj')} @seqs);
		$cls->view($output);
	}
	sub translateSeq {
		my ($cls,$seqId) = @_;
		my $seqObj = $cls->obj->get($seqId,'obj')->translate;
		$cls->newSeqObj($seqObj);
	}
	sub revcomSeq {
		my ($cls,$seqId) = @_;
		my $seqObj = $cls->obj->get($seqId,'obj')->revcom;
		$cls->newSeqObj($seqObj);
	}
	sub writeSeq {
		my ($cls,@seqs) = @_;
		@seqs = $cls->idToObj(@seqs);
		for my $seq (@seqs) {
			$cls->obj->get(seqioout=>'obj')->write_seq($seq);
		}
	}
#Alignments
	##$aln
	#add,remove,indexSeq,slice,remove_gaps,map_chars
	sub displayAln {
		my ($cls,@alns) = @_;
		@alns = $cls->idToObj(@alns);
		my %methods = (score=>'Score',consensus_string=>'Consensus string',
			length=>'Length',no_residues=>'Number of residues',
			no_sequences=>'Number of sequences',percentage_identity=>'Percentage identity');
		my $output;

		for my $aln (@alns) {
			for my $m (sort keys %methods) {
				$output .= $methods{$m}.": ". $aln->$m
				."\n";
			}
		}
		$cls->view($output);
	}
	sub getAlnSeqs {
		my ($cls,$alnId) = @_;
		my $aln = $cls->objObj($alnId);
		my @ids;

		for my $seq ($aln->each_seq) { push(@ids,$cls->newSeqObj($seq)) }
		#for my $seq ($aln->each_seq) {
			#my $obj = Bio::Seq->new(-display_id=>$seq->{display_id},-seq=>$seq->{seq});
			##print "s: ",$seq->{seq},"\n";
			#push(@ids,$cls->newSeqObj($obj))
		#}
		$cls->obj->set($alnId,seqs=>\@ids);
	}
	sub writeAln {
		my ($cls,@alns) = @_;
		@alns = $cls->idToObj(@alns);
		for my $aln (@alns) {
			$cls->obj->get(alnioout=>'obj')->write_aln($aln);
		}
	}
	##misc
	sub seqIds {
		my ($cls,@seqs) = @_;
		my $output;

		if (@seqs > 0) {
			for my $seq  (@seqs) {
				$output .= sprintf "%s\n", $seq->display_id;
			}
		}
		else { $output = "No sequences to display." }
		return $output;
	}
	sub idToObj {
		my ($cls,@idOrObj) = @_;
		my @return;
		if (ref $idOrObj[0]) {
			@return = @idOrObj;
		}
		else { @return = map { $cls->obj->get($_,'obj') } @idOrObj }
		wantarray ? @return : $return[0];
	}
	#d: all Bio::*IO can use this
	sub ripStream { my $fh = $_[1]->fh; my @seqs = <$fh>; return @seqs }
1;
__END__