lib/Bio/Tools/isPcr.pm

   1 #
   2 # BioPerl module for Bio::Tools::isPcr
   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::isPcr - Parse isPcr output and make features
  17
  18 =head1 SYNOPSIS
  19
  20     # A simple annotation pipeline wrapper for isPcr data
  21     # assuming isPcr data is already generated in file seq1.isPcr
  22     # and sequence data is in fasta format in file called seq1.fa
  23
  24     # Note: this parser is meant for the default fasta output from
  25     # isPcr.  bed and psl output formats are not supported.
  26
  27     use Bio::Tools::IsPcr;
  28     use Bio::SeqIO;
  29     my $parser = Bio::Tools::isPcr->new(-file => 'seq1.isPcr');
  30     my $seqio = Bio::SeqIO->new(-format => 'fasta', -file => 'seq1.fa');
  31     my $seq = $seqio->next_seq || die("cannot get a seq object from SeqIO");
  32
  33     while( my $feat = $parser->next_feature ) {
  34         # add isPcr annotation to a sequence
  35         $seq->add_SeqFeature($feat);
  36     }
  37     my $seqout = Bio::SeqIO->new(-format => 'embl');
  38     $seqout->write_seq($seq);
  39
  40
  41 =head1 DESCRIPTION
  42
  43 This object serves as a parser for isPcr data (in the default fasta
  44 format), creating a Bio::SeqFeatureI for each isPcr hit.
  45 These can be processed or added as annotation to an existing
  46 Bio::SeqI object for the purposes of automated annotation.
  47
  48 This module is adapted from the Bio::Tools::EPCR module
  49 written by Jason Stajich (jason-at-bioperl.org).
  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::isPcr;
  97
  98 use strict;
  99
 100 use Bio::SeqIO;
 101 use Bio::SeqFeature::Generic;
 102
 103
 104 use base qw(Bio::Root::Root);
 105
 106
 107 =head2 new
 108
 109  Title   : new
 110  Usage   : my $ispcr = Bio::Tools::isPcr->new( -file => $file,
 111                                               -primary => $fprimary,
 112                                               -source => $fsource,
 113                                               -groupclass => $fgroupclass);
 114
 115  Function: Initializes a new isPcr parser
 116  Returns : Bio::Tools::isPcr
 117  Args    : -fh   => filehandle
 118            OR
 119            -file => filename
 120
 121            -primary => a string to be used as the common value for
 122                        each features '-primary' tag.  Defaults to
 123                        the sequence ontology term 'PCR_product'.
 124                        (This in turn maps to the GFF 'type'
 125                        tag (aka 'method')).
 126
 127             -source => a string to be used as the common value for
 128                        each features '-source' tag.  Defaults to
 129                        'isPcr'. (This in turn maps to the GFF 'source'
 130                        tag)
 131
 132             -groupclass => a string to be used as the name of the tag
 133                            which will hold the sts marker namefirst
 134                            attribute.  Defaults to 'name'.
 135
 136 =cut
 137
 138
 139 sub new {
 140   my($class,@args) = @_;
 141
 142   my $self = $class->SUPER::new(@args);
 143   my ($primary,$source,$groupclass) =
 144       $self->_rearrange([qw/PRIMARY SOURCE GROUPCLASS/],@args);
 145   $self->primary(defined $primary ? $primary : 'PCR_product');
 146   $self->source(defined $source ? $source : 'isPcr');
 147   $self->groupclass(defined $groupclass ? $groupclass : 'name');
 148
 149   # default output for isPcr is fasta format
 150   $self->{io} = Bio::SeqIO->new(-format => 'fasta', @args);
 151
 152   return $self;
 153 }
 154
 155 =head2 next_feature
 156
 157  Title   : next_feature
 158  Usage   : $seqfeature = $obj->next_feature();
 159  Function: Returns the next feature available in the analysis result, or
 160            undef if there are no more features.
 161  Example :
 162  Returns : A Bio::SeqFeatureI implementing object, or undef if there are no
 163            more features.
 164  Args    : none
 165
 166 =cut
 167
 168 sub next_feature {
 169     my ($self) = @_;
 170     my $result = $self->{io}->next_seq;
 171     return unless defined $result;
 172
 173     my ($seqname,$location)    = split ':', $result->primary_id;
 174     my ($pcrname,$left,$right) = split /\s+/, $result->desc;
 175     my ($start,$strand,$end)   = $location =~ /^(\d+)([-+])(\d+)$/;
 176     my $amplicon = $result->seq;
 177
 178     # if there are multiple hits, increment the name for
 179     # the groupclass
 180     if (++$self->{seen}->{$pcrname} > 1) {
 181         $pcrname .= "\.$self->{seen}->{$pcrname}";
 182     }
 183
 184     my $tags = {
 185         $self->groupclass => $pcrname,
 186         amplicon          => $amplicon,
 187         left_primer       => $left,
 188         right_primer       => $right
 189         };
 190
 191     my $markerfeature = Bio::SeqFeature::Generic->new(
 192                                                        '-start'   => $start,
 193                                                        '-end'     => $end,
 194                                                        '-strand'  => $strand,
 195                                                        '-source'  => $self->source,
 196                                                        '-primary' => $self->primary,
 197                                                        '-seq_id'  => $seqname,
 198                                                        '-tag'     => $tags
 199                                                        );
 200
 201     return $markerfeature;
 202 }
 203
 204 =head2 source
 205
 206  Title   : source
 207  Usage   : $obj->source($newval)
 208  Function:
 209  Example :
 210  Returns : value of source (a scalar)
 211  Args    : on set, new value (a scalar or undef, optional)
 212
 213
 214 =cut
 215
 216 sub source{
 217     my $self = shift;
 218     return $self->{'_source'} = shift if @_;
 219     return $self->{'_source'};
 220 }
 221
 222 =head2 primary
 223
 224  Title   : primary
 225  Usage   : $obj->primary($newval)
 226  Function:
 227  Example :
 228  Returns : value of primary (a scalar)
 229  Args    : on set, new value (a scalar or undef, optional)
 230
 231
 232 =cut
 233
 234 sub primary{
 235     my $self = shift;
 236     return $self->{'_primary'} = shift if @_;
 237     return $self->{'_primary'};
 238 }
 239
 240 =head2 groupclass
 241
 242  Title   : groupclass
 243  Usage   : $obj->groupclass($newval)
 244  Function:
 245  Example :
 246  Returns : value of groupclass (a scalar)
 247  Args    : on set, new value (a scalar or undef, optional)
 248
 249
 250 =cut
 251
 252 sub groupclass{
 253     my $self = shift;
 254
 255     return $self->{'_groupclass'} = shift if @_;
 256     return $self->{'_groupclass'};
 257 }
 258
 259 1;