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                         sequence = $1
  80                         number   = $10.to_i
  81                         out_of   = $11.to_i
  82                         env_from = $20.to_i
  83                         env_to   = $21.to_i
  84                         i_e_value  = $13.to_f
  85                         if ( number > max_domain_copy_number_per_protein )
  86                             max_domain_copy_number_sequence    = sequence
  87                             max_domain_copy_number_per_protein = number
  88                         end
  89                         if ( ( ( e_value_threshold.to_f < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
  90                                  ( ( length_threshold.to_f <= 0 )   || ( env_to - env_from + 1 ) >= length_threshold.to_f )  )
  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                             domain_pass_counter += 1
 104                             count_species( sequence, passed_species_counts )
 105                             if !passed_seqs.has?( sequence, true, false )
 106                                 HmmsearchDomainExtractor.add_sequence( sequence, in_msa, passed_seqs )
 107                                 proteins_with_passing_domains += 1
 108                             end
 109                         else
 110                             print( domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)" )
 111                             log << domain_fail_counter.to_s + ": " + sequence.to_s + " did not meet threshold(s)"
 112                             if ( ( e_value_threshold.to_f >= 0.0 ) && ( i_e_value > e_value_threshold ) )
 113                                 print( " iE=" + i_e_value.to_s )
 114                                 log << " iE=" + i_e_value.to_s
 115                             end
 116                             if ( ( length_threshold.to_f > 0 ) && ( env_to - env_from + 1 ) < length_threshold.to_f )
 117                                 le = env_to - env_from + 1
 118                                 print( " l=" + le.to_s )
 119                                 log << " l=" + le.to_s
 120                             end
 121                             print( Constants::LINE_DELIMITER )
 122                             log << Constants::LINE_DELIMITER
 123                             domain_fail_counter  += 1
 124                             count_species( sequence, failed_species_counts )
 125                             if !failed_seqs.has?( sequence, true, false )
 126                                 HmmsearchDomainExtractor.add_sequence( sequence, in_msa, failed_seqs )
 127                                 proteins_with_failing_domains += 1
 128                             end
 129                         end
 130                     end
 131                 end
 132             end
 133
 134             if domain_pass_counter < 1
 135                 error_msg = "no domain sequences were extracted"
 136                 raise StandardError, error_msg
 137             end
 138
 139             log << Constants::LINE_DELIMITER
 140             puts( "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s )
 141             log << "Max domain copy number per protein : " + max_domain_copy_number_per_protein.to_s
 142             log << Constants::LINE_DELIMITER
 143
 144             if ( max_domain_copy_number_per_protein > 1 )
 145                 puts( "First protein with this copy number: " + max_domain_copy_number_sequence )
 146                 log << "First protein with this copy number: " + max_domain_copy_number_sequence
 147                 log << Constants::LINE_DELIMITER
 148             end
 149
 150             io = MsaIO.new()
 151             w = FastaWriter.new()
 152             w.set_line_width( 60 )
 153             w.clean( true )
 154
 155             begin
 156                 io.write_to_file( out_msa, outfile, w )
 157             rescue Exception
 158                 error_msg = "could not write to \"" + outfile + "\""
 159                 raise IOError, error_msg
 160             end
 161
 162             begin
 163                 io.write_to_file( passed_seqs, passed_seqs_outfile, w )
 164             rescue Exception
 165                 error_msg = "could not write to \"" + passed_seqs_outfile + "\""
 166                 raise IOError, error_msg
 167             end
 168
 169             begin
 170                 io.write_to_file( failed_seqs, failed_seqs_outfile, w )
 171             rescue Exception
 172                 error_msg = "could not write to \"" + failed_seqs_outfile + "\""
 173                 raise IOError, error_msg
 174             end
 175
 176             log << ld
 177             log << "passing domains              : " + domain_pass_counter.to_s + ld
 178             log << "failing domains              : " + domain_fail_counter.to_s + ld
 179             log << "proteins with passing domains: " + proteins_with_passing_domains.to_s + ld
 180             log << "proteins with failing domains: " + proteins_with_failing_domains.to_s + ld
 181             log << ld
 182             log << 'passing domains counts per species: ' << ld
 183             passed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 184             log << ld
 185             log << 'failing domains counts per species: ' << ld
 186             failed_species_counts.each_pair { | species, count | log << "#{species}: #{count}" << ld }
 187             log << ld
 188             return domain_pass_counter
 189
 190         end # parse
 191
 192         private
 193
 194
 195         def HmmsearchDomainExtractor.add_sequence( sequence_name, in_msa, add_to_msa )
 196             seqs = in_msa.find_by_name( sequence_name, true, false )
 197             if ( seqs.length < 1 )
 198                 error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
 199                 raise StandardError, error_msg
 200             end
 201             if ( seqs.length > 1 )
 202                 error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
 203                 raise StandardError, error_msg
 204             end
 205             seq = in_msa.get_sequence( seqs[ 0 ] )
 206             add_to_msa.add_sequence( seq )
 207         end
 208
 209         # raises ArgumentError, StandardError
 210         def HmmsearchDomainExtractor.extract_domain( sequence,
 211                 number,
 212                 out_of,
 213                 seq_from,
 214                 seq_to,
 215                 in_msa,
 216                 out_msa,
 217                 add_position,
 218                 add_domain_number,
 219                 add_domain_number_as_digit,
 220                 add_domain_number_as_letter,
 221                 trim_name )
 222             if ( number < 1 || out_of < 1 || number > out_of )
 223                 error_msg = "impossible: number=" + number.to_s + ", out of=" + out_of.to_s
 224                 raise ArgumentError, error_msg
 225             end
 226             if ( seq_from < 1 || seq_to < 1 || seq_from >= seq_to )
 227                 error_msg = "impossible: seq-f=" + seq_from.to_s + ", seq-t=" + seq_to.to_s
 228                 raise ArgumentError, error_msg
 229             end
 230             seqs = in_msa.find_by_name( sequence, true, false )
 231             if ( seqs.length < 1 )
 232                 error_msg = "sequence \"" + sequence + "\" not found in sequence file"
 233                 raise StandardError, error_msg
 234             end
 235             if ( seqs.length > 1 )
 236                 error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
 237                 raise StandardError, error_msg
 238             end
 239             # hmmsearch is 1 based, wheres sequences are 0 bases in this package.
 240             seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
 241             if ( add_position )
 242                 seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
 243             end
 244
 245             if ( trim_name )
 246                 seq.set_name( seq.get_name[ 0, seq.get_name.length - TRIM_BY ] )
 247             end
 248
 249             if ( out_of != 1 )
 250                 if ( add_domain_number_as_digit )
 251                     seq.set_name( seq.get_name + number.to_s )
 252                 elsif ( add_domain_number_as_letter )
 253                     if number > 25
 254                         error_msg = 'too many identical domains per sequence, cannot use letters to distinguish them'
 255                         raise StandardError, error_msg
 256                     end
 257                     seq.set_name( seq.get_name + ( number + 96 ).chr )
 258                 elsif ( add_domain_number )
 259                     seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
 260                 end
 261             end
 262
 263             if ( seq.get_name.length > 10 )
 264                 error_msg = "sequence name [" + seq.get_name + "] is longer than 10 characters"
 265                 raise StandardError, error_msg
 266             end
 267
 268             out_msa.add_sequence( seq )
 269         end
 270
 271         def count_species( sequence, species_counts_map )
 272             species = get_species( sequence )
 273             if species != nil
 274                 if !species_counts_map.has_key?( species )
 275                     species_counts_map[ species ] = 1
 276                 else
 277                     species_counts_map[ species ] = species_counts_map[ species ] + 1
 278                 end
 279             end
 280         end
 281
 282         def get_species( sequence_name )
 283             if sequence_name =~ /^.+_(.+)$/
 284                 return $1
 285             else
 286                 return nil
 287             end
 288         end
 289
 290         def is_ignorable?( line )
 291             return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
 292         end
 293
 294     end # class HmmsearchDomainExtractor
 295
 296 end # module Evoruby
 297