forester/ruby/evoruby/lib/evo/tool/hmmscan_summary.rb

   1 #
   2 # = lib/evo/apps/hmmscan_parser.rb - HmmscanParser class
   3 #
   4 # Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $
   8 #
   9 # last modified: 11/24/2009
  10
  11 require 'lib/evo/util/constants'
  12 require 'lib/evo/util/util'
  13 require 'lib/evo/util/command_line_arguments'
  14 require 'lib/evo/io/parser/hmmscan_parser'
  15
  16 module Evoruby
  17
  18   class HmmscanSummary
  19
  20     PRG_NAME       = "hsp"
  21     PRG_VERSION    = "2.000"
  22     PRG_DESC       = "hmmscan parser"
  23     PRG_DATE       = "2012.10.19"
  24     COPYRIGHT      = "2012 Christian M Zmasek"
  25     CONTACT        = "phylosoft@gmail.com"
  26     WWW            = "www.phylosoft.org"
  27
  28     DELIMITER_OPTION              = "d"
  29     I_E_VALUE_THRESHOLD_OPTION    = "e"
  30     FS_E_VALUE_THRESHOLD_OPTION   = "pe"
  31     HMM_FOR_PROTEIN_OUTPUT        = "m"
  32     IGNORE_DUF_OPTION             = "i"
  33     PARSE_OUT_DESCRIPITION_OPTION = "a"
  34     HELP_OPTION_1                 = "help"
  35     HELP_OPTION_2                 = "h"
  36
  37     def initialize
  38       @domain_counts = Hash.new
  39     end
  40
  41     # raises ArgumentError, IOError
  42     def parse( inpath,
  43         outpath,
  44         column_delimiter,
  45         i_e_value_threshold,
  46         ignore_dufs,
  47         get_descriptions,
  48         fs_e_value_threshold,
  49         hmm_for_protein_output )
  50       Util.check_file_for_readability( inpath )
  51       Util.check_file_for_writability( outpath )
  52
  53       outfile = File.open( outpath, "a" )
  54
  55       query     = ""
  56       desc      = ""
  57       model     = ""
  58       env_from  = ""
  59       env_to    = ""
  60       i_e_value = ""
  61
  62       hmmscan_results_per_protein = []
  63
  64       hmmscan_parser = HmmscanParser.new( inpath )
  65
  66       prev_query = ""
  67
  68       hmmscan_parser.parse.each do | r |
  69         model     = r.model
  70         query     = r.query
  71         i_e_value = r.i_e_value
  72         env_from  = r.env_from
  73         env_to    = r.env_to
  74
  75         if ( ( i_e_value_threshold < 0.0 ) || ( i_e_value <= i_e_value_threshold ) ) &&
  76            ( !ignore_dufs || ( model !~ /^DUF\d+/ ) )
  77           count_model( model )
  78           outfile.print( query +
  79              column_delimiter )
  80           if ( get_descriptions )
  81             outfile.print( desc +
  82                column_delimiter )
  83           end
  84           outfile.print( model +
  85              column_delimiter +
  86              env_from.to_s +
  87              column_delimiter +
  88              env_to.to_s +
  89              column_delimiter +
  90              i_e_value.to_s )
  91           outfile.print( Constants::LINE_DELIMITER )
  92         end
  93
  94         if  !prev_query.empty? && prev_query != query
  95           if !hmmscan_results_per_protein.empty?
  96             process_hmmscan_results_per_protein( hmmscan_results_per_protein,
  97               fs_e_value_threshold,
  98               hmm_for_protein_output  )
  99           end
 100           hmmscan_results_per_protein.clear
 101         end
 102         prev_query = query
 103         hmmscan_results_per_protein << r
 104       end
 105       if !hmmscan_results_per_protein.empty?
 106         process_hmmscan_results_per_protein( hmmscan_results_per_protein,
 107           fs_e_value_threshold,
 108           hmm_for_protein_output  )
 109       end
 110       outfile.flush()
 111       outfile.close()
 112
 113     end # def parse
 114
 115     def count_model( model )
 116       if ( @domain_counts.has_key?( model ) )
 117         count = @domain_counts[ model ].to_i
 118         count += 1
 119         @domain_counts[ model ] = count
 120       else
 121         @domain_counts[ model ] = 1
 122       end
 123     end
 124
 125     def process_hmmscan_results_per_protein( hmmscan_results_per_protein,
 126         fs_e_value_threshold,
 127         hmm_for_protein_output )
 128
 129       fs_e_value = -1
 130       hmmscan_results_per_protein.each do | r |
 131         if r.model ==  hmm_for_protein_output
 132           fs_e_value = r.fs_e_value
 133           if fs_e_value > fs_e_value_threshold
 134             return
 135           end
 136         end
 137       end
 138
 139
 140       first = hmmscan_results_per_protein[ 0 ]
 141       s = ""
 142       s << first.query + "\t"
 143       s << fs_e_value.to_s + "\t"
 144       s << first.qlen.to_s + "\t"
 145       # s << first.fs_e_value.to_s + "\t"
 146       # s << first.out_of.to_s + "\t"
 147       hmmscan_results_per_protein.each do | r |
 148         s <<  r.model + "|"
 149       end
 150       puts s
 151     end
 152
 153
 154     def get_domain_counts()
 155       return @domain_counts
 156     end
 157
 158     def run()
 159
 160       Util.print_program_information( PRG_NAME,
 161         PRG_VERSION,
 162         PRG_DESC,
 163         PRG_DATE,
 164         COPYRIGHT,
 165         CONTACT,
 166         WWW,
 167         STDOUT )
 168
 169       begin
 170         cla = CommandLineArguments.new( ARGV )
 171       rescue ArgumentError => e
 172         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 173       end
 174
 175       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
 176            cla.is_option_set?( HELP_OPTION_2 ) )
 177         print_help
 178         exit( 0 )
 179       end
 180
 181       if ( cla.get_number_of_files != 2 )
 182         print_help
 183         exit( -1 )
 184       end
 185
 186       allowed_opts = Array.new
 187       allowed_opts.push( DELIMITER_OPTION )
 188       allowed_opts.push( I_E_VALUE_THRESHOLD_OPTION )
 189       allowed_opts.push( FS_E_VALUE_THRESHOLD_OPTION )
 190       allowed_opts.push( IGNORE_DUF_OPTION )
 191       allowed_opts.push( PARSE_OUT_DESCRIPITION_OPTION )
 192       allowed_opts.push( HMM_FOR_PROTEIN_OUTPUT )
 193
 194       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
 195       if ( disallowed.length > 0 )
 196         Util.fatal_error( PRG_NAME,
 197           "unknown option(s): " + disallowed,
 198           STDOUT )
 199       end
 200
 201       inpath = cla.get_file_name( 0 )
 202       outpath = cla.get_file_name( 1 )
 203
 204       column_delimiter = "\t"
 205       if ( cla.is_option_set?( DELIMITER_OPTION ) )
 206         begin
 207           column_delimiter = cla.get_option_value( DELIMITER_OPTION )
 208         rescue ArgumentError => e
 209           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 210         end
 211       end
 212
 213       i_e_value_threshold = -1.0
 214       if ( cla.is_option_set?( I_E_VALUE_THRESHOLD_OPTION ) )
 215         begin
 216           i_e_value_threshold = cla.get_option_value_as_float( I_E_VALUE_THRESHOLD_OPTION )
 217         rescue ArgumentError => e
 218           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 219         end
 220         if ( i_e_value_threshold < 0.0 )
 221           Util.fatal_error( PRG_NAME, "attempt to use a negative i-E-value threshold", STDOUT )
 222         end
 223       end
 224
 225       fs_e_value_threshold = -1.0
 226       if ( cla.is_option_set?( FS_E_VALUE_THRESHOLD_OPTION ) )
 227         begin
 228           fs_e_value_threshold = cla.get_option_value_as_float( FS_E_VALUE_THRESHOLD_OPTION )
 229         rescue ArgumentError => e
 230           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 231         end
 232         if ( fs_e_value_threshold < 0.0 )
 233           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 234         end
 235       end
 236
 237
 238       hmm_for_protein_output = ""
 239       if ( cla.is_option_set?( HMM_FOR_PROTEIN_OUTPUT ) )
 240         begin
 241           hmm_for_protein_output = cla.get_option_value( HMM_FOR_PROTEIN_OUTPUT )
 242         rescue ArgumentError => e
 243           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 244         end
 245       end
 246
 247
 248       ignore_dufs = false
 249       if ( cla.is_option_set?( IGNORE_DUF_OPTION ) )
 250         ignore_dufs = true
 251       end
 252
 253       parse_descriptions = false
 254       if ( cla.is_option_set?( PARSE_OUT_DESCRIPITION_OPTION ) )
 255         parse_descriptions = true
 256       end
 257
 258       puts()
 259       puts( "hmmpfam outputfile  : " + inpath )
 260       puts( "outputfile          : " + outpath )
 261       if ( i_e_value_threshold >= 0.0 )
 262         puts( "i-E-value threshold : " + i_e_value_threshold.to_s )
 263       else
 264         puts( "i-E-value threshold : no threshold" )
 265       end
 266       if ( parse_descriptions )
 267         puts( "parse descriptions  : true" )
 268       else
 269         puts( "parse descriptions  : false" )
 270       end
 271       if ( ignore_dufs )
 272         puts( "ignore DUFs         : true" )
 273       else
 274         puts( "ignore DUFs         : false" )
 275       end
 276       if ( column_delimiter == "\t" )
 277         puts( "column delimiter    : TAB" )
 278       else
 279         puts( "column delimiter     : " + column_delimiter )
 280       end
 281       if ( fs_e_value_threshold >= 0.0 )
 282         puts( "E-value threshold   : " + fs_e_value_threshold.to_s )
 283       else
 284         puts( "E-value threshold   : no threshold" )
 285       end
 286       if ( !hmm_for_protein_output.empty? )
 287         puts( "HMM for proteins    : " + hmm_for_protein_output )
 288       end
 289       puts()
 290
 291       begin
 292         parse( inpath,
 293           outpath,
 294           column_delimiter,
 295           i_e_value_threshold,
 296           ignore_dufs,
 297           parse_descriptions,
 298           fs_e_value_threshold,
 299           hmm_for_protein_output )
 300       rescue ArgumentError, IOError => e
 301         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 302       end
 303       domain_counts = get_domain_counts()
 304
 305
 306       puts
 307       puts( "domain counts (considering potential i-E-value threshold and ignoring of DUFs):" )
 308       puts( "(number of different domains: " + domain_counts.length.to_s + ")" )
 309       puts
 310       puts( Util.draw_histogram( domain_counts, "#" ) )
 311       puts
 312       Util.print_message( PRG_NAME, 'OK' )
 313       puts
 314
 315     end # def run()
 316
 317     def print_help()
 318       puts( "Usage:" )
 319       puts()
 320       puts( "  " + PRG_NAME + ".rb [options] <hmmscan outputfile> <outputfile>" )
 321       puts()
 322       puts( "  options: -" + DELIMITER_OPTION + ": column delimiter for outputfile, default is TAB" )
 323       puts( "           -" + I_E_VALUE_THRESHOLD_OPTION  + ": i-E-value threshold, default is no threshold" )
 324       puts( "           -" + PARSE_OUT_DESCRIPITION_OPTION  + ": parse query description (in addition to query name)" )
 325       puts( "           -" + IGNORE_DUF_OPTION  + ": ignore DUFs" )
 326       puts( "           -" + FS_E_VALUE_THRESHOLD_OPTION  + ": E-value threshold for full protein sequences, only for protein summary" )
 327       puts( "           -" + HMM_FOR_PROTEIN_OUTPUT + ": HMM for protein summary" )
 328       puts()
 329     end
 330
 331   end # class HmmscanParser
 332
 333 end # module Evoruby