X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=forester%2Fruby%2Fevoruby%2Flib%2Fevo%2Ftool%2Fhmmscan_summary.rb;h=3b2996849c94fb28f4d88756b595d4a1ff3c4ff0;hb=8a4f118ddf875a7e7c17fbbd902116d5648ce805;hp=db82d4a8bebcde3e6f7253bf9a3809d4cff2347c;hpb=258424df68a8b78b20d2ccf1baedf7697c935c57;p=jalview.git diff --git a/forester/ruby/evoruby/lib/evo/tool/hmmscan_summary.rb b/forester/ruby/evoruby/lib/evo/tool/hmmscan_summary.rb index db82d4a..3b29968 100644 --- a/forester/ruby/evoruby/lib/evo/tool/hmmscan_summary.rb +++ b/forester/ruby/evoruby/lib/evo/tool/hmmscan_summary.rb @@ -6,7 +6,6 @@ # # $Id: hmmscan_parser.rb,v 1.5 2010/12/13 19:00:11 cmzmasek Exp $ # -# last modified: 121003 require 'set' @@ -21,14 +20,15 @@ module Evoruby class HmmscanSummary PRG_NAME = "hsp" - PRG_VERSION = "2.000" + PRG_VERSION = "2.001" PRG_DESC = "hmmscan summary" - PRG_DATE = "2012.10.23" - COPYRIGHT = "2012 Christian M Zmasek" - CONTACT = "phylosoft@gmail.com" - WWW = "www.phylosoft.org" + PRG_DATE = "2013.10.23" + COPYRIGHT = "2013 Christian M Zmasek" + CONTACT = "phyloxml@gmail.com" + WWW = "https://sites.google.com/site/cmzmasek/home/software/forester" DELIMITER_OPTION = "d" + SPECIES_OPTION = "s" I_E_VALUE_THRESHOLD_OPTION = "ie" FS_E_VALUE_THRESHOLD_OPTION = "pe" HMM_FOR_PROTEIN_OUTPUT = "m" @@ -48,8 +48,6 @@ module Evoruby def run - - Util.print_program_information( PRG_NAME, PRG_VERSION, PRG_DESC, @@ -84,6 +82,7 @@ module Evoruby 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 ) @@ -116,6 +115,8 @@ module Evoruby end end + + fs_e_value_threshold = -1.0 if ( cla.is_option_set?( FS_E_VALUE_THRESHOLD_OPTION ) ) begin @@ -146,6 +147,15 @@ module Evoruby 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 @@ -159,6 +169,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 @@ -201,13 +212,13 @@ module Evoruby parse_descriptions, fs_e_value_threshold, hmm_for_protein_output, - uniprot ) + uniprot, + species ) rescue IOError => e Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT ) end domain_counts = get_domain_counts() - puts puts( "domain counts (considering potential i-E-value threshold and ignoring of DUFs):" ) puts( "(number of different domains: " + domain_counts.length.to_s + ")" ) @@ -230,9 +241,8 @@ module Evoruby get_descriptions, fs_e_value_threshold, hmm_for_protein_output, - uniprot ) - - + uniprot, + species ) Util.check_file_for_readability( inpath ) Util.check_file_for_writability( outpath ) @@ -251,8 +261,6 @@ module Evoruby hmmscan_results_per_protein = [] - - prev_query = "" results.each do | r | @@ -288,7 +296,8 @@ module Evoruby fs_e_value_threshold, hmm_for_protein_output, i_e_value_threshold, - true ) + false, + species ) end hmmscan_results_per_protein.clear end @@ -303,17 +312,18 @@ module Evoruby end end 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, - true ) + uniprot, + species ) end outfile.flush() outfile.close() - end # def parse def process_id( id ) @@ -323,8 +333,6 @@ module Evoruby id end - - def count_model( model ) if ( @domain_counts.has_key?( model ) ) count = @domain_counts[ model ].to_i @@ -339,7 +347,8 @@ module Evoruby fs_e_value_threshold, hmm_for_protein_output, i_e_value_threshold, - uniprotkb ) + uniprotkb, + species ) dc = 0 # filter according to i-Evalue threshold @@ -347,12 +356,14 @@ module Evoruby hmmscan_results_per_protein_filtered = [] hmmscan_results_per_protein.each do | r | + + if r.model == hmm_for_protein_output - if r.fs_e_value > fs_e_value_threshold + if i_e_value_threshold > 0.0 && r.fs_e_value > fs_e_value_threshold return end end - if r.i_e_value <= i_e_value_threshold + if i_e_value_threshold <= 0 || r.i_e_value <= i_e_value_threshold hmmscan_results_per_protein_filtered << r if r.model == hmm_for_protein_output dc += 1 @@ -376,7 +387,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" @@ -385,37 +396,19 @@ module Evoruby s << r.model + " " end s << "\t" - e = UniprotKB::get_entry_by_id( process_id( own.query ) ) - # if e != nil && e.de != nil - # e.de.each do |i| - # - # end - # else - # s << "-" - # end - s << "\t" - if e != nil && e.dr != nil - e.dr.each do | dr | - if dr != nil - if dr =~ /PDB;\s+([A-Z0-9]{4});/ - s << $1 - - end - end - end - else - s << "-" - end - s << "\t" - - + if !uniprotkb.empty? + #e = UniprotKB::get_entry_by_id( process_id( own.query ) ) + #if e != nil + # s << uniprot_annotation( e ) + # # s << "\uniprot_annotationt" + #end + end - s << "\t" 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" @@ -439,6 +432,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 = "" @@ -494,7 +499,6 @@ module Evoruby end - def print_help() puts( "Usage:" ) puts() @@ -506,6 +510,7 @@ module Evoruby puts( " -" + IGNORE_DUF_OPTION + ": ignore DUFs" ) puts( " -" + FS_E_VALUE_THRESHOLD_OPTION + ": E-value threshold for full protein sequences, only for protein summary" ) puts( " -" + HMM_FOR_PROTEIN_OUTPUT + ": HMM for protein summary" ) + puts( " -" + SPECIES_OPTION + ": species for protein summary" ) puts() end