forester/ruby/evoruby/lib/evo/apps/domain_sequence_extractor.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     ADD_SPECIES                        = 's'
  34     MIN_LINKER_OPT                        = 'ml'
  35     TRIM_OPTION                        = 't'
  36     LOG_FILE_SUFFIX                    = '_domain_seq_extr.log'
  37     PASSED_SEQS_SUFFIX                 = '_domain_seq_extr_passed'
  38     FAILED_SEQS_SUFFIX                 = '_domain_seq_extr_failed'
  39     HELP_OPTION_1                      = 'help'
  40     HELP_OPTION_2                      = 'h'
  41
  42     def run()
  43
  44       Util.print_program_information( PRG_NAME,
  45         PRG_VERSION,
  46         PRG_DESC ,
  47         PRG_DATE,
  48         COPYRIGHT,
  49         CONTACT,
  50         WWW,
  51         STDOUT )
  52
  53       ld = Constants::LINE_DELIMITER
  54
  55       begin
  56         cla = CommandLineArguments.new( ARGV )
  57       rescue ArgumentError
  58         Util.fatal_error( PRG_NAME, "error: " + $!, STDOUT )
  59       end
  60
  61       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
  62            cla.is_option_set?( HELP_OPTION_2 ) )
  63         print_help
  64         exit( 0 )
  65       end
  66
  67       if ( cla.get_number_of_files != 4 )
  68         print_help
  69         exit( -1 )
  70       end
  71
  72       allowed_opts = Array.new
  73       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
  74       allowed_opts.push( ADD_POSITION_OPTION )
  75       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION )
  76       allowed_opts.push( LENGTH_THRESHOLD_OPTION )
  77       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT )
  78       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION_AS_LETTER )
  79       allowed_opts.push( TRIM_OPTION )
  80       allowed_opts.push( ADD_SPECIES )
  81       allowed_opts.push( MIN_LINKER_OPT )
  82
  83       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
  84       if ( disallowed.length > 0 )
  85         Util.fatal_error( PRG_NAME,
  86           "unknown option(s): " + disallowed,
  87           STDOUT )
  88       end
  89
  90       domain_id           = cla.get_file_name( 0 )
  91       hmmsearch_output    = cla.get_file_name( 1 )
  92       fasta_sequence_file = cla.get_file_name( 2 )
  93       outfile             = cla.get_file_name( 3 )
  94
  95       add_position = false
  96       if ( cla.is_option_set?( ADD_POSITION_OPTION ) )
  97         add_position = true
  98       end
  99
 100       trim = false
 101       if ( cla.is_option_set?( TRIM_OPTION ) )
 102         trim = true
 103       end
 104
 105       add_domain_number           = false
 106       add_domain_number_as_letter = false
 107       add_domain_number_as_digit  = false
 108
 109       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION ) )
 110         add_domain_number = true
 111       end
 112       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION_AS_LETTER ) )
 113         add_domain_number_as_letter = true
 114       end
 115       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT ) )
 116         add_domain_number_as_digit = true
 117       end
 118
 119       add_species = false
 120       if cla.is_option_set? ADD_SPECIES
 121         add_species = true
 122       end
 123
 124       if ( add_domain_number_as_letter && add_domain_number_as_digit )
 125         puts( "attempt to add domain number as letter and digit at the same time" )
 126         print_help
 127         exit( -1 )
 128       end
 129
 130       e_value_threshold = -1.0
 131       if ( cla.is_option_set?( E_VALUE_THRESHOLD_OPTION ) )
 132         begin
 133           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
 134         rescue ArgumentError => e
 135           Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 136         end
 137         if ( e_value_threshold < 0.0 )
 138           Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 139         end
 140       end
 141
 142       length_threshold = -1
 143       if ( cla.is_option_set?( LENGTH_THRESHOLD_OPTION ) )
 144         begin
 145           length_threshold = cla.get_option_value_as_int( LENGTH_THRESHOLD_OPTION )
 146         rescue ArgumentError => e
 147           Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 148         end
 149         if ( length_threshold < 0)
 150           Forester::Util.fatal_error( PRG_NAME, "attempt to use a negative length threshold", STDOUT )
 151         end
 152       end
 153
 154
 155        min_linker = nil
 156        if ( cla.is_option_set?( MIN_LINKER_OPT ) )
 157         begin
 158            min_linker = cla.get_option_value_as_int( MIN_LINKER_OPT )
 159         rescue ArgumentError => e
 160           Forester::Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 161         end
 162          if ( !min_linker || min_linker > 100 || min_linker < -100 )
 163           Forester::Util.fatal_error( PRG_NAME, "unexpected value for min linker " + min_linker.to_s, STDOUT )
 164         end
 165       end
 166
 167
 168       log = String.new
 169
 170       puts()
 171       puts( "Domain                                 : " + domain_id )
 172       log << "Domain                                 : " + domain_id + ld
 173       puts( "Hmmscan outputfile                     : " + hmmsearch_output )
 174       log << "Hmmscan outputfile                     : " + hmmsearch_output + ld
 175       puts( "Fasta sequencefile (complete sequences): " + fasta_sequence_file )
 176       log << "Fasta sequencefile (complete sequences): " + fasta_sequence_file + ld
 177       puts( "Outputfile                             : " + outfile )
 178       log << "Outputfile                             : " + outfile + ld
 179       puts( "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX )
 180       log << "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX + ld
 181       puts( "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX )
 182       log << "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX + ld
 183       puts( "Logfile                                : " + outfile + LOG_FILE_SUFFIX )
 184       log <<  "Logfile                                : " + outfile + LOG_FILE_SUFFIX + ld
 185       if ( e_value_threshold >= 0.0 )
 186         puts( "E-value threshold : " + e_value_threshold.to_s )
 187         log << "E-value threshold : " + e_value_threshold.to_s + ld
 188       else
 189         puts( "E-value threshold : no threshold" )
 190         log << "E-value threshold : no threshold" + ld
 191       end
 192       if ( length_threshold > 0 )
 193         puts( "Length threshold  : " + length_threshold.to_s )
 194         log << "Length threshold  : " + length_threshold.to_s + ld
 195       else
 196         puts( "Length threshold  : no threshold" )
 197         log << "Length threshold  : no threshold" + ld
 198       end
 199
 200       if ( trim )
 201         puts( "Trim last 2 chars : true" )
 202         log << "Trim last 2 chars : true" + ld
 203       else
 204         puts( "Trim names        : false" )
 205         log << "Trim names        : false" + ld
 206       end
 207
 208       if ( min_linker )
 209         puts( "Min linker        : " + min_linker.to_s )
 210         log << "Min linker        :  " + min_linker.to_s +  ld
 211
 212       end
 213
 214
 215       if ( add_position )
 216         puts( "Add positions (rel to complete seq) to extracted domains: true" )
 217         log << "Add positions (rel to complete seq) to extracted domains: true" + ld
 218       else
 219         puts( "Add positions (rel to complete seq) to extracted domains: false" )
 220         log << "Add positions (rel to complete seq) to extracted domains: false" + ld
 221       end
 222
 223       if ( add_domain_number || add_domain_number_as_digit || add_domain_number_as_letter )
 224         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): true" )
 225         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): true" + ld
 226       else
 227         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): false" )
 228         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): false" + ld
 229       end
 230
 231       puts
 232
 233       domain_count = 0
 234       begin
 235         parser = HmmscanDomainExtractor.new()
 236         domain_count = parser.parse( domain_id,
 237           hmmsearch_output,
 238           fasta_sequence_file,
 239           outfile,
 240           outfile + PASSED_SEQS_SUFFIX,
 241           outfile + FAILED_SEQS_SUFFIX,
 242           e_value_threshold,
 243           length_threshold,
 244           add_position,
 245           add_domain_number,
 246           add_domain_number_as_digit,
 247           add_domain_number_as_letter,
 248           trim,
 249           add_species,
 250           min_linker,
 251           log )
 252       rescue ArgumentError, IOError, StandardError => e
 253         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 254       rescue Exception => e
 255         Util.fatal_error( PRG_NAME, "unexpected exception!: " + e.to_s, STDOUT )
 256       end
 257
 258       puts
 259       Util.print_message( PRG_NAME, "extracted a total of " + domain_count.to_s + " domains" )
 260       Util.print_message( PRG_NAME, "wrote;               " + outfile )
 261       Util.print_message( PRG_NAME, "wrote:               " + outfile + LOG_FILE_SUFFIX )
 262       Util.print_message( PRG_NAME, "(wrote:              " + outfile + PASSED_SEQS_SUFFIX + ")" )
 263       Util.print_message( PRG_NAME, "(wrote:              " + outfile + FAILED_SEQS_SUFFIX + ")" )
 264
 265       begin
 266         f = File.open( outfile + LOG_FILE_SUFFIX, 'a' )
 267         f.print( log )
 268         f.close
 269       rescue Exception => e
 270         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 271       end
 272
 273       puts
 274       Util.print_message( PRG_NAME, "OK" )
 275       puts
 276
 277     end
 278
 279     def print_help()
 280       puts()
 281       puts( "Usage:" )
 282       puts()
 283       puts( "  " + PRG_NAME + ".rb [options] <domain> <hmmscan outputfile> <file containing complete sequences in fasta format> <outputfile>" )
 284       puts()
 285       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f>: E-value threshold, default is no threshold" )
 286       puts( "           -" + LENGTH_THRESHOLD_OPTION   + "=<i>: length threshold, default is no threshold" )
 287       puts( "           -" + ADD_POSITION_OPTION  + ": to add positions (rel to complete seq) to extracted domains" )
 288       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\")" )
 289       puts( "           -" + ADD_DOMAIN_NUMBER_OPTION_AS_DIGIT  + ": to add numbers to extracted domains as digit (example \"domain2\")" )
 290       puts( "           -" + ADD_DOMAIN_NUMBER_OPTION_AS_LETTER  + ": to add numbers to extracted domains as letter (example \"domaina\")" )
 291       puts( "           -" + TRIM_OPTION  + ": to remove the last 2 characters from sequence names" )
 292       puts( "           -" + ADD_SPECIES  + ": to add species [in brackets]" )
 293       puts( "           -" + MIN_LINKER_OPT  + "=<i>: to extract pairs of same domains with a distance inbetween shorter than a given value" )
 294       puts()
 295     end
 296
 297   end # class DomainSequenceExtractor
 298
 299 end # module Evoruby