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

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