3 # BioPerl module for Bio::Tools::QRNA
5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
7 # Cared for by Jason Stajich <jason-at-bioperl-dot-org>
9 # Copyright Jason Stajich
11 # You may distribute this module under the same terms as perl itself
13 # POD documentation - main docs before the code
17 Bio::Tools::QRNA - A Parser for qrna output
22 my $parser = Bio::Tools::QRNA->new(-file => $qrnaoutput);
23 while( my $feature = $parser->next_feature ) {
29 Parses QRNA output (E.Rivas;
30 http://www.genetics.wustl.edu/eddy/software/#qrna).
32 This module is not complete, but currently it packs information from
33 each QRNA alignment into a single Bio::SeqFeature::Generic object.
35 Not all options for QRNA output have been tested or tried. It has
36 been tested on sliding window output (-w -x) and shuffled output (-b
39 See t/QRNA.t for example usage.
41 At some point we may have more complicated feature object which will
42 support this data rather than forcing most of the information into
43 tag/value pairs in a SeqFeature::Generic.
45 Running with -verbose =E<gt> 1 will store extra data in the feature. The
46 entire unparsed entry for a particular feature will be stored as a
47 string in the tag 'entry' it is accessible via:
49 my ($entry) = $f->each_tag_value('entry');
51 The winning model for any given alignment test will be the name stored
52 in the primary_tag field of feature. The bit score will stored in the
53 score field. The logoddpost is availble via the a tag/value pair.
54 This example code will show how to print out the score and log odds
57 # assuming you got a feature already
58 print "model score logoddspost\n";
59 foreach my $model ( qw(OTH COD RNA) ) {
60 my ($score) = $f->get_tag_values("$model\_score");
61 my ($logoddspost) = $f->get_tag_values("$model\_logoddspost");
62 print "$model $score $logoddspost\n";
65 The start and end of the alignment for both the query and hit sequence
66 are available through the L<Bio::SeqFeature::FeaturePair> interface,
67 specifically L<Bio::SeqFeature::FeaturePair::feature1> and
68 L<Bio::SeqFeature::FeaturePair::feature2>. Additionally if you have
69 run QRNA with an input file which has the location of the alignment
70 stored in the FASTA filename as in (ID/START-END) which is the default
71 output format from L<Bio::AlignIO::fasta> produced alignment output,
72 this module will re-number start/end for the two sequences so they are
73 in the actual coordinates of the sequence rather than the relative
74 coordinates of the alignment. You may find the bioperl utillity
75 script search2alnblocks useful in creating your input files for QRNA.
77 Some other words of warning, QRNA uses a 0 based numbering system for
78 sequence locations, Bioperl uses a 1 based system. You'll notice that
79 locations will be +1 they are reported in the raw QRNA output.
85 User feedback is an integral part of the evolution of this and other
86 Bioperl modules. Send your comments and suggestions preferably to
87 the Bioperl mailing list. Your participation is much appreciated.
89 bioperl-l@bioperl.org - General discussion
90 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
94 Please direct usage questions or support issues to the mailing list:
96 L<bioperl-l@bioperl.org>
98 rather than to the module maintainer directly. Many experienced and
99 reponsive experts will be able look at the problem and quickly
100 address it. Please include a thorough description of the problem
101 with code and data examples if at all possible.
103 =head2 Reporting Bugs
105 Report bugs to the Bioperl bug tracking system to help us keep track
106 of the bugs and their resolution. Bug reports can be submitted via
109 http://bugzilla.open-bio.org/
111 =head1 AUTHOR - Jason Stajich
113 Email jason-at-bioperl-dot-org
117 The rest of the documentation details each of the object methods.
118 Internal methods are usually preceded with a _
123 # Let the code begin...
126 package Bio
::Tools
::QRNA
;
127 use vars
qw(@Models);
130 use Bio::SeqFeature::Generic;
131 use Bio::SeqFeature::FeaturePair;
133 use base qw(Bio::Root::IO Bio::SeqAnalysisParserI);
134 @Models = qw(OTH COD RNA);
139 Usage : my $obj = Bio::Tools::QRNA->new();
140 Function: Builds a new Bio::Tools::QRNA object
141 Returns : an instance of Bio::Tools::QRNA
142 Args : -fh/-file filehandle/filename standard input for
150 Usage : my $feature = $parser->next_feature
151 Function: Get the next QRNA feature
160 my $f = shift @
{$self->{'_parsed_features'} || []};
161 if( ! defined $f && $self->_parse_pair ) {
162 $f = shift @
{$self->{'_parsed_features'} || []};
168 my ($self,@args) = @_;
169 my (@features,%data);
171 while( defined( $_ = $self->_readline) ) {
172 next if( /^\#\-\-/o );
173 if( /^\#\s+(qrna)\s+(\S+)\s+\(([^\)]+)\)/o ) {
174 $self->program_name($1);
175 $self->program_version($2);
176 $self->program_date($3);
177 } elsif( /^\#\s+(PAM model)\s+\=\s+(.+)\s+$/o ) {
178 $self->PAM_model($2);
179 } elsif( /^\#\s+(RNA model)\s+\=\s+(\S+)/o ) {
180 $self->RNA_model($2);
181 } elsif( /^\#\s+(seq file)\s+\=\s+(.+)\s+$/o ) {
183 } elsif( /^\#\s+(\d+)\s+\[([\-+])\s+strand\]/o ) {
185 if( $data{'alignment_len'} ) {
186 push @features, $self->_make_feature(\
%data);
188 $self->_pushback($_);
194 } elsif( />(\S+)\s+\((\d+)\)/ ) {
195 if( @
{$data{'seqs'} || []} == 2 ) {
196 $self->warn( "already seen seqs ".join(' ', ,map { $_->[0] }
197 @
{$data{'seqs'}}). "\n");
199 push @
{$data{'seqs'}}, [$1,$2];
201 } elsif( /^length alignment:\s+(\d+)\s+\(id\=(\d+(\.\d+)?)\)/o ) {
203 if( $data{'alignment_len'} ) {
204 push @features, $self->_make_feature(\
%data);
205 # reset all the data but the 'seqs' field
206 %data = ( 'seqs' => $data{'seqs'} );
209 if( /\(((sre_)?shuffled)\)/ ) {
210 $data{'shuffled'} = $1;
212 $data{'alignment_len'} = $1;
213 $data{'alignment_pid'} = $2;
214 } elsif ( /^pos([XY
]):\s
+(\d
+)\
-(\d
+)\s
+\
[(\d
+)\
-(\d
+)\
]\
((\d
+)\
)\s
+
215 \
-\
-\s
+\
((\S
+\s
+\S
+\s
+\S
+\s
+\S
+)\
)/ox
) {
216 $data{"seq\_$1"}->{'aln'} = [ $2,$3, $4,$5, $6];
217 @
{$data{"seq\_$1"}->{'base_comp'}} = split(/\s+/,$7);
218 } elsif( /^winner\s+\=\s+(\S{3})/ ) {
219 $data{'winning_model'} = $1;
220 } elsif( /^(\S{3})\s+ends\s+\=\s+(\-?\d+)\s+(\-?\d+)/ ) {
223 $data{'model_location'}->{$1} = [ $2,$3 ];
224 } elsif( /^\s+(logoddspost)?OTH\s+\=\s+/ox ) {
225 while( /(\S+)\s+\=\s+(\-?\d+(\.\d+))/g ) {
226 my ($model,$score)= ($1,$2);
227 if( $model =~ s/^logoddspost// ) {
228 $data{'model_scores'}->{'logoddspost'}->{$model} = $score;
230 $data{'model_scores'}->{'bits'}->{$model} = $score;
234 $data{'entry'} .= $_;
237 push @
{$self->{'_parsed_features'}}, @features;
238 return scalar @features;
246 Usage : $obj->PAM_model($newval)
249 Returns : value of PAM_model (a scalar)
250 Args : on set, new value (a scalar or undef, optional)
257 return $self->{'PAM_model'} = shift if @_;
258 return $self->{'PAM_model'};
264 Usage : $obj->RNA_model($newval)
267 Returns : value of RNA_model (a scalar)
268 Args : on set, new value (a scalar or undef, optional)
276 return $self->{'RNA_model'} = shift if @_;
277 return $self->{'RNA_model'};
283 Usage : $obj->seq_file($newval)
286 Returns : value of seq_file (a scalar)
287 Args : on set, new value (a scalar or undef, optional)
295 return $self->{'seq_file'} = shift if @_;
296 return $self->{'seq_file'};
303 Usage : $obj->program_name($newval)
306 Returns : value of program_name (a scalar)
307 Args : on set, new value (a scalar or undef, optional)
315 return $self->{'program_name'} = shift if @_;
316 return $self->{'program_name'} || 'qrna';
319 =head2 program_version
321 Title : program_version
322 Usage : $obj->program_version($newval)
325 Returns : value of program_version (a scalar)
326 Args : on set, new value (a scalar or undef, optional)
334 return $self->{'program_version'} = shift if @_;
335 return $self->{'program_version'};
341 Usage : $obj->program_date($newval)
344 Returns : value of program_date (a scalar)
345 Args : on set, new value (a scalar or undef, optional)
352 return $self->{'program_date'} = shift if @_;
353 return $self->{'program_date'};
357 my ($self,$data) = @_;
358 my ($qoffset,$hoffset) = (1,1);
359 # when you run qrna and have produced ID/START-END
360 # formatted input strings we can remap the location
363 # name is stored as the first entry in the seq array ref
364 my ($qid,$hid) = ( $data->{'seqs'}->[0]->[0],
365 $data->{'seqs'}->[1]->[0]);
366 if( $qid =~ /(\S+)\/(\d
+)\
-(\d
+)/ ) {
367 ($qid,$qoffset) = ($1,$2);
369 if( $hid =~ /(\S+)\/(\d
+)\
-(\d
+)/ ) {
370 ($hid,$hoffset) = ($1,$2);
373 my $f = Bio
::SeqFeature
::FeaturePair
->new();
375 my ($s,$e) = @
{$data->{'model_location'}->{$data->{'winning_model'}}};
376 my $qf = Bio
::SeqFeature
::Generic
->new
377 ( -primary_tag
=> $data->{'winning_model'},
378 -source_tag
=> $self->program_name,
379 -score
=> $data->{'model_scores'}->{'bits'}->{$data->{'winning_model'}},
380 -start
=> $s+$qoffset,
383 -strand
=> ($s < $e ) ?
1 : -1,
386 my $hf = Bio
::SeqFeature
::Generic
->new
387 ( -primary_tag
=> $qf->primary_tag,
388 -source_tag
=> $qf->source_tag,
389 -score
=> $qf->score,
391 -start
=> $s + $hoffset,
392 -end
=> $e + $hoffset,
393 -strand
=> $qf->strand,
397 $f->add_tag_value('alignment_len', $data->{'alignment_len'});
398 $f->add_tag_value('alignment_pid', $data->{'alignment_pid'});
399 # store the other model scores and data
400 foreach my $model ( @Models ) {
401 $f->add_tag_value("$model\_score", $data->{'model_scores'}->{'bits'}->{$model});
402 $f->add_tag_value("$model\_logoddspost", $data->{'model_scores'}->{'logoddspost'}->{$model});
403 if( ! $data->{'model_location'}->{$model} ) {
404 if( $self->verbose > 0 ) {
405 $self->debug( $data->{'entry'} );
407 $self->throw("no location parsed for $model in ",
408 (map { @
$_ } @
{$data->{'seqs'}}), " ", $f->start, " ", $f->end);
410 $f->add_tag_value("$model\_positions",
411 join("..",@
{$data->{'model_location'}->{$model} }));
414 # probably a better way to store this - as
415 # a seq object perhaps
416 $f->add_tag_value('seq1', @
{$data->{'seqs'}->[0]});
417 $f->add_tag_value('seq2', @
{$data->{'seqs'}->[1]});
418 $f->add_tag_value('entry', $data->{'entry'}) if $self->verbose > 0;
419 if( $data->{'shuffled'} ) {
420 $f->add_tag_value('shuffled', $data->{'shuffled'});