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