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