--- /dev/null
+#!/usr/bin/env ruby
+
+localdb = "sp"
+# database name from which homologues are collected
+# by locally installed blast. Leave this if you do
+# not use the '-l' option.
+
+mafftpath = "/usr/local/bin/mafft"
+# path of mafft. "/usr/local/bin/mafft"
+# if mafft is in your command path, "mafft" is ok.
+
+blastpath = "blastall"
+# path of blastall.
+# if blastall is in your command path, "blastall" is ok.
+
+# mafft-homologs.rb v. 2.1 aligns sequences together with homologues
+# automatically collected from SwissProt via NCBI BLAST.
+#
+# mafft > 5.58 is required
+#
+# Usage:
+# mafft-homologs.rb [options] input > output
+# Options:
+# -a # the number of collected sequences (default: 50)
+# -e # threshold value (default: 1e-10)
+# -o "xxx" options for mafft
+# (default: " --op 1.53 --ep 0.123 --maxiterate 1000")
+# -l locally carries out blast searches instead of NCBI blast
+# (requires locally installed blast and a database)
+# -f outputs collected homologues also (default: off)
+# -w entire sequences are subjected to BLAST search
+# (default: well-aligned region only)
+
+require 'getopts'
+require 'tempfile'
+
+# mktemp
+GC.disable
+temp_vf = Tempfile.new("_vf").path
+temp_if = Tempfile.new("_if").path
+temp_pf = Tempfile.new("_pf").path
+temp_af = Tempfile.new("_af").path
+temp_qf = Tempfile.new("_qf").path
+temp_bf = Tempfile.new("_bf").path
+temp_rid = Tempfile.new("_rid").path
+temp_res = Tempfile.new("_res").path
+
+
+system( mafftpath + " --help > #{temp_vf} 2>&1" )
+pfp = File.open( "#{temp_vf}", 'r' )
+while pfp.gets
+ break if $_ =~ /MAFFT v/
+end
+pfp.close
+if( $_ ) then
+ mafftversion = sub( /^\D*/, "" ).split(" ").slice(0).strip.to_s
+else
+ mafftversion = "0"
+end
+if( mafftversion < "5.58" ) then
+ puts ""
+ puts "======================================================"
+ puts "Install new mafft (v. >= 5.58)"
+ puts "======================================================"
+ puts ""
+ exit
+end
+
+srand ( 0 )
+
+def readfasta( fp, name, seq )
+ nseq = 0
+ tmpseq = ""
+ while fp.gets
+ if $_ =~ /^>/ then
+ name.push( $_.sub(/>/,"").strip )
+ seq.push( tmpseq ) if nseq > 0
+ nseq += 1
+ tmpseq = ""
+ else
+ tmpseq += $_.strip
+ end
+ end
+ seq.push( tmpseq )
+ return nseq
+end
+
+nadd = 50
+eval = 1e-10
+local = 0
+fullout = 0
+entiresearch = 0
+corewin = 50
+corethr = 0.3
+mafftopt = " --op 1.53 --ep 0.123 --localpair --maxiterate 1000 --reorder "
+if getopts( "s", "f", "w", "l", "h", "e:", "a:", "o:", "c:", "d:" ) == nil || ARGV.length == 0 || $OPT_h then
+ puts "Usage: #{$0} [-h -l -e# -a# -o\"[options for mafft]\"] input_file"
+ exit
+end
+
+if $OPT_c then
+ corewin = $OPT_c.to_i
+end
+if $OPT_d then
+ corethr = $OPT_d.to_f
+end
+if $OPT_w
+ entiresearch = 1
+end
+if $OPT_f
+ fullout = 1
+end
+if $OPT_s
+ fullout = 0
+end
+if $OPT_l
+ local = 1
+end
+if $OPT_e then
+ eval = $OPT_e.to_f
+end
+if $OPT_a then
+ nadd = $OPT_a.to_i
+end
+if $OPT_o then
+ mafftopt += " " + $OPT_o + " "
+end
+
+system "cat " + ARGV.to_s + " > #{temp_if}"
+ar = mafftopt.split(" ")
+nar = ar.length
+for i in 0..(nar-1)
+ if ar[i] == "--seed" then
+ system "cat #{ar[i+1]} >> #{temp_if}"
+ end
+end
+
+nseq = 0
+ifp = File.open( "#{temp_if}", 'r' )
+ while ifp.gets
+ nseq += 1 if $_ =~ /^>/
+ end
+ifp.close
+
+if nseq >= 100 then
+ STDERR.puts "The number of input sequences must be <100."
+ exit
+elsif nseq == 1 then
+ system( "cp #{temp_if}" + " #{temp_pf}" )
+else
+ STDERR.puts "Performing preliminary alignment .. "
+ if entiresearch == 1 then
+# system( mafftpath + " --maxiterate 1000 --localpair #{temp_if} > #{temp_pf}" )
+ system( mafftpath + " --maxiterate 0 --retree 2 #{temp_if} > #{temp_pf}" )
+ else
+ system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} #{temp_if} > #{temp_pf}" )
+ end
+end
+
+pfp = File.open( "#{temp_pf}", 'r' )
+inname = []
+inseq = []
+slen = []
+act = []
+nin = 0
+nin = readfasta( pfp, inname, inseq )
+for i in 0..(nin-1)
+ slen.push( inseq[i].gsub(/-/,"").length )
+ act.push( 1 )
+end
+pfp.close
+
+pfp = File.open( "#{temp_if}", 'r' )
+orname = []
+orseq = []
+nin = 0
+nin = readfasta( pfp, orname, orseq )
+pfp.close
+
+allen = inseq[0].length
+for i in 0..(nin-2)
+ for j in (i+1)..(nin-1)
+ next if act[i] == 0
+ next if act[j] == 0
+ pid = 0.0
+ total = 0
+ for a in 0..(allen-1)
+ next if inseq[i][a,1] == "-" || inseq[j][a,1] == "-"
+ total += 1
+ pid += 1.0 if inseq[i][a,1] == inseq[j][a,1]
+ end
+ pid /= total
+# puts "#{i.to_s}, #{j.to_s}, #{pid.to_s}"
+ if pid > 0.5 then
+ if slen[i] < slen[j]
+ act[i] = 0
+ else
+ act[j] = 0
+ end
+ end
+ end
+end
+#p act
+
+
+afp = File.open( "#{temp_af}", 'w' )
+
+STDERR.puts "Searching .. \n"
+ids = []
+add = []
+sco = []
+for i in 0..(nin-1)
+ inseq[i].gsub!(/-/,"")
+ afp.puts ">" + orname[i]
+ afp.puts orseq[i]
+
+# afp.puts ">" + inname[i]
+# afp.puts inseq[i]
+
+ STDERR.puts "Query (#{i+1}/#{nin})\n" + inname[i]
+ if act[i] == 0 then
+ STDERR.puts "Skip.\n\n"
+ next
+ end
+
+ if local == 0 then
+ 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}"
+ system command
+
+ ridp = File.open( "#{temp_rid}", 'r' )
+ while ridp.gets
+ break if $_ =~ / RID = (.*)/
+ end
+ ridp.close
+ rid = $1.strip
+ STDERR.puts "Submitted to NCBI. rid = " + rid
+
+ STDERR.printf "Waiting "
+ while 1
+ STDERR.printf "."
+ sleep 10
+ 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}"
+ system command
+ resp = File.open( "#{temp_res}", 'r' )
+# resp.gets
+# if $_ =~ /WAITING/ then
+# resp.close
+# next
+# end
+ while( resp.gets )
+ break if $_ =~ /QBlastInfoBegin/
+ end
+ resp.gets
+ if $_ =~ /WAITING/ then
+ resp.close
+ next
+ else
+ resp.close
+ break
+ end
+ end
+ else
+# puts "Not supported"
+# exit
+ qfp = File.open( "#{temp_qf}", 'w' )
+ qfp.puts "> "
+ qfp.puts inseq[i]
+ qfp.close
+ command = blastpath + " -p blastp -e #{eval} -b 1000 -m 7 -i #{temp_qf} -d #{localdb} > #{temp_res}"
+ system command
+ resp = File.open( "#{temp_res}", 'r' )
+ end
+ STDERR.puts " Done.\n\n"
+
+ resp = File.open( "#{temp_res}", 'r' )
+ while 1
+ while resp.gets
+ break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ || $_ =~ /(<Iteration_stat>)/
+ end
+ id = $1
+ break if $_ =~ /<Iteration_stat>/
+# p id
+ while resp.gets
+ break if $_ =~ /<Hsp_bit-score>(.*)<\/Hsp_bit-score>/
+ end
+ score = $1.to_f
+# p score
+
+ known = ids.index( id )
+ if known != nil then
+ if sco[known] >= score then
+ next
+ else
+ ids.delete_at( known )
+ add.delete_at( known )
+ sco.delete_at( known )
+ end
+ end
+ while resp.gets
+ break if $_ =~ /<Hsp_hseq>(.*)<\/Hsp_hseq>/
+ end
+# break if $1 == nil
+ target = $1.sub( /-/, "" ).sub( /U/, "X" )
+# p target
+# STDERR.puts "adding 1 seq"
+ ids.push( id )
+ sco.push( score )
+ add.push( target )
+ end
+ resp.close
+end
+
+n = ids.length
+
+outnum = 0
+while n > 0 && outnum < nadd
+ m = rand( n )
+ afp.puts ">_addedbymaffte_" + ids[m]
+ afp.puts add[m]
+ ids.delete_at( m )
+ add.delete_at( m )
+ n -= 1
+ outnum += 1
+end
+afp.close
+
+STDERR.puts "Performing alignment .. "
+system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
+STDERR.puts "done."
+
+bfp = File.open( "#{temp_bf}", 'r' )
+outseq = []
+outnam = []
+readfasta( bfp, outnam, outseq )
+bfp.close
+
+outseq2 = []
+outnam2 = []
+
+len = outseq.length
+for i in 0..(len-1)
+# p outnam[i]
+ if fullout == 0 && outnam[i] =~ /_addedbymaffte_/ then
+ next
+ end
+ outseq2.push( outseq[i] )
+ outnam2.push( outnam[i].sub( /_addedbymaffte_/, "_ho_" ) )
+end
+
+nout = outseq2.length
+len = outseq[0].length
+p = len
+while p>0
+ p -= 1
+ allgap = 1
+ for j in 0..(nout-1)
+ if outseq2[j][p,1] != "-" then
+ allgap = 0
+ break
+ end
+ end
+ if allgap == 1 then
+ for j in 0..(nout-1)
+ outseq2[j][p,1] = ""
+ end
+ end
+end
+for i in 0..(nout-1)
+ puts ">" + outnam2[i]
+ puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
+end
+
+
+system( "rm -rf #{temp_if} #{temp_vf} #{temp_af} #{temp_bf} #{temp_pf} #{temp_qf} #{temp_res} #{temp_rid}" )
git://source.jalview.org / jabaws.git / blobdiff