File Coverage

Bio/Matrix/PSM/IO/meme.pm
Criterion Covered Total %
statement 130 148 87.8
branch 27 40 67.5
condition 8 9 88.8
subroutine 12 13 92.3
pod 3 3 100.0
total 180 213 84.5


line stmt bran cond sub pod time code
1             #---------------------------------------------------------
2              
3             =head1 NAME
4              
5             Bio::Matrix::PSM::IO::meme - PSM meme parser implementation
6              
7             =head1 SYNOPSIS
8              
9             See Bio::Matrix::PSM::IO for detailed documentation on how to use PSM parsers
10              
11             =head1 DESCRIPTION
12              
13             Parser for meme.
14              
15             =head1 FEEDBACK
16              
17             =head2 Mailing Lists
18              
19             User feedback is an integral part of the evolution of this
20             and other Bioperl modules. Send your comments and suggestions preferably
21             to one of the Bioperl mailing lists.
22             Your participation is much appreciated.
23              
24             bioperl-l@bioperl.org - General discussion
25             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
26              
27             =head2 Support
28              
29             Please direct usage questions or support issues to the mailing list:
30              
31             I
32              
33             rather than to the module maintainer directly. Many experienced and
34             reponsive experts will be able look at the problem and quickly
35             address it. Please include a thorough description of the problem
36             with code and data examples if at all possible.
37              
38             =head2 Reporting Bugs
39              
40             Report bugs to the Bioperl bug tracking system to help us keep track
41             the bugs and their resolution. Bug reports can be submitted via the
42             web:
43              
44             https://github.com/bioperl/bioperl-live/issues
45              
46             =head1 AUTHOR - Stefan Kirov
47              
48             Email skirov@utk.edu
49              
50             =head1 APPENDIX
51              
52             =cut
53              
54              
55             # Let the code begin...
56             package Bio::Matrix::PSM::IO::meme;
57 2     2   606 use Bio::Matrix::PSM::InstanceSite;
  2         4  
  2         59  
58 2     2   824 use Bio::Matrix::PSM::SiteMatrix;
  2         4  
  2         56  
59 2     2   768 use Bio::Matrix::PSM::Psm;
  2         4  
  2         56  
60 2     2   8 use vars qw(@HEADER);
  2         3  
  2         86  
61 2     2   6 use strict;
  2         2  
  2         38  
62              
63 2     2   7 use base qw(Bio::Matrix::PSM::PsmHeader Bio::Matrix::PSM::IO);
  2         3  
  2         757  
