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

   1 #
   2 # = lib/evo/io/parser/hmmscan_domain_extractor.rb - HmmscanDomainExtractor class
   3 #
   4 # Copyright::  Copyright (C) 2012 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id:  $
   8
   9
  10 require 'lib/evo/util/constants'
  11 require 'lib/evo/msa/msa_factory'
  12 require 'lib/evo/io/msa_io'
  13 require 'lib/evo/io/writer/fasta_writer'
  14 require 'lib/evo/io/parser/fasta_parser'
  15 require 'lib/evo/io/parser/hmmscan_parser'
  16
  17
  18
  19 module Evoruby
  20
  21   class HmmscanDomainExtractor
  22
  23     ADD_TO_CLOSE_PAIRS = 20
  24
  25     def initialize
  26     end
  27
  28     # raises ArgumentError, IOError, StandardError
  29     def parse( domain_id,
  30         hmmscan_output,
  31         fasta_sequence_file,
  32         outfile,
  33         passed_seqs_outfile,
  34         failed_seqs_outfile,
  35         e_value_threshold,
  36         length_threshold,
  37         add_position,
  38         add_domain_number,
  39         add_species,
  40         min_linker,
  41         log )
  42
  43       Util.check_file_for_readability( hmmscan_output )
  44       Util.check_file_for_readability( fasta_sequence_file )
  45       Util.check_file_for_writability( outfile + ".fasta" )
  46       Util.check_file_for_writability( passed_seqs_outfile )
  47       Util.check_file_for_writability( failed_seqs_outfile )
  48
  49       in_msa = nil
  50       factory = MsaFactory.new()
  51       in_msa = factory.create_msa_from_file( fasta_sequence_file, FastaParser.new() )
  52
  53       if ( in_msa == nil || in_msa.get_number_of_seqs() < 1 )
  54         error_msg = "could not find fasta sequences in " + fasta_sequence_file
  55         raise IOError, error_msg
  56       end
  57
  58       out_msa = Msa.new
  59
  60       failed_seqs = Msa.new
  61       passed_seqs = Msa.new
  62       out_msa_pairs = nil
  63       out_msa_isolated = nil
  64       out_msa_singles = nil
  65       out_msa_single_domains_protein_seqs = nil
  66       out_msa_close_pairs_protein_seqs = nil
  67       out_msa_close_pairs_only_protein_seqs = nil
  68       out_msa_isolated_protein_seqs = nil
  69       out_msa_isolated_only_protein_seqs = nil
  70       out_msa_isolated_and_close_pair_protein_seqs = nil
  71       if min_linker
  72         out_msa_pairs = Msa.new
  73         out_msa_isolated = Msa.new
  74         out_msa_singles = Msa.new
  75         out_msa_single_domains_protein_seqs = Msa.new
  76         out_msa_close_pairs_protein_seqs = Msa.new
  77         out_msa_close_pairs_only_protein_seqs = Msa.new
  78         out_msa_isolated_protein_seqs = Msa.new
  79         out_msa_isolated_only_protein_seqs = Msa.new
  80         out_msa_isolated_and_close_pair_protein_seqs  = Msa.new
  81       end
  82
  83       ld = Constants::LINE_DELIMITER
  84
  85       domain_pass_counter     = 0
  86       domain_fail_counter     = 0
  87       proteins_with_failing_domains = 0
  88       max_domain_copy_number_per_protein = -1
  89       max_domain_copy_number_sequence    = ""
  90
  91       hmmscan_datas = []
  92
  93       hmmscan_parser = HmmscanParser.new( hmmscan_output )
  94       results = hmmscan_parser.parse
  95
  96       results.each do | r |
  97         if domain_id != r.model
  98           next
  99         end
 100
 101         sequence  = r.query
 102         number    = r.number
 103         out_of    = r.out_of
 104         env_from  = r.env_from
 105         env_to    = r.env_to
 106         i_e_value = r.i_e_value
 107
 108         if ( ( ( e_value_threshold < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
 109              ( ( length_threshold <= 0 )   || ( env_to - env_from + 1 ) >= length_threshold.to_f ) )
 110           hmmscan_datas << HmmsearchData.new( sequence, number, out_of, env_from, env_to, i_e_value )
 111           if ( number > max_domain_copy_number_per_protein )
 112             max_domain_copy_number_sequence    = sequence
 113             max_domain_copy_number_per_protein = number
 114           end
 115         else # failed
 116           print( domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)" )
 117           log << domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)"
 118           if ( ( e_value_threshold.to_f >= 0.0 ) && ( i_e_value > e_value_threshold ) )
 119             print( " iE=" + i_e_value.to_s )
 120             log << " iE=" + i_e_value.to_s
 121           end
 122           if ( ( length_threshold.to_f > 0 ) && ( env_to - env_from + 1 ) < length_threshold.to_f )
 123             le = env_to - env_from + 1
 124             print( " l=" + le.to_s )
 125             log << " l=" + le.to_s
 126           end
 127           print( ld )
 128           log << ld
 129           domain_fail_counter  += 1
 130         end
 131
 132         if number > out_of
 133           error_msg = "number > out_of ! (this should not have happened)"
 134           raise StandardError, error_msg
 135         end
 136
 137         if number == out_of
 138           if !hmmscan_datas.empty?
 139             process_hmmscan_datas( hmmscan_datas,
 140               in_msa,
 141               add_position,
 142               add_domain_number,
 143               add_species,
 144               out_msa,
 145               out_msa_singles,
 146               out_msa_pairs,
 147               out_msa_isolated,
 148               min_linker,
 149               out_msa_single_domains_protein_seqs,
 150               out_msa_close_pairs_protein_seqs,
 151               out_msa_close_pairs_only_protein_seqs,
 152               out_msa_isolated_protein_seqs,
 153               out_msa_isolated_only_protein_seqs,
 154               out_msa_isolated_and_close_pair_protein_seqs )
 155             domain_pass_counter += hmmscan_datas.length
 156             if passed_seqs.find_by_name_start( sequence, true ).length < 1
 157               add_sequence( sequence, in_msa, passed_seqs )
 158             else
 159               error_msg = "this should not have happened"
 160               raise StandardError, error_msg
 161             end
 162           else
 163             if failed_seqs.find_by_name_start( sequence, true ).length < 1
 164               add_sequence( sequence, in_msa, failed_seqs )
 165               proteins_with_failing_domains += 1
 166             else
 167               error_msg = "this should not have happened"
 168               raise StandardError, error_msg
 169             end
 170           end
 171           hmmscan_datas.clear
 172         end
 173       end
 174
 175       if domain_pass_counter < 1
 176         error_msg = "no domain sequences were extracted"
 177         raise IOError, error_msg
 178       end
 179
 180       log << ld
 181       puts( "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s )
 182       log << "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s
 183       log << ld
 184
 185       if ( max_domain_copy_number_per_protein > 1 )
 186         puts( "First protein with this copy number: " + max_domain_copy_number_sequence )
 187         log << "First protein with this copy number: " + max_domain_copy_number_sequence
 188         log << ld
 189       end
 190
 191       write_msa( out_msa, outfile + ".fasta"  )
 192       write_msa( passed_seqs, passed_seqs_outfile )
 193       write_msa( failed_seqs, failed_seqs_outfile )
 194
 195       if out_msa_pairs
 196         write_msa( out_msa_pairs, outfile +"_" + min_linker.to_s + ".fasta")
 197       end
 198
 199       if out_msa_singles
 200         write_msa( out_msa_singles, outfile +"_singles.fasta")
 201       end
 202
 203       if out_msa_isolated
 204         write_msa( out_msa_isolated, outfile + "_" + min_linker.to_s + "_isolated.fasta");
 205       end
 206
 207       if out_msa_single_domains_protein_seqs
 208         write_msa( out_msa_single_domains_protein_seqs, outfile +"_proteins_with_singles.fasta" )
 209       end
 210
 211       if out_msa_close_pairs_protein_seqs
 212         write_msa( out_msa_close_pairs_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_close_pairs.fasta" )
 213       end
 214
 215       if out_msa_close_pairs_only_protein_seqs
 216         write_msa( out_msa_close_pairs_only_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_close_pairs_only.fasta" )
 217       end
 218
 219       if  out_msa_isolated_protein_seqs
 220         write_msa(  out_msa_isolated_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_domains.fasta" )
 221       end
 222
 223       if out_msa_isolated_only_protein_seqs
 224         write_msa( out_msa_isolated_only_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_domains_only.fasta" )
 225       end
 226
 227       if out_msa_isolated_and_close_pair_protein_seqs
 228         write_msa( out_msa_isolated_and_close_pair_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_and_close_pairs.fasta" )
 229       end
 230
 231       if min_linker
 232         if ( out_msa_single_domains_protein_seqs.get_number_of_seqs +
 233              out_msa_close_pairs_only_protein_seqs.get_number_of_seqs +
 234              out_msa_isolated_only_protein_seqs.get_number_of_seqs +
 235              out_msa_isolated_and_close_pair_protein_seqs.get_number_of_seqs ) != passed_seqs.get_number_of_seqs
 236           error_msg = "this should not have happened"
 237           raise StandardError, error_msg
 238         end
 239       end
 240
 241       log << ld
 242       log << "passing domains                             : " + domain_pass_counter.to_s + ld
 243       log << "failing domains                             : " + domain_fail_counter.to_s + ld
 244       log << "input proteins                              : " + in_msa.get_number_of_seqs.to_s + ld
 245       log << "proteins with passing domains               : " + passed_seqs.get_number_of_seqs.to_s + ld
 246       log << "proteins with no passing domains            : " + proteins_with_failing_domains.to_s + ld
 247       if min_linker
 248         log << "min linker length                           : " + min_linker.to_s + ld
 249         log << "single domains                              : " + out_msa_singles.get_number_of_seqs.to_s + ld
 250         log << "domains in close pairs                      : " + (2 * out_msa_pairs.get_number_of_seqs).to_s + ld
 251         log << "isolated domains                            : " + out_msa_isolated.get_number_of_seqs.to_s + ld
 252         log << "proteins wih single domains                 : " + out_msa_single_domains_protein_seqs.get_number_of_seqs.to_s + ld
 253         log << "proteins wih close pair domains             : " + out_msa_close_pairs_protein_seqs.get_number_of_seqs.to_s + ld
 254         log << "proteins wih close pair domains only        : " + out_msa_close_pairs_only_protein_seqs.get_number_of_seqs.to_s + ld
 255         log << "proteins wih isolated domains               : " + out_msa_isolated_protein_seqs.get_number_of_seqs.to_s + ld
 256         log << "proteins wih isolated domains only          : " + out_msa_isolated_only_protein_seqs.get_number_of_seqs.to_s + ld
 257         log << "proteins wih close pair and isolated domains: " + out_msa_isolated_and_close_pair_protein_seqs.get_number_of_seqs.to_s + ld
 258       end
 259
 260       log << ld
 261
 262       return domain_pass_counter
 263
 264     end # parse
 265
 266
 267     private
 268
 269     def write_msa( msa, filename )
 270       io = MsaIO.new()
 271       w = FastaWriter.new()
 272       w.set_line_width( 60 )
 273       w.clean( true )
 274       begin
 275         io.write_to_file( msa, filename, w )
 276       rescue Exception
 277         error_msg = "could not write to \"" + filename + "\""
 278         raise IOError, error_msg
 279       end
 280     end
 281
 282
 283     def add_sequence( sequence_name, in_msa, add_to_msa )
 284       seqs = in_msa.find_by_name_start( sequence_name, true )
 285       if ( seqs.length < 1 )
 286         error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 287         raise StandardError, error_msg
 288       end
 289       if ( seqs.length > 1 )
 290         error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 291         raise StandardError, error_msg
 292       end
 293       seq = in_msa.get_sequence( seqs[ 0 ] )
 294       add_to_msa.add_sequence( seq )
 295     end
 296
 297     def process_hmmscan_datas( hmmscan_datas,
 298         in_msa,
 299         add_position,
 300         add_domain_number,
 301         add_species,
 302         out_msa,
 303         out_msa_singles,
 304         out_msa_pairs,
 305         out_msa_isolated,
 306         min_linker,
 307         out_msa_single_domains_protein_seqs,
 308         out_msa_close_pairs_protein_seqs,
 309         out_msa_close_pairs_only_protein_seqs,
 310         out_msa_isolated_protein_seqs,
 311         out_msa_isolated_only_protein_seqs,
 312         out_msa_isolated_and_close_pair_protein_seqs )
 313
 314       actual_out_of = hmmscan_datas.size
 315       saw_close_pair = false
 316       saw_isolated = false
 317
 318       seq_name = ""
 319       prev_seq_name = nil
 320
 321       hmmscan_datas.each_with_index do |hmmscan_data, index|
 322         if hmmscan_data.number < ( index + 1 )
 323           error_msg = "hmmscan_data.number < ( index + 1 ) (this should not have happened)"
 324           raise StandardError, error_msg
 325         end
 326
 327         seq_name =  hmmscan_data.seq_name
 328
 329         extract_domain( seq_name,
 330           index + 1,
 331           actual_out_of,
 332           hmmscan_data.env_from,
 333           hmmscan_data.env_to,
 334           in_msa,
 335           out_msa,
 336           add_position,
 337           add_domain_number,
 338           add_species )
 339
 340         if min_linker
 341           if actual_out_of == 1
 342             extract_domain( seq_name,
 343               1,
 344               1,
 345               hmmscan_data.env_from,
 346               hmmscan_data.env_to,
 347               in_msa,
 348               out_msa_singles,
 349               add_position,
 350               add_domain_number,
 351               add_species )
 352             if out_msa_single_domains_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 353               add_sequence( seq_name, in_msa, out_msa_single_domains_protein_seqs )
 354             else
 355               error_msg = "this should not have happened"
 356               raise StandardError, error_msg
 357             end
 358
 359           else
 360             first = index == 0
 361             last = index == hmmscan_datas.length - 1
 362
 363             if ( ( first && ( ( hmmscan_datas[ index + 1 ].env_from - hmmscan_data.env_to ) > min_linker) )  ||
 364                  ( last && ( ( hmmscan_data.env_from - hmmscan_datas[ index - 1 ].env_to ) > min_linker ) ) ||
 365                  ( !first && !last &&  ( ( hmmscan_datas[ index + 1 ].env_from - hmmscan_data.env_to ) > min_linker ) &&
 366                    ( ( hmmscan_data.env_from - hmmscan_datas[ index - 1 ].env_to ) > min_linker ) ) )
 367
 368               extract_domain( seq_name,
 369                 index + 1,
 370                 actual_out_of,
 371                 hmmscan_data.env_from,
 372                 hmmscan_data.env_to,
 373                 in_msa,
 374                 out_msa_isolated,
 375                 add_position,
 376                 add_domain_number,
 377                 add_species )
 378               saw_isolated = true
 379
 380             elsif !first
 381
 382               from = hmmscan_datas[ index - 1 ].env_from
 383               to = hmmscan_data.env_to
 384
 385               if ADD_TO_CLOSE_PAIRS > 0
 386                 from = from - ADD_TO_CLOSE_PAIRS
 387                 if from < 1
 388                   from = 1
 389                 end
 390                 to = to + ADD_TO_CLOSE_PAIRS
 391                 temp_seqs = in_msa.find_by_name_start( seq_name, true )
 392                 temp_seq = in_msa.get_sequence( temp_seqs[ 0 ] )
 393                 if to >  temp_seq.get_length
 394                   to =  temp_seq.get_length
 395                 end
 396               end
 397
 398               extract_domain( seq_name,
 399                 index.to_s  + "+" + ( index + 1 ).to_s,
 400                 actual_out_of,
 401                 from,
 402                 hmmscan_data.env_to,
 403                 in_msa,
 404                 out_msa_pairs,
 405                 add_position,
 406                 add_domain_number,
 407                 add_species )
 408               saw_close_pair = true
 409             end
 410           end
 411         end
 412         if prev_seq_name && prev_seq_name != seq_name
 413           error_msg = "this should not have happened"
 414           raise StandardError, error_msg
 415         end
 416         prev_seq_name = seq_name
 417       end # each
 418       if saw_isolated
 419         if out_msa_isolated_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 420           add_sequence( seq_name, in_msa, out_msa_isolated_protein_seqs )
 421         else
 422           error_msg = "this should not have happened"
 423           raise StandardError, error_msg
 424         end
 425       end
 426       if saw_close_pair
 427         if out_msa_close_pairs_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 428           add_sequence( seq_name, in_msa, out_msa_close_pairs_protein_seqs )
 429         else
 430           error_msg = "this should not have happened"
 431           raise StandardError, error_msg
 432         end
 433       end
 434       if saw_close_pair && saw_isolated
 435         if out_msa_isolated_and_close_pair_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 436           add_sequence( seq_name, in_msa, out_msa_isolated_and_close_pair_protein_seqs )
 437         else
 438           error_msg = "this should not have happened"
 439           raise StandardError, error_msg
 440         end
 441       elsif saw_close_pair
 442         if out_msa_close_pairs_only_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 443           add_sequence( seq_name, in_msa, out_msa_close_pairs_only_protein_seqs )
 444         else
 445           error_msg = "this should not have happened"
 446           raise StandardError, error_msg
 447         end
 448       elsif saw_isolated
 449         if out_msa_isolated_only_protein_seqs.find_by_name_start( seq_name, true ).length < 1
 450           add_sequence( seq_name, in_msa, out_msa_isolated_only_protein_seqs )
 451         else
 452           error_msg = "this should not have happened"
 453           raise StandardError, error_msg
 454         end
 455       end
 456     end # def process_hmmscan_data
 457
 458     def extract_domain( sequence,
 459         number,
 460         out_of,
 461         seq_from,
 462         seq_to,
 463         in_msa,
 464         out_msa,
 465         add_position,
 466         add_domain_number,
 467         add_species )
 468       if number.is_a?( Fixnum ) && ( number < 1 || out_of < 1 || number > out_of )
 469         error_msg = "number=" + number.to_s + ", out of=" + out_of.to_s
 470         raise StandardError, error_msg
 471       end
 472       if seq_from < 1 || seq_to < 1 || seq_from >= seq_to
 473         error_msg = "impossible: seq-from=" + seq_from.to_s + ", seq-to=" + seq_to.to_s
 474         raise StandardError, error_msg
 475       end
 476       seqs = in_msa.find_by_name_start( sequence, true )
 477       if seqs.length < 1
 478         error_msg = "sequence \"" + sequence + "\" not found in sequence file"
 479         raise IOError, error_msg
 480       end
 481       if seqs.length > 1
 482         error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
 483         raise IOError, error_msg
 484       end
 485       # hmmsearch is 1 based, wheres sequences are 0 bases in this package.
 486       seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 487
 488       orig_name = seq.get_name
 489
 490       seq.set_name( orig_name.split[ 0 ] )
 491
 492       if add_position
 493         seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
 494       end
 495
 496       if out_of != 1 && add_domain_number
 497         seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 498       end
 499
 500       if add_species
 501         a = orig_name.rindex "["
 502         b = orig_name.rindex "]"
 503         unless a && b
 504           error_msg = "species not found in " + orig_name
 505           raise StandardError, error_msg
 506         end
 507         species = orig_name[ a .. b ]
 508         seq.set_name( seq.get_name + " " + species )
 509       end
 510       out_msa.add_sequence( seq )
 511     end
 512
 513     def is_ignorable?( line )
 514       return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
 515     end
 516
 517   end # class HmmscanDomainExtractor
 518
 519   class HmmsearchData
 520     def initialize( seq_name, number, out_of, env_from, env_to, i_e_value )
 521       @seq_name = seq_name
 522       @number = number
 523       @out_of = out_of
 524       @env_from = env_from
 525       @env_to = env_to
 526       @i_e_value = i_e_value
 527     end
 528     attr_reader :seq_name, :number, :out_of, :env_from, :env_to, :i_e_value
 529   end
 530
 531 end # module Evoruby
 532