From 1bdb72cf6c17a75d5ce0c5187cff502bcc0c2dce Mon Sep 17 00:00:00 2001 From: martinahansen Date: Thu, 14 Oct 2010 14:04:08 +0000 Subject: [PATCH] ported analyze_seq to ruby git-svn-id: http://biopieces.googlecode.com/svn/trunk@1131 74ccb610-7750-0410-82ae-013aeee3265d --- bp_bin/analyze_seq | 121 +++++++++------------------------------------ 1 file changed, 23 insertions(+), 98 deletions(-) diff --git a/bp_bin/analyze_seq b/bp_bin/analyze_seq index 1418889..d668c99 100755 --- a/bp_bin/analyze_seq +++ b/bp_bin/analyze_seq @@ -1,6 +1,6 @@ -#!/usr/bin/env perl +#!/usr/bin/env ruby -# Copyright (C) 2007-2009 Martin A. Hansen. +# Copyright (C) 2007-2010 Martin A. Hansen. # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License @@ -18,117 +18,42 @@ # http://www.gnu.org/copyleft/gpl.html - -# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - - -# Analyze sequences in the stream. - - # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< +# This program is part of the Biopieces framework (www.biopieces.org). -use warnings; -use strict; -use Maasha::Biopieces; -use Maasha::Common; -use Maasha::Seq; - - -# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - - -my ( $options, $in, $out, $record, $analysis ); - -$options = Maasha::Biopieces::parse_options(); - -$in = Maasha::Biopieces::read_stream( $options->{ "stream_in" } ); -$out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } ); - -while ( $record = Maasha::Biopieces::get_record( $in ) ) -{ - seq_analyze( $record ) if $record->{ "SEQ" }; - - Maasha::Biopieces::put_record( $record, $out ); -} - -Maasha::Biopieces::close_stream( $in ); -Maasha::Biopieces::close_stream( $out ); +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< +# Analyze sequences in the stream. # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< -sub seq_analyze -{ - # Martin A. Hansen, July 2009. +require 'biopieces' +require 'seq' - # Analyzes the composition of the string in the record and appends - # the analysis to the record. +casts = [] - my ( $record, # Biopiece record with a SEQ entry. - ) = @_; +bp = Biopieces.new - # Returns nothing. - - my ( %char_hash, @indels, @alph, $char, $gc, $at, $lc, $max, $indels ); - - %char_hash = Maasha::Common::str_analyze( $record->{ 'SEQ' } ); - - $record->{ 'SEQ_TYPE' } = Maasha::Seq::seq_guess_type( $record->{ 'SEQ' } ); - $record->{ 'SEQ_LEN' } = length $record->{ 'SEQ' }; - - @alph = Maasha::Seq::seq_alph( $record->{ 'SEQ_TYPE' } . "_AMBI" ); - @indels = qw( - ~ . _ ); - - $max = 0; - - foreach $char ( @alph ) - { - $char_hash{ $char } += $char_hash{ lc $char } || 0; - - $record->{ "RES[$char]" } = $char_hash{ $char }; - - $max = $char_hash{ $char } if $char_hash{ $char } > $max; - - $record->{ "RES_SUM" } += $char_hash{ $char }; - } - - $indels = 0; - - map { $record->{ "RES[$_]" } = $char_hash{ $_ }; $indels += $char_hash{ $_ } } @indels; - - if ( $record->{ "SEQ_TYPE" } =~ /DNA|RNA/i ) - { - $gc = $char_hash{ "G" } + $char_hash{ "C" }; - $at = $char_hash{ "A" } + $char_hash{ "T" } + $char_hash{ "U" }; - - $lc = 0; - - map { $lc += $char_hash{ lc $_ } || 0 } @alph; - - $record->{ "MIX_INDEX" } = sprintf( "%.2f", $max / ( $record->{ "SEQ_LEN" } - $indels ) ); - $record->{ "GC%" } = sprintf( "%.2f", 100 * $gc / ( $record->{ "SEQ_LEN" } - $indels ) ); - $record->{ "SOFT_MASK%" } = sprintf( "%.2f", 100 * $lc / ( $record->{ "SEQ_LEN" } - $indels ) ); - $record->{ "HARD_MASK%" } = sprintf( "%.2f", 100 * ( $char_hash{ "n" } + $char_hash{ "N" } ) / ( $record->{ "SEQ_LEN" } - $indels ) ); - $record->{ "MELT_TEMP" } = sprintf( "%.2f", 4 * $gc + 2 * $at ); - } -} - - -# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< +options = bp.parse(ARGV, casts) +bp.each_record do |record| + if record.has_key? :SEQ + seq = Seq.new(record[:SEQ_NAME], record[:SEQ], record[:SEQ_TYPE], record[:SCORE]) + comp = seq.composition -BEGIN -{ - Maasha::Biopieces::status_set(); -} + comp.each_pair do |key,val| + record["RES[#{key}]".to_sym] = val + end + record["HARD_MASK%".to_sym] = seq.hard_mask + record["SOFT_MASK%".to_sym] = seq.soft_mask + record["GC%".to_sym] = ((comp["G"] + comp["C"]).to_f / (seq.len - seq.indels).to_f * 100.0).round(2) + end -END -{ - Maasha::Biopieces::status_log(); -} + bp.puts record +end # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< -- 2.39.2