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

   1 #
   2 # = lib/evo/io/parser/hmmsearch_domain_extractor.rb - HmmsearchDomainExtractor class
   3 #
   4 # Copyright::  Copyright (C) 2006-2008 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: hmmsearch_domain_extractor.rb,v 1.24 2009/11/25 06:30:24 cmzmasek Exp $
   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 HmmsearchDomainExtractor
  20
  21         TRIM_BY = 2
  22
  23         def initialize
  24         end
  25
  26         # raises ArgumentError, IOError, StandardError
  27         def parse( hmmsearch_output,
  28                 fasta_sequence_file,
  29                 outfile,
  30                 passed_seqs_outfile,
  31                 failed_seqs_outfile,
  32                 e_value_threshold,
  33                 length_threshold,
  34                 add_position,
  35                 add_domain_number,
  36                 add_domain_number_as_digit,
  37                 add_domain_number_as_letter,
  38                 trim_name,
  39                 log )
  40
  41             Util.check_file_for_readability( hmmsearch_output )
  42             Util.check_file_for_readability( fasta_sequence_file )
  43             Util.check_file_for_writability( outfile )
  44             Util.check_file_for_writability( passed_seqs_outfile )
  45             Util.check_file_for_writability( failed_seqs_outfile )
  46
  47             in_msa = nil
  48             factory = MsaFactory.new()
  49             in_msa = factory.create_msa_from_file( fasta_sequence_file, FastaParser.new() )
  50
  51             if ( in_msa == nil || in_msa.get_number_of_seqs() < 1 )
  52                 error_msg = "could not find fasta sequences in " + fasta_sequence_file
  53                 raise IOError, error_msg
  54             end
  55
  56             out_msa = Msa.new
  57             failed_seqs = Msa.new
  58             passed_seqs = Msa.new
  59
  60             ld = Constants::LINE_DELIMITER
  61
  62             domain_pass_counter     = 0
  63             domain_fail_counter     = 0
  64             proteins_with_passing_domains = 0
  65             proteins_with_failing_domains = 0
  66             max_domain_copy_number_per_protein = -1
  67             max_domain_copy_number_sequence    = ""
  68
  69
  70             File.open( hmmsearch_output ) do | file |
  71                 while line = file.gets
  72                     if !is_ignorable?( line ) && line =~ /^\S+\s+/
  73
  74                         #         tn      acc     tlen    query   acc     qlen    Evalue  score   bias    #       of      c-E     i-E     score   bias    hf      ht      af      at      ef      et      acc     desc
  75                         #         1       2       3       4       5       6       7       8       9       10      11      12      13      14      15      16      17      18      19      20      21      22      23
  76                         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+(.*)/
  77
  78                         sequence = $1
  79                         number   = $10.to_i
  80                         out_of   = $11.to_i
  81                         env_from = $20.to_i
  82                         env_to   = $21.to_i
  83                         i_e_value  = $13.to_f
  84                         if ( number > max_domain_copy_number_per_protein )
  85                             max_domain_copy_number_sequence    = sequence
  86                             max_domain_copy_number_per_protein = number
  87                         end
  88                         if ( ( ( e_value_threshold.to_f < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
  89                                  ( ( length_threshold.to_f <= 0 )   || ( env_to - env_from + 1 ) >= length_threshold.to_f )  )
  90
  91                             HmmsearchDomainExtractor.extract_domain( sequence,
  92                                 number,
  93                                 out_of,
  94                                 env_from,
  95                                 env_to,
  96                                 in_msa,
  97                                 out_msa,
  98                                 add_position,
  99                                 add_domain_number,
 100                                 add_domain_number_as_digit,
 101                                 add_domain_number_as_letter,
 102                                 trim_name )
 103
 104                             domain_pass_counter += 1
 105
 106                             if !passed_seqs.has?( sequence, true, false )
 107                                 HmmsearchDomainExtractor.add_sequence( sequence, in_msa, passed_seqs )
 108                                 proteins_with_passing_domains += 1
 109                             end
 110                         else
 111                             print( domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)" )
 112                             log << domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)"
 113                             if ( ( e_value_threshold.to_f >= 0.0 ) && ( i_e_value > e_value_threshold ) )
 114                                 print( " iE=" + i_e_value.to_s )
 115                                 log << " iE=" + i_e_value.to_s
 116                             end
 117                             if ( ( length_threshold.to_f > 0 ) && ( env_to - env_from + 1 ) < length_threshold.to_f )
 118                                 le = env_to - env_from + 1
 119                                 print( " l=" + le.to_s )
 120                                 log << " l=" + le.to_s
 121                             end
 122                             print( Constants::LINE_DELIMITER )
 123                             log << Constants::LINE_DELIMITER
 124                             domain_fail_counter  += 1
 125
 126                             if !failed_seqs.has?( sequence, true, false )
 127                                 HmmsearchDomainExtractor.add_sequence( sequence, in_msa, failed_seqs )
 128                                 proteins_with_failing_domains += 1
 129                             end
 130                         end
 131                     end
 132                 end
 133             end
 134
 135             if domain_pass_counter < 1
 136                 error_msg = "no domain sequences were extracted"
 137                 raise StandardError, error_msg
 138             end
 139
 140             log << Constants::LINE_DELIMITER
 141             puts( "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s )
 142             log << "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s
 143             log << Constants::LINE_DELIMITER
 144
 145             if ( max_domain_copy_number_per_protein > 1 )
 146                 puts( "First protein with this copy number: " + max_domain_copy_number_sequence )
 147                 log << "First protein with this copy number: " + max_domain_copy_number_sequence
 148                 log << Constants::LINE_DELIMITER
 149             end
 150
 151             io = MsaIO.new()
 152             w = FastaWriter.new()
 153             w.set_line_width( 60 )
 154             w.clean( true )
 155
 156             begin
 157                 io.write_to_file( out_msa, outfile, w )
 158             rescue Exception
 159                 error_msg = "could not write to \"" + outfile + "\""
 160                 raise IOError, error_msg
 161             end
 162
 163             begin
 164                 io.write_to_file( passed_seqs, passed_seqs_outfile, w )
 165             rescue Exception
 166                 error_msg = "could not write to \"" + passed_seqs_outfile + "\""
 167                 raise IOError, error_msg
 168             end
 169
 170             begin
 171                 io.write_to_file( failed_seqs, failed_seqs_outfile, w )
 172             rescue Exception
 173                 error_msg = "could not write to \"" + failed_seqs_outfile + "\""
 174                 raise IOError, error_msg
 175             end
 176
 177             log << ld
 178             log << "passing domains              : " + domain_pass_counter.to_s + ld
 179             log << "failing domains              : " + domain_fail_counter.to_s + ld
 180             log << "proteins with passing domains: " + proteins_with_passing_domains.to_s + ld
 181             log << "proteins with failing domains: " + proteins_with_failing_domains.to_s + ld
 182             log << ld
 183
 184             return domain_pass_counter
 185
 186         end # parse
 187
 188         private
 189
 190         def HmmsearchDomainExtractor.add_sequence( sequence_name, in_msa, add_to_msa )
 191             seqs = in_msa.find_by_name( sequence_name, true, false )
 192             if ( seqs.length < 1 )
 193                 error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 194                 raise StandardError, error_msg
 195             end
 196             if ( seqs.length > 1 )
 197                 error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 198                 raise StandardError, error_msg
 199             end
 200             seq = in_msa.get_sequence( seqs[ 0 ] )
 201             add_to_msa.add_sequence( seq )
 202         end
 203
 204         # raises ArgumentError, StandardError
 205         def HmmsearchDomainExtractor.extract_domain( sequence,
 206                 number,
 207                 out_of,
 208                 seq_from,
 209                 seq_to,
 210                 in_msa,
 211                 out_msa,
 212                 add_position,
 213                 add_domain_number,
 214                 add_domain_number_as_digit,
 215                 add_domain_number_as_letter,
 216                 trim_name )
 217             if ( number < 1 || out_of < 1 || number > out_of )
 218                 error_msg = "impossible: number=" + number.to_s + ", out of=" + out_of.to_s
 219                 raise ArgumentError, error_msg
 220             end
 221             if ( seq_from < 1 || seq_to < 1 || seq_from >= seq_to )
 222                 error_msg = "impossible: seq-f=" + seq_from.to_s + ", seq-t=" + seq_to.to_s
 223                 raise ArgumentError, error_msg
 224             end
 225             seqs = in_msa.find_by_name( sequence, true, false )
 226             if ( seqs.length < 1 )
 227                 error_msg = "sequence \"" + sequence + "\" not found in sequence file"
 228                 raise StandardError, error_msg
 229             end
 230             if ( seqs.length > 1 )
 231                 error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
 232                 raise StandardError, error_msg
 233             end
 234             # hmmsearch is 1 based, wheres sequences are 0 bases in this package.
 235             seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 236             if ( add_position )
 237                 seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
 238             end
 239
 240             if ( trim_name )
 241                 seq.set_name( seq.get_name[ 0, seq.get_name.length - TRIM_BY ] )
 242             end
 243
 244             if ( out_of != 1 )
 245                 if ( add_domain_number_as_digit )
 246                     seq.set_name( seq.get_name + number.to_s )
 247                 elsif ( add_domain_number_as_letter )
 248                     if number > 25
 249                         error_msg = 'too many identical domains per sequence, cannot use letters to distinguish them'
 250                         raise StandardError, error_msg
 251                     end
 252                     seq.set_name( seq.get_name + ( number + 96 ).chr )
 253                 elsif ( add_domain_number )
 254                     seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 255                 end
 256             end
 257
 258             if ( seq.get_name.length > 10 )
 259                 error_msg = "sequence name [" + seq.get_name + "] is longer than 10 characters"
 260                 raise StandardError, error_msg
 261             end
 262
 263             out_msa.add_sequence( seq )
 264         end
 265
 266
 267
 268
 269
 270         def is_ignorable?( line )
 271             return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
 272         end
 273
 274     end # class HmmsearchDomainExtractor
 275
 276 end # module Evoruby
 277