Bio/Tools/ipcress.pm

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