2 # = lib/evo/tool/hmmscan_summary.rb - HmmscanSummary class
4 # Copyright:: Copyright (C) 2012 Christian M. Zmasek
5 # License:: GNU Lesser General Public License (LGPL)
7 # $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $
12 require 'lib/evo/util/constants'
13 require 'lib/evo/util/util'
14 require 'lib/evo/util/command_line_arguments'
15 require 'lib/evo/io/parser/hmmscan_parser'
16 require 'lib/evo/io/web/uniprotkb'
24 PRG_DESC = "hmmscan summary"
25 PRG_DATE = "2013.10.23"
26 COPYRIGHT = "2013 Christian M Zmasek"
27 CONTACT = "phyloxml@gmail.com"
28 WWW = "https://sites.google.com/site/cmzmasek/home/software/forester"
30 DELIMITER_OPTION = "d"
32 I_E_VALUE_THRESHOLD_OPTION = "ie"
33 FS_E_VALUE_THRESHOLD_OPTION = "pe"
34 HMM_FOR_PROTEIN_OUTPUT = "m"
35 IGNORE_DUF_OPTION = "i"
36 PARSE_OUT_DESCRIPITION_OPTION = "a"
38 HELP_OPTION_1 = "help"
42 AVOID_HHMS = [ "RRM_1", "RRM_2", "RRM_3", "RRM_4", "RRM_5", "RRM_6" ]
43 LIMIT_FOR_CLOSE_DOMAINS = 20
46 @domain_counts = Hash.new
51 Util.print_program_information( PRG_NAME,
61 cla = CommandLineArguments.new( ARGV )
62 rescue ArgumentError => e
63 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
66 if ( cla.is_option_set?( HELP_OPTION_1 ) ||
67 cla.is_option_set?( HELP_OPTION_2 ) )
72 if ( cla.get_number_of_files != 2 )
77 allowed_opts = Array.new
78 allowed_opts.push( DELIMITER_OPTION )
79 allowed_opts.push( I_E_VALUE_THRESHOLD_OPTION )
80 allowed_opts.push( FS_E_VALUE_THRESHOLD_OPTION )
81 allowed_opts.push( IGNORE_DUF_OPTION )
82 allowed_opts.push( PARSE_OUT_DESCRIPITION_OPTION )
83 allowed_opts.push( HMM_FOR_PROTEIN_OUTPUT )
84 allowed_opts.push( UNIPROT )
85 allowed_opts.push( SPECIES_OPTION )
87 disallowed = cla.validate_allowed_options_as_str( allowed_opts )
88 if ( disallowed.length > 0 )
89 Util.fatal_error( PRG_NAME,
90 "unknown option(s): " + disallowed,
94 inpath = cla.get_file_name( 0 )
95 outpath = cla.get_file_name( 1 )
97 column_delimiter = "\t"
98 if ( cla.is_option_set?( DELIMITER_OPTION ) )
100 column_delimiter = cla.get_option_value( DELIMITER_OPTION )
101 rescue ArgumentError => e
102 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
106 i_e_value_threshold = -1.0
107 if ( cla.is_option_set?( I_E_VALUE_THRESHOLD_OPTION ) )
109 i_e_value_threshold = cla.get_option_value_as_float( I_E_VALUE_THRESHOLD_OPTION )
110 rescue ArgumentError => e
111 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
113 if ( i_e_value_threshold < 0.0 )
114 Util.fatal_error( PRG_NAME, "attempt to use a negative i-E-value threshold", STDOUT )
120 fs_e_value_threshold = -1.0
121 if ( cla.is_option_set?( FS_E_VALUE_THRESHOLD_OPTION ) )
123 fs_e_value_threshold = cla.get_option_value_as_float( FS_E_VALUE_THRESHOLD_OPTION )
124 rescue ArgumentError => e
125 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
127 if ( fs_e_value_threshold < 0.0 )
128 Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
132 hmm_for_protein_output = ""
133 if ( cla.is_option_set?( HMM_FOR_PROTEIN_OUTPUT ) )
135 hmm_for_protein_output = cla.get_option_value( HMM_FOR_PROTEIN_OUTPUT )
136 rescue ArgumentError => e
137 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
142 if ( cla.is_option_set?( UNIPROT ) )
144 uniprot = cla.get_option_value( UNIPROT )
145 rescue ArgumentError => e
146 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
151 if ( cla.is_option_set?( SPECIES_OPTION ) )
153 species = cla.get_option_value( SPECIES_OPTION )
154 rescue ArgumentError => e
155 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
160 if ( cla.is_option_set?( IGNORE_DUF_OPTION ) )
164 parse_descriptions = false
165 if ( cla.is_option_set?( PARSE_OUT_DESCRIPITION_OPTION ) )
166 parse_descriptions = true
170 puts( "hmmpfam outputfile : " + inpath )
171 puts( "outputfile : " + outpath )
172 puts( "species : " + species )
173 if ( i_e_value_threshold >= 0.0 )
174 puts( "i-E-value threshold : " + i_e_value_threshold.to_s )
176 puts( "i-E-value threshold : no threshold" )
178 if ( parse_descriptions )
179 puts( "parse descriptions : true" )
181 puts( "parse descriptions : false" )
184 puts( "ignore DUFs : true" )
186 puts( "ignore DUFs : false" )
188 if ( column_delimiter == "\t" )
189 puts( "column delimiter : TAB" )
191 puts( "column delimiter : " + column_delimiter )
193 if fs_e_value_threshold >= 0.0
194 puts( "E-value threshold : " + fs_e_value_threshold.to_s )
196 puts( "E-value threshold : no threshold" )
198 if !hmm_for_protein_output.empty?
199 puts( "HMM for proteins : " + hmm_for_protein_output )
202 puts( "Uniprot : " + uniprot )
213 fs_e_value_threshold,
214 hmm_for_protein_output,
218 Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
220 domain_counts = get_domain_counts()
223 puts( "domain counts (considering potential i-E-value threshold and ignoring of DUFs):" )
224 puts( "(number of different domains: " + domain_counts.length.to_s + ")" )
226 puts( Util.draw_histogram( domain_counts, "#" ) )
228 Util.print_message( PRG_NAME, 'OK' )
235 # raises ArgumentError, IOError
242 fs_e_value_threshold,
243 hmm_for_protein_output,
247 Util.check_file_for_readability( inpath )
248 Util.check_file_for_writability( outpath )
250 hmmscan_parser = HmmscanParser.new( inpath )
251 results = hmmscan_parser.parse
253 outfile = File.open( outpath, "a" )
262 hmmscan_results_per_protein = []
266 results.each do | r |
269 i_e_value = r.i_e_value
270 env_from = r.env_from
273 if ( ( i_e_value_threshold < 0.0 ) || ( i_e_value <= i_e_value_threshold ) ) &&
274 ( !ignore_dufs || ( model !~ /^DUF\d+/ ) )
276 outfile.print( query +
278 if ( get_descriptions )
279 outfile.print( desc +
282 outfile.print( model +
289 outfile.print( Constants::LINE_DELIMITER )
292 if !hmm_for_protein_output.empty?
293 if !prev_query.empty? && prev_query != query
294 if !hmmscan_results_per_protein.empty?
295 process_hmmscan_results_per_protein( hmmscan_results_per_protein,
296 fs_e_value_threshold,
297 hmm_for_protein_output,
302 hmmscan_results_per_protein.clear
307 if !AVOID_HHMS.include? r.model
308 hmmscan_results_per_protein << r
311 hmmscan_results_per_protein << r
316 if !hmm_for_protein_output.empty? && !hmmscan_results_per_protein.empty?
317 process_hmmscan_results_per_protein( hmmscan_results_per_protein,
318 fs_e_value_threshold,
319 hmm_for_protein_output,
330 if id =~ /(sp|tr)\|\S+\|(\S+)/
336 def count_model( model )
337 if ( @domain_counts.has_key?( model ) )
338 count = @domain_counts[ model ].to_i
340 @domain_counts[ model ] = count
342 @domain_counts[ model ] = 1
346 def process_hmmscan_results_per_protein( hmmscan_results_per_protein,
347 fs_e_value_threshold,
348 hmm_for_protein_output,
354 # filter according to i-Evalue threshold
355 # abort if fs Evalue too high
356 hmmscan_results_per_protein_filtered = []
358 hmmscan_results_per_protein.each do | r |
361 if r.model == hmm_for_protein_output
362 if fs_e_value_threshold > 0.0 && r.fs_e_value > fs_e_value_threshold
366 if i_e_value_threshold <= 0 || r.i_e_value <= i_e_value_threshold
367 hmmscan_results_per_protein_filtered << r
368 if r.model == hmm_for_protein_output
375 # passed on protein E-value, failed in per domain E-values
379 hmmscan_results_per_protein_filtered.sort! { |r1,r2| r1.env_from <=> r2.env_from }
382 hmmscan_results_per_protein_filtered.each do | r |
383 if r.model == hmm_for_protein_output
389 s << own.query + "\t"
391 s << own.fs_e_value.to_s + "\t"
392 s << own.qlen.to_s + "\t"
394 s << hmmscan_results_per_protein_filtered.length.to_s + "\t"
395 hmmscan_results_per_protein_filtered.each do | r |
401 #e = UniprotKB::get_entry_by_id( process_id( own.query ) )
404 # s << uniprot_annotation( e )
405 # # s << "\uniprot_annotationt"
409 overview = make_overview( hmmscan_results_per_protein_filtered, hmm_for_protein_output )
413 s << calc_linkers( hmmscan_results_per_protein_filtered, hmm_for_protein_output ) + "\t"
416 hmmscan_results_per_protein_filtered.each do | r |
419 s << make_interdomain_sequence( r.env_from - prev_r.env_to - 1 )
421 s << make_interdomain_sequence( r.env_from, false )
425 s << r.env_from.to_s << "-" << r.env_to.to_s
426 s << "|ie=" << r.i_e_value.to_s
427 s << "|ce=" << r.c_e_value.to_s
431 s << make_interdomain_sequence( own.qlen - prev_r.env_from, false )
435 def uniprot_annotation( e )
437 pdb_ids = e.get_pdb_ids
439 pdb_ids.each do | pdb |
448 def calc_linkers( hmmscan_results_per_protein_filtered, hmm_for_protein_output )
451 hmmscan_results_per_protein_filtered.each do | r |
452 if r.model == hmm_for_protein_output
454 linkers << ( r.env_from - prev_r.env_to - 1 ).to_s + " "
462 def get_domain_counts()
463 return @domain_counts
466 def make_overview( hmmscan_results_per_protein_filtered, hmm_for_protein_output )
469 hmmscan_results_per_protein_filtered.each do | r |
470 if r.model == hmm_for_protein_output
472 overview << hmm_for_protein_output
474 if ( r.env_from - prev_r.env_to - 1 ) <= LIMIT_FOR_CLOSE_DOMAINS
475 overview << "~" << hmm_for_protein_output
477 overview << "----" << hmm_for_protein_output
486 def make_interdomain_sequence( d, mark_short = true )
505 puts( " " + PRG_NAME + ".rb [options] <hmmscan outputfile> <outputfile>" )
507 puts( " options: -" + DELIMITER_OPTION + ": column delimiter for outputfile, default is TAB" )
508 puts( " -" + I_E_VALUE_THRESHOLD_OPTION + ": i-E-value threshold, default is no threshold" )
509 puts( " -" + PARSE_OUT_DESCRIPITION_OPTION + ": parse query description (in addition to query name)" )
510 puts( " -" + IGNORE_DUF_OPTION + ": ignore DUFs" )
511 puts( " -" + FS_E_VALUE_THRESHOLD_OPTION + ": E-value threshold for full protein sequences, only for protein summary" )
512 puts( " -" + HMM_FOR_PROTEIN_OUTPUT + ": HMM for protein summary" )
513 puts( " -" + SPECIES_OPTION + ": species for protein summary" )