File Coverage

Bio/Matrix/PSM/IO/masta.pm

Criterion	Covered	Total	%
statement	139	142	97.8
branch	36	46	78.2
condition	11	20	55.0
subroutine	12	12	100.0
pod	4	4	100.0
total	202	224	90.1

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

=head1 NAME

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

=head1 SYNOPSIS

MASTA is a position frequency matrix format similar to
fasta. It contains one ID row just like fasta and then the actual
data, which is tab delimited:

  0.1	0.62	.017	0.11
  0.22	0.13	0.54	0.11

Or A,C,G and T could be horizontally positioned (positioning is
automatically detected).  Please note masta will parse only DNA at the
moment.

It will also convert a set of aligned sequences:
ACATGCAT
ACAGGGAT
ACAGGCAT
ACCGGCAT

to a PFM (SiteMatrix object). When writing if you supply SEQ it will
write 10 random instances, which represent correctly the frequency and
can be used as an input for weblogo creation purposes.

See Bio::Matrix::PSM::IO for detailed documentation on how to use masta parser

=head1 DESCRIPTION

Parser for meme.

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this
and other Bioperl modules. Send your comments and suggestions preferably
 to one of the Bioperl mailing lists.
Your participation is much appreciated.

  [e-mail]                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

=head2 Support

Please direct usage questions or support issues to the mailing list:

I

rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution.  Bug reports can be submitted via the
web:

  https://github.com/bioperl/bioperl-live/issues

=head1 AUTHOR - Stefan Kirov

Email [e-mail]

=head1 APPENDIX

=cut


# 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);



=head2 new

 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);

=cut

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;
}

=head2 write_psm

 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 :

=cut

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);
	}
    }
}

=head2 next_matrix

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

=cut

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

=head2 next_psm

 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

=cut

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;