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

   1 #
   2 # = lib/evo/apps/taxonomy_processor - TaxonomyProcessor class
   3 #
   4 # Copyright::    Copyright (C) 2017 Christian M. Zmasek
   5 # License::      GNU Lesser General Public License (LGPL)
   6
   7 require 'lib/evo/util/constants'
   8 require 'lib/evo/util/util'
   9 require 'lib/evo/util/command_line_arguments'
  10 require 'lib/evo/io/parser/hmmscan_domain_extractor'
  11
  12 module Evoruby
  13   class DomainSequenceExtractor
  14
  15     PRG_NAME       = "dsx"
  16     PRG_VERSION    = "2.002"
  17     PRG_DESC       = "Extraction of domain sequences from hmmscan output"
  18     PRG_DATE       = "20170214"
  19     WWW            = "https://sites.google.com/site/cmzmasek/home/software/forester"
  20
  21     E_VALUE_THRESHOLD_OPTION           = 'e'
  22     LENGTH_THRESHOLD_OPTION            = 'l'
  23     ADD_POSITION_OPTION                = 'p'
  24     ADD_DOMAIN_NUMBER_OPTION           = 'd'
  25     ADD_SPECIES                        = 's'
  26     LOG_FILE_SUFFIX                    = '_domain_seq_extr.log'
  27     PASSED_SEQS_SUFFIX                 = '_with_passing_domains.fasta'
  28     FAILED_SEQS_SUFFIX                 = '_with_no_passing_domains.fasta'
  29     HELP_OPTION_1                      = 'help'
  30     HELP_OPTION_2                      = 'h'
  31     def run()
  32
  33       Util.print_program_information( PRG_NAME,
  34       PRG_VERSION,
  35       PRG_DESC ,
  36       PRG_DATE,
  37       WWW,
  38       STDOUT )
  39
  40       if ( ARGV == nil || ( ARGV.length < 1 )  )
  41         print_help
  42         exit( -1 )
  43       end
  44
  45       begin
  46         cla = CommandLineArguments.new( ARGV )
  47       rescue ArgumentError
  48         Util.fatal_error( PRG_NAME, "error: " + $!, STDOUT )
  49       end
  50
  51       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
  52       cla.is_option_set?( HELP_OPTION_2 ) )
  53         print_help
  54         exit( 0 )
  55       end
  56
  57       unless ( cla.get_number_of_files == 2 || cla.get_number_of_files == 3 || cla.get_number_of_files == 4 )
  58         print_help
  59         exit( -1 )
  60       end
  61
  62       allowed_opts = Array.new
  63       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
  64       allowed_opts.push( ADD_POSITION_OPTION )
  65       allowed_opts.push( ADD_DOMAIN_NUMBER_OPTION )
  66       allowed_opts.push( LENGTH_THRESHOLD_OPTION )
  67       allowed_opts.push( ADD_SPECIES )
  68
  69       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
  70       if ( disallowed.length > 0 )
  71         Util.fatal_error( PRG_NAME,
  72         "unknown option(s): " + disallowed,
  73         STDOUT )
  74       end
  75
  76       e_value_threshold = 0.1
  77       if ( cla.is_option_set?( E_VALUE_THRESHOLD_OPTION ) )
  78         begin
  79           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
  80         rescue ArgumentError => e
  81           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
  82         end
  83         if ( e_value_threshold < 0.0 )
  84           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
  85         end
  86       end
  87
  88       length_threshold = 50
  89       if ( cla.is_option_set?( LENGTH_THRESHOLD_OPTION ) )
  90         begin
  91           length_threshold = cla.get_option_value_as_int( LENGTH_THRESHOLD_OPTION )
  92         rescue ArgumentError => e
  93           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
  94         end
  95         if ( length_threshold < 0)
  96           Util.fatal_error( PRG_NAME, "attempt to use a negative length threshold", STDOUT )
  97         end
  98       end
  99
 100       domain_id           = cla.get_file_name( 0 )
 101       hmmscan_output      = cla.get_file_name( 1 )
 102       fasta_sequence_file = ""
 103       outfile             = ""
 104
 105       if (cla.get_number_of_files == 4)
 106         fasta_sequence_file = cla.get_file_name( 2 )
 107         outfile             = cla.get_file_name( 3 )
 108       elsif (cla.get_number_of_files == 2 || cla.get_number_of_files == 3)
 109         if (cla.get_number_of_files == 3)
 110           fasta_sequence_file = cla.get_file_name( 2 )
 111         else
 112           hmmscan_index = hmmscan_output.index(Constants::HMMSCAN)
 113           if ( hmmscan_index != nil )
 114             prefix = hmmscan_output[0 .. hmmscan_index-1 ]
 115             suffix = Constants::ID_NORMALIZED_FASTA_FILE_SUFFIX
 116             files = Dir.entries( "." )
 117             matching_files = Util.get_matching_files( files, prefix, suffix)
 118             if matching_files.length < 1
 119               Util.fatal_error( PRG_NAME, 'no file matching [' + prefix +
 120               '...' + suffix + '] present in current directory: need to indicate <file containing complete sequences in fasta format> as second argument' )
 121             end
 122             if matching_files.length > 1
 123               Util.fatal_error( PRG_NAME, 'more than one file matching [' +
 124               prefix  + '...' + suffix + '] present in current directory: need to indicate <file containing complete sequences in fasta format> as second argument' )
 125             end
 126             fasta_sequence_file = matching_files[ 0 ]
 127           else
 128             Util.fatal_error( PRG_NAME, 'input files do not seem in format for standard analysis pipeline, need to explicitly indicate all' )
 129           end
 130         end
 131         hmmscan_index = hmmscan_output.index(Constants::HMMSCAN)
 132         if ( hmmscan_index != nil )
 133           outfile = hmmscan_output.sub(Constants::HMMSCAN, "_") + "_" + domain_id
 134           e = e_value_threshold >= 1 ? e_value_threshold.to_i : e_value_threshold.to_s.sub!('.','_')
 135           outfile += "_E" + e.to_s
 136           outfile += "_L" + length_threshold.to_i.to_s
 137         else
 138           Util.fatal_error( PRG_NAME, 'input files do not seem in format for standard analysis pipeline, need to explicitly indicate all' )
 139         end
 140       end
 141
 142       if outfile.downcase.end_with?( ".fasta" )
 143         outfile = outfile[ 0 .. outfile.length - 7 ]
 144       elsif outfile.downcase.end_with?( ".fsa" )
 145         outfile = outfile[ 0 .. outfile.length - 5 ]
 146       end
 147
 148       add_position = false
 149       if ( cla.is_option_set?( ADD_POSITION_OPTION ) )
 150         add_position = true
 151       end
 152
 153       add_domain_number = false
 154       if ( cla.is_option_set?( ADD_DOMAIN_NUMBER_OPTION ) )
 155         add_domain_number = true
 156       end
 157
 158       add_species = false
 159       if cla.is_option_set?( ADD_SPECIES )
 160         add_species = true
 161       end
 162
 163       log = String.new
 164       ld = Constants::LINE_DELIMITER
 165
 166       puts()
 167       puts( "Domain                                 : " + domain_id )
 168       log << "Domain                                 : " + domain_id + ld
 169       puts( "Hmmscan outputfile                     : " + hmmscan_output )
 170       log << "Hmmscan outputfile                     : " + hmmscan_output + ld
 171       puts( "Fasta sequencefile (complete sequences): " + fasta_sequence_file )
 172       log << "Fasta sequencefile (complete sequences): " + fasta_sequence_file + ld
 173       puts( "Outputfile                             : " + outfile + ".fasta" )
 174       log << "Outputfile                             : " + outfile + ld
 175       puts( "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX )
 176       log << "Passed sequences outfile (fasta)       : " + outfile + PASSED_SEQS_SUFFIX + ld
 177       puts( "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX )
 178       log << "Failed sequences outfile (fasta)       : " + outfile + FAILED_SEQS_SUFFIX + ld
 179       puts( "Logfile                                : " + outfile + LOG_FILE_SUFFIX )
 180       log <<  "Logfile                                : " + outfile + LOG_FILE_SUFFIX + ld
 181       if ( e_value_threshold >= 0.0 )
 182         puts( "E-value threshold : " + e_value_threshold.to_s )
 183         log << "E-value threshold : " + e_value_threshold.to_s + ld
 184       else
 185         puts( "E-value threshold : no threshold" )
 186         log << "E-value threshold : no threshold" + ld
 187       end
 188       if ( length_threshold > 0 )
 189         puts( "Length threshold  : " + length_threshold.to_s )
 190         log << "Length threshold  : " + length_threshold.to_s + ld
 191       else
 192         puts( "Length threshold  : no threshold" )
 193         log << "Length threshold  : no threshold" + ld
 194       end
 195
 196       if ( add_position )
 197         puts( "Add positions (rel to complete seq) to extracted domains: true" )
 198         log << "Add positions (rel to complete seq) to extracted domains: true" + ld
 199       else
 200         puts( "Add positions (rel to complete seq) to extracted domains: false" )
 201         log << "Add positions (rel to complete seq) to extracted domains: false" + ld
 202       end
 203
 204       if ( add_domain_number )
 205         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): true" )
 206         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): true" + ld
 207       else
 208         puts( "Add numbers to extracted domains (in case of more than one domain per complete seq): false" )
 209         log << "Add numbers to extracted domains (in case of more than one domain per complete seq): false" + ld
 210       end
 211
 212       puts
 213
 214       domain_count = 0
 215       begin
 216         parser = HmmscanDomainExtractor.new()
 217         domain_count = parser.parse( domain_id,
 218         hmmscan_output,
 219         fasta_sequence_file,
 220         outfile,
 221         outfile + PASSED_SEQS_SUFFIX,
 222         outfile + FAILED_SEQS_SUFFIX,
 223         e_value_threshold,
 224         length_threshold,
 225         add_position,
 226         add_domain_number,
 227         add_species,
 228         false,
 229         log )
 230       rescue ArgumentError, IOError => e
 231         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 232
 233       rescue Exception => e
 234         puts e.backtrace
 235         Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
 236
 237       end
 238
 239       puts
 240       Util.print_message( PRG_NAME, "extracted a total of " + domain_count.to_s + " domains" )
 241       Util.print_message( PRG_NAME, "wrote;               " + outfile + ".fasta")
 242       Util.print_message( PRG_NAME, "wrote:               " + outfile + LOG_FILE_SUFFIX )
 243       Util.print_message( PRG_NAME, "wrote:               " + outfile + PASSED_SEQS_SUFFIX )
 244       Util.print_message( PRG_NAME, "wrote:               " + outfile + FAILED_SEQS_SUFFIX )
 245
 246       begin
 247         f = File.open( outfile + LOG_FILE_SUFFIX, 'a' )
 248         f.print( log )
 249         f.close
 250       rescue Exception => e
 251         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 252       end
 253
 254       puts
 255       Util.print_message( PRG_NAME, "OK" )
 256       puts
 257
 258     end
 259
 260     def print_help()
 261       puts()
 262       puts( "Usage:" )
 263       puts()
 264       puts( "  " + PRG_NAME + ".rb [options] <domain> <hmmscan outputfile> [file containing complete sequences in fasta format] [outputfile]" )
 265       puts()
 266       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f> : iE-value threshold, default is 0.1" )
 267       puts( "           -" + LENGTH_THRESHOLD_OPTION   + "=<i> : length threshold, default is 50" )
 268       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\")" )
 269       puts( "           -" + ADD_POSITION_OPTION  + "     : to add positions (rel to complete seq) to extracted domains" )
 270       puts( "           -" + ADD_SPECIES  + "     : to add species [in brackets]" )
 271       puts()
 272       puts( "Examples:" )
 273       puts
 274       puts( "  " + PRG_NAME + ".rb -d -e=1e-6 -l=50 Pkinase P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10 P53_ni.fasta P53_#{Constants::PFAM_V_FOR_EX}_10_Pkinase_E1_0e-06_L50" )
 275       puts
 276       puts( "  " + PRG_NAME + ".rb -d -e=1e-6 -l=50 Pkinase P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10 P53_ni.fasta" )
 277       puts
 278       puts( "  " + PRG_NAME + ".rb -d -e=1e-6 -l=50 Pkinase P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10" )
 279       puts()
 280     end
 281
 282   end # class DomainSequenceExtractor
 283
 284 end # module Evoruby