← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/Matrix/PSM/IO/meme.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Matrix/PSM/IO/meme.pm
1
1
#---------------------------------------------------------
2
 
# $Id: meme.pm,v 1.7 2003/12/12 17:01:42 skirov Exp $
 
2
# $Id: meme.pm,v 1.20.4.1 2006/10/02 23:10:22 sendu Exp $
3
3
 
4
4
=head1 NAME
5
5
 
6
 
Bio::Matrix::PSM::meme - PSM meme parser implementation
 
6
Bio::Matrix::PSM::IO::meme - PSM meme parser implementation
7
7
 
8
8
=head1 SYNOPSIS
9
9
 
22
22
 to one of the Bioperl mailing lists.
23
23
Your participation is much appreciated.
24
24
 
25
 
  bioperl-l@bioperl.org                 - General discussion
26
 
  http://bio.perl.org/MailList.html             - About the mailing lists
 
25
  bioperl-l@bioperl.org                  - General discussion
 
26
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
27
27
 
28
28
=head2 Reporting Bugs
29
29
 
30
30
Report bugs to the Bioperl bug tracking system to help us keep track
31
 
 the bugs and their resolution.
32
 
 Bug reports can be submitted via email or the web:
 
31
the bugs and their resolution.  Bug reports can be submitted via the
 
32
web:
33
33
 
34
 
  bioperl-bugs@bio.perl.org
35
 
  http://bugzilla.bioperl.org/
 
34
  http://bugzilla.open-bio.org/
36
35
 
37
36
=head1 AUTHOR - Stefan Kirov
38
37
 
45
44
 
46
45
# Let the code begin...
47
46
package Bio::Matrix::PSM::IO::meme;
48
 
use Bio::Matrix::PSM::IO;
49
47
use Bio::Matrix::PSM::InstanceSite;
50
48
use Bio::Matrix::PSM::SiteMatrix;
51
49
use Bio::Matrix::PSM::Psm;
52
 
use Bio::Matrix::PSM::PsmHeader;
53
 
use vars qw(@ISA @HEADER);
 
50
use vars qw(@HEADER);
54
51
use strict;
55
52
 
56
 
@ISA=qw(Bio::Matrix::PSM::PsmHeader Bio::Matrix::PSM::IO Bio::Root::Root);
 
53
use base qw(Bio::Matrix::PSM::PsmHeader Bio::Matrix::PSM::IO);
57
54
 
58
55
@Bio::Matrix::PSM::IO::meme::HEADER = qw(e_val sites IC width);
59
56
 
79
76
    $self->{file} = $file;
80
77
    $self->{query}= $query;
81
78
    $self->{end}  = 0;
 
79
    $self->{_strand}=0; #This we'll need to see if revcom option is used
82
80
    $self->_initialize_io(@args) || warn "Did you intend to use STDIN?"; #Read only for now
83
81
    #Skip header
84
82
    my $line;
114
112
    my $line=$self->_readline;
115
113
    while ($line !~ /^\*{10,}/ ) {
116
114
        chomp $line;
117
 
        $line=~s/\s+/,/g;
 
115
        $line =~ s/\s+/,/g;
118
116
        my ($id1,$w1,$l1,$id2,$w2,$l2)=split(/,/,$line);
119
117
        push @{$self->{hid}},$id1;
120
118
        $self->{weight}->{$id1}=$w1;
178
176
sub next_psm {
179
177
    #Parses the next prediction and returns a psm objects
180
178
    my $self=shift;
181
 
    return undef if ($self->{end});
182
 
    my ($endm,$line,$instances,$tr,$width,$motif_id,$sites,$e_val,$id,$ic);
 
179
    return if ($self->{end});
 
180
    my ($endm,$line,$instances,$tr,$width,$motif_id,$sites,$e_val,$id,$ic,$lA,$lC,$lG,$lT);
183
181
    while (defined( $line = $self->_readline) ) {
 
182
#Check if revcom is enabled, not very original check....
 
183
  $self->{_strand}=1 if (($line=~/^Sequence name/) && ($line=~/Strand/));
184
184
        if ($line=~ m/\sSite\s/) {
185
 
            $instances=$self->_parseInstance;
 
185
            $instances= $self->_parseInstance;
186
186
        }
187
187
        #Here starts the next motif
188
188
        if ( ($line=~/width/) && ($line=~/sites/)) {
201
201
            ($ic)=split(/\s/,$line);
202
202
        }
203
203
        #Last info-prob matrix data
204
 
        if ($line=~/letter\-probability\s+matrix/) {
 
204
        if ($line=~/position-specific\s+scoring matrix/) {
 
205
                ($lA,$lC,$lG,$lT)=_parse_logs($self);
 
206
        }
 
207
        if ($line=~/^letter-probability\smatrix/) {
205
208
            my %matrix_dat=$self->_parseMatrix($motif_id);
206
209
            my $psm= new Bio::Matrix::PSM::Psm(%matrix_dat, 
207
210
                                               -instances=>$instances, 
208
211
                                               -e_val=>$e_val,
209
212
                                               -IC=>$ic, 
210
213
                                               -width=>$width, 
211
 
                                               -sites=>$sites);
 
214
                                               -sites=>$sites,
 
215
                                                   -lA=>$lA,
 
216
                                                   -lC=>$lC,
 
217
                                                   -lG=>$lG,
 
218
                                                   -lT=>$lT,
 
219
                                                   );
212
220
            return $psm;
213
221
        }
214
222
        if ($line=~"SUMMARY OF MOTIFS") {
215
223
            $self->{end}=1;
216
 
            return undef;
 
224
            return;
217
225
        }
218
226
        $endm=1 if ($line=~/^Time\s/); 
219
227
    }
220
228
        if ($endm) { #End of file found, end of current motif too, but not all predictions were made as requested (No summary)
221
229
            $self->{end}=1;
222
 
                warn "This MEME analysis was terminated prematurely, you may have less motifs than you requested\n";
223
 
            return undef;
 
230
            warn "This MEME analysis was terminated prematurely, you may have less motifs than you requested\n";
 
231
            return;
224
232
        }
225
233
    $self->throw("Wrong format\n"); # Multiple keywords not found, probably wrong format
226
234
}
246
254
    do {
247
255
        chomp $line;
248
256
        last if ($line eq '');
249
 
        $line=~s/\s\s/,/g;
250
 
        $line=~s/\s//g;
 
257
  $line=~s/^\s+//;
 
258
        $line=~s/\s+/,/g;
251
259
        ($pA[$i],$pC[$i],$pG[$i],$pT[$i])=split(/,/,$line);
252
260
        $i++;
253
261
        $line=$self->_readline;
254
262
    } until $line =~ /\-{10,}/;
 
263
    return (-pA=>\@pA,-pC=>\@pC,-pG=>\@pG,-pT=>\@pT,-id=>$id);
 
264
}
 
