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) 2006-2007 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: Exp $
   8 #
   9 # last modified: 06/11/2007
  10
  11 require 'lib/evo/util/constants'
  12 require 'lib/evo/util/util'
  13 require 'lib/evo/util/command_line_arguments'
  14 require 'lib/evo/msa/msa_factory'
  15 require 'lib/evo/io/parser/fasta_parser'
  16 require 'lib/evo/sequence/protein_domain'
  17 require 'lib/evo/sequence/domain_structure'
  18
  19 module Evoruby
  20
  21   class DomainsToForester
  22
  23     PRG_NAME       = "d2f"
  24     PRG_DESC       = "parsed hmmpfam output to forester format"
  25     PRG_VERSION    = "1.001"
  26     PRG_DATE       = "20120807"
  27     COPYRIGHT      = "2012 Christian M Zmasek"
  28     CONTACT        = "phylosoft@gmail.com"
  29     WWW            = "www.phylosoft.org"
  30
  31     E_VALUE_THRESHOLD_OPTION         = "e"
  32     OVERWRITE_IF_SAME_FROM_TO_OPTION = "o"
  33     HELP_OPTION_1                    = "help"
  34     HELP_OPTION_2                    = "h"
  35
  36     def parse( domains_list_file,
  37         original_seqs_file,
  38         outfile,
  39         column_delimiter,
  40         e_value_threshold,
  41         overwrite_if_same_from_to )
  42       Util.check_file_for_readability( domains_list_file )
  43       Util.check_file_for_readability( original_seqs_file )
  44       Util.check_file_for_writability( outfile )
  45
  46       domain_structures = Hash.new() # protein name is key, domain structure is value
  47
  48       f = MsaFactory.new
  49
  50       original_seqs = f.create_msa_from_file( original_seqs_file, FastaParser.new )
  51       if ( original_seqs.get_number_of_seqs < 1 )
  52         error_msg = "\"" + original_seqs_file + "\" appears devoid of sequences in fasta-format"
  53         raise ArgumentError, error_msg
  54       end
  55
  56       File.open( domains_list_file ) do | file |
  57         while line = file.gets
  58           if !is_ignorable?( line )
  59
  60             a = line.split( column_delimiter )
  61             l = a.length
  62             if ( ( l < 4 ) || ( e_value_threshold >= 0.0 && l < 5 ) )
  63               error_msg = "unexpected format at line: " + line
  64               raise IOError, error_msg
  65             end
  66             protein_name = a[ 0 ]
  67             domain_name  = a[ 1 ]
  68             seq_from     = -1
  69             seq_to       = -1
  70             ##########################################
  71             if domain_name =~ /RRM_\d/
  72               puts "ignoring " + line
  73               next
  74             end
  75             ##########################################
  76
  77
  78             begin
  79               seq_from = a[ 2 ].to_i
  80             rescue Exception
  81               error_msg = "failed to parse seq from from \"" + a[ 2 ] + "\" [line: " + line + "]"
  82               raise IOError, error_msg
  83             end
  84             begin
  85               seq_to = a[ 3 ].to_i
  86             rescue Exception
  87               error_msg = "failed to parse seq to from \"" + a[ 3 ] + "\" [line: " + line + "]"
  88               raise IOError, error_msg
  89             end
  90
  91             e_value = -1
  92             if l > 4
  93               begin
  94                 e_value = a[ 4 ].to_f
  95               rescue Exception
  96                 error_msg = "failed to parse E-value from \"" + a[ 4 ] + "\" [line: " + line + "]"
  97                 raise IOError, error_msg
  98               end
  99             end
 100
 101             seq = original_seqs.get_by_name_start( protein_name )
 102
 103             total_length = seq.get_length
 104
 105             if ( ( ( e_value_threshold < 0.0 ) || ( e_value <= e_value_threshold ) )  )
 106               pd = ProteinDomain.new( domain_name, seq_from, seq_to, "", e_value )
 107               ds = nil
 108               if ( domain_structures.has_key?( protein_name ) )
 109                 ds = domain_structures[ protein_name ]
 110               else
 111                 ds = DomainStructure.new( total_length )
 112                 domain_structures[ protein_name ] = ds
 113               end
 114               ds.add_domain( pd, overwrite_if_same_from_to )
 115             end
 116
 117           end
 118         end
 119       end
 120
 121       out = File.open( outfile, "a" )
 122       ds = domain_structures.sort
 123       for d in ds
 124         protein_name     = d[ 0 ]
 125         domain_structure = d[ 1 ]
 126         out.print( protein_name.to_s )
 127         out.print( "\t" )
 128         out.print( domain_structure.to_NHX )
 129         out.print( Constants::LINE_DELIMITER  )
 130       end
 131
 132       out.flush()
 133       out.close()
 134
 135     end # parse
 136
 137
 138
 139
 140     def run()
 141
 142       Util.print_program_information( PRG_NAME,
 143         PRG_VERSION,
 144         PRG_DESC,
 145         PRG_DATE,
 146         COPYRIGHT,
 147         CONTACT,
 148         WWW,
 149         STDOUT )
 150
 151       begin
 152         cla = CommandLineArguments.new( ARGV )
 153       rescue ArgumentError => e
 154         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 155       end
 156
 157       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
 158            cla.is_option_set?( HELP_OPTION_2 ) )
 159         print_help
 160         exit( 0 )
 161       end
 162
 163       if cla.get_number_of_files != 3
 164         print_help
 165         exit( -1 )
 166       end
 167
 168       allowed_opts = Array.new
 169       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
 170       allowed_opts.push( OVERWRITE_IF_SAME_FROM_TO_OPTION )
 171
 172       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
 173       if ( disallowed.length > 0 )
 174         Util.fatal_error( PRG_NAME,
 175           "unknown option(s): " + disallowed,
 176           STDOUT )
 177       end
 178
 179       domains_list_file       = cla.get_file_name( 0 )
 180       original_sequences_file = cla.get_file_name( 1 )
 181       outfile                 = cla.get_file_name( 2 )
 182
 183
 184       e_value_threshold = -1.0
 185       if cla.is_option_set?( E_VALUE_THRESHOLD_OPTION )
 186         begin
 187           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
 188         rescue ArgumentError => e
 189           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 190         end
 191         if ( e_value_threshold < 0.0 )
 192           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 193         end
 194       end
 195       overwrite_if_same_from_to = false
 196       if ( cla.is_option_set?( OVERWRITE_IF_SAME_FROM_TO_OPTION ) )
 197         overwrite_if_same_from_to = true
 198       end
 199
 200       puts
 201       puts( "Domains list file                      : " + domains_list_file )
 202       puts( "Fasta sequencefile (complete sequences): " + original_sequences_file )
 203       puts( "Outputfile                             : " + outfile )
 204       if ( e_value_threshold >= 0.0 )
 205         puts( "E-value threshold                      : " + e_value_threshold.to_s )
 206       else
 207         puts( "E-value threshold                      : no threshold" )
 208       end
 209       if ( overwrite_if_same_from_to )
 210         puts( "Overwrite if same from and to          : true" )
 211       else
 212         puts( "Overwrite if same from and to          : false" )
 213       end
 214
 215       puts
 216
 217       begin
 218         parse( domains_list_file,
 219           original_sequences_file,
 220           outfile,
 221           " ",
 222           e_value_threshold,
 223           overwrite_if_same_from_to )
 224
 225       rescue ArgumentError, IOError, StandardError => e
 226         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 227       rescue Exception => e
 228         Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
 229       end
 230
 231
 232       puts
 233       Util.print_message( PRG_NAME, 'OK' )
 234       puts
 235
 236     end
 237
 238     private
 239
 240     def print_help()
 241       puts
 242       puts( "Usage:" )
 243       puts
 244       puts( "  " + PRG_NAME + ".rb [options] <domains list file (parsed hmmpfam output)> <file containing complete sequences in fasta format> <outputfile>" )
 245       puts()
 246       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f> : E-value threshold, default is no threshold" )
 247       puts( "               -" + OVERWRITE_IF_SAME_FROM_TO_OPTION  + " : overwrite domain with same start and end with domain with better E-value" )
 248       puts
 249     end
 250
 251
 252
 253     def is_ignorable?( line )
 254       return ( line !~ /[A-Za-z0-9-]/ || line =~ /^\s*#/)
 255     end
 256
 257
 258   end # class DomainsToForester
 259
 260
 261 end # module Evoruby