64              
65             @Bio::Matrix::PSM::IO::meme::HEADER = qw(e_val sites IC width);
66              
67             =head2 new
68              
69             Title : new
70             Usage : my $psmIO = Bio::Matrix::PSM::IO->new(-format=>'meme',
71             -file=>$file);
72             Function: Associates a file with the appropriate parser
73             Throws : Throws if the file passed is in HTML format or
74             if the MEME header cannot be found.
75             Example :
76             Args : hash
77             Returns : "Bio::Matrix::PSM::$format"->new(@args);
78              
79             =cut
80              
81             sub new {
82 2     2 1 6 my($class, @args)=@_;
83 2         10 my $self = $class->SUPER::new(@args);
84 2         12 my ($file)=$self->_rearrange(['FILE'], @args);
85 2         8 my ($query,$tr1)=split(/\./,$file,2);
86 2         4 $self->{file} = $file;
87 2         2 $self->{query}= $query;
88 2         3 $self->{end} = 0;
89 2         4 $self->{_strand}=0; #This we'll need to see if revcom option is used
90 2 50       14 $self->_initialize_io(@args) || warn "Did you intend to use STDIN?"; #Read only for now
91             #Skip header
92 2         3 my $line;
93 2         12 while (my $line=$self->_readline) {
94 64 50       89 $self->throw('Cannot parse HTML, please use text output\n') if ($line=~//); #Should start parsing HTML output, not a bug deal
95 64         42 chomp($line);
96 64 100       72 if ($line=~"^ALPHABET") {
97 2         5 $self=_parse_coordinates($self);
98 2         5 last;
99             }
100 62 100 100     203 push @{$self->{unstructured}},$line unless (($line=~/\*{10,}/) || ($line eq ''));
  32         67  
101             }
102 2         10 $self->_initialize;
103 2         11 return $self;
104             }
105              
106             =head2 _parse_coordinates
107              
108             Title : _parse_coordinates
109             Usage :
110             Function:
111             Throws :
112             Example : Internal stuff
113             Returns :
114             Args :
115              
116             =cut
117              
118             sub _parse_coordinates {
119 2     2   2 my $self=shift;
120 2         5 $self->_readline;
121 2         4 $self->_readline;
122 2         4 my $line=$self->_readline;
123 2         10 while ($line !~ /^\*{10,}/ ) {
124 4         4 chomp $line;
125 4         23 $line =~ s/\s+/,/g;
126 4         15 my ($id1,$w1,$l1,$id2,$w2,$l2)=split(/,/,$line);
127 4         5 push @{$self->{hid}},$id1;
  4         8  
128 4         8 $self->{weight}->{$id1}=$w1;
129 4         7 $self->{length}->{$id1}=$l1;
130 4 50       9 if ($id2) {
131 4         3 push @{$self->{hid}},$id2;
  4         4  
132 4         8 $self->{weight}->{$id2}=$w2;
133 4         4 $self->{length}->{$id2}=$l2;
134             }
135 4         6 $line=$self->_readline;
136             }
137 2         3 return $self;
138             }
139              
140             =head2 header
141              
142             Title : header
143             Usage : my %header=$psmIO->header;
144             Function: Returns the header for the MEME file
145             Throws :
146             Example : Fetching all the sequences included in the MEME analysis,
147             being parsed
148             my %header=$psmIO->header;
149             foreach my $seqid (@{$header{instances}}) {
150             my $seq=$db->get_Seq_by_acc($id);
151             #Do something with the sequence
152             }
153             where $db might be Bio::DB:GenBank object, see
154             Returns : Hash with three keys: instances, weights and lengths, which
155             should be self-explenatory. Each value is an array
156             reference. Each array element corresponds to the same
157             element in the other two arrays. So $header{instances}->[$i]
158             will refer to the same sequence in the motif file as
159             $header{weights}->[$i] and $header{lengths}->[$i]
160             Args : none
161             Notes : OBSOLETE!
162              
163             =cut
164              
165             sub header {
166 0     0 1 0 my $self=shift;
167 0         0 my @instances=@{$self->{_inst_name}};
  0         0  
168 0         0 my @weights=@{$self->{_inst_weight}};
  0         0  
169 0         0 my @lengths=@{$self->{_inst_coord}};
  0         0  
170 0         0 return (instances=>\@instances,weights=>\@weights,lengths=>\@lengths);
171             }
172              
173             =head2 next_psm
174              
175             Title : next_psm
176             Usage : my $psm=$psmIO->next_psm();
177             Function: Reads the next PSM from the input file, associated with this object
178             Throws : Throws if the format is inconsistent with the rules for MEME 3.0.4:
179             no SUMMARY Section present or some keywords are missing/altered.
180             Example :
181             Returns : Bio::Matrix::PSM::Psm object
182             Args : none
183              
184             =cut
185              
186             sub next_psm {
187             #Parses the next prediction and returns a psm objects
188 2     2 1 3 my $self=shift;
189 2 50       6 return if ($self->{end});
190 2         2 my ($endm,$line,$instances,$tr,$width,$motif_id,$sites,$e_val,$id,$ic,$lA,$lC,$lG,$lT);
191 2         6 while (defined( $line = $self->_readline) ) {
192             #Check if revcom is enabled, not very original check....
193 182 50 66     264 $self->{_strand}=1 if (($line=~/^Sequence name/) && ($line=~/Strand/));
194 182 100       221 if ($line=~ m/\sSite\s/) {
195 2         6 $instances= $self->_parseInstance;
196             }
197             #Here starts the next motif
198 182 100 100     265 if ( ($line=~/width/) && ($line=~/sites/)) {
199 2         4 chomp($line);
200 2         15 $line=~s/[\t\s=]+/,/g;
201 2         4 $line=~s/\t/,/g;
202             #Parsing the general information for this prediction
203 2         9 ($tr,$motif_id,$tr,$width,$tr,$sites,
204             $tr,$tr,$tr,$e_val)=split(/,/,$line);
205 2         7 $self->{id}=$self->{query} . $motif_id;
206             }
207 182 100       252 if ($line =~ /content/i) {
208 2         4 $line=$self->_readline;
209 2         4 chomp($line);
210 2         7 $line=~s/[\)\(]//g;
211 2         7 ($ic)=split(/\s/,$line);
212             }
213             #Last info-prob matrix data
214 182 100       232 if ($line=~/position-specific\s+scoring matrix/) {
215 2         5 ($lA,$lC,$lG,$lT)=_parse_logs($self);
216             }
217 182 100       192 if ($line=~/^letter-probability\smatrix/) {
218 2         6 my %matrix_dat=$self->_parseMatrix($motif_id);
219 2         27 my $psm= Bio::Matrix::PSM::Psm->new(%matrix_dat,
220             -instances=>$instances,
221             -e_val=>$e_val,
222             -IC=>$ic,
223             -width=>$width,
224             -sites=>$sites,
225             -lA=>$lA,
226             -lC=>$lC,
227             -lG=>$lG,
228             -lT=>$lT,
229             );
230 2         10 return $psm;
231             }
232 180 50       187 if ($line=~"SUMMARY OF MOTIFS") {
233 0         0 $self->{end}=1;
234 0         0 return;
235             }
236 180 50       307 $endm=1 if ($line=~/^Time\s/);
237             }
238 0 0       0 if ($endm) { #End of file found, end of current motif too, but not all predictions were made as requested (No summary)
239 0         0 $self->{end}=1;
240 0         0 warn "This MEME analysis was terminated prematurely, you may have less motifs than you requested\n";
241 0         0 return;
242             }
243 0         0 $self->throw("Wrong format\n"); # Multiple keywords not found, probably wrong format
244             }
245              
246             =head2 _parseMatrix
247              
248             Title : _parseMatrix
249             Usage :
250             Function: Parses the next site matrix information in the meme file
251             Throws :
252             Example : Internal stuff
253             Returns : hash as for constructing a SiteMatrix object (see SiteMatrixI)
254             Args : string
255              
256             =cut
257              
258             sub _parseMatrix {
259 2     2   3 my ($self,$id)=@_;
260 2         3 my (@pA,@pC,@pG,@pT);
261 2         6 my $i=0;
262 2         4 my $line = $self->_readline;
263             #Most important part- the probability matrix
264 2         3 do {
265 50         36 chomp $line;
266 50 50       59 last if ($line eq '');
267 50         85 $line=~s/^\s+//;
268 50         127 $line=~s/\s+/,/g;
269 50         145 ($pA[$i],$pC[$i],$pG[$i],$pT[$i])=split(/,/,$line);
270 50         40 $i++;
271 50         66 $line=$self->_readline;
272             } until $line =~ /\-{10,}/;
273 2         16 return (-pA=>\@pA,-pC=>\@pC,-pG=>\@pG,-pT=>\@pT,-id=>$id);
274             }
275              
276             =head2 _parse_logs
277              
278             Title : _parse_logs
279             Usage :
280             Function: Parses the next site matrix log values in the meme file
281             Throws :
282             Example : Internal stuff
283             Returns : array of array refs
284             Args : string
285              
286             =cut
287              
288             sub _parse_logs {
289 2     2   2 my $self=shift;
290 2         3 my (@lA,@lC,@lG,@lT);
291 2         2 my $i=0;
292 2         4 $self->_readline; $self->_readline;
  2         4  
293 2         5 my $line = $self->_readline;
294             #Most important part- the probability matrix
295 2         3 do {
296 50         33 chomp $line;
297 50 50       57 last if ($line eq '');
298 50         88 $line=~s/^\s+//;
299 50         111 $line=~s/\s+/,/g;
300 50         150 ($lA[$i],$lC[$i],$lG[$i],$lT[$i])=split(/,/,$line);
301 50         46 $i++;
302 50         60 $line=$self->_readline;
303             } until $line =~ /\-{10,}/;
304            
305 2         5 return (\@lA,\@lC,\@lG,\@lT);
306             }
307              
308             =head2 _parseInstance
309              
310             Title : _parseInstance
311             Usage :
312             Function: Parses the next sites instances from the meme file
313             Throws :
314             Example : Internal stuff
315             Returns : Bio::Matrix::PSM::InstanceSite object
316             Args : none
317              
318             =cut
319              
320             sub _parseInstance {
321 2     2   9 my $self = shift;
322 2         2 my $i=0;
323 2         22 $self->_readline;
324 2         2 my ($line,@instance);
325 2         12 while (defined($line=$self->_readline) ) {
326 10 100       23 last if ($line =~ /\-{5}/ );
327 8         7 chomp($line);
328 8         41 my @comp=split(/\s+/,$line);
329 8         5 my ($id,$start,$score,$strand,$s1,$s2,$s3);
330 8 50       13 if ( $self->{_strand}) {
331 0         0 ($id,$strand,$start,$score,$s1,$s2,$s3)=@comp;
332             } else {
333 8         12 ($id,$start,$score,$s1,$s2,$s3)=@comp;
334 8         11 $strand=1;
335             }
336 8         14 my $seq= $s1.$s2.$s3;
337 8 50       21 if ($seq =~ /[^ACGTacgtNnXx\-\.]/) {
338 0         0 my $col=$#comp;
339 0         0 $self->throw("I have not been able to parse the correct instance sequence: $seq, $col columns\n");
340             }
341 8         11 my $sid = $self->{id} . '@' . $id;
342             $instance[$i] = Bio::Matrix::PSM::InstanceSite->new
343             (-mid => $self->{id},
344 8         47 -start => $start,
345             -score => $score,
346             -seq => $seq,
347             -strand => $strand,
348             -accession_number => $id,
349             -primary_id => $sid,
350             -desc => 'Bioperl MEME parser object' );
351 8         29 $i++;
352             }
353 2         4 $self->{instances} = \@instance;
354 2         4 return \@instance;
355             }
356              
357            
358            
359              
360             1;