Bio::Matrix::PSM::IO::masta - motif fasta format parser

  0.1	0.62	.017	0.11
  0.22	0.13	0.54	0.11

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://redmine.open-bio.org/projects/bioperl/

 Title   : new
 Usage   : my $psmIO =  Bio::Matrix::PSM::IO->new(-format=> 'masta',
						 -file  => $file, 
                                                 -mtype => 'PWM');
 Function: Associates a file with the appropriate parser
 Throws  :
 Example :
 Args    : hash
 Returns : "Bio::Matrix::PSM::$format"->new(@args);

 Title   : write_psm
 Usage   : 
 Function: writes a pfm/pwm/raw sequence in a simple masta format
 Throws  :
 Example : 
 Args    : SiteMatrix object, type (optional string: PWM, SEQ or PFM)
 Returns : 

  Title   : next_matrix
  Usage   : my $matrix = $psmio->next_matrix;
  Function: Alias of next_psm function

 Title   : next_psm
 Usage   : my $matrix=$psmio->next_psm;
 Function: returns the next matrix in the stream
 Throws  : If there is you mix different types, for example weights and
           frequencies occur in the same entry You can mix weights, but these
           should be designated by different ID lines
 Example :
 Args    :
 Returns : Bio::Matrix::PSM::SiteMatrix

#---------------------------------------------------------

# Let the code begin...
package Bio::Matrix::PSM::IO::masta;
use Bio::Matrix::PSM::SiteMatrix;
use vars qw(@HEADER);
use strict;

use base qw(Bio::Matrix::PSM::IO Bio::Root::Root);

sub new {
    my($class, @args)=@_;
    my $self = $class->SUPER::new(@args);
    my ($file)=$self->_rearrange(['FILE'], @args);
    my ($query,$tr1)=split(/\./,$file,2);
    $self->{file}  = $file;
    $self->{_end}  = 0;
    $self->{mtype} = uc($self->_rearrange(['MTYPE'], @args) || "PFM");
    $self->_initialize_io(@args) || $self->warn("Did you intend to use STDIN?"); #Read only for now
    return $self;
}

sub write_psm {
    my ($self,$matrix,$type)=@_;
    $self->{mtype} = uc($type) if ($type);
    my $idline=">". $matrix->id . "\n";
    $self->_print($idline);
    unless ($self->{mtype} eq 'SEQ') {
	while (my %h=$matrix->next_pos) {
	    my $row=$self->{mtype} eq 'PWM' ? join("\t",$h{lA},$h{lC},$h{lG},$h{lT},"\n"):join("\t",$h{pA},$h{pC},$h{pG},$h{pT},"\n");
	    $self->_print ($row);
	}
    } else {
	my @seq;
	while (my %h=$matrix->next_pos) {
	    my ($a,$c,$g,$t)=_freq_to_count(\%h);
	    $self->throw("Could not convert from frequency to count\n") if (($a+$c+$g+$t) !=10);
	    for my $i (0..$a-1) {$seq[$i].='A';}
	    my $m=$a+$c;
	    for my $i ($a..$m-1) {$seq[$i].='C';}
	    my $n=$a+$c+$g;
	    for my $i ($m..$n-1) {$seq[$i].='G';}
	    for my $i ($n..9) {$seq[$i].='T';}
	}	
	foreach my $s (@seq) {
	    $s.="\n";
	    $self->_print ($s);
	}
    }
}

sub next_matrix { 
    shift->next_psm(@_);
}

