in progress...
[jalview.git] / forester / ruby / evoruby / lib / evo / io / parser / hmmscan_domain_extractor.rb
1 #
2 # = lib/evo/io/parser/hmmscan_domain_extractor.rb - HmmscanDomainExtractor class
3 #
4 # Copyright::    Copyright (C) 2017 Christian M. Zmasek
5 # License::      GNU Lesser General Public License (LGPL)
6 #
7 # Last modified: 2017/02/16
8
9 require 'lib/evo/util/constants'
10 require 'lib/evo/msa/msa_factory'
11 require 'lib/evo/io/msa_io'
12 require 'lib/evo/io/writer/fasta_writer'
13 require 'lib/evo/io/parser/fasta_parser'
14 require 'lib/evo/io/parser/hmmscan_parser'
15
16 module Evoruby
17   class HmmscanDomainExtractor
18
19     ADD_TO_CLOSE_PAIRS = 0
20     def initialize
21     end
22
23     # raises ArgumentError, IOError, StandardError
24     def parse( domain_id,
25       hmmscan_output,
26       fasta_sequence_file,
27       outfile,
28       passed_seqs_outfile,
29       failed_seqs_outfile,
30       e_value_threshold,
31       length_threshold,
32       add_position,
33       add_domain_number,
34       add_species,
35       min_linker,
36       log )
37
38       Util.check_file_for_readability( hmmscan_output )
39       Util.check_file_for_readability( fasta_sequence_file )
40       Util.check_file_for_writability( outfile + ".fasta" )
41       Util.check_file_for_writability( passed_seqs_outfile )
42       Util.check_file_for_writability( failed_seqs_outfile )
43
44       in_msa = nil
45       factory = MsaFactory.new()
46       in_msa = factory.create_msa_from_file( fasta_sequence_file, FastaParser.new() )
47
48       if ( in_msa == nil || in_msa.get_number_of_seqs() < 1 )
49         error_msg = "could not find fasta sequences in " + fasta_sequence_file
50         raise IOError, error_msg
51       end
52
53       out_msa = Msa.new
54
55       failed_seqs = Msa.new
56       passed_seqs = Msa.new
57       out_msa_pairs = nil
58       out_msa_isolated = nil
59       out_msa_singles = nil
60       out_msa_single_domains_protein_seqs = nil
61       out_msa_close_pairs_protein_seqs = nil
62       out_msa_close_pairs_only_protein_seqs = nil
63       out_msa_isolated_protein_seqs = nil
64       out_msa_isolated_only_protein_seqs = nil
65       out_msa_isolated_and_close_pair_protein_seqs = nil
66       if min_linker
67         out_msa_pairs = Msa.new
68         out_msa_isolated = Msa.new
69         out_msa_singles = Msa.new
70         out_msa_single_domains_protein_seqs = Msa.new
71         out_msa_close_pairs_protein_seqs = Msa.new
72         out_msa_close_pairs_only_protein_seqs = Msa.new
73         out_msa_isolated_protein_seqs = Msa.new
74         out_msa_isolated_only_protein_seqs = Msa.new
75         out_msa_isolated_and_close_pair_protein_seqs  = Msa.new
76       end
77
78       ld = Constants::LINE_DELIMITER
79
80       domain_pass_counter                = 0
81       domain_fail_counter                = 0
82       proteins_with_failing_domains      = 0
83       domain_not_present_counter         = 0
84       protein_counter                    = 1
85       max_domain_copy_number_per_protein = -1
86       max_domain_copy_number_sequence    = ""
87       passing_target_length_sum          = 0
88       overall_target_length_sum          = 0
89       overall_target_length_min          = 10000000
90       overall_target_length_max          = 0
91       passing_target_length_min          = 10000000
92       passing_target_length_max          = 0
93
94       hmmscan_datas = []
95
96       hmmscan_parser = HmmscanParser.new( hmmscan_output )
97       results = hmmscan_parser.parse
98
99       prev_query = nil
100       saw_target = false
101
102       results.each do | r |
103
104         if ( prev_query != nil ) && ( r.query != prev_query )
105           protein_counter += 1
106           if !saw_target
107             log << domain_not_present_counter.to_s + ": " + prev_query.to_s + " lacks target domain" + ld
108             domain_not_present_counter += 1
109           end
110           saw_target = false
111         end
112
113         prev_query = r.query
114
115         if domain_id != r.model
116           next
117         end
118
119         saw_target = true
120
121         sequence  = r.query
122         number    = r.number
123         out_of    = r.out_of
124         env_from  = r.env_from
125         env_to    = r.env_to
126         i_e_value = r.i_e_value
127         prev_query = r.query
128
129         length = env_to - env_from + 1
130
131         overall_target_length_sum += length
132         if length > overall_target_length_max
133           overall_target_length_max = length
134         end
135         if length < overall_target_length_min
136           overall_target_length_min = length
137         end
138
139         if ( ( ( e_value_threshold < 0.0 ) || ( i_e_value <= e_value_threshold ) ) &&
140         ( ( length_threshold <= 0 ) || ( length >= length_threshold.to_f ) ) )
141           hmmscan_datas << HmmsearchData.new( sequence, number, out_of, env_from, env_to, i_e_value )
142           passing_target_length_sum += length
143           if length > passing_target_length_max
144             passing_target_length_max = length
145           end
146           if length < passing_target_length_min
147             passing_target_length_min = length
148           end
149           if ( number > max_domain_copy_number_per_protein )
150             max_domain_copy_number_sequence    = sequence
151             max_domain_copy_number_per_protein = number
152           end
153         else # no pass
154           log << domain_fail_counter.to_s + ": " + sequence.to_s + " fails threshold(s)"
155           if ( ( e_value_threshold.to_f >= 0.0 ) && ( i_e_value > e_value_threshold ) )
156             log << " iE=" + i_e_value.to_s
157           end
158           if ( ( length_threshold.to_f > 0 ) && ( env_to - env_from + 1 ) < length_threshold.to_f )
159             le = env_to - env_from + 1
160             log << " l=" + le.to_s
161           end
162           log << ld
163           domain_fail_counter += 1
164         end
165
166         if number > out_of
167           error_msg = "number > out_of (this should not have happened)"
168           raise StandardError, error_msg
169         end
170
171         if number == out_of
172           if !hmmscan_datas.empty?
173             process_hmmscan_datas( hmmscan_datas,
174             in_msa,
175             add_position,
176             add_domain_number,
177             add_species,
178             out_msa,
179             out_msa_singles,
180             out_msa_pairs,
181             out_msa_isolated,
182             min_linker,
183             out_msa_single_domains_protein_seqs,
184             out_msa_close_pairs_protein_seqs,
185             out_msa_close_pairs_only_protein_seqs,
186             out_msa_isolated_protein_seqs,
187             out_msa_isolated_only_protein_seqs,
188             out_msa_isolated_and_close_pair_protein_seqs )
189             domain_pass_counter += hmmscan_datas.length
190             if passed_seqs.find_by_name_start( sequence, true ).length < 1
191               add_sequence( sequence, in_msa, passed_seqs )
192             else
193               error_msg = "this should not have happened"
194               raise StandardError, error_msg
195             end
196           else # no pass
197             if failed_seqs.find_by_name_start( sequence, true ).length < 1
198               add_sequence( sequence, in_msa, failed_seqs )
199               proteins_with_failing_domains += 1
200             else
201               error_msg = "this should not have happened"
202               raise StandardError, error_msg
203             end
204           end
205           hmmscan_datas.clear
206         end
207
208       end # results.each do | r |
209
210       if (prev_query != nil) && (!saw_target)
211         log << domain_not_present_counter.to_s + ": " + prev_query.to_s + " lacks target domain" + ld
212         domain_not_present_counter += 1
213       end
214
215       if domain_pass_counter < 1
216         error_msg = "no domain sequences were extracted"
217         raise IOError, error_msg
218       end
219
220       if ( domain_not_present_counter + passed_seqs.get_number_of_seqs + proteins_with_failing_domains ) != protein_counter
221         error_msg = "not present + passing + not passing != proteins in sequence (fasta) file (this should not have happened)"
222         raise StandardError, error_msg
223       end
224
225       puts
226       log << ld
227
228       log << ld
229       avg_pass = ( passing_target_length_sum / domain_pass_counter )
230       puts( "Passing target domain lengths: average: " + avg_pass.to_s  )
231       log << "Passing target domain lengths: average: " + avg_pass.to_s
232       log << ld
233       puts( "Passing target domain lengths: min-max: " + passing_target_length_min.to_s + "-"  + passing_target_length_max.to_s)
234       log << "Passing target domain lengths: min-max: " + passing_target_length_min.to_s + "-"  + passing_target_length_max.to_s
235       log << ld
236       puts( "Passing target domain lengths:     sum: " + domain_pass_counter.to_s  )
237       log << "Passing target domain lengths:     sum: " + domain_pass_counter.to_s
238       log << ld
239       log << ld
240       puts
241       sum = domain_pass_counter + domain_fail_counter
242       avg_all = overall_target_length_sum / sum
243       puts( "All target domain lengths:     average: " + avg_all.to_s  )
244       log << "All target domain lengths:     average: " +avg_all.to_s
245       log << ld
246       puts( "All target domain lengths:     min-max: " + overall_target_length_min.to_s + "-"  + overall_target_length_max.to_s)
247       log << "All target domain lengths:     min-max: " + overall_target_length_min.to_s + "-"  + overall_target_length_max.to_s
248       log << ld
249       puts( "All target domain lengths:         sum: " + sum.to_s  )
250       log << "All target domain lengths:         sum: " + sum.to_s
251
252       puts
253       puts( "Proteins with passing target domain(s): " + passed_seqs.get_number_of_seqs.to_s )
254       puts( "Proteins with no passing target domain: " + proteins_with_failing_domains.to_s )
255       puts( "Proteins with no target domain        : " + domain_not_present_counter.to_s )
256
257       log << ld
258       log << ld
259       puts
260       puts( "Max target domain copy number per protein (includes non-passing): " + max_domain_copy_number_per_protein.to_s )
261       log << "Max target domain copy number per protein (includes non-passing): " + max_domain_copy_number_per_protein.to_s
262       log << ld
263
264       if ( max_domain_copy_number_per_protein > 1 )
265         puts( "First target protein with this copy number: " + max_domain_copy_number_sequence )
266         log << "First target protein with this copy number: " + max_domain_copy_number_sequence
267         log << ld
268       end
269
270       write_msa( out_msa, outfile + ".fasta"  )
271       write_msa( passed_seqs, passed_seqs_outfile )
272       write_msa( failed_seqs, failed_seqs_outfile )
273
274       if out_msa_pairs
275         write_msa( out_msa_pairs, outfile + "_" + min_linker.to_s + ".fasta")
276       end
277
278       if out_msa_singles
279         write_msa( out_msa_singles, outfile + "_singles.fasta")
280       end
281
282       if out_msa_isolated
283         write_msa( out_msa_isolated, outfile + "_" + min_linker.to_s + "_isolated.fasta");
284       end
285
286       if out_msa_single_domains_protein_seqs
287         write_msa( out_msa_single_domains_protein_seqs, outfile + "_proteins_with_singles.fasta" )
288       end
289
290       if out_msa_close_pairs_protein_seqs
291         write_msa( out_msa_close_pairs_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_close_pairs.fasta" )
292       end
293
294       if out_msa_close_pairs_only_protein_seqs
295         write_msa( out_msa_close_pairs_only_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_close_pairs_only.fasta" )
296       end
297
298       if  out_msa_isolated_protein_seqs
299         write_msa(  out_msa_isolated_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_domains.fasta" )
300       end
301
302       if out_msa_isolated_only_protein_seqs
303         write_msa( out_msa_isolated_only_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_domains_only.fasta" )
304       end
305
306       if out_msa_isolated_and_close_pair_protein_seqs
307         write_msa( out_msa_isolated_and_close_pair_protein_seqs, outfile + "_" + min_linker.to_s + "_proteins_with_isolated_and_close_pairs.fasta" )
308       end
309
310       if min_linker
311         if ( out_msa_single_domains_protein_seqs.get_number_of_seqs +
312         out_msa_close_pairs_only_protein_seqs.get_number_of_seqs +
313         out_msa_isolated_only_protein_seqs.get_number_of_seqs +
314         out_msa_isolated_and_close_pair_protein_seqs.get_number_of_seqs ) != passed_seqs.get_number_of_seqs
315           error_msg = "this should not have happened"
316           raise StandardError, error_msg
317         end
318       end
319
320       log << ld
321       log << "passing target domains                       : " + domain_pass_counter.to_s + ld
322       log << "failing target domains                       : " + domain_fail_counter.to_s + ld
323       log << "proteins in sequence (fasta) file            : " + in_msa.get_number_of_seqs.to_s + ld
324       log << "proteins in hmmscan outputfile               : " + protein_counter.to_s + ld
325       log << "proteins with passing target domain(s)       : " + passed_seqs.get_number_of_seqs.to_s + ld
326       log << "proteins with no passing target domain       : " + proteins_with_failing_domains.to_s + ld
327       log << "proteins with no target domain               : " + domain_not_present_counter.to_s + ld
328       if min_linker
329         log << "min linker length                            : " + min_linker.to_s + ld
330         log << "single domains                               : " + out_msa_singles.get_number_of_seqs.to_s + ld
331         log << "domains in close pairs                       : " + (2 * out_msa_pairs.get_number_of_seqs).to_s + ld
332         log << "isolated domains                             : " + out_msa_isolated.get_number_of_seqs.to_s + ld
333         log << "proteins with single domains                 : " + out_msa_single_domains_protein_seqs.get_number_of_seqs.to_s + ld
334         log << "proteins with close pair domains             : " + out_msa_close_pairs_protein_seqs.get_number_of_seqs.to_s + ld
335         log << "proteins with close pair domains only        : " + out_msa_close_pairs_only_protein_seqs.get_number_of_seqs.to_s + ld
336         log << "proteins with isolated domains               : " + out_msa_isolated_protein_seqs.get_number_of_seqs.to_s + ld
337         log << "proteins with isolated domains only          : " + out_msa_isolated_only_protein_seqs.get_number_of_seqs.to_s + ld
338         log << "proteins with close pair and isolated domains: " + out_msa_isolated_and_close_pair_protein_seqs.get_number_of_seqs.to_s + ld
339       end
340
341       log << ld
342
343       return domain_pass_counter
344
345     end # parse
346
347     private
348
349     def write_msa( msa, filename )
350       io = MsaIO.new()
351       w = FastaWriter.new()
352       w.set_line_width( 60 )
353       w.clean( true )
354       begin
355         io.write_to_file( msa, filename, w )
356       rescue Exception
357         error_msg = "could not write to \"" + filename + "\""
358         raise IOError, error_msg
359       end
360     end
361
362     def add_sequence( sequence_name, in_msa, add_to_msa )
363       seqs = in_msa.find_by_name_start( sequence_name, true )
364       if ( seqs.length < 1 )
365         error_msg = "sequence \"" + sequence_name + "\" not found in sequence file"
366         raise StandardError, error_msg
367       end
368       if ( seqs.length > 1 )
369         error_msg = "sequence \"" + sequence_name + "\" not unique in sequence file"
370         raise StandardError, error_msg
371       end
372       seq = in_msa.get_sequence( seqs[ 0 ] )
373       add_to_msa.add_sequence( seq )
374     end
375
376     def process_hmmscan_datas( hmmscan_datas,
377       in_msa,
378       add_position,
379       add_domain_number,
380       add_species,
381       out_msa,
382       out_msa_singles,
383       out_msa_pairs,
384       out_msa_isolated,
385       min_linker,
386       out_msa_single_domains_protein_seqs,
387       out_msa_close_pairs_protein_seqs,
388       out_msa_close_pairs_only_protein_seqs,
389       out_msa_isolated_protein_seqs,
390       out_msa_isolated_only_protein_seqs,
391       out_msa_isolated_and_close_pair_protein_seqs )
392
393       actual_out_of = hmmscan_datas.size
394       saw_close_pair = false
395       saw_isolated = false
396
397       seq_name = ""
398       prev_seq_name = nil
399
400       hmmscan_datas.each_with_index do |hmmscan_data, index|
401         if hmmscan_data.number < ( index + 1 )
402           error_msg = "hmmscan_data.number < ( index + 1 ) (this should not have happened)"
403           raise StandardError, error_msg
404         end
405
406         seq_name =  hmmscan_data.seq_name
407
408         extract_domain( seq_name,
409         index + 1,
410         actual_out_of,
411         hmmscan_data.env_from,
412         hmmscan_data.env_to,
413         in_msa,
414         out_msa,
415         add_position,
416         add_domain_number,
417         add_species )
418
419         if min_linker
420           if actual_out_of == 1
421             extract_domain( seq_name,
422             1,
423             1,
424             hmmscan_data.env_from,
425             hmmscan_data.env_to,
426             in_msa,
427             out_msa_singles,
428             add_position,
429             add_domain_number,
430             add_species )
431             if out_msa_single_domains_protein_seqs.find_by_name_start( seq_name, true ).length < 1
432               add_sequence( seq_name, in_msa, out_msa_single_domains_protein_seqs )
433             else
434               error_msg = "this should not have happened"
435               raise StandardError, error_msg
436             end
437
438           else
439             first = index == 0
440             last = index == hmmscan_datas.length - 1
441
442             if ( ( first && ( ( hmmscan_datas[ index + 1 ].env_from - hmmscan_data.env_to ) > min_linker) )  ||
443             ( last && ( ( hmmscan_data.env_from - hmmscan_datas[ index - 1 ].env_to ) > min_linker ) ) ||
444             ( !first && !last &&  ( ( hmmscan_datas[ index + 1 ].env_from - hmmscan_data.env_to ) > min_linker ) &&
445             ( ( hmmscan_data.env_from - hmmscan_datas[ index - 1 ].env_to ) > min_linker ) ) )
446
447               extract_domain( seq_name,
448               index + 1,
449               actual_out_of,
450               hmmscan_data.env_from,
451               hmmscan_data.env_to,
452               in_msa,
453               out_msa_isolated,
454               add_position,
455               add_domain_number,
456               add_species )
457               saw_isolated = true
458
459             elsif !first
460
461               from = hmmscan_datas[ index - 1 ].env_from
462               to = hmmscan_data.env_to
463
464               if ADD_TO_CLOSE_PAIRS > 0
465                 from = from - ADD_TO_CLOSE_PAIRS
466                 if from < 1
467                   from = 1
468                 end
469                 to = to + ADD_TO_CLOSE_PAIRS
470                 temp_seqs = in_msa.find_by_name_start( seq_name, true )
471                 temp_seq = in_msa.get_sequence( temp_seqs[ 0 ] )
472                 if to >  temp_seq.get_length
473                   to =  temp_seq.get_length
474                 end
475               end
476
477               extract_domain( seq_name,
478               index.to_s  + "+" + ( index + 1 ).to_s,
479               actual_out_of,
480               from,
481               to,
482               in_msa,
483               out_msa_pairs,
484               add_position,
485               add_domain_number,
486               add_species )
487               saw_close_pair = true
488             end
489           end
490         end
491         if prev_seq_name && prev_seq_name != seq_name
492           error_msg = "this should not have happened"
493           raise StandardError, error_msg
494         end
495         prev_seq_name = seq_name
496       end # each
497       if saw_isolated
498         if out_msa_isolated_protein_seqs.find_by_name_start( seq_name, true ).length < 1
499           add_sequence( seq_name, in_msa, out_msa_isolated_protein_seqs )
500         else
501           error_msg = "this should not have happened"
502           raise StandardError, error_msg
503         end
504       end
505       if saw_close_pair
506         if out_msa_close_pairs_protein_seqs.find_by_name_start( seq_name, true ).length < 1
507           add_sequence( seq_name, in_msa, out_msa_close_pairs_protein_seqs )
508         else
509           error_msg = "this should not have happened"
510           raise StandardError, error_msg
511         end
512       end
513       if saw_close_pair && saw_isolated
514         if out_msa_isolated_and_close_pair_protein_seqs.find_by_name_start( seq_name, true ).length < 1
515           add_sequence( seq_name, in_msa, out_msa_isolated_and_close_pair_protein_seqs )
516         else
517           error_msg = "this should not have happened"
518           raise StandardError, error_msg
519         end
520       elsif saw_close_pair
521         if out_msa_close_pairs_only_protein_seqs.find_by_name_start( seq_name, true ).length < 1
522           add_sequence( seq_name, in_msa, out_msa_close_pairs_only_protein_seqs )
523         else
524           error_msg = "this should not have happened"
525           raise StandardError, error_msg
526         end
527       elsif saw_isolated
528         if out_msa_isolated_only_protein_seqs.find_by_name_start( seq_name, true ).length < 1
529           add_sequence( seq_name, in_msa, out_msa_isolated_only_protein_seqs )
530         else
531           error_msg = "this should not have happened"
532           raise StandardError, error_msg
533         end
534       end
535     end # def process_hmmscan_data
536
537     def extract_domain( sequence,
538       number,
539       out_of,
540       seq_from,
541       seq_to,
542       in_msa,
543       out_msa,
544       add_position,
545       add_domain_number,
546       add_species )
547       if number.is_a?( Fixnum ) && ( number < 1 || out_of < 1 || number > out_of )
548         error_msg = "number=" + number.to_s + ", out of=" + out_of.to_s
549         raise StandardError, error_msg
550       end
551       if seq_from < 1 || seq_to < 1 || seq_from >= seq_to
552         error_msg = "impossible: seq-from=" + seq_from.to_s + ", seq-to=" + seq_to.to_s
553         raise StandardError, error_msg
554       end
555       seqs = in_msa.find_by_name_start( sequence, true )
556       if seqs.length < 1
557         error_msg = "sequence \"" + sequence + "\" not found in sequence file"
558         raise IOError, error_msg
559       end
560       if seqs.length > 1
561         error_msg = "sequence \"" + sequence + "\" not unique in sequence file"
562         raise IOError, error_msg
563       end
564       # hmmscan is 1 based, whereas sequences are 0 bases in this package.
565       seq = in_msa.get_sequence( seqs[ 0 ] ).get_subsequence( seq_from - 1, seq_to - 1 )
566
567       orig_name = seq.get_name
568
569       seq.set_name( orig_name.split[ 0 ] )
570
571       if add_position
572         seq.set_name( seq.get_name + "_" + seq_from.to_s + "-" + seq_to.to_s )
573       end
574
575       if out_of != 1 && add_domain_number
576         seq.set_name( seq.get_name + "~" + number.to_s + "-" + out_of.to_s )
577       end
578
579       if add_species
580         a = orig_name.rindex "["
581         b = orig_name.rindex "]"
582         unless a && b
583           error_msg = "species not found in " + orig_name
584           raise StandardError, error_msg
585         end
586         species = orig_name[ a .. b ]
587         seq.set_name( seq.get_name + " " + species )
588       end
589       out_msa.add_sequence( seq )
590     end
591
592     def is_ignorable?( line )
593       return ( line !~ /[A-Za-z0-9-]/ || line =~/^#/ )
594     end
595
596   end # class HmmscanDomainExtractor
597
598   class HmmsearchData
599     def initialize( seq_name, number, out_of, env_from, env_to, i_e_value )
600       @seq_name = seq_name
601       @number = number
602       @out_of = out_of
603       @env_from = env_from
604       @env_to = env_to
605       @i_e_value = i_e_value
606     end
607     attr_reader :seq_name, :number, :out_of, :env_from, :env_to, :i_e_value
608   end
609
610 end # module Evoruby
611