← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to Bio/AlignIO/meme.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/AlignIO/meme.pm
35
35
and end coordinates are derived from the "Start" field. See
36
36
L<Bio::SimpleAlign> and L<Bio::LocatableSeq> for more information.
37
37
 
38
 
This module can only parse MEME version 3.0 and greater.  Previous
 
38
This module can only parse MEME version 3 and 4.  Previous
39
39
versions have output formats that are more difficult to parse
40
 
correctly.  If the meme output file is not version 3.0 or greater,
 
40
correctly.  If the meme output file is not version 3.0 or greater
41
41
we signal an error.
42
42
 
43
43
=head1 FEEDBACK
84
84
 
85
85
# Constants
86
86
my $MEME_VERS_ERR =
87
 
"MEME output file must be generated by version 3.0 or higher";
 
87
  "MEME output file must be generated by version 3.0 or higher";
88
88
my $MEME_NO_HEADER_ERR =
89
 
"MEME output file contains no header line (ex: MEME version 3.0)";
90
 
my $HTML_VERS_ERR =
91
 
"MEME output file must be generated with the -text option";
 
89
  "MEME output file contains no header line (ex: MEME version 3.0)";
 
90
my $HTML_VERS_ERR = "MEME output file must be generated with the -text option";
92
91
 
93
92
=head2 next_aln
94
93
 
102
101
=cut
103
102
 
