bp_bin/analyze_seq

   1 #!/usr/bin/env perl
   2
   3 # Copyright (C) 2007-2009 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Analyze sequences in the stream.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use Maasha::Biopieces;
  34 use Maasha::Common;
  35 use Maasha::Seq;
  36
  37
  38 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  39
  40
  41 my ( $options, $in, $out, $record, $analysis );
  42
  43 $options = Maasha::Biopieces::parse_options();
  44
  45 $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
  46 $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
  47
  48 while ( $record = Maasha::Biopieces::get_record( $in ) )
  49 {
  50     seq_analyze( $record ) if $record->{ "SEQ" };
  51
  52     Maasha::Biopieces::put_record( $record, $out );
  53 }
  54
  55 Maasha::Biopieces::close_stream( $in );
  56 Maasha::Biopieces::close_stream( $out );
  57
  58
  59 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  60
  61
  62 sub seq_analyze
  63 {
  64     # Martin A. Hansen, July 2009.
  65
  66     # Analyzes the composition of the string in the record and appends
  67     # the analysis to the record.
  68
  69     my ( $record,   # Biopiece record with a SEQ entry.
  70        ) = @_;
  71
  72     # Returns nothing.
  73
  74     my ( %char_hash, @indels, @alph, $char, $gc, $at, $lc, $max, $indels );
  75
  76     %char_hash = Maasha::Common::str_analyze( $record->{ 'SEQ' } );
  77
  78     $record->{ 'SEQ_TYPE' } = Maasha::Seq::seq_guess_type( $record->{ 'SEQ' } );
  79     $record->{ 'SEQ_LEN' }  = length $record->{ 'SEQ' };
  80
  81     @alph   = Maasha::Seq::seq_alph( $record->{ 'SEQ_TYPE' } . "_AMBI" );
  82     @indels = qw( - ~ . _ );
  83
  84     $max = 0;
  85
  86     foreach $char ( @alph )
  87     {
  88         $char_hash{ $char } += $char_hash{ lc $char } || 0;
  89
  90         $record->{ "RES[$char]" } = $char_hash{ $char };
  91
  92         $max = $char_hash{ $char } if $char_hash{ $char } > $max;
  93
  94         $record->{ "RES_SUM" } += $char_hash{ $char };
  95     }
  96
  97     $indels = 0;
  98
  99     map { $record->{ "RES[$_]" } = $char_hash{ $_ }; $indels += $char_hash{ $_ } } @indels;
 100
 101     if ( $record->{ "SEQ_TYPE" } =~ /DNA|RNA/i )
 102     {
 103         $gc = $char_hash{ "G" } + $char_hash{ "C" };
 104         $at = $char_hash{ "A" } + $char_hash{ "T" } + $char_hash{ "U" };
 105
 106         $lc = 0;
 107
 108         map { $lc += $char_hash{ lc $_ } || 0 } @alph;
 109
 110         $record->{ "MIX_INDEX" }  = sprintf( "%.2f", $max / ( $record->{ "SEQ_LEN" } - $indels ) );
 111         $record->{ "GC%" }        = sprintf( "%.2f", 100 * $gc / ( $record->{ "SEQ_LEN" } - $indels ) );
 112         $record->{ "SOFT_MASK%" } = sprintf( "%.2f", 100 * $lc / ( $record->{ "SEQ_LEN" } - $indels ) );
 113         $record->{ "HARD_MASK%" } = sprintf( "%.2f", 100 * ( $char_hash{ "n" } + $char_hash{ "N" } ) / ( $record->{ "SEQ_LEN" } - $indels ) );
 114         $record->{ "MELT_TEMP" }  = sprintf( "%.2f", 4 * $gc + 2 * $at );
 115     }
 116 }
 117
 118
 119 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 120
 121
 122 BEGIN
 123 {
 124     Maasha::Biopieces::status_set();
 125 }
 126
 127
 128 END
 129 {
 130     Maasha::Biopieces::status_log();
 131 }
 132
 133
 134 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 135
 136
 137 __END__