Summary

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

Package variables

Included modules

Inherit

Bio::Matrix::PSM::IO Bio::Root::Root

Synopsis

No synopsis!

Description

Parser for meme.

Methods

Methods description

 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

Methods code

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

}

General documentation

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

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.

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

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:

  http://bugzilla.open-bio.org/

Email skirov@utk.edu