forester/ruby/evoruby/lib/evo/apps/domain_sequence_extractor_new.rb

   1 #
   2 # = lib/evo/apps/domain_sequence_extractor.rb - DomainSequenceExtractor class
   3 #
   4 # Copyright::  Copyright (C) 2012 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id:Exp $
   8
   9
  10 require 'lib/evo/util/constants'
  11 require 'lib/evo/util/util'
  12 require 'lib/evo/util/command_line_arguments'
  13 require 'lib/evo/io/parser/hmmscan_domain_extractor'
  14
  15 module Evoruby
  16
  17   class DomainSequenceExtractor
  18
  19     PRG_NAME       = "dsx"
  20     PRG_VERSION    = "1.000"
  21     PRG_DESC       = "extraction of domain sequences from hmmscan output"
  22     PRG_DATE       = "20120906"
  23     COPYRIGHT      = "2012 Christian M Zmasek"
  24     CONTACT        = "phylosoft@gmail.com"
  25     WWW            = "www.phylosoft.org"
  26
  27     E_VALUE_THRESHOLD_OPTION           = 'e'
  28     LENGTH_THRESHOLD_OPTION            = 'l'
  29     ADD_POSITION_OPTION                = 'p'
  30     ADD_DOMAIN_NUMBER_OPTION           = 'd'
  31     ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT  = 'dd'
  32     ADD_DOMAIN_NUMBER_OPTION_AS_LETTER = 'dl'
  33     TRIM_OPTION                        = 't'
  34     LOG_FILE_SUFFIX                    = '_domain_seq_extr.log'
  35     PASSED_SEQS_SUFFIX                 = '_domain_seq_extr_passed'
  36     FAILED_SEQS_SUFFIX                 = '_domain_seq_extr_failed'
  37     HELP_OPTION_1                      = 'help'
  38     HELP_OPTION_2                      = 'h'
  39
  40     def run()
  41
  42       Util.print_program_information( PRG_NAME,
  43         PRG_VERSION,
  44         PRG_DESC ,
  45         PRG_DATE,
  46         COPYRIGHT,
  47         CONTACT,
  48         WWW,
  49         STDOUT )
  50
  51       ld = Constants::LINE_DELIMITER
  52
  53       begin
  54         cla = CommandLineArguments.new( ARGV )
  55       rescue ArgumentError
  56         Util.fatal_error( PRG_NAME, "error: " + $!, STDOUT )
  57       end
  58
  59       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
  60            cla.is_option_set?( HELP_OPTION_2 ) )
  61         print_help
  62         exit( 0 )
  63       end
  64
  65       if ( cla.get_number_of_files != 4 )
  66         print_help
  67         exit( -1 )
  68       end
  69
  70       allowed_opts = Array.new
  71       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
  72       allowed_opts.push( ADD_POSITION_OPTION )
  73       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION )
  74       allowed_opts.push( LENGTH_THRESHOLD_OPTION )
  75       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT )
  76       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION_AS_LETTER )
  77       allowed_opts.push( TRIM_OPTION )
  78
  79       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
  80       if ( disallowed.length > 0 )
  81         Util.fatal_error( PRG_NAME,
  82           "unknown option(s): " + disallowed,
  83           STDOUT )
  84       end
  85
  86       domain_id           = cla.get_file_name( 0 )
  87       hmmsearch_output    = cla.get_file_name( 1 )
  88       fasta_sequence_file = cla.get_file_name( 2 )
  89       outfile             = cla.get_file_name( 3 )
  90
  91       add_position = false
  92       if ( cla.is_option_set?( ADD_POSITION_OPTION ) )
  93         add_position = true
  94       end
  95
  96       trim = false
  97       if ( cla.is_option_set?( TRIM_OPTION ) )
  98         trim = true
  99       end
 100
 101       add_domain_number           = false
 102       add_domain_number_as_letter = false
 103       add_domain_number_as_digit  = false
 104
 105       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION ) )
 106         add_domain_number = true
 107       end
 108       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION_AS_LETTER ) )
 109         add_domain_number_as_letter = true
 110       end
 111       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT ) )
 112         add_domain_number_as_digit = true
 113       end
 114
 115       if ( add_domain_number_as_letter && add_domain_number_as_digit )
 116         puts( "attempt to add domain number as letter and digit at the same time" )
 117         print_help
 118         exit( -1 )
 119       end
 120
 121       e_value_threshold = -1.0
 122       if ( cla.is_option_set?( E_VALUE_THRESHOLD_OPTION ) )
 123         begin
 124           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
 125         rescue ArgumentError => e
 126           Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 127         end
 128         if ( e_value_threshold < 0.0 )
 129           Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 130         end
 131       end
 132
 133       length_threshold = -1
 134       if ( cla.is_option_set?( LENGTH_THRESHOLD_OPTION ) )
 135         begin
 136           length_threshold = cla.get_option_value_as_int( LENGTH_THRESHOLD_OPTION )
 137         rescue ArgumentError => e
 138           Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 139         end
 140         if ( length_threshold < 0)
 141           Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative length threshold", STDOUT )
 142         end
 143       end
 144
 145       log = String.new
 146
 147       puts()
 148       puts( "Domain                                 : " + domain_id )
 149       log << "Domain                                 : " + domain_id + ld
 150       puts( "Hmmscan outputfile                     : " + hmmsearch_output )
 151       log << "Hmmscan outputfile                     : " + hmmsearch_output + ld
 152       puts( "Fasta sequencefile (complete sequences): " + fasta_sequence_file )
 153       log << "Fasta sequencefile (complete sequences): " + fasta_sequence_file + ld
 154       puts( "Outputfile                             : " + outfile )
 155       log << "Outputfile                             : " + outfile + ld
 156       puts( "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX )
 157       log << "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX + ld
 158       puts( "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX )
 159       log << "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX + ld
 160       puts( "Logfile                                : " + outfile + LOG_FILE_SUFFIX )
 161       log <<  "Logfile                                : " + outfile + LOG_FILE_SUFFIX + ld
 162       if ( e_value_threshold >= 0.0 )
 163         puts( "E-value threshold : " + e_value_threshold.to_s )
 164         log << "E-value threshold : " + e_value_threshold.to_s + ld
 165       else
 166         puts( "E-value threshold : no threshold" )
 167         log << "E-value threshold : no threshold" + ld
 168       end
 169       if ( length_threshold > 0 )
 170         puts( "Length threshold  : " + length_threshold.to_s )
 171         log << "Length threshold  : " + length_threshold.to_s + ld
 172       else
 173         puts( "Length threshold  : no threshold" )
 174         log << "Length threshold  : no threshold" + ld
 175       end
 176
 177       if ( trim )
 178         puts( "Trim last 2 chars : true" )
 179         log << "Trim last 2 chars : true" + ld
 180       else
 181         puts( "Trim names        : false" )
 182         log << "Trim names        : false" + ld
 183       end
 184
 185
 186       if ( add_position )
 187         puts( "Add positions (rel to complete seq) to extracted domains: true" )
 188         log << "Add positions (rel to complete seq) to extracted domains: true" + ld
 189       else
 190         puts( "Add positions (rel to complete seq) to extracted domains: false" )
 191         log << "Add positions (rel to complete seq) to extracted domains: false" + ld
 192       end
 193
 194       if ( add_domain_number || add_domain_number_as_digit || add_domain_number_as_letter )
 195         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): true" )
 196         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): true" + ld
 197       else
 198         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): false" )
 199         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): false" + ld
 200       end
 201
 202       puts
 203
 204       domain_count = 0
 205       begin
 206         parser = HmmscanDomainExtractor.new()
 207         domain_count = parser.parse( domain_id,
 208           hmmsearch_output,
 209           fasta_sequence_file,
 210           outfile,
 211           outfile + PASSED_SEQS_SUFFIX,
 212           outfile + FAILED_SEQS_SUFFIX,
 213           e_value_threshold,
 214           length_threshold,
 215           add_position,
 216           add_domain_number,
 217           add_domain_number_as_digit,
 218           add_domain_number_as_letter,
 219           trim,
 220           log )
 221       rescue ArgumentError, IOError, StandardError => e
 222         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 223       rescue Exception => e
 224         Util.fatal_error( PRG_NAME, "unexpected exception!: " + e.to_s, STDOUT )
 225       end
 226
 227       puts
 228       Util.print_message( PRG_NAME, "extracted a total of " + domain_count.to_s + " domains" )
 229       Util.print_message( PRG_NAME, "wrote;               " + outfile )
 230       Util.print_message( PRG_NAME, "wrote:               " + outfile + LOG_FILE_SUFFIX )
 231       Util.print_message( PRG_NAME, "(wrote:              " + outfile + PASSED_SEQS_SUFFIX + ")" )
 232       Util.print_message( PRG_NAME, "(wrote:              " + outfile + FAILED_SEQS_SUFFIX + ")" )
 233
 234       begin
 235         f = File.open( outfile + LOG_FILE_SUFFIX, 'a' )
 236         f.print( log )
 237         f.close
 238       rescue Exception => e
 239         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 240       end
 241
 242       puts
 243       Util.print_message( PRG_NAME, "OK" )
 244       puts
 245
 246     end
 247
 248     def print_help()
 249       puts()
 250       puts( "Usage:" )
 251       puts()
 252       puts( "  " + PRG_NAME + ".rb [options] <domain> <hmmscan outputfile> <file containing complete sequences in fasta format> <outputfile>" )
 253       puts()
 254       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f>: E-value threshold, default is no threshold" )
 255       puts( "           -" + LENGTH_THRESHOLD_OPTION   + "=<i>: length threshold, default is no threshold" )
 256       puts( "           -" + ADD_POSITION_OPTION  + ": to add positions (rel to complete seq) to extracted domains" )
 257       puts( "           -" + ADD_DOMAIN_NUMBER_OPTION  + ": to add numbers to extracted domains (in case of more than one domain per complete seq) (example \"domain~2-3\")" )
 258       puts( "           -" + ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT  + ": to add numbers to extracted domains as digit (example \"domain2\")" )
 259       puts( "           -" + ADD_DOMAIN_NUMBER_OPTION_AS_LETTER  + ": to add numbers to extracted domains as letter (example \"domaina\")" )
 260       puts( "           -" + TRIM_OPTION  + ": to remove the last 2 characters from sequence names" )
 261       puts()
 262     end
 263
 264   end # class DomainSequenceExtractor
 265
 266 end # module Evoruby