inprogress
[jalview.git] / forester / ruby / evoruby / lib / evo / tool / hmmscan_summary.rb
index 9c5e3a9..92bceb2 100644 (file)
@@ -1,17 +1,20 @@
 #
-# = lib/evo/apps/hmmscan_parser.rb - HmmscanParser class
+# = lib/evo/tool/hmmscan_summary.rb - HmmscanSummary class
 #
-# Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
+# Copyright::  Copyright (C) 2012 Christian M. Zmasek
 # License::    GNU Lesser General Public License (LGPL)
 #
 # $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $
 #
-# last modified: 11/24/2009
+# last modified: 121003
+
+require 'set'
 
 require 'lib/evo/util/constants'
 require 'lib/evo/util/util'
 require 'lib/evo/util/command_line_arguments'
 require 'lib/evo/io/parser/hmmscan_parser'
+require 'lib/evo/io/web/uniprotkb'
 
 module Evoruby
 
@@ -19,18 +22,20 @@ module Evoruby
 
     PRG_NAME       = "hsp"
     PRG_VERSION    = "2.000"
-    PRG_DESC       = "hmmscan parser"
-    PRG_DATE       = "2012.10.19"
+    PRG_DESC       = "hmmscan summary"
+    PRG_DATE       = "2012.10.23"
     COPYRIGHT      = "2012 Christian M Zmasek"
     CONTACT        = "phylosoft@gmail.com"
     WWW            = "www.phylosoft.org"
 
     DELIMITER_OPTION              = "d"
+    SPECIES_OPTION                = "s"
     I_E_VALUE_THRESHOLD_OPTION    = "ie"
     FS_E_VALUE_THRESHOLD_OPTION   = "pe"
     HMM_FOR_PROTEIN_OUTPUT        = "m"
     IGNORE_DUF_OPTION             = "i"
     PARSE_OUT_DESCRIPITION_OPTION = "a"
+    UNIPROT                       = "u"
     HELP_OPTION_1                 = "help"
     HELP_OPTION_2                 = "h"
 
@@ -44,6 +49,8 @@ module Evoruby
 
     def run
 
+
+
       Util.print_program_information( PRG_NAME,
         PRG_VERSION,
         PRG_DESC,
@@ -60,7 +67,7 @@ module Evoruby
       end
 
       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
-           cla.is_option_set?( HELP_OPTION_2 ) )
+            cla.is_option_set?( HELP_OPTION_2 ) )
         print_help
         exit( 0 )
       end
@@ -77,6 +84,8 @@ module Evoruby
       allowed_opts.push( IGNORE_DUF_OPTION )
       allowed_opts.push( PARSE_OUT_DESCRIPITION_OPTION )
       allowed_opts.push( HMM_FOR_PROTEIN_OUTPUT )
+      allowed_opts.push( UNIPROT )
+      allowed_opts.push( SPECIES_OPTION )
 
       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
       if ( disallowed.length > 0 )
@@ -109,6 +118,8 @@ module Evoruby
         end
       end
 
+      
+
       fs_e_value_threshold = -1.0
       if ( cla.is_option_set?( FS_E_VALUE_THRESHOLD_OPTION ) )
         begin
@@ -130,6 +141,24 @@ module Evoruby
         end
       end
 
+      uniprot = ""
+      if ( cla.is_option_set?( UNIPROT ) )
+        begin
+          uniprot = cla.get_option_value( UNIPROT )
+        rescue ArgumentError => e
+          Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
+        end
+      end
+
+      species = "HUMAN"
+       if ( cla.is_option_set?( SPECIES_OPTION ) )
+        begin
+          species = cla.get_option_value( SPECIES_OPTION )
+        rescue ArgumentError => e
+          Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
+        end
+      end
+
       ignore_dufs = false
       if ( cla.is_option_set?( IGNORE_DUF_OPTION ) )
         ignore_dufs = true
@@ -143,6 +172,7 @@ module Evoruby
       puts()
       puts( "hmmpfam outputfile  : " + inpath )
       puts( "outputfile          : " + outpath )
+      puts( "species             : " + species )
       if ( i_e_value_threshold >= 0.0 )
         puts( "i-E-value threshold : " + i_e_value_threshold.to_s )
       else
@@ -163,14 +193,17 @@ module Evoruby
       else
         puts( "column delimiter     : " + column_delimiter )
       end
-      if ( fs_e_value_threshold >= 0.0 )
+      if fs_e_value_threshold >= 0.0
         puts( "E-value threshold   : " + fs_e_value_threshold.to_s )
       else
         puts( "E-value threshold   : no threshold" )
       end
-      if ( !hmm_for_protein_output.empty? )
+      if !hmm_for_protein_output.empty?
         puts( "HMM for proteins    : " + hmm_for_protein_output )
       end
+      if !uniprot.empty?
+        puts( "Uniprot             : " + uniprot )
+      end
       puts()
 
       begin
@@ -181,8 +214,10 @@ module Evoruby
           ignore_dufs,
           parse_descriptions,
           fs_e_value_threshold,
-          hmm_for_protein_output )
-      rescue ArgumentError, IOError => e
+          hmm_for_protein_output,
+          uniprot,
+          species )
+      rescue IOError => e
         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
       end
       domain_counts = get_domain_counts()
@@ -209,10 +244,18 @@ module Evoruby
         ignore_dufs,
         get_descriptions,
         fs_e_value_threshold,
