inprogress
[jalview.git] / forester / ruby / evoruby / exe / run_phylo_pipeline_x.rb
1 #!/usr/local/bin/ruby -w
2 #
3 # = run_phylo_pipeline
4 #
5 # Copyright::  Copyright (C) 2010 Christian M. Zmasek
6 # License::    GNU Lesser General Public License (LGPL)
7 #
8 # $Id Exp $
9 #
10 #
11
12
13 module Evoruby
14
15   class RunPhyloPipeline
16
17     PFAM      = "/home/czmasek/DATA/PFAM/PFAM270X/"
18     HMMSCAN  = "/home/czmasek/SOFTWARE/HMMER/hmmer-3.0/src/hmmscan"
19     HSP       = "/home/czmasek/SOFTWARE/FORESTER/DEV/forester/forester/ruby/evoruby/exe/hsp.rb"
20     D2F       = "/home/czmasek/SOFTWARE/FORESTER/DEV/forester/forester/ruby/evoruby/exe/d2f.rb"
21     DSX       = "/home/czmasek/SOFTWARE/FORESTER/DEV/forester/forester/ruby/evoruby/exe/dsx.rb"
22     TAP       = "/home/czmasek/SOFTWARE/FORESTER/DEV/forester/forester/ruby/evoruby/exe/tap.rb"
23
24     def run
25       unless ARGV.length >= 2 && ARGV.length <= 4
26         error "arguments are:  <min-length> " +
27          "<neg E-value exponent for domain extraction> [E-value for hmmscan, default is 10] [hmmscan option, default is --nobias, --max for no heuristics]"
28       end
29       
30       length      = ARGV[ 0 ].to_i
31       e_value_exp = ARGV[ 1 ].to_i
32
33       e_for_hmmscan = 10
34       hmmscan_option = "--nobias"
35
36       if ARGV.length == 4
37         hmmscan_option = ARGV[ 3 ]
38       end
39       if ARGV.length == 3 || ARGV.length == 4
40         e_for_hmmscan = ARGV[ 2 ].to_i
41       end
42
43       if e_value_exp < 0
44         error "E-value exponent for domain extraction cannot be negative"
45       end
46       if length <= 1
47         error "length cannot be smaller than or equal to 1"
48       end
49       if e_for_hmmscan < 1
50         error "E-value for hmmscan cannot be smaller than 1"
51       end
52
53       input_files = Dir.entries(".").select { |f| !File.directory?( f ) && f.downcase.end_with?( ".fasta" ) }
54
55       puts "Input files:"
56       input_files.each do | input |
57         puts input
58       end
59       puts
60
61       counter = 1
62       input_files.each do | input |
63
64         puts counter.to_s + "/" +  input_files.size.to_s + " " + input + ": "
65         
66         counter += 1
67         
68         hmm_name = ""
69
70         if input.downcase.end_with?( "_ni.fasta" )
71           hmm_name = input[ 0 .. input.length - 10 ]
72         elsif input.downcase.end_with?( ".fasta" )
73           hmm_name = input[ 0 .. input.length - 7 ]
74           puts
75           puts "a. identifier normalization:"
76           cmd = "#{TAP} #{input} #{hmm_name}_ni.fasta #{hmm_name}.nim"
77           run_command( cmd )
78           input = hmm_name + "_ni.fasta"
79         else
80           error "illegal name: " + input
81         end
82
83         puts
84         puts "b. hmmscan:"
85         cmd = "#{HMMSCAN} #{hmmscan_option} --domtblout #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s} -E #{e_for_hmmscan.to_s} #{PFAM}Pfam-A.hmm #{input}"
86         run_command( cmd )
87         puts
88
89         puts "c. hmmscan to simple domain table:"
90         cmd = "#{HSP} #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s} #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s}_domain_table"
91         run_command( cmd )
92         puts
93
94         puts "d. domain table to forester format:"
95         cmd = "#{D2F} -e=10 #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s}_domain_table #{input} #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s}.dff"
96         run_command( cmd )
97         puts
98
99         puts "e. dsx:"
100         cmd = "#{DSX} -d -e=1e-#{e_value_exp.to_s} -l=#{length} #{hmm_name} #{hmm_name}_hmmscan_#{e_for_hmmscan.to_s} #{input} #{hmm_name}__#{hmm_name}__ee#{e_value_exp.to_s}_#{length}"
101         run_command( cmd )
102         puts
103
104       end
105
106     end
107
108     def run_command cmd
109       puts cmd
110       `#{cmd}`
111     end
112
113     def get_base_name n
114       if n.downcase.end_with?( "_ni.fasta" )
115         n[ 0 .. n.length - 10 ]
116       elsif n.downcase.end_with?( ".fasta" )
117         n[ 0 .. n.length - 7 ]
118       else
119         error "illegal name: " + n
120       end
121     end
122
123     def error msg
124       puts
125       puts msg
126       puts
127       exit
128     end
129
130   end
131
132   p = RunPhyloPipeline.new()
133
134   p.run()
135
136 end
137
138
139