in progress
[jalview.git] / forester / ruby / evoruby / lib / evo / tool / domains_to_forester.rb
1 #
2 # = lib/evo/apps/domains_to_forester - DomainsToForester class
3 #
4 # Copyright::  Copyright (C) 2006-2007 Christian M. Zmasek
5 # License::    GNU Lesser General Public License (LGPL)
6 #
7 # $Id: Exp $
8 #
9 # last modified: 06/11/2007
10
11 require 'lib/evo/util/constants'
12 require 'lib/evo/util/util'
13 require 'lib/evo/util/command_line_arguments'
14 require 'lib/evo/msa/msa_factory'
15 require 'lib/evo/io/parser/fasta_parser'
16 require 'lib/evo/sequence/protein_domain'
17 require 'lib/evo/sequence/domain_structure'
18
19 module Evoruby
20
21   class DomainsToForester
22
23     PRG_NAME       = "d2f"
24     PRG_DESC       = "parsed hmmpfam output to forester format"
25     PRG_VERSION    = "1.001"
26     PRG_DATE       = "20120807"
27     COPYRIGHT      = "2012 Christian M Zmasek"
28     CONTACT        = "phylosoft@gmail.com"
29     WWW            = "www.phylosoft.org"
30
31     E_VALUE_THRESHOLD_OPTION         = "e"
32     OVERWRITE_IF_SAME_FROM_TO_OPTION = "o"
33     HELP_OPTION_1                    = "help"
34     HELP_OPTION_2                    = "h"
35
36     def parse( domains_list_file,
37         original_seqs_file,
38         outfile,
39         column_delimiter,
40         e_value_threshold,
41         overwrite_if_same_from_to )
42       Util.check_file_for_readability( domains_list_file )
43       Util.check_file_for_readability( original_seqs_file )
44       Util.check_file_for_writability( outfile )
45
46       domain_structures = Hash.new() # protein name is key, domain structure is value
47
48       f = MsaFactory.new
49
50       original_seqs = f.create_msa_from_file( original_seqs_file, FastaParser.new )
51       if ( original_seqs.get_number_of_seqs < 1 )
52         error_msg = "\"" + original_seqs_file + "\" appears devoid of sequences in fasta-format"
53         raise ArgumentError, error_msg
54       end
55
56       File.open( domains_list_file ) do | file |
57         while line = file.gets
58           if !is_ignorable?( line )
59             
60             a = line.split( column_delimiter )
61             l = a.length
62             if ( ( l < 4 ) || ( e_value_threshold >= 0.0 && l < 5 ) )
63               error_msg = "unexpected format at line: " + line
64               raise IOError, error_msg
65             end
66             protein_name = a[ 0 ]
67             domain_name  = a[ 1 ]
68             seq_from     = -1
69             seq_to       = -1
70             ##########################################
71             if domain_name =~ /RRM_\d/
72               puts "ignoring " + line 
73               next
74             end
75             ##########################################
76             
77             
78             begin
79               seq_from = a[ 2 ].to_i
80             rescue Exception
81               error_msg = "failed to parse seq from from \"" + a[ 2 ] + "\" [line: " + line + "]"
82               raise IOError, error_msg
83             end
84             begin
85               seq_to = a[ 3 ].to_i
86             rescue Exception
87               error_msg = "failed to parse seq to from \"" + a[ 3 ] + "\" [line: " + line + "]"
88               raise IOError, error_msg
89             end
90
91             e_value = -1
92             if l > 4
93               begin
94                 e_value = a[ 4 ].to_f
95               rescue Exception
96                 error_msg = "failed to parse E-value from \"" + a[ 4 ] + "\" [line: " + line + "]"
97                 raise IOError, error_msg
98               end
99             end
100
101             seq = original_seqs.get_by_name_start( protein_name )
102
103             total_length = seq.get_length
104
105             if ( ( ( e_value_threshold < 0.0 ) || ( e_value <= e_value_threshold ) )  )
106               pd = ProteinDomain.new( domain_name, seq_from, seq_to, "", e_value )
107               ds = nil
108               if ( domain_structures.has_key?( protein_name ) )
109                 ds = domain_structures[ protein_name ]
110               else
111                 ds = DomainStructure.new( total_length )
112                 domain_structures[ protein_name ] = ds
113               end
114               ds.add_domain( pd, overwrite_if_same_from_to )
115             end
116
117           end
118         end
119       end
120
121       out = File.open( outfile, "a" )
122       ds = domain_structures.sort
123       for d in ds
124         protein_name     = d[ 0 ]
125         domain_structure = d[ 1 ]
126         out.print( protein_name.to_s )
127         out.print( "\t" )
128         out.print( domain_structure.to_NHX )
129         out.print( Constants::LINE_DELIMITER  )
130       end
131
132       out.flush()
133       out.close()
134
135     end # parse
136
137
138
139
140     def run()
141
142       Util.print_program_information( PRG_NAME,
143         PRG_VERSION,
144         PRG_DESC,
145         PRG_DATE,
146         COPYRIGHT,
147         CONTACT,
148         WWW,
149         STDOUT )
150
151       begin
152         cla = CommandLineArguments.new( ARGV )
153       rescue ArgumentError => e
154         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
155       end
156
157       if ( cla.is_option_set?( HELP_OPTION_1 ) ||
158            cla.is_option_set?( HELP_OPTION_2 ) )
159         print_help
160         exit( 0 )
161       end
162
163       if cla.get_number_of_files != 3
164         print_help
165         exit( -1 )
166       end
167
168       allowed_opts = Array.new
169       allowed_opts.push( E_VALUE_THRESHOLD_OPTION )
170       allowed_opts.push( OVERWRITE_IF_SAME_FROM_TO_OPTION )
171
172       disallowed = cla.validate_allowed_options_as_str( allowed_opts )
173       if ( disallowed.length > 0 )
174         Util.fatal_error( PRG_NAME,
175           "unknown option(s): " + disallowed,
176           STDOUT )
177       end
178
179       domains_list_file       = cla.get_file_name( 0 )
180       original_sequences_file = cla.get_file_name( 1 )
181       outfile                 = cla.get_file_name( 2 )
182
183
184       e_value_threshold = -1.0
185       if cla.is_option_set?( E_VALUE_THRESHOLD_OPTION )
186         begin
187           e_value_threshold = cla.get_option_value_as_float( E_VALUE_THRESHOLD_OPTION )
188         rescue ArgumentError => e
189           Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
190         end
191         if ( e_value_threshold < 0.0 )
192           Util.fatal_error( PRG_NAME, "attempt to use a negative E-value threshold", STDOUT )
193         end
194       end
195       overwrite_if_same_from_to = false
196       if ( cla.is_option_set?( OVERWRITE_IF_SAME_FROM_TO_OPTION ) )
197         overwrite_if_same_from_to = true
198       end
199
200       puts
201       puts( "Domains list file                      : " + domains_list_file )
202       puts( "Fasta sequencefile (complete sequences): " + original_sequences_file )
203       puts( "Outputfile                             : " + outfile )
204       if ( e_value_threshold >= 0.0 )
205         puts( "E-value threshold                      : " + e_value_threshold.to_s )
206       else
207         puts( "E-value threshold                      : no threshold" )
208       end
209       if ( overwrite_if_same_from_to )
210         puts( "Overwrite if same from and to          : true" )
211       else
212         puts( "Overwrite if same from and to          : false" )
213       end
214
215       puts
216
217       begin
218         parse( domains_list_file,
219           original_sequences_file,
220           outfile,
221           " ",
222           e_value_threshold,
223           overwrite_if_same_from_to )
224
225       rescue ArgumentError, IOError, StandardError => e
226         Util.fatal_error( PRG_NAME, "error: " + e.to_s, STDOUT )
227       rescue Exception => e
228         Util.fatal_error( PRG_NAME, "unexpected exception: " + e.to_s, STDOUT )
229       end
230
231
232       puts
233       Util.print_message( PRG_NAME, 'OK' )
234       puts
235
236     end
237
238     private
239
240     def print_help()
241       puts
242       puts( "Usage:" )
243       puts
244       puts( "  " + PRG_NAME + ".rb [options] <domains list file (parsed hmmpfam output)> <file containing complete sequences in fasta format> <outputfile>" )
245       puts()
246       puts( "  options: -" + E_VALUE_THRESHOLD_OPTION  + "=<f> : E-value threshold, default is no threshold" )
247       puts( "               -" + OVERWRITE_IF_SAME_FROM_TO_OPTION  + " : overwrite domain with same start and end with domain with better E-value" )
248       puts
249     end
250
251
252
253     def is_ignorable?( line )
254       return ( line !~ /[A-Za-z0-9-]/ || line =~ /^\s*#/)
255     end
256
257
258   end # class DomainsToForester
259
260
261 end # module Evoruby