265
 
 
266
=head2 _parse_logs
 
267
 
 
268
 Title   : _parse_logs
 
269
 Usage   :
 
270
 Function: Parses the next site matrix log values in the meme file
 
271
 Throws  :
 
272
 Example :  Internal stuff
 
273
 Returns :  array of array refs
 
274
 Args    :  string
 
275
 
 
276
=cut
 
277
 
 
278
sub _parse_logs {
 
279
    my $self=shift;
 
280
    my (@lA,@lC,@lG,@lT);
 
281
    my $i=0;
 
282
    $self->_readline;   $self->_readline;
 
283
    my $line = $self->_readline;
 
284
    #Most important part- the probability matrix
 
285
    do {
 
286
        chomp $line;
 
287
        last if ($line eq '');
 
288
  $line=~s/^\s+//;
 
289
        $line=~s/\s+/,/g;
 
290
        ($lA[$i],$lC[$i],$lG[$i],$lT[$i])=split(/,/,$line);
 
291
        $i++;
 
292
        $line=$self->_readline;
 
293
    } until $line =~ /\-{10,}/;
255
294
    
256
 
    return (-pA=>\@pA,-pC=>\@pC,-pG=>\@pG,-pT=>\@pT,-id=>$id);
 
295
    return (\@lA,\@lC,\@lG,\@lT);
257
296
}
258
297
 
259
 
 
260
298
=head2 _parseInstance
261
299
 
262
300
 Title   : _parseInstance
277
315
    while (defined($line=$self->_readline) ) {
278
316
        last if ($line =~ /\-{5}/ );
279
317
        chomp($line);
280
 
        my $seq = $line;
281
 
        $seq  =~ s/[^AGCTXN-]//g;
282
 
        $line =~ s/\s[AGCTXN-]+\s[AGCTXN-]+\s[AGCTXN-]+//g;
283
 
        $line=~s/[\s\+]+/,/g;
284
 
        my ($id,$start,$score)=split(/,/,$line);
 
318
        my @comp=split(/\s+/,$line);
 
319
        my ($id,$start,$score,$strand,$s1,$s2,$s3);
 
320
        if ( $self->{_strand}) {
 
321
            ($id,$strand,$start,$score,$s1,$s2,$s3)=@comp;
 
322
        } else {
 
323
            ($id,$start,$score,$s1,$s2,$s3)=@comp;
 
324
            $strand=1;
 
325
        }
 
326
        my $seq= $s1.$s2.$s3;
 
327
        if ($seq =~ /[^ACGTacgtNnXx\-\.]/) {
 
328
            my $col=$#comp;
 
329
            $self->throw("I have not been able to parse the correct instance sequence: $seq, $col columns\n");
 
330
        }
285
331
        my $sid = $self->{id} . '@' . $id;
286
332
        $instance[$i] = new Bio::Matrix::PSM::InstanceSite
287
 
            (-mid   => $self->{id}, 
288
 
             -start => $start, 
289
 
             -score => $score,
290
 
             -seq   => $seq, 
 
333
            (-mid      => $self->{id}, 
 
334
             -start    => $start, 
 
335
             -score    => $score,
 
336
             -seq      => $seq, 
 
337
             -strand   => $strand,
291
338
             -accession_number => $id, 
292
339
             -primary_id => $sid, 
293
340
             -desc => 'Bioperl MEME parser object' );
297
344
    return \@instance;
298
345
}
299
346
 
 
347
                                
 
348
                        
 
349
 
300
350
1;