scripts/seq/extract_feature_seq.PLS

   1 #!perl -w
   2 use strict;
   3 # $Id: extract_feature_seq.PLS,v 1.5 2006-07-04 22:23:29 mauricio Exp $
   4 # Author Jason Stajich <jason@bioperl.org>
   5
   6 =head1 NAME
   7
   8 extract_feature_seq - extract the corresponding sequence for a specified feature type
   9
  10 =head1 SYNOPSIS
  11
  12 extract_feature_seq.PLS -i file --format genbank --feature=CDS -o output.fa
  13
  14 =head1 DESCRIPTION
  15
  16 This script will extract the sequence for all the features you specify.
  17
  18 =head1 FEEDBACK
  19
  20 =head2 Mailing Lists
  21
  22 User feedback is an integral part of the evolution of this and other
  23 Bioperl modules. Send your comments and suggestions preferably to
  24 the Bioperl mailing list.  Your participation is much appreciated.
  25
  26   bioperl-l@bioperl.org                  - General discussion
  27   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  28
  29 =head2 Reporting Bugs
  30
  31 Report bugs to the Bioperl bug tracking system to help us keep track
  32 of the bugs and their resolution. Bug reports can be submitted via
  33 email or the web:
  34
  35   http://bugzilla.open-bio.org/
  36
  37 =head1 AUTHOR
  38
  39  Jason Stajich E<lt>jason-at-bioperl-dot-orgE<gt>
  40
  41 =cut
  42
  43 use Bio::SeqIO;
  44 use Getopt::Long;
  45
  46 my ($input,$format,$featuretype,$output);
  47 $featuretype ='CDS';
  48 GetOptions(
  49            'i|input:s' => \$input,
  50            'format:s'  => \$format,
  51            'feature:s' => \$featuretype,
  52            'o|output:s'=> \$output);
  53
  54 $input || shift if @ARGV;
  55
  56 my $in = new Bio::SeqIO(-file => $input,
  57                         -format => $format);
  58 my $out;
  59 if ($output ) {
  60     $out = new Bio::SeqIO(-file => ">$output");
  61 } else {
  62     $out = new Bio::SeqIO(); # use STDOUT for output
  63 }
  64
  65 my $count = 1;
  66 while( my $seq = $in->next_seq ) {
  67     foreach my $f ( grep { $_->primary_tag =~ /$featuretype/i }
  68                     $seq->get_SeqFeatures ) {
  69         my $s = $f->spliced_seq;
  70         if( $featuretype =~ /gene|CDS/ ) {
  71             $s->display_id($f->has_tag('gene') ? join(',',sort $f->each_tag_value('gene')) :
  72                            $f->has_tag('label') ? join(',',$f->each_tag_value('label')):
  73                            $s->display_id);
  74         } else {
  75             $s->display_id(sprintf("%s_%s_%d",
  76                                    $seq->display_id,
  77                                    $f->primary_tag,
  78                                    $count++));
  79         }
  80         $out->write_seq($s);
  81     }
  82 }
  83
  84 __END__