source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/mafft-homologs.tmpl

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)

File size: 8.1 KB
Line 
1#!/usr/bin/env ruby
2
3localdb = "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
8mafftpath = "_BINDIR/mafft"   
9# path of mafft. "/usr/local/bin/mafft"
10# if mafft is in your command path, "mafft" is ok.
11
12blastpath = "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
34require 'getopts'
35require 'tempfile'
36
37# mktemp
38GC.disable
39temp_vf = Tempfile.new("_vf").path
40temp_if = Tempfile.new("_if").path
41temp_pf = Tempfile.new("_pf").path
42temp_af = Tempfile.new("_af").path
43temp_qf = Tempfile.new("_qf").path
44temp_bf = Tempfile.new("_bf").path
45temp_rid = Tempfile.new("_rid").path
46temp_res = Tempfile.new("_res").path
47
48
49system( mafftpath + " --help > #{temp_vf} 2>&1" )
50pfp = File.open( "#{temp_vf}", 'r' )
51while pfp.gets
52        break if $_ =~ /MAFFT v/
53end
54pfp.close
55if( $_ ) then
56        mafftversion = sub( /^\D*/, "" ).split(" ").slice(0).strip.to_s
57else
58        mafftversion = "0"
59end
60if( mafftversion < "5.58" ) then
61        puts ""
62        puts "======================================================"
63        puts "Install new mafft (v. >= 5.58)"
64        puts "======================================================"
65        puts ""
66        exit
67end
68
69srand ( 0 )
70
71def 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
86end
87
88nadd = 50
89eval = 1e-10
90local = 0
91fullout = 0
92entiresearch = 0
93corewin = 50
94corethr = 0.3
95mafftopt = " --op 1.53 --ep 0.123 --localpair --maxiterate 1000 --reorder "
96if 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
99end
100
101if $OPT_c then
102        corewin = $OPT_c.to_i
103end
104if $OPT_d then
105        corethr = $OPT_d.to_f
106end
107if $OPT_w
108        entiresearch = 1
109end
110if $OPT_f
111        fullout = 1
112end
113if $OPT_s
114        fullout = 0
115end
116if $OPT_l
117        local = 1
118end
119if $OPT_e then
120        eval = $OPT_e.to_f
121end
122if $OPT_a then
123        nadd = $OPT_a.to_i
124end
125if $OPT_o then
126        mafftopt += " " + $OPT_o + " "
127end
128
129system "cat " + ARGV.to_s + " > #{temp_if}"
130ar = mafftopt.split(" ")
131nar = ar.length
132for i in 0..(nar-1)
133        if ar[i] == "--seed" then
134                system "cat #{ar[i+1]} >> #{temp_if}"
135        end
136end
137
138nseq = 0
139ifp = File.open( "#{temp_if}", 'r' )
140        while ifp.gets
141                nseq += 1 if $_ =~ /^>/
142        end
143ifp.close
144
145if nseq >= 100 then
146        STDERR.puts "The number of input sequences must be <100."
147        exit
148elsif nseq == 1 then
149        system( "cp #{temp_if}"  + " #{temp_pf}" )
150else
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
158end
159
160pfp = File.open( "#{temp_pf}", 'r' )
161inname = []
162inseq = []
163slen = []
164act = []
165nin = 0
166nin = readfasta( pfp, inname, inseq )
167for i in 0..(nin-1)
168        slen.push( inseq[i].gsub(/-/,"").length )
169        act.push( 1 )
170end
171pfp.close
172
173pfp = File.open( "#{temp_if}", 'r' )
174orname = []
175orseq = []
176nin = 0
177nin = readfasta( pfp, orname, orseq )
178pfp.close
179
180allen = inseq[0].length
181for 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
202end
203#p act
204
205
206afp = File.open( "#{temp_af}", 'w' )
207
208STDERR.puts "Searching .. \n"
209ids = []
210add = []
211sco = []
212for 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( /-/, "" ).sub( /U/, "X" )
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
311end
312
313n = ids.length
314
315outnum = 0
316while 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
324end
325afp.close
326
327STDERR.puts "Performing alignment .. "
328system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
329STDERR.puts "done."
330
331bfp = File.open( "#{temp_bf}", 'r' )
332outseq = []
333outnam = []
334readfasta( bfp, outnam, outseq )
335bfp.close
336
337outseq2 = []
338outnam2 = []
339
340len = outseq.length
341for 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_" ) )
348end
349
350nout = outseq2.length
351len = outseq[0].length
352p = len
353while 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
367end
368for i in 0..(nout-1)
369        puts ">" + outnam2[i]
370        puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
371end
372
373
374system( "rm -rf #{temp_if} #{temp_vf} #{temp_af} #{temp_bf} #{temp_pf} #{temp_qf} #{temp_res} #{temp_rid}" )
Note: See TracBrowser for help on using the repository browser.