lib/Bio/Tools/Phylo/Gerp.pm

   1 # $Id: Gumby.pm,v 1.2 2007/06/14 18:01:52 nathan Exp $
   2 #
   3 # BioPerl module for Bio::Tools::Phylo::Gerp
   4 #
   5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   6 #
   7 # Cared for by Sendu Bala <bix@sendu.me.uk>
   8 #
   9 # Copyright Sendu Bala
  10 #
  11 # You may distribute this module under the same terms as perl itself
  12
  13 # POD documentation - main docs before the code
  14
  15 =head1 NAME
  16
  17 Bio::Tools::Phylo::Gerp - Parses output from GERP
  18
  19 =head1 SYNOPSIS
  20
  21   use strict;
  22
  23   use Bio::Tools::Phylo::Gerp;
  24
  25   my $parser = Bio::Tools::Phylo::Gerp->new(-file => "alignment.rates.elems");
  26
  27   while (my $feat = $parser->next_result) {
  28     my $start = $feat->start;
  29     my $end = $feat->end;
  30     my $rs_score = $feat->score;
  31     my $p_value = ($feat->annotation->get_Annotations('p-value'))[0]->value;
  32   }
  33
  34 =head1 DESCRIPTION
  35
  36 This module is used to parse the output from 'GERP' (v2) by Eugene Davydov
  37 (originally Gregory M. Cooper et al.). You can get details here:
  38 http://mendel.stanford.edu/sidowlab/
  39
  40 It works on the .elems files produced by gerpelem.
  41
  42 Each result is a Bio::SeqFeature::Annotated representing a single constrained
  43 element.
  44
  45 =head1 FEEDBACK
  46
  47 =head2 Mailing Lists
  48
  49 User feedback is an integral part of the evolution of this and other
  50 Bioperl modules. Send your comments and suggestions preferably to
  51 the Bioperl mailing list. Your participation is much appreciated.
  52
  53   bioperl-l@bioperl.org                  - General discussion
  54   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  55
  56 =head2 Support
  57
  58 Please direct usage questions or support issues to the mailing list:
  59
  60 I<bioperl-l@bioperl.org>
  61
  62 rather than to the module maintainer directly. Many experienced and
  63 reponsive experts will be able look at the problem and quickly
  64 address it. Please include a thorough description of the problem
  65 with code and data examples if at all possible.
  66
  67 =head2 Reporting Bugs
  68
  69 Report bugs to the Bioperl bug tracking system to help us keep track
  70 of the bugs and their resolution. Bug reports can be submitted via the
  71 web:
  72
  73   https://github.com/bioperl/bioperl-live/issues
  74
  75 =head1 AUTHOR - Sendu Bala
  76
  77 Email bix@sendu.me.uk
  78
  79 =head1 APPENDIX
  80
  81 The rest of the documentation details each of the object methods.
  82 Internal methods are usually preceded with a _
  83
  84 =cut
  85
  86 # Let the code begin...
  87
  88 package Bio::Tools::Phylo::Gerp;
  89 use strict;
  90
  91 use Bio::SeqFeature::Generic;
  92 use Bio::Annotation::SimpleValue;
  93
  94 use base qw(Bio::Root::Root Bio::Root::IO);
  95
  96
  97 =head2 new
  98
  99  Title   : new
 100  Usage   : my $obj = Bio::Tools::Phylo::Gerp->new();
 101  Function: Builds a new Bio::Tools::Phylo::Gerp object
 102  Returns : Bio::Tools::Phylo::Gerp
 103  Args    : -file (or -fh) should contain the contents of a gerpelem .elems file
 104
 105 =cut
 106
 107 sub new {
 108     my ($class, @args) = @_;
 109     my $self = $class->SUPER::new(@args);
 110
 111     $self->_initialize_io(@args);
 112
 113     return $self;
 114 }
 115
 116 =head2 next_result
 117
 118  Title   : next_result
 119  Usage   : $result = $obj->next_result();
 120  Function: Returns the next result available from the input, or undef if there
 121            are no more results.
 122  Returns : Bio::SeqFeature::Annotated object. Features are annotated with a tag
 123            for 'pvalue', and a 'predicted' tag. They have no sequence id unless
 124            the input GERP file is non-standard, with the seq id as the 6th
 125            column.
 126
 127            NB: feature coordinates are alignment columns of the alignment
 128            used to create the result file.
 129  Args    : none
 130
 131 =cut
 132
 133 sub next_result {
 134     my ($self) = @_;
 135
 136     my $line = $self->_readline || return;
 137
 138     while ($line !~ /^\d+\s+\d+\s+\d+\s+\S+\s+\S+\s*(?:\S+\s*)?$/) {
 139         $line = $self->_readline || return;
 140     }
 141
 142     #start   end     length  RS-score        p-value
 143     # code elsewhere adds seq_id on the end (not valid GERP), so we capture that
 144     # if present
 145     my ($start, $end, undef, $rs_score, $p_value, $seq_id) = split(/\s+/, $line);
 146     my $feat = Bio::SeqFeature::Generic->new(
 147         $seq_id ? (-seq_id => $seq_id) : (),
 148         -start        => $start,
 149         -end          => $end,
 150         -strand       => 1,
 151         -score        => $rs_score,
 152         #-type         => 'conserved_region', ***causes 740x increase in SeqFeatureDB storage requirements!
 153         -source       => 'GERP');
 154
 155     my $sv = Bio::Annotation::SimpleValue->new(-tagname => 'predicted', -value => 1);
 156     $feat->annotation->add_Annotation($sv);
 157     $sv = Bio::Annotation::SimpleValue->new(-tagname => 'pvalue', -value => $p_value);
 158     $feat->annotation->add_Annotation($sv);
 159
 160     return $feat;
 161 }
 162
 163 1;