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