104
103
sub next_aln {
105
 
        my ($self) = @_;
106
 
        my $aln = Bio::SimpleAlign->new(-source => 'meme');
107
 
        my $line;
108
 
        my $good_align_sec = 0;
109
 
        my $in_align_sec = 0;
110
 
        my $evalue;
111
 
        while (!$good_align_sec && defined($line = $self->_readline())) {
112
 
                if (!$in_align_sec) {
113
 
                        # Check for the meme header
114
 
                        if ($line =~ /^\s*MEME\s+version\s+(\S+)/ ){
115
 
                                $self->{'meme_vers'} = $1;
 
104
    my ($self) = @_;
 
105
    my $aln = Bio::SimpleAlign->new( -source => 'meme' );
 
106
    my $line;
 
107
    my $good_align_sec = 0;
 
108
    my $in_align_sec   = 0;
 
109
    my $evalue;
 
110
    while ( !$good_align_sec && defined( $line = $self->_readline() ) ) {
 
111
        if ( !$in_align_sec ) {
 
112
 
 
113
            # Check for the meme header
 
114
            if ( $line =~ /^\s*MEME\s+version\s+(\S+)/ ) {
 
115
                $self->{'meme_vers'} = $1;
116
116
                my ($vers) = $self->{'meme_vers'} =~ /^(\d)/;
117
 
                                $self->throw($MEME_VERS_ERR) unless ($vers >= 3);
118
 
                                $self->{'seen_header'} = 1;
119
 
              }
120
 
                        # Check if they've output the HTML version
121
 
                        if ($line =~ /\<TITLE\>/i){
122
 
                                $self->throw($HTML_VERS_ERR);
123
 
              }
124
 
                        # Grab the evalue
125
 
                        if ($line =~ /MOTIF\s+\d+\s+width.+E-value = (\S+)/) {
126
 
                                $self->throw($MEME_NO_HEADER_ERR) unless ($self->{'seen_header'});
127
 
                                $evalue = $1;
128
 
                        }
129
 
                        # Check if we're going into an alignment section
130
 
                        if ($line =~ /sites sorted by position/) {
131
 
                                $self->throw($MEME_NO_HEADER_ERR) unless ($self->{'seen_header'});
132
 
                                $in_align_sec = 1;
133
 
                        }
134
 
                } elsif ($line =~ /^(\S+)\s+([+-]?)\s+(\d+)\s+
135
 
                       (\S+)\s+([.ACTGNX\-]*)\s+([ACTGNX\-]+)\s+
136
 
                       ([.ACTGNX\-]*)/xi ) {
137
 
                        # Got a sequence line
138
 
                        my $seq_name = $1;
139
 
                        my $strand = ($2 eq '-') ? -1 : 1;
140
 
                        my $start_pos = $3;
141
 
                        my $central = uc($6);
142
 
 
143
 
                        # my $p_val = $4;
144
 
                        # my $left_flank = uc($5);
145
 
                        # my $right_flank = uc($7);
146
 
 
147
 
                        # Info about the flanking sequence
148
 
                        # my $start_len = ($strand > 0) ? length($left_flank) :
149
 
                        # length($right_flank);
150
 
                        # my $end_len = ($strand > 0) ? length($right_flank) :
151
 
                        # length($left_flank);
152
 
 
153
 
                        # Make the sequence.  Meme gives the start coordinate at the left
154
 
                        # hand side of the motif relative to the INPUT sequence.
155
 
                        my $end_pos = $start_pos + length($central) - 1;
156
 
                        my $seq = Bio::LocatableSeq->new
157
 
                            ('-seq'            => $central,
158
 
                             '-display_id'     => $seq_name,
159
 
                             '-start'          => $start_pos,
160
 
                             '-end'            => $end_pos,
161
 
                             '-strand'         => $strand,
162
 
                             '-alphabet'       => $self->alphabet,
163
 
                             );
164
 
                        # Add the sequence motif to the alignment
165
 
                        $aln->add_seq($seq);
166
 
                } elsif (($line =~ /^\-/) || ($line =~ /Sequence name/)){
167
 
                        # These are acceptable things to be in the site section
168
 
                } elsif ($line =~ /^\s*$/){
169
 
                        # This ends the site section
170
 
                        $in_align_sec = 0;
171
 
                        $good_align_sec = 1;
172
 
                } else{
173
 
                        $self->warn("Unrecognized format:\n$line");
174
 
                        return 0;
175
 
                }
176
 
        }
177
 
        # Signal an error if we didn't find a header section
178
 
        $self->throw($MEME_NO_HEADER_ERR) unless ($self->{'seen_header'});
179
 
        
180
 
        if ($good_align_sec) {
181
 
                $aln->score($evalue);
182
 
                return $aln;
183
 
        }
184
 
        
185
 
        return;
 
117
 
 
118
                $self->throw($MEME_VERS_ERR) unless ( $vers >= 3 );
 
119
                $self->{'seen_header'} = 1;
 
120
            }
 
121
 
 
122
            # Check if they've output the HTML version
 
123
            if ( $line =~ /\<TITLE\>/i ) {
 
124
                $self->throw($HTML_VERS_ERR);
 
125
            }
 
126
 
 
127
            # Grab the evalue
 
128
            if ( $line =~ /MOTIF\s+\d+\s+width.+E-value = (\S+)/ ) {
 
129
                $self->throw($MEME_NO_HEADER_ERR)
 
130
                  unless ( $self->{'seen_header'} );
 
131
                $evalue = $1;
 
132
            }
 
133
 
 
134
            # Check if we're going into an alignment section
 
135
            if ( $line =~ /sites sorted by position/ ) {
 
136
                $self->throw($MEME_NO_HEADER_ERR)
 
137
                  unless ( $self->{'seen_header'} );
 
138
                $in_align_sec = 1;
 
139
            }
 
140
        }
 
141
        # The first regexp is for version 3, the second is for version 4
 
142
        elsif ( $line =~ /^(\S+)\s+([+-]?)\s+(\d+)\s+
 
143
                           \S+\s+[.A-Z\-]*\s+([A-Z\-]+)\s+
 
144
                           ([.A-Z\-]*)/xi
 
145
                ||
 
146
                $line =~ /^(\S+)\s+([+-]?)\s+(\d+)\s+
 
147
                           \S+\s+\.\s+([A-Z\-]+)/xi
 
148
          )
 
149
        {
 
150
            # Got a sequence line
 
151
            my $seq_name  = $1;
 
152
            my $strand    = ( $2 eq '-' ) ? -1 : 1;
 
153
            my $start_pos = $3;
 
154
            my $central   = uc($4);
 
155
 
 
156
            # my $p_val = $4;
 
157
            # my $left_flank = uc($5);
 
158
            # my $right_flank = uc($7);
 
159
 
 
160
            # Info about the flanking sequence
 
161
            # my $start_len = ($strand > 0) ? length($left_flank) :
 
162
            # length($right_flank);
 
163
            # my $end_len = ($strand > 0) ? length($right_flank) :
 
164
            # length($left_flank);
 
165
 
 
166
            # Make the sequence.  Meme gives the start coordinate at the left
 
167
            # hand side of the motif relative to the INPUT sequence.
 
168
            my $end_pos = $start_pos + length($central) - 1;
 
169
            my $seq     = Bio::LocatableSeq->new(
 
170
                -seq        => $central,
 
171
                -display_id => $seq_name,
 
172
                -start      => $start_pos,
 
173
                -end        => $end_pos,
 
174
                -strand     => $strand,
 
175
                -alphabet   => $self->alphabet,
 
176
            );
 
177
 
 
178
            # Add the sequence motif to the alignment
 
179
            $aln->add_seq($seq);
 
180
        }
 
181
        elsif ( ( $line =~ /^\-/ ) || ( $line =~ /Sequence name/ ) ) {
 
182
 
 
183
            # These are acceptable things to be in the site section
 
184
        }
 
185
        elsif ( $line =~ /^\s*$/ ) {
 
186
 
 
187
            # This ends the site section
 
188
            $in_align_sec   = 0;
 
189
            $good_align_sec = 1;
 
190
        }
 
191
        else {
 
192
            $self->warn("Unrecognized format:\n$line");
 
193
            return 0;
 
194
        }
 
195
    }
 
196
 
 
197
    # Signal an error if we didn't find a header section
 
198
    $self->throw($MEME_NO_HEADER_ERR) unless ( $self->{'seen_header'} );
 
199
 
 
200
    if ($good_align_sec) {
 
201
        $aln->score($evalue);
 
202
        return $aln;
 
203
    }
 
204
 
 
205
    return;
186
206
}
187
207
 
188
208
=head2 write_aln
196
216
=cut
197
217
 
198
218
sub write_aln {
199
 
   my ($self,@aln) = @_;
200
 
   $self->throw_not_implemented();
 
219
    my ( $self, @aln ) = @_;
 
220
    $self->throw_not_implemented();
201
221
}
202
222
 
203
223
# ----------------------------------------
205
225
# ----------------------------------------
206
226
 
207
227
sub _initialize {
208
 
  my($self,@args) = @_;
209
 
 
210
 
  # Call into our base version
211
 
  $self->SUPER::_initialize(@args);
212
 
 
213
 
  # Then initialize our data variables
214
 
  $self->{'seen_header'} = 0;
 
228
    my ( $self, @args ) = @_;
 
229
 
 
230
    # Call into our base version
 
231
    $self->SUPER::_initialize(@args);
 
232
 
 
233
    # Then initialize our data variables
 
234
    $self->{'seen_header'} = 0;
215
235
}
216
236
 
217
237
1;