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     # raises ArgumentError, IOError, StandardError
  27     def parse( domain_id,
  28         hmmsearch_output,
  29         fasta_sequence_file,
  30         outfile,
  31         passed_seqs_outfile,
  32         failed_seqs_outfile,
  33         e_value_threshold,
  34         length_threshold,
  35         add_position,
  36         add_domain_number,
  37         add_domain_number_as_digit,
  38         add_domain_number_as_letter,
  39         trim_name,
  40         add_species,
  41         log )
  42
  43       Util.check_file_for_readability( hmmsearch_output )
  44       Util.check_file_for_readability( fasta_sequence_file )
  45       Util.check_file_for_writability( outfile )
  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       failed_seqs = Msa.new
  60       passed_seqs = Msa.new
  61
  62       ld = Constants::LINE_DELIMITER
  63
  64       domain_pass_counter     = 0
  65       domain_fail_counter     = 0
  66       proteins_with_passing_domains = 0
  67       proteins_with_failing_domains = 0
  68       max_domain_copy_number_per_protein = -1
  69       max_domain_copy_number_sequence    = ''
  70       failed_species_counts         = Hash.new
  71       passed_species_counts         = Hash.new
  72
  73       File.open( hmmsearch_output ) do | file |
  74         while line = file.gets
  75           if !is_ignorable?( line ) && line =~ /^\S+\s+/
  76
  77             #         tn      acc     tlen    query   acc     qlen    Evalue  score   bias    #       of      c-E     i-E     score   bias    hf      ht      af      at      ef      et      acc     desc
  78             #         1       2       3       4       5       6       7       8       9       10      11      12      13      14      15      16      17      18      19      20      21      22      23
  79             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+(.*)/
  80
  81             target_name = $1
  82             if domain_id != target_name
  83               next
  84             end
  85
  86
  87             sequence = $4
  88             number   = $10.to_i
  89             out_of   = $11.to_i
  90             env_from = $20.to_i
  91             env_to   = $21.to_i
  92             i_e_value  = $13.to_f
  93             if ( number > max_domain_copy_number_per_protein )
  94               max_domain_copy_number_sequence    = sequence
  95               max_domain_copy_number_per_protein = number
  96             end
  97             if ( ( ( e_value_threshold.to_f < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
  98                  ( ( length_threshold.to_f <= 0 )   || ( env_to - env_from + 1 ) >= length_threshold.to_f )  )
  99               extract_domain( sequence,
 100                 number,
 101                 out_of,
 102                 env_from,
 103                 env_to,
 104                 in_msa,
 105                 out_msa,
 106                 add_position,
 107                 add_domain_number,
 108                 add_domain_number_as_digit,
 109                 add_domain_number_as_letter,
 110                 trim_name ,
 111                 add_species )
 112               domain_pass_counter += 1
 113               count_species( sequence, passed_species_counts )
 114               if passed_seqs.find_by_name_start( sequence, true ).length < 1
 115                 add_sequence( sequence, in_msa, passed_seqs )
 116                 proteins_with_passing_domains += 1
 117               end
 118             else
 119               print( domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)" )
 120               log << domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)"
 121               if ( ( e_value_threshold.to_f >= 0.0 ) && ( i_e_value > e_value_threshold ) )
 122                 print( " iE=" + i_e_value.to_s )
 123                 log << " iE=" + i_e_value.to_s
 124               end
 125               if ( ( length_threshold.to_f > 0 ) && ( env_to - env_from + 1 ) < length_threshold.to_f )
 126                 le = env_to - env_from + 1
 127                 print( " l=" + le.to_s )
 128                 log << " l=" + le.to_s
 129               end
 130               print( Constants::LINE_DELIMITER )
 131               log << Constants::LINE_DELIMITER
 132               domain_fail_counter  += 1
 133               count_species( sequence, failed_species_counts )
 134               if failed_seqs.find_by_name_start( sequence, true ).length < 1
 135                 add_sequence( sequence, in_msa, failed_seqs )
 136                 proteins_with_failing_domains += 1
 137               end
 138             end
 139           end
 140         end
 141       end
 142
 143       if domain_pass_counter < 1
 144         error_msg = "no domain sequences were extracted"
 145         raise StandardError, error_msg
 146       end
 147
 148       log << Constants::LINE_DELIMITER
 149       puts( "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s )
 150       log << "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s
 151       log << Constants::LINE_DELIMITER
 152
 153       if ( max_domain_copy_number_per_protein > 1 )
 154         puts( "First protein with this copy number: " + max_domain_copy_number_sequence )
 155         log << "First protein with this copy number: " + max_domain_copy_number_sequence
 156         log << Constants::LINE_DELIMITER
 157       end
 158
 159       io = MsaIO.new()
 160       w = FastaWriter.new()
 161       w.set_line_width( 60 )
 162       w.clean( true )
 163
 164       begin
 165         io.write_to_file( out_msa, outfile, w )
 166       rescue Exception
 167         error_msg = "could not write to \"" + outfile + "\""
 168         raise IOError, error_msg
 169       end
 170
 171       begin
 172         io.write_to_file( passed_seqs, passed_seqs_outfile, w )
 173       rescue Exception
 174         error_msg = "could not write to \"" + passed_seqs_outfile + "\""
 175         raise IOError, error_msg
 176       end
 177
 178       begin
 179         io.write_to_file( failed_seqs, failed_seqs_outfile, w )
 180       rescue Exception
 181         error_msg = "could not write to \"" + failed_seqs_outfile + "\""
 182         raise IOError, error_msg
 183       end
 184
 185       log << ld
 186       log << "passing domains              : " + domain_pass_counter.to_s + ld
 187       log << "failing domains              : " + domain_fail_counter.to_s + ld
 188       log << "proteins with passing domains: " + proteins_with_passing_domains.to_s + ld
 189       log << "proteins with failing domains: " + proteins_with_failing_domains.to_s + ld
 190       log << ld
 191       log << 'passing domains counts per species: ' << ld
 192       passed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 193       log << ld
 194       log << 'failing domains counts per species: ' << ld
 195       failed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 196       log << ld
 197       return domain_pass_counter
 198
 199     end # parse
 200
 201     private
 202
 203
 204     def add_sequence( sequence_name, in_msa, add_to_msa )
 205       seqs = in_msa.find_by_name_start( sequence_name, true )
 206       if ( seqs.length < 1 )
 207         error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 208         raise StandardError, error_msg
 209       end
 210       if ( seqs.length > 1 )
 211         error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 212         raise StandardError, error_msg
 213       end
 214       seq = in_msa.get_sequence( seqs[ 0 ] )
 215       add_to_msa.add_sequence( seq )
 216     end
 217
 218     # raises ArgumentError, StandardError
 219     def extract_domain( sequence,
 220         number,
 221         out_of,
 222         seq_from,
 223         seq_to,
 224         in_msa,
 225         out_msa,
 226         add_position,
 227         add_domain_number,
 228         add_domain_number_as_digit,
 229         add_domain_number_as_letter,
 230         trim_name,
 231         add_species )
 232       if ( number < 1 || out_of < 1 || number > out_of )
 233         error_msg = "impossible: number=" + number.to_s + ", out of=" + out_of.to_s
 234         raise ArgumentError, error_msg
 235       end
 236       if ( seq_from < 1 || seq_to < 1 || seq_from >= seq_to )
 237         error_msg = "impossible: seq-f=" + seq_from.to_s + ", seq-t=" + seq_to.to_s
 238         raise ArgumentError, error_msg
 239       end
 240       seqs = in_msa.find_by_name_start( sequence, true )
 241       if seqs.length < 1
 242         error_msg = "sequence \"" + sequence + "\" not found in sequence file"
 243         raise StandardError, error_msg
 244       end
 245       if seqs.length > 1
 246         error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
 247         raise StandardError, error_msg
 248       end
 249       # hmmsearch is 1 based, wheres sequences are 0 bases in this package.
 250       seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 251
 252       orig_name = seq.get_name
 253
 254       seq.set_name( orig_name.split[ 0 ] )
 255
 256       if add_position
 257         seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
 258       end
 259
 260       if trim_name
 261         seq.set_name( seq.get_name[ 0, seq.get_name.length - TRIM_BY ] )
 262       end
 263
 264       if out_of != 1
 265         if ( add_domain_number_as_digit )
 266           seq.set_name( seq.get_name + number.to_s )
 267         elsif ( add_domain_number_as_letter )
 268           if number > 25
 269             error_msg = 'too many identical domains per sequence, cannot use letters to distinguish them'
 270             raise StandardError, error_msg
 271           end
 272           seq.set_name( seq.get_name + ( number + 96 ).chr )
 273         elsif ( add_domain_number )
 274           seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 275         end
 276
 277         if add_species
 278           a = orig_name.rindex "["
 279           b = orig_name.rindex "]"
 280           unless a && b
 281             error_msg = "species not found in " + orig_name
 282             raise StandardError, error_msg
 283           end
 284           species = orig_name[ a .. b ]
 285           seq.set_name( seq.get_name + " [" + species + "]" )
 286         end
 287
 288       end
 289
 290       # if ( seq.get_name.length > 10 )
 291       #   error_msg = "sequence name [" + seq.get_name + "] is longer than 10 characters"
 292       #   raise StandardError, error_msg
 293       # end
 294
 295       out_msa.add_sequence( seq )
 296     end
 297
 298     def count_species( sequence, species_counts_map )
 299       species = get_species( sequence )
 300       if species != nil
 301         if !species_counts_map.has_key?( species )
 302           species_counts_map[ species ] = 1
 303         else
 304           species_counts_map[ species ] = species_counts_map[ species ] + 1
 305         end
 306       end
 307     end
 308
 309     def get_species( sequence_name )
 310       if sequence_name =~ /^.+_(.+)$/
 311         return $1
 312       else
 313         return nil
 314       end
 315     end
 316
 317     def is_ignorable?( line )
 318       return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
 319     end
 320
 321   end # class HmmscanDomainExtractor
 322
 323 end # module Evoruby
 324