Bio::Matrix::PSM::IO masta
SummaryIncluded librariesPackage variablesDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::Matrix::PSM::IO::masta - motif fasta format parser
Package variables
Globals (from "use vars" definitions)
@HEADER
Included modules
Bio::Matrix::PSM::SiteMatrix
strict
Inherit
Bio::Matrix::PSM::IO Bio::Root::Root
Synopsis
No synopsis!
Description
Parser for meme.
Methods
newDescriptionCode
write_psmDescriptionCode
next_matrixDescriptionCode
next_psmDescriptionCode
_make_matrix
No description
Code
_rearrange_matrix
No description
Code
_count_positions
No description
Code
_freq_to_count
No description
Code
Methods description
newcode    nextTop
 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);
write_psmcodeprevnextTop
 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 :
next_matrixcodeprevnextTop
  Title   : next_matrix
Usage : my $matrix = $psmio->next_matrix;
Function: Alias of next_psm function
next_psmcodeprevnextTop
 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
newdescriptionprevnextTop
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;
}
write_psmdescriptionprevnextTop
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);
	}
    }
}
next_matrixdescriptionprevnextTop
sub next_matrix {
     shift->next_psm(@_);
}
next_psmdescriptionprevnextTop
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);
}
_make_matrixdescriptionprevnextTop
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);
}
_rearrange_matrixdescriptionprevnextTop
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;
}
_count_positionsdescriptionprevnextTop
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};
}
_freq_to_countdescriptionprevnextTop
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;
}
General documentation
SYNOPSIS Top
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
FEEDBACKTop
Mailing ListsTop
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
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
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.
Reporting BugsTop
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://redmine.open-bio.org/projects/bioperl/
AUTHOR - Stefan KirovTop
Email skirov@utk.edu
APPENDIXTop