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             failed_species_counts         = Hash.new
  69             passed_species_counts         = Hash.new
  70
  71             File.open( hmmsearch_output ) do | file |
  72                 while line = file.gets
  73                     if !is_ignorable?( line ) && line =~ /^\S+\s+/
  74
  75                         #         tn      acc     tlen    query   acc     qlen    Evalue  score   bias    #       of      c-E     i-E     score   bias    hf      ht      af      at      ef      et      acc     desc
  76                         #         1       2       3       4       5       6       7       8       9       10      11      12      13      14      15      16      17      18      19      20      21      22      23
  77                         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+(.*)/
  78
  79                         # line =~ /^(\S+)\s+(\d+)\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+)/
  80                         sequence = $1
  81                         number   = $10.to_i
  82                         out_of   = $11.to_i
  83                         env_from = $20.to_i
  84                         env_to   = $21.to_i
  85                         i_e_value  = $13.to_f
  86                         if ( number > max_domain_copy_number_per_protein )
  87                             max_domain_copy_number_sequence    = sequence
  88                             max_domain_copy_number_per_protein = number
  89                         end
  90                         if ( ( ( e_value_threshold.to_f < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
  91                                  ( ( length_threshold.to_f <= 0 )   || ( env_to - env_from + 1 ) >= length_threshold.to_f )  )
  92                             HmmsearchDomainExtractor.extract_domain( sequence,
  93                                 number,
  94                                 out_of,
  95                                 env_from,
  96                                 env_to,
  97                                 in_msa,
  98                                 out_msa,
  99                                 add_position,
 100                                 add_domain_number,
 101                                 add_domain_number_as_digit,
 102                                 add_domain_number_as_letter,
 103                                 trim_name )
 104                             domain_pass_counter += 1
 105                             count_species( sequence, passed_species_counts )
 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                             count_species( sequence, failed_species_counts )
 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             log << 'passing domains counts per species: ' << ld
 184             passed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 185             log << ld
 186             log << 'failing domains counts per species: ' << ld
 187             failed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 188             log << ld
 189             return domain_pass_counter
 190
 191         end # parse
 192
 193         private
 194
 195
 196         def HmmsearchDomainExtractor.add_sequence( sequence_name, in_msa, add_to_msa )
 197             seqs = in_msa.find_by_name( sequence_name, true, false )
 198             if ( seqs.length < 1 )
 199                 error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 200                 raise StandardError, error_msg
 201             end
 202             if ( seqs.length > 1 )
 203                 error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 204                 raise StandardError, error_msg
 205             end
 206             seq = in_msa.get_sequence( seqs[ 0 ] )
 207             add_to_msa.add_sequence( seq )
 208         end
 209
 210         # raises ArgumentError, StandardError
 211         def HmmsearchDomainExtractor.extract_domain( sequence,
 212                 number,
 213                 out_of,
 214                 seq_from,
 215                 seq_to,
 216                 in_msa,
 217                 out_msa,
 218                 add_position,
 219                 add_domain_number,
 220                 add_domain_number_as_digit,
 221                 add_domain_number_as_letter,
 222                 trim_name )
 223             if ( number < 1 || out_of < 1 || number > out_of )
 224                 error_msg = "impossible: number=" + number.to_s + ", out of=" + out_of.to_s
 225                 raise ArgumentError, error_msg
 226             end
 227             if ( seq_from < 1 || seq_to < 1 || seq_from >= seq_to )
 228                 error_msg = "impossible: seq-f=" + seq_from.to_s + ", seq-t=" + seq_to.to_s
 229                 raise ArgumentError, error_msg
 230             end
 231             seqs = in_msa.find_by_name( sequence, true, false )
 232             if ( seqs.length < 1 )
 233                 error_msg = "sequence \"" + sequence + "\" not found in sequence file"
 234                 raise StandardError, error_msg
 235             end
 236             if ( seqs.length > 1 )
 237                 error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
 238                 raise StandardError, error_msg
 239             end
 240             # hmmsearch is 1 based, wheres sequences are 0 bases in this package.
 241             seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 242             if ( add_position )
 243                 seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
 244             end
 245
 246             if ( trim_name )
 247                 seq.set_name( seq.get_name[ 0, seq.get_name.length - TRIM_BY ] )
 248             end
 249
 250             if ( out_of != 1 )
 251                 if ( add_domain_number_as_digit )
 252                     seq.set_name( seq.get_name + number.to_s )
 253                 elsif ( add_domain_number_as_letter )
 254                     if number > 25
 255                         error_msg = 'too many identical domains per sequence, cannot use letters to distinguish them'
 256                         raise StandardError, error_msg
 257                     end
 258                     seq.set_name( seq.get_name + ( number + 96 ).chr )
 259                 elsif ( add_domain_number )
 260                     seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 261                 end
 262             end
 263
 264             if ( seq.get_name.length > 10 )
 265                 error_msg = "sequence name [" + seq.get_name + "] is longer than 10 characters"
 266                 raise StandardError, error_msg
 267             end
 268
 269             out_msa.add_sequence( seq )
 270         end
 271
 272         def count_species( sequence, species_counts_map )
 273             species = get_species( sequence )
 274             if species != nil
 275                 if !species_counts_map.has_key?( species )
 276                     species_counts_map[ species ] = 1
 277                 else
 278                     species_counts_map[ species ] = species_counts_map[ species ] + 1
 279                 end
 280             end
 281         end
 282
 283         def get_species( sequence_name )
 284             if sequence_name =~ /^.+_(.+)$/
 285                 return $1
 286             else
 287                 return nil
 288             end
 289         end
 290
 291         def is_ignorable?( line )
 292             return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
 293         end
 294
 295     end # class HmmsearchDomainExtractor
 296
 297 end # module Evoruby
 298