sub next_psm {
    my $self=shift;
    return if ($self->{_end});
    my $line=$self->_readline;
    $self->throw("No ID line- wrong format\n") unless ($line=~/^>/);
    my ($id,$desc)=split(/[\t\s]+/,$line,2);
    $id=~s/>//;
    my ($mtype,$format,@mdata,$len);
    $self->{_mtype} = 0;
    while ($line=$self->_readline) {
	next if $line =~ /^\s+$/;# There should not be empty lines, but just in case...
	chomp $line;
	if ($line =~ /^>/) {
	    $self->_pushback($line);
	    last;
	}

	if ($line !~ /[^ACGTacgt]/g) {
	    # This is a set of aligned sequences
	    $self->throw("Mixing between types is not allowed or a parsing error occured\n") 
		if (($self->{_mtype} != 3) && ($mtype)) ;
	    $self->throw("Bad sequence- different length: $line\n") 
		if (($len) && ($len!=length($line)));
	    $len=length($line) unless ($len);
	    push @mdata,$line;
	    $self->{_mtype}=3;
	} else {
		# do not strip 'e's since they are part of number notation for small/big numbers
	    $line=~s/[a-df-zA-DF-Z]//g; #Well we may wanna do a hash and auto check for letter order if there is a really boring talk...
	    $line=~s/^[\s\t]+//;
	    $line=~s/[\s\t]+/\t/g;
	    my @data=split(/[\s\t]+/,$line);
	    if ($#data==3) {
		$self->throw("Mixing between types is not allowed or a parsing error occured\n") if (($mtype)&&($self->{_mtype} !=1)) ;
		$self->{_mtype}=1;
		$mtype=1;
	    }
	    else   {
		$self->throw("Mixing between types is not allowedor a parsing error occured\n") if (($mtype)&&($self->{_mtype} !=2)) ;
		$self->{_mtype}=2;
		$mtype=1;
	    }
	    push @mdata,\@data;
	}
    }
    $self->{_end} = 1 if (!defined $line || $line !~ /^>/);
    return _make_matrix(\@mdata,$self->{_mtype},$id,$desc);
}

sub _make_matrix {
    my ($mdata,$type,$id,$desc)=@_;
    if ($type==1) {
	my @rearr=_rearrange_matrix($mdata); 
	$mdata=\@rearr;
    }
#Auto recognition for what type is this entry (PFM, PWM or simple count)
#A bit dangerous, I hate too much auto stuff, but I want to be able to mix different
#types in a single file
    my $mformat='count';
    my ($a,$c,$g,$t);
    if ($type == 3 ) {
	($a,$c,$g,$t)= &_count_positions($mdata);
    } else {
	($a,$c,$g,$t)=@{$mdata};	
	my $k=$a->[0]+$c->[0]+$g->[0]+$t->[0];
	my $l= ($a->[0]+$c->[0]+$g->[0]+$t->[0]) - 
	    (abs($a->[0])+abs($c->[0])+abs($g->[0])+abs($t->[0]));
	$mformat='freq' if (($k==1) && ($l==0));
	$mformat='pwm' if ($l!=0);
    }
    my (@fa,@fc,@fg,@ft,%mparam);

    if ($mformat eq 'pwm') {
	foreach my $i (0..$#{$a}) {
	    my $ca=exp $a->[$i];
	    my $cc=exp $c->[$i];
	    my $cg=exp $g->[$i];
	    my $ct=exp $t->[$i];
	    my $all=$ca+$cc+$cg+$ct;
	    push @fa,($ca/$all)*100;
	    push @fc,($cc/$all)*100;
	    push @fg,($cg/$all)*100;
	    push @ft,($ct/$all)*100;
	}
    }
    $desc.=", source is $mformat";
    if ($mformat eq 'pwm') {
	$desc=~s/^pwm//;
	%mparam=(-pA=>\@fa,-pC=>\@fc,-pG=>\@fg,-pT=>\@ft,-id=>$id,-desc=>$desc,
		 -lA=>$a,-lC=>$c,-lG=>$g,-lT=>$t);
    }
    else {
	%mparam=(-pA=>$a,-pC=>$c,-pG=>$g,-pT=>$t,-id=>$id,-desc=>$desc);
    }
    return new Bio::Matrix::PSM::SiteMatrix(%mparam);
}

sub _rearrange_matrix {
    my $mdata=shift;
    my (@a,@c,@g,@t);
    foreach my $entry (@{$mdata}) {
	my ($a,$c,$g,$t)=@$entry;
	push @a,$a;
	push @c,$c;
	push @g,$g;
	push @t,$t;
    }
    return \@a,\@c,\@g,\@t;
}


sub _count_positions {
    my $seq=shift;
    my %pos;
    my $l=length($seq->[0])-1;
    for( my $i = 0; $i <= $l; $i++ ) {
	for ( qw(A C G T) ) {
	    $pos{$_}->[$i] = 0;
	}
    }
    foreach my $sequence (@{$seq}) {
	my @let= split(//,$sequence);
	for my $i (0..$#let) {
	    $pos{uc($let[$i])}->[$i]++;
	}
    }
    return $pos{A},$pos{C},$pos{G},$pos{T};
}


sub _freq_to_count {
    my $h=shift;
    my $a=int(10*$h->{pA}+0.5);
    my $c=int(10*$h->{pC}+0.5);
    my $g=int(10*$h->{pG}+0.5);
    my $t=int(10*$h->{pT}+0.5);
    return ($a,$c,$g,$t);
}

1;