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 use strict;
  98
  99 use Bio::SeqIO;
 100 use Bio::SeqFeature::Generic;
 101
 102
 103 use base qw(Bio::Root::Root);
 104
 105
 106 =head2 new
 107
 108  Title   : new
 109  Usage   : my $ispcr = Bio::Tools::isPcr->new( -file => $file,
 110                                               -primary => $fprimary,
 111                                               -source => $fsource,
 112                                               -groupclass => $fgroupclass);
 113
 114  Function: Initializes a new isPcr parser
 115  Returns : Bio::Tools::isPcr
 116  Args    : -fh   => filehandle
 117            OR
 118            -file => filename
 119
 120            -primary => a string to be used as the common value for
 121                        each features '-primary' tag.  Defaults to
 122                        the sequence ontology term 'PCR_product'.
 123                        (This in turn maps to the GFF 'type'
 124                        tag (aka 'method')).
 125
 126             -source => a string to be used as the common value for
 127                        each features '-source' tag.  Defaults to
 128                        'isPcr'. (This in turn maps to the GFF 'source'
 129                        tag)
 130
 131             -groupclass => a string to be used as the name of the tag
 132                            which will hold the sts marker namefirst
 133                            attribute.  Defaults to 'name'.
 134
 135 =cut
 136
 137
 138 sub new {
 139   my($class,@args) = @_;
 140
 141   my $self = $class->SUPER::new(@args);
 142   my ($primary,$source,$groupclass) =
 143       $self->_rearrange([qw/PRIMARY SOURCE GROUPCLASS/],@args);
 144   $self->primary(defined $primary ? $primary : 'PCR_product');
 145   $self->source(defined $source ? $source : 'isPcr');
 146   $self->groupclass(defined $groupclass ? $groupclass : 'name');
 147
 148   # default output for isPcr is fasta format
 149   $self->{io} = Bio::SeqIO->new(-format => 'fasta', @args);
 150
 151   return $self;
 152 }
 153
 154 =head2 next_feature
 155
 156  Title   : next_feature
 157  Usage   : $seqfeature = $obj->next_feature();
 158  Function: Returns the next feature available in the analysis result, or
 159            undef if there are no more features.
 160  Example :
 161  Returns : A Bio::SeqFeatureI implementing object, or undef if there are no
 162            more features.
 163  Args    : none
 164
 165 =cut
 166
 167 sub next_feature {
 168     my ($self) = @_;
 169     my $result = $self->{io}->next_seq;
 170     return unless defined $result;
 171
 172     my ($seqname,$location)    = split ':', $result->primary_id;
 173     my ($pcrname,$left,$right) = split /\s+/, $result->desc;
 174     my ($start,$strand,$end)   = $location =~ /^(\d+)([-+])(\d+)$/;
 175     my $amplicon = $result->seq;
 176
 177     # if there are multiple hits, increment the name for
 178     # the groupclass
 179     if (++$self->{seen}->{$pcrname} > 1) {
 180         $pcrname .= "\.$self->{seen}->{$pcrname}";
 181     }
 182
 183     my $tags = {
 184         $self->groupclass => $pcrname,
 185         amplicon          => $amplicon,
 186         left_primer       => $left,
 187         right_primer       => $right
 188         };
 189
 190     my $markerfeature = Bio::SeqFeature::Generic->new(
 191                                                        '-start'   => $start,
 192                                                        '-end'     => $end,
 193                                                        '-strand'  => $strand,
 194                                                        '-source'  => $self->source,
 195                                                        '-primary' => $self->primary,
 196                                                        '-seq_id'  => $seqname,
 197                                                        '-tag'     => $tags
 198                                                        );
 199
 200     return $markerfeature;
 201 }
 202
 203 =head2 source
 204
 205  Title   : source
 206  Usage   : $obj->source($newval)
 207  Function:
 208  Example :
 209  Returns : value of source (a scalar)
 210  Args    : on set, new value (a scalar or undef, optional)
 211
 212
 213 =cut
 214
 215 sub source{
 216     my $self = shift;
 217     return $self->{'_source'} = shift if @_;
 218     return $self->{'_source'};
 219 }
 220
 221 =head2 primary
 222
 223  Title   : primary
 224  Usage   : $obj->primary($newval)
 225  Function:
 226  Example :
 227  Returns : value of primary (a scalar)
 228  Args    : on set, new value (a scalar or undef, optional)
 229
 230
 231 =cut
 232
 233 sub primary{
 234     my $self = shift;
 235     return $self->{'_primary'} = shift if @_;
 236     return $self->{'_primary'};
 237 }
 238
 239 =head2 groupclass
 240
 241  Title   : groupclass
 242  Usage   : $obj->groupclass($newval)
 243  Function:
 244  Example :
 245  Returns : value of groupclass (a scalar)
 246  Args    : on set, new value (a scalar or undef, optional)
 247
 248
 249 =cut
 250
 251 sub groupclass{
 252     my $self = shift;
 253
 254     return $self->{'_groupclass'} = shift if @_;
 255     return $self->{'_groupclass'};
 256 }
 257
 258 1;