Bio/Tools/ipcress.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tools::ipcress
   4 #
   5 # Cared for by Sheldon McKay <mckays@cshl.edu>
   6 #
   7 # Copyright Sheldon McKay
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::Tools::ipcress - Parse ipcress output and make features
  16
  17 =head1 SYNOPSIS
  18
  19     # A simple annotation pipeline wrapper for ipcress data
  20     # assuming ipcress data is already generated in file seq1.ipcress
  21     # and sequence data is in fasta format in file called seq1.fa
  22
  23     use Bio::Tools::ipcress;
  24     use Bio::SeqIO;
  25     my $parser = Bio::Tools::ipcress->new(-file => 'seq1.ipcress');
  26     my $seqio = Bio::SeqIO->new(-format => 'fasta', -file => 'seq1.fa');
  27     my $seq = $seqio->next_seq || die("cannot get a seq object from SeqIO");
  28
  29     while( my $feat = $parser->next_feature ) {
  30         # add ipcress annotation to a sequence
  31         $seq->add_SeqFeature($feat);
  32     }
  33     my $seqout = Bio::SeqIO->new(-format => 'embl');
  34     $seqout->write_seq($seq);
  35
  36
  37 =head1 DESCRIPTION
  38
  39 This object serves as a parser for ipcress data, creating a
  40 Bio::SeqFeatureI for each ipcress hit.  These can be processed or added
  41 as annotation to an existing Bio::SeqI object for the purposes of
  42 automated annotation.
  43
  44 This module is adapted from the Bio::Tools::EPCR module
  45 written by Jason Stajich (jason-at-bioperl.org).
  46
  47 Ipcress is available through Guy Slater's Exonerate package
  48 http://www.ebi.ac.uk/~guy/exonerate/
  49
  50 =head1 FEEDBACK
  51
  52 =head2 Mailing Lists
  53
  54 User feedback is an integral part of the evolution of this and other
  55 Bioperl modules. Send your comments and suggestions preferably to
  56 the Bioperl mailing list.  Your participation is much appreciated.
  57
  58   bioperl-l@bioperl.org                  - General discussion
  59   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  60
  61 =head2 Reporting Bugs
  62
  63 Report bugs to the Bioperl bug tracking system to help us keep track
  64 of the bugs and their resolution. Bug reports can be submitted via the
  65 web:
  66
  67   http://bugzilla.open-bio.org/
  68
  69 =head1 AUTHOR - Sheldon McKay
  70
  71 Email mckays@cshl.edu
  72
  73 =head1 APPENDIX
  74
  75 The rest of the documentation details each of the object methods.
  76 Internal methods are usually preceded with a _
  77
  78 =cut
  79
  80
  81 # Let the code begin...
  82
  83
  84 package Bio::Tools::ipcress;
  85 use strict;
  86
  87 use Bio::SeqFeature::Generic;
  88
  89 use base qw(Bio::Root::Root);
  90
  91 =head2 new
  92
  93  Title   : new
  94  Usage   : my $ipcress = Bio::Tools::ipcress->new(-file => $file,
  95                                            -primary => $fprimary,
  96                                            -source => $fsource,
  97                                            -groupclass => $fgroupclass);
  98  Function: Initializes a new ipcress parser
  99  Returns : Bio::Tools::ipcress
 100  Args    : -fh   => filehandle
 101            OR
 102            -file => filename
 103
 104            -primary => a string to be used as the common value for
 105                        each features '-primary' tag.  Defaults to
 106                        the sequence ontology term 'PCR_product'.
 107                        (This in turn maps to the GFF 'type'
 108                        tag (aka 'method')).
 109
 110             -source => a string to be used as the common value for
 111                        each features '-source' tag.  Defaults to
 112                        'ipcress'. (This in turn maps to the GFF 'source'
 113                        tag)
 114
 115              -groupclass => a string to be used as the name of the tag
 116                            which will hold the sts marker namefirst
 117                            attribute.  Defaults to 'name'.
 118
 119 =cut
 120
 121 sub new {
 122   my($class,@args) = @_;
 123
 124   my $self = $class->SUPER::new(@args);
 125   my ($primary, $source,
 126       $groupclass, $file, $fh) = $self->_rearrange([qw(PRIMARY
 127                                                        SOURCE
 128                                                        GROUPCLASS
 129                                                        FILE FH)],@args);
 130   $self->primary(defined $primary ? $primary : 'PCR_product');
 131   $self->source(defined $source ? $source : 'ipcress');
 132   $self->groupclass(defined $groupclass ? $groupclass : 'name');
 133
 134   local $/ = 'Ipcress result';
 135   my @result;
 136
 137   if ($file) {
 138       open FH, $file;
 139       @result = (<FH>);
 140       close FH;
 141   }
 142   elsif ($fh) {
 143       @result = (<$fh>);
 144   }
 145   else {
 146       $self->throw("Bio::Tools::ipcress: no input file");
 147   }
 148
 149
 150   shift @result;
 151
 152   $self->{result} = \@result;
 153
 154   return $self;
 155 }
 156
 157 =head2 next_feature
 158
 159  Title   : next_feature
 160  Usage   : $seqfeature = $obj->next_feature();
 161  Function: Returns the next feature available in the analysis result, or
 162            undef if there are no more features.
 163  Example :
 164  Returns : A Bio::SeqFeatureI implementing object, or undef if there are no
 165            more features.
 166  Args    : none
 167
 168 =cut
 169
 170 sub next_feature {
 171     my ($self) = @_;
 172     my $result = shift @{$self->{result}};
 173     return unless defined($result);
 174
 175     chomp $result;
 176     my @lines = split "\n", $result;
 177     my ($ipcress) = grep /ipcress: /, @lines;
 178
 179     my (undef,$seqname,$mkrname,$length,undef,$start,$mismatchL,
 180        undef,undef,$mismatchR,$desc) = split /\s+/, $ipcress;
 181
 182     my $end = $start + $length;
 183     $start += 1;
 184
 185     my $strand = $desc eq 'forward' ? '+' : $desc eq 'revcomp' ? '-' : 0;
 186
 187     my ($left)  = grep /\# forward/, @lines;
 188     $left  =~ s/[^A-Z]+//g;
 189     my ($right) = grep /\# revcomp/, @lines;
 190     $right =~ s/[^A-Z]+//g;
 191     $right = reverse $right;
 192
 193     # if there are multiple hits, increment the name for
 194     # the groupclass
 195     if (++$self->{seen}->{$mkrname} > 1) {
 196         $mkrname .= "\.$self->{seen}->{$mkrname}";
 197     }
 198
 199
 200     my $markerfeature = Bio::SeqFeature::Generic->new
 201         ( '-start'   => $start,
 202           '-end'     => $end,
 203           '-strand'  => $strand,
 204           '-source'  => $self->source,
 205           '-primary' => $self->primary,
 206           '-seq_id'  => $seqname,
 207           '-tag'     => {
 208               $self->groupclass => $mkrname,
 209           });
 210
 211     if (!$strand) {
 212         $markerfeature->add_tag_value('Note' => "bad product: single primer amplification");
 213     }
 214
 215     $markerfeature->add_tag_value('left_primer' => $left);
 216     $markerfeature->add_tag_value('right_primer' => $right);
 217     $markerfeature->add_tag_value('left_mismatches' => $mismatchL) if $mismatchL;
 218     $markerfeature->add_tag_value('right_mismatches' => $mismatchR) if $mismatchR;
 219
 220     return $markerfeature;
 221 }
 222
 223 =head2 source
 224
 225  Title   : source
 226  Usage   : $obj->source($newval)
 227  Function:
 228  Example :
 229  Returns : value of source (a scalar)
 230  Args    : on set, new value (a scalar or undef, optional)
 231
 232
 233 =cut
 234
 235 sub source{
 236     my $self = shift;
 237     return $self->{'_source'} = shift if @_;
 238     return $self->{'_source'};
 239 }
 240
 241 =head2 primary
 242
 243  Title   : primary
 244  Usage   : $obj->primary($newval)
 245  Function:
 246  Example :
 247  Returns : value of primary (a scalar)
 248  Args    : on set, new value (a scalar or undef, optional)
 249
 250
 251 =cut
 252
 253 sub primary{
 254     my $self = shift;
 255     return $self->{'_primary'} = shift if @_;
 256     return $self->{'_primary'};
 257 }
 258
 259 =head2 groupclass
 260
 261  Title   : groupclass
 262  Usage   : $obj->groupclass($newval)
 263  Function:
 264  Example :
 265  Returns : value of groupclass (a scalar)
 266  Args    : on set, new value (a scalar or undef, optional)
 267
 268
 269 =cut
 270
 271 sub groupclass{
 272     my $self = shift;
 273
 274     return $self->{'_groupclass'} = shift if @_;
 275     return $self->{'_groupclass'};
 276 }
 277
 278 1;