QRNA documentation.

Parses QRNA output (E.Rivas;
http://www.genetics.wustl.edu/eddy/software/#qrna).
This module is not complete, but currently it packs information from
each QRNA alignment into a single Bio::SeqFeature::Generic object.
Not all options for QRNA output have been tested or tried. It has
been tested on sliding window output (-w -x) and shuffled output (-b
or -B).
See t/QRNA.t for example usage.
At some point we may have more complicated feature object which will
support this data rather than forcing most of the information into
tag/value pairs in a SeqFeature::Generic.
Running with -verbose => 1 will store extra data in the feature. The
entire unparsed entry for a particular feature will be stored as a
string in the tag 'entry' it is accessible via:

  my ($entry) = $f->each_tag_value('entry');

The winning model for any given alignment test will be the name stored
in the primary_tag field of feature. The bit score will stored in the
score field. The logoddpost is availble via the a tag/value pair.
This example code will show how to print out the score and log odds
post for each model.

  # assuming you got a feature already
  print "model score logoddspost\n";
  foreach my $model ( qw(OTH COD RNA) ) {
    my ($score)       = $f->get_tag_values("$model\_score");
    my ($logoddspost) = $f->get_tag_values("$model\_logoddspost");
    print "$model $score $logoddspost\n";
  }

The start and end of the alignment for both the query and hit sequence
are available through the Bio::SeqFeature::FeaturePair interface,
specifically Bio::SeqFeature::FeaturePair::feature1 and
Bio::SeqFeature::FeaturePair::feature2. Additionally if you have
run QRNA with an input file which has the location of the alignment
stored in the FASTA filename as in (ID/START-END) which is the default
output format from Bio::AlignIO::fasta produced alignment output,
this module will re-number start/end for the two sequences so they are
in the actual coordinates of the sequence rather than the relative
coordinates of the alignment. You may find the bioperl utillity
script search2alnblocks useful in creating your input files for QRNA.
Some other words of warning, QRNA uses a 0 based numbering system for
sequence locations, Bioperl uses a 1 based system. You'll notice that
locations will be +1 they are reported in the raw QRNA output.

sub _parse_pair {

   my ($self,@args) = @_;
   my (@features,%data);
   my $seenstart = 0;
   while( defined( $_ = $self->_readline) ) {
       next if( /^\#\-\-/o );
       if( /^\#\s+(qrna)\s+(\S+)\s+\(([^\)]+)\)/o ) {
	   $self->program_name($1);
	   $self->program_version($2);
	   $self->program_date($3);
       } elsif( /^\#\s+(PAM model)\s+\=\s+(.+)\s+$/o ) {
	   $self->PAM_model($2);
       } elsif( /^\#\s+(RNA model)\s+\=\s+(\S+)/o ) {
	   $self->RNA_model($2);
       } elsif( /^\#\s+(seq file)\s+\=\s+(.+)\s+$/o ) {
	   $self->seq_file($2);	   
       } elsif( /^\#\s+(\d+)\s+\[([\-+])\s+strand\]/o ) {
	   if( $seenstart ) { 
	       if( $data{'alignment_len'} ) {
		   push @features, $self->_make_feature(\%data);
	       }
	       $self->_pushback($_);
	       last;
	   }
	   $seenstart = 1;
       } elsif( /^\#/ ) {
	   next;
       } elsif( />(\S+)\s+\((\d+)\)/ ) {
	   if( @{$data{'seqs'} || []} == 2 ) { 
	       $self->warn( "already seen seqs ".join(' ', ,map { $_->[0] } 
						      @{$data{'seqs'}}). "\n");
	   } else { 
	       push @{$data{'seqs'}}, [$1,$2];
	   }
       } elsif( /^length alignment:\s+(\d+)\s+\(id\=(\d+(\.\d+)?)\)/o ) {
	   
	   if( $data{'alignment_len'} ) {
	       push @features, $self->_make_feature(\%data);	
	       # reset all the data but the 'seqs' field
	       %data  = ( 'seqs' => $data{'seqs'} );
	   }
	   
	   if( /\(((sre_)?shuffled)\)/ ) { 
	       $data{'shuffled'} = $1;
	   }
	   $data{'alignment_len'} = $1;
	   $data{'alignment_pid'} = $2;
       } elsif ( /^pos([XY]):\s+(\d+)\-(\d+)\s+\[(\d+)\-(\d+)\]\((\d+)\)\s+
		 \-\-\s+\((\S+\s+\S+\s+\S+\s+\S+)\)/ox ) {
	   $data{"seq\_$1"}->{'aln'} = [ $2,$3, $4,$5, $6];
	   @{$data{"seq\_$1"}->{'base_comp'}} = split(/\s+/,$7);
       } elsif( /^winner\s+\=\s+(\S{3})/ ) {
	   $data{'winning_model'} = $1;
       } elsif( /^(\S{3})\s+ends\s+\=\s+(\-?\d+)\s+(\-?\d+)/ ) {
	   # QRNA is 0-based
	   # Bioperl is 1 based
	   $data{'model_location'}->{$1} = [ $2,$3 ];
       }  elsif( /^\s+(logoddspost)?OTH\s+\=\s+/ox ) {
	   while( /(\S+)\s+\=\s+(\-?\d+(\.\d+))/g ) {
	       my ($model,$score)= ($1,$2);
	       if( $model =~ s/^logoddspost// ) {
		   $data{'model_scores'}->{'logoddspost'}->{$model} = $score;
	       } else {
		   $data{'model_scores'}->{'bits'}->{$model} = $score;
	       }
	   }
       }
       $data{'entry'} .= $_;

}

sub _make_feature {

     my ($self,$data) = @_; 
    my ($qoffset,$hoffset) = (1,1);
    # when you run qrna and have produced ID/START-END
    # formatted input strings we can remap the location
    # to the original

    # name is stored as the first entry in the seq array ref
    my ($qid,$hid) = ( $data->{'seqs'}->[0]->[0],
		       $data->{'seqs'}->[1]->[0]);
    if( $qid =~ /(\S+)\/(\d+)\-(\d+)/ ) {
	($qid,$qoffset) = ($1,$2);
    }
    if( $hid =~ /(\S+)\/(\d+)\-(\d+)/ ) {
	($hid,$hoffset) = ($1,$2);
    }

    my $f = new Bio::SeqFeature::FeaturePair;

    my ($s,$e) = @{$data->{'model_location'}->{$data->{'winning_model'}}};
    my $qf = new Bio::SeqFeature::Generic
	( -primary_tag => $data->{'winning_model'},
	  -source_tag  => $self->program_name,
	  -score       => $data->{'model_scores'}->{'bits'}->{$data->{'winning_model'}},
	  -start       => $s+$qoffset,
	  -end         => $e+$qoffset,
	  -seq_id      => $qid,
	  -strand      => ($s < $e ) ? 1 : -1,
	  );

    my $hf = new Bio::SeqFeature::Generic
	( -primary_tag => $qf->primary_tag,
	  -source_tag  => $qf->source_tag,
	  -score       => $qf->score,
	  -seq_id      => $hid,
	  -start       => $s + $hoffset,
	  -end         => $e + $hoffset,
	  -strand      => $qf->strand,
	  );
    $f->feature1($qf);
    $f->feature2($hf);
    $f->add_tag_value('alignment_len', $data->{'alignment_len'});
    $f->add_tag_value('alignment_pid', $data->{'alignment_pid'});
    # store the other model scores and data
    foreach my $model ( @Models ) {
	$f->add_tag_value("$model\_score", $data->{'model_scores'}->{'bits'}->{$model});
	$f->add_tag_value("$model\_logoddspost", $data->{'model_scores'}->{'logoddspost'}->{$model});
	if( ! $data->{'model_location'}->{$model} ) {
	    if( $self->verbose > 0 ) {
		$self->debug( $data->{'entry'} );
	    }
	    die "no location parsed for $model in ",
	    (map { @$_ } @{$data->{'seqs'}}), " ", $f->start, " ", $f->end, "\n";
	} else { 
	    $f->add_tag_value("$model\_positions", 
			      join("..",@{$data->{'model_location'}->{$model} }));
	}
    }
    # probably a better way to store this - as 
    # a seq object perhaps
    $f->add_tag_value('seq1', @{$data->{'seqs'}->[0]});
    $f->add_tag_value('seq2', @{$data->{'seqs'}->[1]});
    $f->add_tag_value('entry', $data->{'entry'}) if $self->verbose > 0;
    if( $data->{'shuffled'} ) {
	$f->add_tag_value('shuffled', $data->{'shuffled'});
    }
    return $f;

}

next_feature	Description	Code
_parse_pair	No description	Code
PAM_model	Description	Code
RNA_model	Description	Code
seq_file	Description	Code
program_name	Description	Code
program_version	Description	Code
program_date	Description	Code
_make_feature	No description	Code