scripts/seq/bp_seqcut.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 =head1 NAME
   6
   7 I<bp_seqcut.pl>
   8
   9 =head1 USAGE
  10
  11  bp_seqcut.pl [options -h,-s,-e,-f,-w] <FILES>...
  12  bp_seqcut.pl [options -h,-f,-w] s-e <FILES>...
  13
  14    -h this help message
  15    -s which residue to start cutting on
  16    -e which residue to finish cutting on
  17    -f format of the files, defaults to FASTA but you can specify anything supported by SeqIO from BioPerl
  18    -w hard wrap width, this might not be supported depending on which format you are using
  19
  20 =head1 Description
  21
  22 A script to cut FASTA sequences with a given range `fastacut -s 1 -e 10 *.fasta` or `fastacut 1-10 *.fasta`.
  23 This is just a convenience wrapper around the Bio::SeqIO module. Useful if you wish to trim out a section of an
  24 alignment to build a profile of a specific region of sequence.
  25
  26 =head1 AUTHOR
  27
  28 B<Matt Oates> - I<Matt.Oates@bristol.ac.uk>
  29
  30 =head1 FEEDBACK
  31
  32 =head2 Mailing Lists
  33
  34 User feedback is an integral part of the evolution of this and other
  35 Bioperl modules. Send your comments and suggestions preferably to
  36 the Bioperl mailing list.  Your participation is much appreciated.
  37
  38   bioperl-l@bioperl.org                  - General discussion
  39   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  40
  41 =head2 Reporting Bugs
  42
  43 Report bugs to the Bioperl bug tracking system to help us keep track
  44 of the bugs and their resolution. Bug reports can be submitted via
  45 email or the web:
  46
  47   https://github.com/bioperl/bioperl-live/issues
  48
  49 =head1 EDIT HISTORY
  50
  51 2010-11-22 - Matt Oates
  52         First features added.
  53 =cut
  54
  55
  56
  57 # Includes
  58 =head1 DEPENDANCY
  59 B<Getopt::Long> Used to parse command line options.
  60 B<Pod::Usage> Used for usage and help output.
  61 B<Bio::SeqIO> Used to cut up sequences and parse FASTA.
  62 =cut
  63 use Getopt::Long;                     #Deal with command line options
  64 use Pod::Usage;                       #Print a usage man page from the POD comments after __END__
  65 use Bio::SeqIO;
  66
  67 # Command Line Options
  68 my $help;    #Same again but this time should we output the POD man page defined after __END__
  69 my $format = "Fasta";
  70 my $start;
  71 my $end;
  72 my $width = 72; #Default for Jalview output
  73 my $outfile = '/dev/stdout';
  74
  75 #Set command line flags and parameters.
  76 GetOptions("help|h!" => \$help,
  77            "start|s=s" => \$start,
  78            "format|f=s" => \$format,
  79            "end|e=s" => \$end,
  80            "width|w=s" => \$width,
  81            "outfile|o=s" => \$outfile,
  82         ) or die "Fatal Error: Problem parsing command-line ".$!;
  83
  84 #Get other command line arguments that weren't optional flags.
  85 ($start,$end) = split (/-/, shift) unless ($start and $end);
  86 my @files = @ARGV;
  87
  88 #Print out some help if it was asked for or if no arguments were given.
  89 pod2usage(-exitstatus => 0, -verbose => 2) if $help;
  90
  91 pod2usage(-exitstatus => 0, -verbose => 1, -msg => 'Please specify the sequence files you wish to cut.')
  92     unless scalar @files;
  93
  94 pod2usage(-exitstatus => 0, -verbose => 1, -msg => 'Please specify the region you wish to cut -s 1 -e 10 or 1-10.')
  95     unless defined $end;
  96
  97 my $out = Bio::SeqIO->newFh(-file => ">$outfile", -format => $format) or die "Couldn't open selected output sequence file.";
  98
  99 #Open and iterate over all sequence in all files
 100 foreach my $file (@files) {
 101         my $in  = Bio::SeqIO->new(-file => $file, -format => $format);
 102         while ( my $seq = $in->next_seq() ) {
 103             #Alter the ID to be postfixed with '/s-e'
 104             $seq->display_id($seq->display_id."/$start-$end");
 105             #Edit the sequence we have cut out
 106                 my $sequence = $seq->subseq($start,$end);
 107                 $sequence =~ s/([^\n]{0,$width})/$1\n/gi;
 108                 chomp $sequence;
 109             $seq->seq($sequence);
 110             #Print the sequence back out
 111                 print $out $seq;
 112         }
 113 }
 114
 115 1;
 116 __END__
 117