forester/ruby/evoruby/lib/evo/tool/hmmscan_analysis.rb

   1 #
   2 # = lib/evo/tool/hmmscan_summary.rb - HmmscanSummary class
   3 #
   4 # Copyright::  Copyright (C) 2012 Christian M. Zmasek
   5 # License::    GNU Lesser General Public License (LGPL)
   6 #
   7 # $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $
   8 #
   9
  10 require 'lib/evo/util/constants'
  11 require 'lib/evo/util/util'
  12 require 'lib/evo/util/command_line_arguments'
  13 require 'lib/evo/io/parser/hmmscan_parser'
  14 require 'lib/evo/msa/msa'
  15 require 'lib/evo/msa/msa_factory'
  16 require 'lib/evo/io/msa_io'
  17 require 'lib/evo/io/parser/fasta_parser'
  18
  19 module Evoruby
  20
  21   class HmmscanAnalysis
  22
  23     PRG_NAME       = "hsa"
  24     PRG_VERSION    = "1.001"
  25     PRG_DESC       = "hmmscan analysis"
  26     PRG_DATE       = "140127"
  27     COPYRIGHT      = "2014 Christian Zmasek"
  28     CONTACT        = "phyloxml@gmail.com"
  29     WWW            = "https://sites.google.com/site/cmzmasek/home/software/forester"
  30
  31     I_E_VALUE_THRESHOLD_OPTION    = "ie"
  32     FS_E_VALUE_THRESHOLD_OPTION   = "pe"
  33     TARGET_MODELS                 = "m"
  34     EXTRACTION                    = "x"
  35     HELP_OPTION_1                 = "help"
  36     HELP_OPTION_2                 = "h"
  37
  38     USE_AVOID_HMMS = true
  39     AVOID_HHMS = [ "RRM_1", "RRM_2", "RRM_3", "RRM_4", "RRM_5", "RRM_6" ]
  40     LIMIT_FOR_CLOSE_DOMAINS = 20
  41
  42
  43     def run
  44
  45       Util.print_program_information( PRG_NAME,
  46         PRG_VERSION,
  47         PRG_DESC,
  48         PRG_DATE,
  49         COPYRIGHT,
  50         CONTACT,
  51         WWW,
  52         STDOUT )
  53
  54       begin
  55         cla = CommandLineArguments.new( ARGV )
  56       rescue ArgumentError => e
  57         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
  58       end
  59
  60       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
  61            cla.is_option_set?( HELP_OPTION_2 ) )
  62         print_help
  63         exit( 0 )
  64       end
  65
  66       if ( cla.get_number_of_files != 2 && cla.get_number_of_files != 3 )
  67         print_help
  68         exit( -1 )
  69       end
  70
  71       allowed_opts = Array.new
  72       allowed_opts.push( I_E_VALUE_THRESHOLD_OPTION )
  73       allowed_opts.push( FS_E_VALUE_THRESHOLD_OPTION )
  74       allowed_opts.push( TARGET_MODELS )
  75       allowed_opts.push( EXTRACTION )
  76
  77       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
  78       if ( disallowed.length > 0 )
  79         Util.fatal_error( PRG_NAME,
  80           "unknown option(s): " + disallowed,
  81           STDOUT )
  82       end
  83
  84       inpath = cla.get_file_name( 0 )
  85       Util.check_file_for_readability( inpath )
  86
  87       seq_file_path = cla.get_file_name( 1 )
  88       Util.check_file_for_readability(  seq_file_path  )
  89
  90       extraction_output = nil
  91       if ( cla.get_number_of_files == 3 )
  92         extraction_output = cla.get_file_name( 2 )
  93         if File.exist?( extraction_output )
  94           Util.fatal_error( PRG_NAME, "error: [#{extraction_output}] already exists" )
  95         end
  96
  97       end
  98
  99       msa = nil
 100       if seq_file_path != nil
 101         msa = read_fasta_file( seq_file_path )
 102       end
 103
 104       i_e_value_threshold = -1
 105       if ( cla.is_option_set?( I_E_VALUE_THRESHOLD_OPTION ) )
 106         begin
 107           i_e_value_threshold = cla.get_option_value_as_float( I_E_VALUE_THRESHOLD_OPTION )
 108         rescue ArgumentError => e
 109           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 110         end
 111         if ( i_e_value_threshold < 0.0 )
 112           Util.fatal_error( PRG_NAME, "attempt to use a negative i-E-value threshold", STDOUT )
 113         end
 114       end
 115
 116       fs_e_value_threshold = -1
 117       if ( cla.is_option_set?( FS_E_VALUE_THRESHOLD_OPTION ) )
 118         begin
 119           fs_e_value_threshold = cla.get_option_value_as_float( FS_E_VALUE_THRESHOLD_OPTION )
 120         rescue ArgumentError => e
 121           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 122         end
 123         if ( fs_e_value_threshold < 0.0 )
 124           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
 125         end
 126       end
 127
 128       target_models = []
 129       if ( cla.is_option_set?( TARGET_MODELS ) )
 130         begin
 131           hmm_for_protein_output = cla.get_option_value( TARGET_MODELS )
 132           target_models = hmm_for_protein_output.split( "/" )
 133         rescue ArgumentError => e
 134           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 135         end
 136       end
 137
 138       x_models = []
 139       if ( cla.is_option_set?( EXTRACTION ) )
 140         begin
 141           x_models = cla.get_option_value( EXTRACTION ).split( "~" )
 142           puts x_models
 143         rescue ArgumentError => e
 144           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 145         end
 146       end
 147
 148
 149       if i_e_value_threshold > 0
 150         puts "i E-value threshold:" + "\t" + i_e_value_threshold.to_s
 151       else
 152         puts "i E-value threshold:" + "\t" + "n/a"
 153       end
 154
 155       if fs_e_value_threshold > 0
 156         puts "fs E-value threshold:" + "\t" + fs_e_value_threshold.to_s
 157       else
 158         puts "fs E-value threshold:" + "\t" + "n/a"
 159       end
 160
 161       puts
 162
 163       title = "ID" + "\t"
 164       title << "SPECIES" + "\t"
 165       target_models.each do | t |
 166         title << t + "\t"
 167       end
 168       title << "LENGTH" + "\t"
 169       title << "#DOMAINS" + "\t"
 170       title << "DOMAINS" + "\t"
 171       title << "DA" + "\t"
 172       title << "DETAILED DA"
 173       puts title
 174
 175       begin
 176         parse( inpath,
 177           i_e_value_threshold,
 178           fs_e_value_threshold,
 179           target_models,
 180           x_models,
 181           msa,
 182           extraction_output )
 183       rescue IOError => e
 184         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
 185       end
 186
 187     end # def run
 188
 189
 190     private
 191
 192     def read_fasta_file( input )
 193       f = MsaFactory.new()
 194       msa = nil
 195       begin
 196         msa = f.create_msa_from_file( input, FastaParser.new() )
 197       rescue Exception => e
 198         Util.fatal_error( PRG_NAME, "error: " + e.to_s )
 199       end
 200       msa
 201     end
 202
 203
 204     # raises ArgumentError, IOError
 205     def parse( inpath,
 206         i_e_value_threshold,
 207         fs_e_value_threshold,
 208         target_models,
 209         x_models,
 210         msa,
 211         extraction_output )
 212
 213
 214
 215       extraction_output_file = nil
 216       if extraction_output != nil
 217         extraction_output_file = File.open( extraction_output, "a" )
 218       end
 219
 220       hmmscan_parser = HmmscanParser.new( inpath )
 221
 222       results = hmmscan_parser.parse
 223
 224       hmmscan_results_per_protein = []
 225
 226       query     = ""
 227       prev_query = ""
 228       results.each do | r |
 229         query = r.query
 230         if !prev_query.empty? && prev_query != query
 231           if !hmmscan_results_per_protein.empty?
 232             process_hmmscan_results_per_protein( hmmscan_results_per_protein,
 233               fs_e_value_threshold,
 234               target_models,
 235               x_models,
 236               i_e_value_threshold,
 237               msa,
 238               extraction_output_file )
 239           end
 240           hmmscan_results_per_protein.clear
 241         end
 242         prev_query = query
 243
 244         if USE_AVOID_HMMS
 245           if !AVOID_HHMS.include? r.model
 246             hmmscan_results_per_protein << r
 247           end
 248         else
 249           hmmscan_results_per_protein << r
 250         end
 251
 252       end
 253
 254       if !hmmscan_results_per_protein.empty?
 255         process_hmmscan_results_per_protein( hmmscan_results_per_protein,
 256           fs_e_value_threshold,
 257           target_models,
 258           x_models,
 259           i_e_value_threshold,
 260           msa,
 261           extraction_output_file )
 262       end
 263       if extraction_output_file != nil
 264         extraction_output_file.close
 265       end
 266     end # def parse
 267
 268
 269     def process_hmmscan_results_per_protein( hmmscan_results_per_protein,
 270         fs_e_value_threshold,
 271         target_hmms,
 272         x_models,
 273         i_e_value_threshold,
 274         msa,
 275         extraction_output_file )
 276
 277       raise StandardError, "target hmms is empty" if target_hmms.length < 1
 278       raise StandardError, "results is empty" if hmmscan_results_per_protein.length < 1
 279       raise StandardError, "msa is empty" if ( msa == nil || msa.get_number_of_seqs < 1 )
 280
 281
 282       # filter according to i-Evalue threshold
 283       # abort if fs Evalue too high
 284
 285       if fs_e_value_threshold >= 0
 286         hmmscan_results_per_protein.each do | r |
 287           target_hmms.each do | hmm |
 288             if r.model == hmm && r.fs_e_value > fs_e_value_threshold
 289               return
 290             end
 291           end
 292         end
 293       end
 294
 295       hmmscan_results_per_protein_filtered = []
 296       matched = Set.new
 297
 298       hmmscan_results_per_protein.each do | r |
 299         if i_e_value_threshold < 0 || r.i_e_value <= i_e_value_threshold
 300           hmmscan_results_per_protein_filtered << r
 301           target_hmms.each do | hmm |
 302             if r.model == hmm
 303               matched << hmm
 304               break
 305             end
 306           end
 307
 308         end
 309       end
 310
 311       if matched.length < target_hmms.length
 312         return
 313       end
 314       if  hmmscan_results_per_protein_filtered.length < 1
 315         return
 316       end
 317
 318       hmmscan_results_per_protein_filtered.sort! { |r1,r2| r1.env_from <=> r2.env_from }
 319
 320       owns = []
 321       target_hmms.each do | hmm |
 322         hmmscan_results_per_protein_filtered.each do | r |
 323           if r.model == hmm
 324             owns << r
 325             break
 326           end
 327         end
 328       end
 329
 330       query = nil
 331       qlen = nil
 332       owns.each do | own |
 333         if query == nil
 334           query = own.query
 335           qlen = own.qlen
 336         else
 337           raise StandardError, "failed sanity check" if query != own.query || qlen != own.qlen
 338           raise StandardError, "failed sanity check: qlen != own.qlen" if qlen != own.qlen
 339         end
 340       end
 341
 342       species = nil
 343       ll = nil
 344       seq = get_sequence( query, msa  )
 345       species = get_species( seq )
 346       raise StandardError, "could not get species" if species == nil || species.empty?
 347       if x_models != nil &&  x_models.length > 0
 348         ll = extract_linker( hmmscan_results_per_protein_filtered, x_models, seq,  extraction_output_file )
 349       end
 350
 351       s = query + "\t"
 352       s << species + "\t"
 353       owns.each do | own |
 354         s << own.fs_e_value.to_s + "\t"
 355       end
 356
 357       s << qlen.to_s + "\t"
 358
 359       s << hmmscan_results_per_protein_filtered.length.to_s + "\t"
 360       hmmscan_results_per_protein_filtered.each do | r |
 361         s << r.model + " "
 362       end
 363       s << "\t"
 364       if msa != nil
 365         if ll != nil
 366           s << ll.to_s
 367         end
 368         s << "\t"
 369       end
 370       s << make_overview_da( hmmscan_results_per_protein_filtered )
 371       s << "\t"
 372       s << make_detailed_da( hmmscan_results_per_protein_filtered, qlen )
 373       puts s
 374
 375     end
 376
 377
 378     def extract_linker( hmmscan_results_per_protein_filtered, x_models, seq, extraction_output_file )
 379       raise StandardError, "extraction output file is nil" if extraction_output_file == nil
 380       prev_r = nil
 381       hmmscan_results_per_protein_filtered.each do | r |
 382         if  prev_r != nil
 383           if ( x_models.length == 2 && prev_r.model == x_models[ 0 ] && r.model == x_models[ 1 ] )
 384             ll = output_linker( prev_r.env_to, r.env_from, seq, extraction_output_file )
 385             return ll
 386           end
 387         end
 388         prev_r = r
 389       end
 390       return nil
 391     end
 392
 393
 394     def get_sequence( query, msa  )
 395       seq = nil
 396       indices = msa.find_by_name_start( query , true )
 397       if indices.length != 1
 398         if query[ -1, 1 ] == "|"
 399           query.chop!
 400         end
 401         seq = msa.get_by_name_pattern( /\b#{Regexp.quote(query)}\b/ )
 402       else
 403         seq = msa.get_sequence( indices[ 0 ] )
 404       end
 405       seq
 406     end
 407
 408
 409     def get_species( seq )
 410       species = nil
 411       if seq.get_name =~ /\[([A-Z0-9]{3,5})\]/
 412         species = $1
 413       end
 414       species
 415     end
 416
 417
 418     def output_linker( first, last , seq,  output_file )
 419       if ( last - first >= 1 )
 420         output_file.print( ">" + seq.get_name + " [" +  first.to_s + "-" + last.to_s +  "]" + "\n")
 421         output_file.print( seq.get_subsequence( first - 1 , last - 1 ).get_sequence_as_string + "\n" )
 422       end
 423       last - first + 1
 424
 425     end
 426
 427
 428     def make_detailed_da( hmmscan_results_per_protein_filtered , qlen )
 429       s = ""
 430       prev_r = nil
 431       hmmscan_results_per_protein_filtered.each do | r |
 432         if  prev_r != nil
 433           s << make_interdomain_sequence( r.env_from - prev_r.env_to - 1 )
 434         else
 435           s << make_interdomain_sequence( r.env_from, false )
 436         end
 437         s << r.model
 438         s << "["
 439         s << r.env_from.to_s << "-" << r.env_to.to_s
 440         s << " " << r.i_e_value.to_s
 441         s << "]"
 442         prev_r = r
 443       end
 444       s << make_interdomain_sequence( qlen - prev_r.env_to, false )
 445       s
 446     end
 447
 448
 449     def make_overview_da( hmmscan_results_per_protein_filtered )
 450       overview = ""
 451       prev_r = nil
 452       hmmscan_results_per_protein_filtered.each do | r |
 453         if prev_r == nil
 454           overview << r.model
 455         else
 456           if  ( r.env_from - prev_r.env_to - 1 ) <= LIMIT_FOR_CLOSE_DOMAINS
 457             overview << "~" << r.model
 458           else
 459             overview << "----" << r.model
 460           end
 461         end
 462         prev_r = r
 463
 464       end
 465       overview
 466     end
 467
 468
 469     def make_interdomain_sequence( d, mark_short = true )
 470       s = ""
 471       d /= 20
 472       if d >= 10
 473         s << "----//----"
 474       elsif d >= 1
 475         d.times do
 476           s << "-"
 477         end
 478       elsif mark_short
 479         s << "~"
 480       end
 481       s
 482     end
 483
 484
 485     def print_help()
 486       puts "Usage:"
 487       puts
 488       puts "  " + PRG_NAME + ".rb [options] <hmmscan outputfile> <sequences in fasta format> [outputfile]"
 489       puts
 490       puts "  options: -" + I_E_VALUE_THRESHOLD_OPTION  + ": i-E-value threshold, default is no threshold"
 491       puts "           -" + FS_E_VALUE_THRESHOLD_OPTION  + ": E-value threshold for full protein sequences, only for protein summary"
 492       puts "           -" + TARGET_MODELS + ": target HMMs"
 493       puts "           -" + EXTRACTION + ": to extract matching sequences to [outputfile]"
 494       puts
 495     end
 496
 497   end # class
 498
 499 end # module Evoruby