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

   1 #
   2 # = lib/evo/apps/domains_to_forester - DomainsToForester 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/msa/msa_factory'
  11 require 'lib/evo/io/parser/fasta_parser'
  12 require 'lib/evo/sequence/protein_domain'
  13 require 'lib/evo/sequence/domain_structure'
  14
  15 module Evoruby
  16   class DomainsToForester
  17
  18     PRG_NAME       = "d2f"
  19     PRG_DESC       = "converting of parsed hmmpfam output to forester format"
  20     PRG_VERSION    = "1.002"
  21     PRG_DATE       = "20170213"
  22     WWW            = "https://sites.google.com/site/cmzmasek/home/software/forester"
  23
  24     E_VALUE_THRESHOLD_OPTION         = "e"
  25     OVERWRITE_IF_SAME_FROM_TO_OPTION = "o"
  26     HELP_OPTION_1                    = "help"
  27     HELP_OPTION_2                    = "h"
  28     def parse( domains_list_file,
  29       original_seqs_file,
  30       outfile,
  31       column_delimiter,
  32       e_value_threshold,
  33       overwrite_if_same_from_to )
  34       Util.check_file_for_readability( domains_list_file )
  35       Util.check_file_for_readability( original_seqs_file )
  36       Util.check_file_for_writability( outfile )
  37
  38       domain_structures = Hash.new() # protein name is key, domain structure is value
  39
  40       f = MsaFactory.new
  41
  42       original_seqs = f.create_msa_from_file( original_seqs_file, FastaParser.new )
  43       if ( original_seqs.get_number_of_seqs < 1 )
  44         error_msg = "\"" + original_seqs_file + "\" appears devoid of sequences in fasta-format"
  45         raise ArgumentError, error_msg
  46       end
  47
  48       File.open( domains_list_file ) do | file |
  49         while line = file.gets
  50           if !is_ignorable?( line )
  51
  52             a = line.split( column_delimiter )
  53             l = a.length
  54             if ( ( l < 4 ) || ( e_value_threshold >= 0.0 && l < 5 ) )
  55               error_msg = "unexpected format at line: " + line
  56               raise IOError, error_msg
  57             end
  58             protein_name = a[ 0 ]
  59             domain_name  = a[ 1 ]
  60             seq_from     = -1
  61             seq_to       = -1
  62
  63             begin
  64               seq_from = a[ 2 ].to_i
  65             rescue Exception
  66               error_msg = "failed to parse seq from from \"" + a[ 2 ] + "\" [line: " + line + "]"
  67               raise IOError, error_msg
  68             end
  69             begin
  70               seq_to = a[ 3 ].to_i
  71             rescue Exception
  72               error_msg = "failed to parse seq to from \"" + a[ 3 ] + "\" [line: " + line + "]"
  73               raise IOError, error_msg
  74             end
  75
  76             e_value = -1
  77             if l > 4
  78               begin
  79                 e_value = a[ 4 ].to_f
  80               rescue Exception
  81                 error_msg = "failed to parse E-value from \"" + a[ 4 ] + "\" [line: " + line + "]"
  82                 raise IOError, error_msg
  83               end
  84             end
  85
  86             seq = original_seqs.get_by_name_start( protein_name )
  87
  88             total_length = seq.get_length
  89
  90             if ( ( ( e_value_threshold < 0.0 ) || ( e_value <= e_value_threshold ) )  )
  91               pd = ProteinDomain.new( domain_name, seq_from, seq_to, "", e_value )
  92               ds = nil
  93               if ( domain_structures.has_key?( protein_name ) )
  94                 ds = domain_structures[ protein_name ]
  95               else
  96                 ds = DomainStructure.new( total_length )
  97                 domain_structures[ protein_name ] = ds
  98               end
  99               ds.add_domain( pd, overwrite_if_same_from_to )
 100             end
 101
 102           end
 103         end
 104       end
 105
 106       out = File.open( outfile, "a" )
 107       ds = domain_structures.sort
 108       for d in ds
 109         protein_name     = d[ 0 ]
 110         domain_structure = d[ 1 ]
 111         out.print( protein_name.to_s )
 112         out.print( "\t" )
 113         out.print( domain_structure.to_NHX )
 114         out.print( Constants::LINE_DELIMITER  )
 115       end
 116
 117       out.flush()
 118       out.close()
 119
 120     end # parse
 121
 122     def run()
 123
 124       Util.print_program_information( PRG_NAME,
 125       PRG_VERSION,
 126       PRG_DESC,
 127       PRG_DATE,
 128       WWW,
 129       STDOUT )
 130
 131       if ( ARGV == nil || ( ARGV.length < 1 )  )
 132         print_help
 133         exit( -1 )
 134       end
 135
 136       begin
 137         cla = CommandLineArguments.new( ARGV )
 138       rescue ArgumentError => e
 139         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 140       end
 141
 142       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
 143       cla.is_option_set?( HELP_OPTION_2 ) )
 144         print_help
 145         exit( 0 )
 146       end
 147
 148       unless ( cla.get_number_of_files == 1 || cla.get_number_of_files == 2 || cla.get_number_of_files == 3 )
 149         print_help
 150         exit( -1 )
 151       end
 152
 153       allowed_opts = Array.new
 154       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
 155       allowed_opts.push( OVERWRITE_IF_SAME_FROM_TO_OPTION )
 156
 157       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
 158       if ( disallowed.length > 0 )
 159         Util.fatal_error( PRG_NAME,
 160         "unknown option(s): " + disallowed,
 161         STDOUT )
 162       end
 163
 164       e_value_threshold = -1.0
 165       if cla.is_option_set?( E_VALUE_THRESHOLD_OPTION )
 166         begin
 167           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
 168         rescue ArgumentError => e
 169           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 170         end
 171         if ( e_value_threshold < 0.0 )
 172           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 173         end
 174       end
 175
 176       domains_list_file = cla.get_file_name( 0 )
 177       original_sequences_file = ""
 178       outfile = ""
 179       if (cla.get_number_of_files == 3)
 180         original_sequences_file = cla.get_file_name( 1 )
 181         outfile = cla.get_file_name( 2 )
 182       elsif (cla.get_number_of_files == 1 || cla.get_number_of_files == 2 )
 183         if ( cla.get_number_of_files == 2 )
 184           original_sequences_file = cla.get_file_name( 1 )
 185         else
 186           hmmscan_index = domains_list_file.index(Constants::HMMSCAN)
 187           if ( hmmscan_index != nil )
 188             prefix = domains_list_file[0 .. hmmscan_index-1 ]
 189             suffix = Constants::ID_NORMALIZED_FASTA_FILE_SUFFIX
 190             files = Dir.entries( "." )
 191             matching_files = Util.get_matching_files( files, prefix, suffix)
 192             if matching_files.length < 1
 193               Util.fatal_error( PRG_NAME, 'no file matching [' + prefix +
 194               '...' + suffix + '] present in current directory: need to indicate <file containing complete sequences in fasta format> as second argument' )
 195             end
 196             if matching_files.length > 1
 197               Util.fatal_error( PRG_NAME, 'more than one file matching [' +
 198               prefix  + '...' + suffix + '] present in current directory: need to indicate <file containing complete sequences in fasta format> as second argument' )
 199             end
 200             original_sequences_file = matching_files[ 0 ]
 201           else
 202             Util.fatal_error( PRG_NAME, 'input files do not seem in format for standard analysis pipeline, need to explicitly indicate all' )
 203           end
 204
 205         end
 206         outfile = domains_list_file
 207         if (outfile.end_with?(Constants::DOMAIN_TABLE_SUFFIX) )
 208           outfile = outfile.chomp(Constants::DOMAIN_TABLE_SUFFIX)
 209         end
 210         if ( e_value_threshold >= 0.0 )
 211           e = e_value_threshold >= 1 ? e_value_threshold.to_i : e_value_threshold.to_s.sub!('.','_')
 212           outfile = outfile + Constants::DOMAINS_TO_FORESTER_EVALUE_CUTOFF_SUFFIX + e.to_s
 213         end
 214         outfile = outfile + Constants::DOMAINS_TO_FORESTER_OUTFILE_SUFFIX
 215       end
 216
 217       overwrite_if_same_from_to = false
 218       if ( cla.is_option_set?( OVERWRITE_IF_SAME_FROM_TO_OPTION ) )
 219         overwrite_if_same_from_to = true
 220       end
 221
 222       puts
 223       puts( "Domain table                            : " + domains_list_file )
 224       puts( "Fasta sequence file (complete sequences): " + original_sequences_file )
 225       puts( "Outputfile                              : " + outfile )
 226       if ( e_value_threshold >= 0.0 )
 227         puts( "E-value threshold                       : " + e_value_threshold.to_s )
 228       else
 229         puts( "E-value threshold                       : no threshold" )
 230       end
 231       if ( overwrite_if_same_from_to )
 232         puts( "Overwrite if same from and to           : true" )
 233       else
 234         puts( "Overwrite if same from and to           : false" )
 235       end
 236
 237       puts
 238
 239       begin
 240         parse( domains_list_file,
 241         original_sequences_file,
 242         outfile,
 243         " ",
 244         e_value_threshold,
 245         overwrite_if_same_from_to )
 246
 247       rescue ArgumentError, IOError, StandardError => e
 248         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 249       rescue Exception => e
 250         Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
 251       end
 252
 253       puts
 254       Util.print_message( PRG_NAME, "wrote: " + outfile )
 255       Util.print_message( PRG_NAME, "next step in standard analysis pipeline: dsx.rb")
 256       Util.print_message( PRG_NAME, 'OK' )
 257       puts
 258
 259     end
 260
 261     private
 262
 263     def print_help()
 264       puts
 265       puts( "Usage:" )
 266       puts
 267       puts( "  " + PRG_NAME + ".rb [options] <domain table (parsed hmmpfam output)> [file containing complete sequences in fasta format] [outputfile]" )
 268       puts()
 269       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f>: E-value threshold, default is no threshold" )
 270       puts( "           -" + OVERWRITE_IF_SAME_FROM_TO_OPTION  + "    : overwrite domain with same start and end with domain with better E-value" )
 271       puts
 272       puts( "  [next step in standard analysis pipeline: dsx.rb]")
 273       puts()
 274       puts( "Examples:" )
 275       puts
 276       puts( "  " + PRG_NAME + ".rb -o P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10_domain_table P53_ni.fasta P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10.dff" )
 277       puts
 278       puts( "  " + PRG_NAME + ".rb -o P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10_domain_table P53_ni.fasta" )
 279       puts
 280       puts( "  " + PRG_NAME + ".rb -o P53_hmmscan_#{Constants::PFAM_V_FOR_EX}_10_domain_table" )
 281       puts()
 282     end
 283
 284     def is_ignorable?( line )
 285       return ( line !~ /[A-Za-z0-9-]/ || line =~ /^\s*#/)
 286     end
 287
 288   end # class DomainsToForester
 289
 290 end # module Evoruby