forester/ruby/evoruby/lib/evo/io/parser/hmmscan_multi_domain_extractor.rb

   1 #
   2 # = lib/evo/io/parser/hmmscan_domain_extractor.rb - HmmscanMultiDomainExtractor class
   3 #
   4 # Copyright::    Copyright (C) 2017 Christian M. Zmasek
   5 # License::      GNU Lesser General Public License (LGPL)
   6 #
   7 # Last modified: 2017/02/20
   8
   9 require 'lib/evo/util/constants'
  10 require 'lib/evo/msa/msa_factory'
  11 require 'lib/evo/io/msa_io'
  12 require 'lib/evo/io/writer/fasta_writer'
  13 require 'lib/evo/io/parser/fasta_parser'
  14 require 'lib/evo/io/parser/hmmscan_parser'
  15
  16 module Evoruby
  17   class HmmscanMultiDomainExtractor
  18
  19     TESTING = false
  20     OUTPUT_ID = 'mdsx'
  21     DOMAIN_DELIMITER = ' -- '
  22
  23     PASSING_FL_SEQS_SUFFIX    = "_#{OUTPUT_ID}_passing_full_length_seqs.fasta"
  24     FAILING_FL_SEQS_SUFFIX    = "_#{OUTPUT_ID}_failing_full_length_seqs.fasta"
  25     TARGET_DA_SUFFIX          = "_#{OUTPUT_ID}_target_da.fasta"
  26     CONCAT_TARGET_DOM_SUFFIX  = "_#{OUTPUT_ID}_concat_target_doms.fasta"
  27     TARGET_DOM_OUTPUT_MIDPART = "_#{OUTPUT_ID}_target_dom_"
  28     LOG_FILE_SUFFIX           = "_#{OUTPUT_ID}.log"
  29     def initialize
  30       @passing_domains_data = nil
  31       @failing_domains_data = nil
  32       @failing_domain_architectures = nil
  33       @passsing_domain_architectures = nil
  34       @failing_proteins_bc_not_all_target_doms_present = nil
  35       @failing_proteins_bc_missing_cutoffs = nil
  36       @log = nil
  37     end
  38
  39     # raises ArgumentError, IOError, StandardError
  40     def parse( target_da,
  41       hmmscan_output,
  42       fasta_sequence_file,
  43       outfile_base,
  44       log_str )
  45
  46       passing_fl_seqs_outfile   = outfile_base + PASSING_FL_SEQS_SUFFIX
  47       failing_fl_seqs_outfile   = outfile_base + FAILING_FL_SEQS_SUFFIX
  48       target_da_outfile         = outfile_base + TARGET_DA_SUFFIX
  49       concat_target_dom_outfile = outfile_base + CONCAT_TARGET_DOM_SUFFIX
  50       logfile                   = outfile_base + LOG_FILE_SUFFIX
  51
  52       Util.check_file_for_readability( hmmscan_output )
  53       Util.check_file_for_readability( fasta_sequence_file )
  54       Util.check_file_for_writability( passing_fl_seqs_outfile )
  55       Util.check_file_for_writability( failing_fl_seqs_outfile )
  56       Util.check_file_for_writability( target_da_outfile )
  57       Util.check_file_for_writability( concat_target_dom_outfile )
  58       Util.check_file_for_writability( logfile )
  59
  60       @log = log_str
  61
  62       in_msa = nil
  63       factory = MsaFactory.new()
  64       in_msa = factory.create_msa_from_file( fasta_sequence_file, FastaParser.new() )
  65
  66       if ( in_msa == nil || in_msa.get_number_of_seqs() < 1 )
  67         error_msg = "could not find fasta sequences in " + fasta_sequence_file
  68         raise IOError, error_msg
  69       end
  70
  71       failed_seqs_msa = Msa.new
  72       passed_seqs_msa = Msa.new
  73
  74       hmmscan_parser = HmmscanParser.new( hmmscan_output )
  75       results = hmmscan_parser.parse
  76
  77       ####
  78       # Import: if multiple copies of same domain, thresholds need to be same!
  79
  80       #target_domain_ary.push(TargetDomain.new('DNA_pol_B_exo1', 1e-6, -1, 0.6, 0 ))
  81       #target_domain_ary.push(TargetDomain.new('DNA_pol_B', 1e-6, -1, 0.6, 1 ))
  82
  83       #target_domain_ary.push(TargetDomain.new('Hexokinase_1', 1e-6, -1, 0.9, 0 ))
  84       #target_domain_ary.push(TargetDomain.new('Hexokinase_2', 1e-6, -1, 0.9, 1 ))
  85       # target_domain_ary.push(TargetDomain.new('Hexokinase_2', 0.1, -1, 0.5, 1 ))
  86       # target_domain_ary.push(TargetDomain.new('Hexokinase_1', 0.1, -1, 0.5, 1 ))
  87
  88       target_domain_ary = parse_da target_da
  89
  90       # target_domain_ary.push(TargetDomain.new('BH4', 1, -1, 0.2, 0 ))
  91       # target_domain_ary.push(TargetDomain.new('Bcl-2', 1, -1, 0.2, 1 ))
  92       # target_domain_ary.push(TargetDomain.new('Bcl-2', 1, -1, 0.2, 2 ))
  93       # target_domain_ary.push(TargetDomain.new('Bcl-2_3', 0.01, -1, 0.5, 2 ))
  94
  95       #  target_domain_ary.push(TargetDomain.new('Nitrate_red_del', 1000, -1, 0.1, 0 ))
  96       #  target_domain_ary.push(TargetDomain.new('Nitrate_red_del', 1000, -1, 0.1, 1 ))
  97
  98       #target_domain_ary.push(TargetDomain.new('Chordopox_A33R', 1000, -1, 0.1 ))
  99
 100       target_domains = Hash.new
 101
 102       target_domain_architecure = ""
 103
 104       target_domain_ary.each do |target_domain|
 105         target_domains[target_domain.name] = target_domain
 106         if target_domain_architecure.length > 0
 107           target_domain_architecure += DOMAIN_DELIMITER
 108         end
 109         target_domain_architecure += target_domain.name
 110       end
 111
 112       target_domain_architecure.freeze
 113
 114       log 'Hmmscan outputfile             : ' + hmmscan_output
 115       log 'Full length fasta sequence file: ' + fasta_sequence_file
 116       log 'Target domain architecture     : ' + target_domain_architecure
 117       target_domain_ary.each do |x|
 118         log x.to_str
 119       end
 120
 121       target_domain_names = Set.new
 122
 123       target_domains.each_key {|key| target_domain_names.add( key ) }
 124
 125       prev_query_seq_name = nil
 126       domains_in_query_seq = Array.new
 127       passing_sequences = Array.new # This will be an Array of Array of HmmscanResult for the same query seq.
 128       failing_sequences = Array.new # This will be an Array of Array of HmmscanResult for the same query seq.
 129       total_sequences = 0
 130
 131       @failing_domain_architectures = Hash.new
 132       @passsing_domain_architectures = Hash.new
 133       @passing_domains_data = Hash.new
 134       @failing_domains_data = Hash.new
 135       @failing_proteins_bc_not_all_target_doms_present = 0
 136       @failing_proteins_bc_missing_cutoffs = 0
 137       out_domain_msas = Hash.new
 138       out_domain_architecture_msa = Msa.new
 139       out_concatenated_domains_msa = Msa.new
 140       results.each do |hmmscan_result|
 141         if ( prev_query_seq_name != nil ) && ( hmmscan_result.query != prev_query_seq_name )
 142           if compare(domains_in_query_seq, target_domain_names, target_domains, in_msa, out_domain_msas, out_domain_architecture_msa, out_concatenated_domains_msa, target_domain_architecure)
 143             passing_sequences.push(domains_in_query_seq)
 144           else
 145             failing_sequences.push(domains_in_query_seq)
 146           end
 147           domains_in_query_seq = Array.new
 148           total_sequences += 1
 149         end
 150         prev_query_seq_name = hmmscan_result.query
 151         domains_in_query_seq.push(hmmscan_result)
 152       end # each result
 153
 154       if prev_query_seq_name != nil
 155         if compare(domains_in_query_seq, target_domain_names, target_domains, in_msa, out_domain_msas, out_domain_architecture_msa, out_concatenated_domains_msa, target_domain_architecure)
 156           passing_sequences.push(domains_in_query_seq)
 157         else
 158           failing_sequences.push(domains_in_query_seq)
 159         end
 160         total_sequences += 1
 161       end
 162
 163       if in_msa.get_number_of_seqs < total_sequences
 164         error_msg = "hmmscan output contains more protein sequences than fasta sequence file"
 165         raise IOError, error_msg
 166       end
 167
 168       log
 169       log 'Passing domain architectures containing target domain(s):'
 170       @passsing_domain_architectures = @passsing_domain_architectures.sort{|a, b|a<=>b}.to_h
 171       passing_da_sum = 0
 172       @passsing_domain_architectures.each do |da, count|
 173         passing_da_sum += count
 174         log count.to_s.rjust(4) + ': ' + da
 175       end
 176       log
 177       log 'Failing domain architectures containing one or more target domain(s):'
 178       @failing_domain_architectures = @failing_domain_architectures .sort{|a, b|a<=>b}.to_h
 179       failing_da_sum = 0
 180       @failing_domain_architectures .each do |da, count|
 181         failing_da_sum += count
 182         log count.to_s.rjust(4) + ': ' + da
 183       end
 184       log
 185       log 'Passing target domain(s):'
 186       @passing_domains_data = @passing_domains_data.sort{|a, b|a<=>b}.to_h
 187       @passing_domains_data.each do |n, d|
 188         log d.to_str
 189       end
 190       log
 191       log'Failing target domain(s):'
 192       @failing_domains_data = @failing_domains_data.sort{|a, b|a<=>b}.to_h
 193       @failing_domains_data.each do |n, d|
 194         log d.to_str
 195       end
 196
 197       unless total_sequences == (passing_sequences.size + failing_sequences.size)
 198         error_msg = "this should not have happened: total seqs not equal to passing plus failing seqs"
 199         raise StandardError, error_msg
 200       end
 201
 202       unless failing_sequences.size == (@failing_proteins_bc_not_all_target_doms_present + @failing_proteins_bc_missing_cutoffs)
 203         error_msg = "this should not have happened: failing seqs sums not consistent"
 204         raise StandardError, error_msg
 205       end
 206
 207       unless @failing_proteins_bc_missing_cutoffs == failing_da_sum
 208         error_msg = "this should not have happened: failing seqs not equal to failing da sum"
 209         raise StandardError, error_msg
 210       end
 211
 212       unless passing_sequences.size == passing_da_sum
 213         error_msg = "this should not have happened: passing seqs not equal to passing da sum"
 214         raise StandardError, error_msg
 215       end
 216
 217       log
 218       log "Protein sequences in sequence (fasta) file: " + in_msa.get_number_of_seqs.to_s.rjust(5)
 219       log "Protein sequences in hmmscan output file  : " + total_sequences.to_s.rjust(5)
 220       log "  Passing protein sequences               : " + passing_sequences.size.to_s.rjust(5)
 221       log "  Failing protein sequences               : " + failing_sequences.size.to_s.rjust(5)
 222       log "    Not all target domain present         : " + @failing_proteins_bc_not_all_target_doms_present.to_s.rjust(5)
 223       log "    Target domain(s) failing cutoffs      : " + @failing_proteins_bc_missing_cutoffs.to_s.rjust(5)
 224       log
 225
 226       out_domain_msas.keys.sort.each do |domain_name|
 227         file_name = outfile_base + TARGET_DOM_OUTPUT_MIDPART + domain_name + '.fasta'
 228         write_msa out_domain_msas[domain_name], file_name
 229         log "Wrote passing target domain sequence for " +  domain_name.ljust(16) + ': ' + file_name
 230       end
 231
 232       write_msa out_domain_architecture_msa, target_da_outfile
 233       log 'Wrote target domain architecture                         : ' + target_da_outfile
 234
 235       write_msa out_concatenated_domains_msa, concat_target_dom_outfile
 236       log 'Wrote concatenated target domain(s)                      : ' + concat_target_dom_outfile
 237
 238       passing_sequences.each do | domains |
 239         query_name = domains[0].query
 240         if (!TESTING) || (passed_seqs_msa.find_by_name_start( query_name, true ).length < 1)
 241           add_sequence( query_name, in_msa, passed_seqs_msa )
 242         else
 243           error_msg = 'this should not have happened'
 244           raise StandardError, error_msg
 245         end
 246       end
 247
 248       failing_sequences.each do | domains |
 249         query_name = domains[0].query
 250         if (!TESTING) || (failed_seqs_msa.find_by_name_start( query_name, true ).length < 1)
 251           add_sequence( query_name, in_msa, failed_seqs_msa )
 252         else
 253           error_msg = 'this should not have happened'
 254           raise StandardError, error_msg
 255         end
 256       end
 257
 258       write_msa passed_seqs_msa, passing_fl_seqs_outfile
 259       log 'Wrote passing full length protein sequences              : ' + passing_fl_seqs_outfile
 260       write_msa failed_seqs_msa, failing_fl_seqs_outfile
 261       log 'Wrote failing full length protein sequences              : ' + failing_fl_seqs_outfile
 262
 263       begin
 264         f = File.open( logfile, 'w' )
 265         f.print( @log )
 266         f.close
 267       rescue Exception => e
 268         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 269       end
 270       log 'Wrote log file                                           : ' + logfile
 271
 272     end # parse
 273
 274     private
 275
 276     # domains_in_query_seq: Array of HmmscanResult
 277     # target_domain_names: Set of String
 278     # target_domains: Hash String->TargetDomain
 279     # target_domain_architecture: String
 280     def compare(domains_in_query_seq,
 281       target_domain_names,
 282       target_domains,
 283       in_msa,
 284       out_domain_msas,
 285       out_domain_architecture_msa,
 286       out_contactenated_domains_msa,
 287       target_domain_architecture)
 288
 289       domain_names_in_query_seq = Set.new
 290       domains_in_query_seq.each do |domain|
 291         domain_names_in_query_seq.add(domain.model)
 292       end
 293       if target_domain_names.subset?(domain_names_in_query_seq)
 294
 295         passed_domains = Array.new
 296         passed_domains_counts = Hash.new
 297
 298         domains_in_query_seq.each do |domain|
 299           if target_domains.has_key?(domain.model)
 300             target_domain = target_domains[domain.model]
 301
 302             if (target_domain.i_e_value != nil) && (target_domain.i_e_value >= 0)
 303               if domain.i_e_value > target_domain.i_e_value
 304                 addToFailingDomainData(domain)
 305                 next
 306               end
 307             end
 308             if (target_domain.abs_len != nil) && (target_domain.abs_len > 0)
 309               length = 1 + domain.env_to - domain.env_from
 310               if length < target_domain.abs_len
 311                 addToFailingDomainData(domain)
 312                 next
 313               end
 314             end
 315             if (target_domain.rel_len != nil) && (target_domain.rel_len > 0)
 316               length = 1 + domain.env_to - domain.env_from
 317               if length < (target_domain.rel_len * domain.tlen)
 318                 addToFailingDomainData(domain)
 319                 next
 320               end
 321             end
 322
 323             passed_domains.push(domain)
 324
 325             if (passed_domains_counts.key?(domain.model))
 326               passed_domains_counts[domain.model] = passed_domains_counts[domain.model] + 1
 327             else
 328               passed_domains_counts[domain.model] = 1
 329             end
 330
 331             addToPassingDomainData(domain)
 332
 333           end # if target_domains.has_key?(domain.model)
 334         end # domains_in_query_seq.each do |domain|
 335       else
 336         @failing_proteins_bc_not_all_target_doms_present += 1
 337         return false
 338       end # target_domain_names.subset?(domain_names_in_query_seq)
 339
 340       passed_domains.sort! { |a,b| a.ali_from <=> b.ali_from }
 341       # Compare architectures
 342       if !compareArchitectures(target_domain_architecture, passed_domains)
 343         @failing_proteins_bc_missing_cutoffs += 1
 344         return false
 345       end
 346
 347       domain_counter = Hash.new
 348
 349       min_env_from = 10000000
 350       max_env_to = 0
 351
 352       query_seq = nil
 353
 354       concatenated_domains = ''
 355
 356       passed_domains.each do |domain|
 357         domain_name = domain.model
 358
 359         unless out_domain_msas.has_key? domain_name
 360           out_domain_msas[ domain_name ] = Msa.new
 361         end
 362
 363         if domain_counter.key?(domain_name)
 364           domain_counter[domain_name] = domain_counter[domain_name] + 1
 365         else
 366           domain_counter[domain_name] = 1
 367         end
 368
 369         if query_seq == nil
 370           query_seq = domain.query
 371           query_seq.freeze
 372         elsif query_seq != domain.query
 373           error_msg = "seq names do not match: this should not have happened"
 374           raise StandardError, error_msg
 375         end
 376
 377         extracted_domain_seq = extract_domain( query_seq,
 378         domain_counter[domain_name],
 379         passed_domains_counts[domain_name],
 380         domain.env_from,
 381         domain.env_to,
 382         in_msa,
 383         out_domain_msas[ domain_name ] )
 384
 385         concatenated_domains += extracted_domain_seq
 386
 387         if domain.env_from < min_env_from
 388           min_env_from = domain.env_from
 389         end
 390         if domain.env_to > max_env_to
 391           max_env_to = domain.env_to
 392         end
 393       end
 394
 395       extract_domain( query_seq,
 396       1,
 397       1,
 398       min_env_from,
 399       max_env_to,
 400       in_msa,
 401       out_domain_architecture_msa )
 402
 403       out_contactenated_domains_msa.add( query_seq, concatenated_domains)
 404
 405       true
 406     end
 407
 408     def addToPassingDomainData(domain)
 409       unless ( @passing_domains_data.key?(domain.model))
 410         @passing_domains_data[domain.model] = DomainData.new(domain.model)
 411       end
 412       @passing_domains_data[domain.model].add( domain.model, (1 + domain.env_to - domain.env_from), domain.i_e_value)
 413     end
 414
 415     def addToFailingDomainData(domain)
 416       unless ( @failing_domains_data.key?(domain.model))
 417         @failing_domains_data[domain.model] = DomainData.new(domain.model)
 418       end
 419       @failing_domains_data[domain.model].add( domain.model, (1 + domain.env_to - domain.env_from), domain.i_e_value)
 420     end
 421
 422     # passed_domains needs to be sorted!
 423     def compareArchitectures(target_domain_architecture, passed_domains, strict = false)
 424       domain_architecture = ''
 425       passed_domains.each do |domain|
 426         if domain_architecture.length > 0
 427           domain_architecture += DOMAIN_DELIMITER
 428         end
 429         domain_architecture += domain.model
 430       end
 431       pass = false
 432       if strict
 433         pass = (target_domain_architecture == domain_architecture)
 434       else
 435         pass = (domain_architecture.index(target_domain_architecture) != nil)
 436       end
 437       if pass
 438         if (@passsing_domain_architectures.key?(domain_architecture))
 439           @passsing_domain_architectures[domain_architecture] = @passsing_domain_architectures[domain_architecture] + 1
 440         else
 441           @passsing_domain_architectures[domain_architecture] = 1
 442         end
 443       else
 444         if ( @failing_domain_architectures .key?(domain_architecture))
 445           @failing_domain_architectures [domain_architecture] =  @failing_domain_architectures [domain_architecture] + 1
 446         else
 447           @failing_domain_architectures [domain_architecture] = 1
 448         end
 449       end
 450       return pass
 451     end
 452
 453     def write_msa( msa, filename )
 454       io = MsaIO.new()
 455       w = FastaWriter.new()
 456       w.set_line_width( 60 )
 457       w.clean( true )
 458       File.delete(filename) if File.exist?(filename)
 459       begin
 460         io.write_to_file( msa, filename, w )
 461       rescue Exception
 462         error_msg = "could not write to \"" + filename + "\""
 463         raise IOError, error_msg
 464       end
 465     end
 466
 467     def add_sequence( sequence_name, in_msa, add_to_msa )
 468       seqs = in_msa.find_by_name_start( sequence_name, true )
 469       if ( seqs.length < 1 )
 470         error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 471         raise StandardError, error_msg
 472       end
 473       if ( seqs.length > 1 )
 474         error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 475         raise StandardError, error_msg
 476       end
 477       seq = in_msa.get_sequence( seqs[ 0 ] )
 478       add_to_msa.add_sequence( seq )
 479     end
 480
 481     def extract_domain( seq_name,
 482       number,
 483       out_of,
 484       seq_from,
 485       seq_to,
 486       in_msa,
 487       out_msa)
 488       if number < 1 || out_of < 1 || number > out_of
 489         error_msg = "number=" + number.to_s + ", out of=" + out_of.to_s
 490         raise StandardError, error_msg
 491       end
 492       if seq_from < 1 || seq_to < 1 || seq_from >= seq_to
 493         error_msg = "impossible: seq-from=" + seq_from.to_s + ", seq-to=" + seq_to.to_s
 494         raise StandardError, error_msg
 495       end
 496       seqs = in_msa.find_by_name_start(seq_name, true)
 497       if seqs.length < 1
 498         error_msg = "sequence \"" + seq_name + "\" not found in sequence file"
 499         raise IOError, error_msg
 500       end
 501       if seqs.length > 1
 502         error_msg = "sequence \"" + seq_name + "\" not unique in sequence file"
 503         raise IOError, error_msg
 504       end
 505       # hmmscan is 1 based, whereas sequences are 0 bases in this package.
 506       seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 507
 508       orig_name = seq.get_name
 509
 510       seq.set_name( orig_name.split[ 0 ] )
 511
 512       if out_of != 1
 513         seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 514       end
 515
 516       out_msa.add_sequence( seq )
 517
 518       seq.get_sequence_as_string
 519     end
 520
 521     def parse_da( target_da_str )
 522       target_domain_hash = Hash.new
 523       target_domain_ary = Array.new
 524       target_das = target_da_str.split '--'
 525       target_das.each do |x|
 526         inds = x.split '='
 527         unless inds.size == 4
 528           raise IOError, 'domain architecture is ill formatted: ' + x
 529         end
 530         target_domain_name = inds[0]
 531         ie_cutoff = Float(inds[1])
 532         abs_len_cutoff = Integer(inds[2])
 533         rel_len_cutoff = Float(inds[3])
 534         if target_domain_hash.has_key? target_domain_name
 535           target_domain_ary.push target_domain_hash[target_domain_name]
 536         else
 537           td = TargetDomain.new(target_domain_name, ie_cutoff, abs_len_cutoff, rel_len_cutoff)
 538           target_domain_hash[target_domain_name] = td
 539           target_domain_ary.push td
 540         end
 541       end
 542       target_domain_ary
 543     end
 544
 545     def log(str = '')
 546       puts str
 547       @log << str << Constants::LINE_DELIMITER
 548     end
 549
 550   end # class HmmscanMultiDomainExtractor
 551
 552   class DomainData
 553     def initialize( name )
 554       if (name == nil) || name.size < 1
 555         error_msg = "domain name nil or empty"
 556         raise IOError, error_msg
 557       end
 558       @name = name
 559       @count = 0
 560       @i_e_value_min = 10000000.0
 561       @i_e_value_max = -1.0
 562       @i_e_value_sum = 0.0
 563       @len_min = 10000000
 564       @len_max = -1
 565       @len_sum = 0.0
 566     end
 567
 568     def add( name, length, i_e_value)
 569       if name != @name
 570         error_msg = "domain names do not match"
 571         raise IOError, error_msg
 572       end
 573
 574       if length < 0
 575         error_msg = "length cannot me negative"
 576         raise IOError, error_msg
 577       end
 578       if i_e_value < 0
 579         error_msg = "iE-value cannot me negative"
 580         raise IOError, error_msg
 581       end
 582       @count += 1
 583       @i_e_value_sum += i_e_value
 584       @len_sum += length
 585       if i_e_value > @i_e_value_max
 586         @i_e_value_max = i_e_value
 587       end
 588       if i_e_value < @i_e_value_min
 589         @i_e_value_min = i_e_value
 590       end
 591       if length > @len_max
 592         @len_max = length
 593       end
 594       if length < @len_min
 595         @len_min = length
 596       end
 597     end
 598
 599     def avg_length
 600       if @count == 0
 601         return 0
 602       end
 603       @len_sum / @count
 604     end
 605
 606     def avg_i_e_value
 607       if @count == 0
 608         return 0
 609       end
 610       @i_e_value_sum / @count
 611     end
 612
 613     def to_str
 614       s = ''
 615       s << @name.rjust(16) + ': '
 616       s << @count.to_s.rjust(4) + '  '
 617       s << avg_length.round(1).to_s.rjust(6) + ' '
 618       s << @len_min.to_s.rjust(4) + ' -'
 619       s << @len_max.to_s.rjust(4) + '  '
 620       s << ("%.2E" % avg_i_e_value).rjust(9) + ' '
 621       s << ("%.2E" % @i_e_value_min).rjust(9) + ' -'
 622       s << ("%.2E" % @i_e_value_max).rjust(9) + ' '
 623       s
 624     end
 625
 626     attr_reader :name, :count, :i_e_value_min, :i_e_value_max, :length_min, :length_max
 627
 628   end
 629
 630   class TargetDomain
 631     def initialize(name, i_e_value, abs_len, rel_len)
 632       if (name == nil) || name.size < 1
 633         error_msg = "target domain name nil or empty"
 634         raise IOError, error_msg
 635       end
 636       if rel_len > 1
 637         error_msg = name + ": target domain relative length is greater than 1"
 638         raise IOError, error_msg
 639       end
 640       if (abs_len <= 0) && (rel_len <= 0)
 641         error_msg = name + ": need to have either absolute length or relative length cutoff"
 642         raise IOError, error_msg
 643       end
 644       if (abs_len > 0) && (rel_len > 0)
 645         error_msg = name + ": cannot have both absolute length and relative length cutoff"
 646         raise IOError, error_msg
 647       end
 648       @name = name
 649       @i_e_value = i_e_value
 650       @abs_len = abs_len
 651       @rel_len = rel_len
 652     end
 653
 654     def to_str
 655       s = @name.rjust(16) + ':'
 656       s << ' iE-cutoff: ' + ("%.2E" % @i_e_value).rjust(9)
 657       if @abs_len > 0
 658         s << ', abs len-cutoff: ' + @abs_len.to_s.rjust(4)
 659       end
 660       if @rel_len > 0
 661         s << ', rel len-cutoff: ' + @rel_len.to_s.rjust(4)
 662       end
 663       s
 664     end
 665     attr_reader :name, :i_e_value, :abs_len, :rel_len
 666   end
 667
 668 end # module Evoruby
 669