#--------------------------------------------------------- # $Id: masta.pm 16123 2009-09-17 12:57:27Z cjfields $ =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. bioperl-l@bioperl.org - 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: http://bugzilla.open-bio.org/ =head1 AUTHOR - Stefan Kirov Email skirov@utk.edu =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;