-        hmm_for_protein_output )
+        hmm_for_protein_output,
+        uniprot,
+        species )
+
+
+
       Util.check_file_for_readability( inpath )
       Util.check_file_for_writability( outpath )
 
+      hmmscan_parser = HmmscanParser.new( inpath )
+      results = hmmscan_parser.parse
+
       outfile = File.open( outpath, "a" )
 
       query     = ""
@@ -224,11 +267,11 @@ module Evoruby
 
       hmmscan_results_per_protein = []
 
-      hmmscan_parser = HmmscanParser.new( inpath )
+
 
       prev_query = ""
 
-      hmmscan_parser.parse.each do | r |
+      results.each do | r |
         model     = r.model
         query     = r.query
         i_e_value = r.i_e_value
@@ -236,21 +279,21 @@ module Evoruby
         env_to    = r.env_to
 
         if ( ( i_e_value_threshold < 0.0 ) || ( i_e_value <= i_e_value_threshold ) ) &&
-           ( !ignore_dufs || ( model !~ /^DUF\d+/ ) )
+            ( !ignore_dufs || ( model !~ /^DUF\d+/ ) )
           count_model( model )
           outfile.print( query +
-             column_delimiter )
+              column_delimiter )
           if ( get_descriptions )
             outfile.print( desc +
-               column_delimiter )
+                column_delimiter )
           end
           outfile.print( model +
-             column_delimiter +
-             env_from.to_s +
-             column_delimiter +
-             env_to.to_s +
-             column_delimiter +
-             i_e_value.to_s )
+              column_delimiter +
+              env_from.to_s +
+              column_delimiter +
+              env_to.to_s +
+              column_delimiter +
+              i_e_value.to_s )
           outfile.print( Constants::LINE_DELIMITER )
         end
 
@@ -260,7 +303,9 @@ module Evoruby
               process_hmmscan_results_per_protein( hmmscan_results_per_protein,
                 fs_e_value_threshold,
                 hmm_for_protein_output,
-                i_e_value_threshold )
+                i_e_value_threshold,
+                false,
+                species )
             end
             hmmscan_results_per_protein.clear
           end
@@ -275,19 +320,30 @@ module Evoruby
           end
         end
       end
-      if !hmm_for_protein_output.empty?
-        if !hmmscan_results_per_protein.empty?
-          process_hmmscan_results_per_protein( hmmscan_results_per_protein,
-            fs_e_value_threshold,
-            hmm_for_protein_output,
-            i_e_value_threshold )
-        end
+    
+      if !hmm_for_protein_output.empty? && !hmmscan_results_per_protein.empty?
+        process_hmmscan_results_per_protein( hmmscan_results_per_protein,
+          fs_e_value_threshold,
+          hmm_for_protein_output,
+          i_e_value_threshold,
+          false,
+          species )
       end
+
       outfile.flush()
       outfile.close()
 
     end # def parse
 
+    def process_id( id )
+      if  id =~ /(sp|tr)\|\S+\|(\S+)/
+        id = $2
+      end
+      id
+    end
+
+
+
     def count_model( model )
       if ( @domain_counts.has_key?( model ) )
         count = @domain_counts[ model ].to_i
@@ -301,7 +357,11 @@ module Evoruby
     def process_hmmscan_results_per_protein( hmmscan_results_per_protein,
         fs_e_value_threshold,
         hmm_for_protein_output,
-        i_e_value_threshold )
+        i_e_value_threshold,
+        uniprotkb,
+        species      )
+      
+    
 
       dc = 0
       # filter according to i-Evalue threshold
@@ -309,6 +369,9 @@ module Evoruby
       hmmscan_results_per_protein_filtered = []
 
       hmmscan_results_per_protein.each do | r |
+        
+        puts r
+        
         if r.model == hmm_for_protein_output
           if r.fs_e_value > fs_e_value_threshold
             return
@@ -338,7 +401,7 @@ module Evoruby
 
       s = ""
       s << own.query + "\t"
-      s << "HUMAN" + "\t"
+      s << species + "\t"
       s << own.fs_e_value.to_s + "\t"
       s << own.qlen.to_s + "\t"
       s << dc.to_s + "\t"
@@ -347,10 +410,17 @@ module Evoruby
         s << r.model + " "
       end
       s << "\t"
+      puts s
+      #e = UniprotKB::get_entry_by_id( process_id( own.query ) )
+
+      #if e != nil
+      #  s << uniprot_annotation( e )
+      # # s << "\uniprot_annotationt"
+      #end
 
       overview = make_overview( hmmscan_results_per_protein_filtered, hmm_for_protein_output )
 
-      s << overview   + "\t"
+      s << overview  + "\t"
 
       s << calc_linkers(  hmmscan_results_per_protein_filtered, hmm_for_protein_output )  + "\t"
 
@@ -374,6 +444,18 @@ module Evoruby
       puts s
     end
 
+    def uniprot_annotation( e )
+      s = ""
+      pdb_ids = e.get_pdb_ids
+      if !pdb_ids.empty?
+        pdb_ids.each do | pdb |
+          s << pdb << ", "
+        end
+      else
+        s << "-"
+      end
+      s
+    end
 
     def calc_linkers(  hmmscan_results_per_protein_filtered, hmm_for_protein_output )
       linkers = ""