lib/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
  98 use strict;
  99
 100 use Bio::SeqFeature::Generic;
 101
 102 use base qw(Bio::Root::Root);
 103
 104 =head2 new
 105
 106  Title   : new
 107  Usage   : my $ipcress = Bio::Tools::ipcress->new(-file => $file,
 108                                            -primary => $fprimary,
 109                                            -source => $fsource,
 110                                            -groupclass => $fgroupclass);
 111  Function: Initializes a new ipcress parser
 112  Returns : Bio::Tools::ipcress
 113  Args    : -fh   => filehandle
 114            OR
 115            -file => filename
 116
 117            -primary => a string to be used as the common value for
 118                        each features '-primary' tag.  Defaults to
 119                        the sequence ontology term 'PCR_product'.
 120                        (This in turn maps to the GFF 'type'
 121                        tag (aka 'method')).
 122
 123             -source => a string to be used as the common value for
 124                        each features '-source' tag.  Defaults to
 125                        'ipcress'. (This in turn maps to the GFF 'source'
 126                        tag)
 127
 128              -groupclass => a string to be used as the name of the tag
 129                            which will hold the sts marker namefirst
 130                            attribute.  Defaults to 'name'.
 131
 132 =cut
 133
 134 sub new {
 135   my($class,@args) = @_;
 136
 137   my $self = $class->SUPER::new(@args);
 138   my ($primary, $source,
 139       $groupclass, $file, $fh) = $self->_rearrange([qw(PRIMARY
 140                                                        SOURCE
 141                                                        GROUPCLASS
 142                                                        FILE FH)],@args);
 143   $self->primary(defined $primary ? $primary : 'PCR_product');
 144   $self->source(defined $source ? $source : 'ipcress');
 145   $self->groupclass(defined $groupclass ? $groupclass : 'name');
 146
 147   local $/ = 'Ipcress result';
 148   my @result;
 149
 150   if ($file) {
 151       open my $FH, '<', $file or $self->throw("Could not read file '$file': $!");
 152       @result = (<$FH>);
 153       close $FH;
 154   }
 155   elsif ($fh) {
 156       @result = (<$fh>);
 157   }
 158   else {
 159       $self->throw("Bio::Tools::ipcress: no input file");
 160   }
 161
 162
 163   shift @result;
 164
 165   $self->{result} = \@result;
 166
 167   return $self;
 168 }
 169
 170 =head2 next_feature
 171
 172  Title   : next_feature
 173  Usage   : $seqfeature = $obj->next_feature();
 174  Function: Returns the next feature available in the analysis result, or
 175            undef if there are no more features.
 176  Example :
 177  Returns : A Bio::SeqFeatureI implementing object, or undef if there are no
 178            more features.
 179  Args    : none
 180
 181 =cut
 182
 183 sub next_feature {
 184     my ($self) = @_;
 185     my $result = shift @{$self->{result}};
 186     return unless defined($result);
 187
 188     chomp $result;
 189     my @lines = split "\n", $result;
 190     my ($ipcress) = grep /ipcress: /, @lines;
 191
 192     my (undef,$seqname,$mkrname,$length,undef,$start,$mismatchL,
 193        undef,undef,$mismatchR,$desc) = split /\s+/, $ipcress;
 194
 195     my $end = $start + $length;
 196     $start += 1;
 197
 198     my $strand = $desc eq 'forward' ? '+' : $desc eq 'revcomp' ? '-' : 0;
 199
 200     my ($left)  = grep /\# forward/, @lines;
 201     $left  =~ s/[^A-Z]+//g;
 202     my ($right) = grep /\# revcomp/, @lines;
 203     $right =~ s/[^A-Z]+//g;
 204     $right = reverse $right;
 205
 206     # if there are multiple hits, increment the name for
 207     # the groupclass
 208     if (++$self->{seen}->{$mkrname} > 1) {
 209         $mkrname .= "\.$self->{seen}->{$mkrname}";
 210     }
 211
 212
 213     my $markerfeature = Bio::SeqFeature::Generic->new
 214         ( '-start'   => $start,
 215           '-end'     => $end,
 216           '-strand'  => $strand,
 217           '-source'  => $self->source,
 218           '-primary' => $self->primary,
 219           '-seq_id'  => $seqname,
 220           '-tag'     => {
 221               $self->groupclass => $mkrname,
 222           });
 223
 224     if (!$strand) {
 225         $markerfeature->add_tag_value('Note' => "bad product: single primer amplification");
 226     }
 227
 228     $markerfeature->add_tag_value('left_primer' => $left);
 229     $markerfeature->add_tag_value('right_primer' => $right);
 230     $markerfeature->add_tag_value('left_mismatches' => $mismatchL) if $mismatchL;
 231     $markerfeature->add_tag_value('right_mismatches' => $mismatchR) if $mismatchR;
 232
 233     return $markerfeature;
 234 }
 235
 236 =head2 source
 237
 238  Title   : source
 239  Usage   : $obj->source($newval)
 240  Function:
 241  Example :
 242  Returns : value of source (a scalar)
 243  Args    : on set, new value (a scalar or undef, optional)
 244
 245
 246 =cut
 247
 248 sub source{
 249     my $self = shift;
 250     return $self->{'_source'} = shift if @_;
 251     return $self->{'_source'};
 252 }
 253
 254 =head2 primary
 255
 256  Title   : primary
 257  Usage   : $obj->primary($newval)
 258  Function:
 259  Example :
 260  Returns : value of primary (a scalar)
 261  Args    : on set, new value (a scalar or undef, optional)
 262
 263
 264 =cut
 265
 266 sub primary{
 267     my $self = shift;
 268     return $self->{'_primary'} = shift if @_;
 269     return $self->{'_primary'};
 270 }
 271
 272 =head2 groupclass
 273
 274  Title   : groupclass
 275  Usage   : $obj->groupclass($newval)
 276  Function:
 277  Example :
 278  Returns : value of groupclass (a scalar)
 279  Args    : on set, new value (a scalar or undef, optional)
 280
 281
 282 =cut
 283
 284 sub groupclass{
 285     my $self = shift;
 286
 287     return $self->{'_groupclass'} = shift if @_;
 288     return $self->{'_groupclass'};
 289 }
 290
 291 1;