Next version of JABA
[jabaws.git] / binaries / src / mafft / core / mafft-homologs.tmpl
1 #!/usr/bin/env ruby
2
3 localdb = "sp"        
4 # database name from which homologues are collected 
5 # by locally installed blast. Leave this if you do 
6 # not use the '-l' option.
7
8 mafftpath = "_BINDIR/mafft"   
9 # path of mafft. "/usr/local/bin/mafft"
10 # if mafft is in your command path, "mafft" is ok.
11
12 blastpath = "blastall"   
13 # path of blastall. 
14 # if blastall is in your command path, "blastall" is ok.
15
16 # mafft-homologs.rb  v. 2.1 aligns sequences together with homologues 
17 # automatically collected from SwissProt via NCBI BLAST.
18 #
19 # mafft > 5.58 is required
20 #
21 # Usage:
22 #   mafft-homologs.rb [options] input > output
23 # Options:
24 #   -a #      the number of collected sequences (default: 50)
25 #   -e #      threshold value (default: 1e-10)
26 #   -o "xxx"  options for mafft 
27 #             (default: " --op 1.53 --ep 0.123 --maxiterate 1000")
28 #   -l        locally carries out blast searches instead of NCBI blast
29 #             (requires locally installed blast and a database)
30 #   -f        outputs collected homologues also (default: off)
31 #   -w        entire sequences are subjected to BLAST search 
32 #             (default: well-aligned region only)
33
34 require 'getopts'
35 require 'tempfile'
36
37 # mktemp
38 GC.disable
39 temp_vf = Tempfile.new("_vf").path
40 temp_if = Tempfile.new("_if").path
41 temp_pf = Tempfile.new("_pf").path
42 temp_af = Tempfile.new("_af").path
43 temp_qf = Tempfile.new("_qf").path
44 temp_bf = Tempfile.new("_bf").path
45 temp_rid = Tempfile.new("_rid").path
46 temp_res = Tempfile.new("_res").path
47
48
49 system( mafftpath + " --help > #{temp_vf} 2>&1" )
50 pfp = File.open( "#{temp_vf}", 'r' )
51 while pfp.gets
52         break if $_ =~ /MAFFT v/
53 end
54 pfp.close
55 if( $_ ) then
56         mafftversion = sub( /^\D*/, "" ).split(" ").slice(0).strip.to_s
57 else
58         mafftversion = "0"
59 end
60 if( mafftversion < "5.58" ) then
61         puts ""
62         puts "======================================================"
63         puts "Install new mafft (v. >= 5.58)"
64         puts "======================================================"
65         puts ""
66         exit
67 end
68
69 srand ( 0 )
70
71 def readfasta( fp, name, seq )
72         nseq = 0
73         tmpseq = ""
74         while fp.gets
75                 if $_ =~ /^>/ then
76                         name.push( $_.sub(/>/,"").strip )
77                         seq.push( tmpseq ) if nseq > 0
78                         nseq += 1
79                         tmpseq = ""
80                 else
81                         tmpseq += $_.strip
82                 end
83         end
84         seq.push( tmpseq )
85         return nseq
86 end
87
88 nadd = 50
89 eval = 1e-10
90 local = 0
91 fullout = 0
92 entiresearch = 0
93 corewin = 50
94 corethr = 0.3
95 mafftopt = " --op 1.53 --ep 0.123 --localpair --maxiterate 1000 --reorder "
96 if getopts( "s", "f", "w", "l", "h", "e:", "a:", "o:", "c:", "d:" ) == nil ||  ARGV.length == 0 || $OPT_h then
97         puts "Usage: #{$0} [-h -l -e# -a# -o\"[options for mafft]\"] input_file"
98         exit
99 end
100
101 if $OPT_c then
102         corewin = $OPT_c.to_i
103 end
104 if $OPT_d then
105         corethr = $OPT_d.to_f
106 end
107 if $OPT_w
108         entiresearch = 1
109 end
110 if $OPT_f
111         fullout = 1
112 end
113 if $OPT_s
114         fullout = 0
115 end
116 if $OPT_l
117         local = 1
118 end
119 if $OPT_e then
120         eval = $OPT_e.to_f
121 end
122 if $OPT_a then
123         nadd = $OPT_a.to_i
124 end
125 if $OPT_o then
126         mafftopt += " " + $OPT_o + " "
127 end
128
129 system "cat " + ARGV.to_s + " > #{temp_if}"
130 ar = mafftopt.split(" ")
131 nar = ar.length
132 for i in 0..(nar-1)
133         if ar[i] == "--seed" then
134                 system "cat #{ar[i+1]} >> #{temp_if}"
135         end
136 end
137
138 nseq = 0
139 ifp = File.open( "#{temp_if}", 'r' )
140         while ifp.gets
141                 nseq += 1 if $_ =~ /^>/
142         end
143 ifp.close
144
145 if nseq >= 100 then
146         STDERR.puts "The number of input sequences must be <100."
147         exit
148 elsif nseq == 1 then
149         system( "cp #{temp_if}"  + " #{temp_pf}" )
150 else
151         STDERR.puts "Performing preliminary alignment .. "
152         if entiresearch == 1 then
153 #               system( mafftpath + " --maxiterate 1000 --localpair #{temp_if} > #{temp_pf}" )
154                 system( mafftpath + " --maxiterate 0 --retree 2 #{temp_if} > #{temp_pf}" )
155         else
156                 system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} #{temp_if} > #{temp_pf}" )
157         end
158 end
159
160 pfp = File.open( "#{temp_pf}", 'r' )
161 inname = []
162 inseq = []
163 slen = []
164 act = []
165 nin = 0
166 nin = readfasta( pfp, inname, inseq )
167 for i in 0..(nin-1)
168         slen.push( inseq[i].gsub(/-/,"").length )
169         act.push( 1 )
170 end
171 pfp.close
172
173 pfp = File.open( "#{temp_if}", 'r' )
174 orname = []
175 orseq = []
176 nin = 0
177 nin = readfasta( pfp, orname, orseq )
178 pfp.close
179
180 allen = inseq[0].length
181 for i in 0..(nin-2)
182         for j in (i+1)..(nin-1)
183                 next if act[i] == 0
184                 next if act[j] == 0
185                 pid = 0.0
186                 total = 0
187                 for a in 0..(allen-1)
188                         next if inseq[i][a,1] == "-" || inseq[j][a,1] == "-"
189                         total += 1
190                         pid += 1.0 if inseq[i][a,1] == inseq[j][a,1]
191                 end
192                 pid /= total
193 #               puts "#{i.to_s}, #{j.to_s}, #{pid.to_s}"
194                 if pid > 0.5 then
195                         if slen[i] < slen[j]
196                                 act[i] = 0 
197                         else
198                                 act[j] = 0 
199                         end
200                 end
201         end
202 end
203 #p act
204
205
206 afp = File.open( "#{temp_af}", 'w' )
207
208 STDERR.puts "Searching .. \n"
209 ids = []
210 add = []
211 sco = []
212 for i in 0..(nin-1)
213         inseq[i].gsub!(/-/,"")
214         afp.puts ">" + orname[i]
215         afp.puts orseq[i]
216
217 #       afp.puts ">" + inname[i]
218 #       afp.puts inseq[i]
219
220         STDERR.puts "Query (#{i+1}/#{nin})\n" + inname[i]
221         if act[i] == 0 then
222                 STDERR.puts "Skip.\n\n"
223                 next 
224         end
225
226         if local == 0 then
227                 command = "lynx -source 'http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?QUERY=" + inseq[i] + "&DATABASE=swissprot&HITLIST_SIZE=" + nadd.to_s + "&FILTER=L&EXPECT='" + eval.to_s + "'&FORMAT_TYPE=TEXT&PROGRAM=blastp&SERVICE=plain&NCBI_GI=on&PAGE=Proteins&CMD=Put' > #{temp_rid}"
228                 system command
229         
230                 ridp = File.open( "#{temp_rid}", 'r' )
231                 while ridp.gets
232                         break if $_ =~ / RID = (.*)/
233                 end
234                 ridp.close
235                 rid = $1.strip
236                 STDERR.puts "Submitted to NCBI. rid = " + rid
237         
238                 STDERR.printf "Waiting "
239                 while 1 
240                         STDERR.printf "."
241                         sleep 10
242                         command = "lynx -source 'http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?RID=" + rid + "&DESCRIPTIONS=500&ALIGNMENTS=" + nadd.to_s + "&ALIGNMENT_TYPE=Pairwise&OVERVIEW=no&CMD=Get&FORMAT_TYPE=XML' > #{temp_res}"
243                         system command
244                         resp = File.open( "#{temp_res}", 'r' )
245 #                       resp.gets
246 #                       if $_ =~ /WAITING/ then
247 #                               resp.close
248 #                               next
249 #                       end
250                         while( resp.gets )
251                                 break if $_ =~ /QBlastInfoBegin/
252                         end
253                         resp.gets
254                         if $_ =~ /WAITING/ then
255                                 resp.close
256                                 next
257                         else
258                                 resp.close
259                                 break
260                         end
261                 end
262         else
263 #               puts "Not supported"
264 #               exit
265                 qfp = File.open( "#{temp_qf}", 'w' )
266                         qfp.puts "> "
267                         qfp.puts inseq[i]
268                 qfp.close
269                 command = blastpath + "  -p blastp  -e #{eval} -b 1000 -m 7 -i #{temp_qf} -d #{localdb} > #{temp_res}"
270                 system command
271                 resp = File.open( "#{temp_res}", 'r' )
272         end
273         STDERR.puts " Done.\n\n"
274
275         resp = File.open( "#{temp_res}", 'r' )
276         while 1
277                 while resp.gets
278                         break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ || $_ =~ /(<Iteration_stat>)/
279                 end
280                 id = $1
281                 break if $_ =~ /<Iteration_stat>/
282 #               p id
283                 while resp.gets
284                         break if $_ =~ /<Hsp_bit-score>(.*)<\/Hsp_bit-score>/
285                 end
286                 score = $1.to_f
287 #               p score
288
289                 known = ids.index( id )
290                 if known != nil then
291                         if sco[known] >= score then
292                                 next
293                         else
294                                 ids.delete_at( known )
295                                 add.delete_at( known )
296                                 sco.delete_at( known )
297                         end
298                 end
299                 while resp.gets
300                         break if $_ =~ /<Hsp_hseq>(.*)<\/Hsp_hseq>/
301                 end
302 #               break if $1 == nil
303                 target = $1.sub( /-/, "" )
304 #               p target
305 #               STDERR.puts "adding 1 seq"
306                 ids.push( id )
307                 sco.push( score )
308                 add.push( target )
309         end
310         resp.close
311 end
312
313 n = ids.length
314
315 outnum = 0
316 while n > 0 && outnum < nadd
317         m = rand( n )
318         afp.puts ">_addedbymaffte_" + ids[m]
319         afp.puts add[m]
320         ids.delete_at( m )
321         add.delete_at( m )
322         n -= 1
323         outnum += 1
324 end
325 afp.close
326
327 STDERR.puts "Performing alignment .. "
328 system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
329 STDERR.puts "done."
330
331 bfp = File.open( "#{temp_bf}", 'r' )
332 outseq = []
333 outnam = []
334 readfasta( bfp, outnam, outseq )
335 bfp.close
336
337 outseq2 = []
338 outnam2 = []
339
340 len = outseq.length
341 for i in 0..(len-1)
342 #       p outnam[i]
343         if fullout == 0 && outnam[i] =~ /^_addedbymaffte_/ then
344                 next
345         end
346         outseq2.push( outseq[i] )
347         outnam2.push( outnam[i].sub( /^_addedbymaffte_/, "_ho_" ) )
348 end
349
350 nout = outseq2.length
351 len = outseq[0].length
352 p = len
353 while p>0
354         p -= 1
355     allgap = 1
356     for j in 0..(nout-1)
357                 if outseq2[j][p,1] != "-" then
358                         allgap = 0
359                         break
360                 end
361     end
362     if allgap == 1 then
363         for j in 0..(nout-1)
364             outseq2[j][p,1] = ""
365         end
366     end
367 end
368 for i in 0..(nout-1)
369         puts ">" + outnam2[i]
370         puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
371 end
372
373
374 system( "rm -rf #{temp_if} #{temp_vf} #{temp_af} #{temp_bf} #{temp_pf} #{temp_qf} #{temp_res} #{temp_rid}" )