scripts/seqstats/aacomp.PLS

   1 #!perl -w
   2 use strict;
   3 use Carp;
   4
   5 use Bio::SeqIO;
   6 use Getopt::Long;
   7 use Bio::SeqUtils;
   8 use Bio::Tools::IUPAC;
   9
  10 my $table = new Bio::SeqUtils;
  11 my @BASES = $table->valid_aa(0);
  12 my %all = $table->valid_aa(2);
  13 my ($file,$format,$help) = ( undef, 'fasta');
  14 GetOptions(
  15            'i|in:s'  => \$file,
  16            'f|format:s' => \$format,
  17            'h|help|?'  => \$help,
  18            );
  19
  20 my $USAGE = "usage: aacomp [-f format] filename\n\tdefault format is fasta\n";
  21 $file = shift unless $file;
  22
  23 die $USAGE if $help;
  24
  25 my $seqin;
  26 if( defined $file ) {
  27     print "Could not open file [$file]\n$USAGE" and exit unless -e $file;
  28     $seqin = new Bio::SeqIO(-format => $format,
  29                             -file   => $file);
  30 } else {
  31     $seqin = new Bio::SeqIO(-format => $format,
  32                             -fh     => \*STDIN);
  33 }
  34
  35 my %composition;
  36 my $total;
  37 foreach my $base ( @BASES ) {
  38     $composition{$base} = 0;
  39 }
  40 while ( my $seq = $seqin->next_seq ) {
  41     if( $seq->alphabet ne 'protein' ) {
  42         confess("Must only provide amino acid sequences to aacomp...skipping this seq");
  43         next;
  44     }
  45     foreach my $base ( split(//,$seq->seq()) ) {
  46         $composition{uc $base}++;
  47         $total++;
  48     }
  49 }
  50
  51 printf("%d aa\n",$total);
  52 printf("%5s %4s\n", 'aa', '#' );
  53 my $ct = 0;
  54 foreach my $base ( @BASES ) {
  55     printf(" %s %s %3d\n", $base, $all{$base}, $composition{$base} );
  56     $ct += $composition{$base};
  57 }
  58 printf( "%6s %s\n", '','-'x5);
  59 printf( "%6s %3d\n", '',$ct);
  60
  61
  62 __END__
  63
  64
  65 =head1 NAME
  66
  67 aacomp - amino acid composition of protein sequences
  68
  69 =head1 SYNOPSIS
  70
  71   aacomp [-f/--format FORMAT] [-h/--help] filename
  72   or
  73   aacomp [-f/--format FORMAT] < filename
  74   or
  75   aacomp [-f/--format FORMAT] -i filename
  76
  77 =head1 DESCRIPTION
  78
  79 This scripts prints out the count of amino acids over all protein
  80 sequences from the input file.
  81
  82 =head1 OPTIONS
  83
  84 The default sequence format is fasta.
  85
  86 The sequence input can be provided using any of the three methods:
  87
  88 =over 3
  89
  90 =item unnamed argument
  91
  92   aacomp filename
  93
  94 =item named argument
  95
  96   aacomp -i filename
  97
  98 =item standard input
  99
 100   aacomp < filename
 101
 102 =back
 103
 104 =head1 FEEDBACK
 105
 106 =head2 Mailing Lists
 107
 108 User feedback is an integral part of the evolution of this and other
 109 Bioperl modules. Send your comments and suggestions preferably to
 110 the Bioperl mailing list.  Your participation is much appreciated.
 111
 112   bioperl-l@bioperl.org                  - General discussion
 113   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 114
 115 =head2 Reporting Bugs
 116
 117 Report bugs to the Bioperl bug tracking system to help us keep track
 118 of the bugs and their resolution. Bug reports can be submitted via the
 119 web:
 120
 121   https://redmine.open-bio.org/projects/bioperl/
 122
 123 =head1 AUTHOR - Jason Stajich
 124
 125 Email jason@bioperl.org
 126
 127 =head1 HISTORY
 128
 129 Based on aacomp.c from an old version of EMBOSS
 130
 131 =cut