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