Attachment #164660 for bug #235804

View | Details | Raw Unified | Return to bug 235804
Collapse All | Expand All




#   -w        entire sequences are subjected to BLAST search 
#             (default: well-aligned region only)


require 'getopts'
require 'tempfile'

# create temporary files
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

	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 = []

end
pfp.close

pfp = File.open( "#{temp_if}", 'r' )
orname = []
orseq = []
nin = 0

#p act


afp = File.open( "#{temp_af}", 'w' )

STDERR.puts "Searching .. \n"
ids = []

	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

		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

	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>)/

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 = []


	puts ">" + outnam2[i]
	puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
end

system( "rm -rf /tmp/_if#{$$} /tmp/_vf#{$$} /tmp/_af#{$$} /tmp/_bf#{$$} /tmp/_pf#{$$} /tmp/_q#{$$} /tmp/_res#{$$} /tmp/_rid#{$$}" )

Lines 240-250 Link Here

(-)mafft-6.240/src/mafft.tmpl (-4 / +7 lines)
240	shift	240	shift
241	done;	241	done;
242		242
243	# TMPFILE=/tmp/`basename $0`.`whoami`.$$.`date +%y%m%d%H%M%S`	243	# create temporary directory; terminate script in case of failure
244	TMPFILE=/tmp/$progname.$$	244	if ! TMPFILE=`mktemp -dt $progname.XXXXXXXXXX`; then
		245	echo "Failed to create temporary directory"
		246	exit 1
		247	fi
		248
245	umask 077	249	umask 077
246	mkdir $TMPFILE \|\| er=1	250	trap "rm -rf $TMPFILE " 0
247	trap "rm -r $TMPFILE " 0
248	if [ $# -eq 1 ]; then	251	if [ $# -eq 1 ]; then
249	if [ -r "$1" -o "$1" = - ]; then	252	if [ -r "$1" -o "$1" = - ]; then
250	cat "$1" \| tr "\r" "\n" > $TMPFILE/infile	253	cat "$1" \| tr "\r" "\n" > $TMPFILE/infile

Return to bug 235804

Lines 31-41 Link Here

(-)mafft-6.240/src/mafft-homologs.tmpl (-29 / +39 lines)
31	# -w entire sequences are subjected to BLAST search	31	# -w entire sequences are subjected to BLAST search
32	# (default: well-aligned region only)	32	# (default: well-aligned region only)
33		33
34
35	require 'getopts'	34	require 'getopts'
		35	require 'tempfile'
		36
		37	# create temporary files
		38	temp_vf = Tempfile.new("_vf").path
		39	temp_if = Tempfile.new("_if").path
		40	temp_pf = Tempfile.new("_pf").path
		41	temp_af = Tempfile.new("_af").path
		42	temp_qf = Tempfile.new("_qf").path
		43	temp_bf = Tempfile.new("_bf").path
		44	temp_rid = Tempfile.new("_rid").path
		45	temp_res = Tempfile.new("_res").path
36		46
37	system( mafftpath + " --help > /tmp/_vf#{$$} 2>&1" )	47
38	pfp = File.open( "/tmp/_vf#{$$}", 'r' )	48	system( mafftpath + " --help > #{temp_vf} 2>&1" )
		49	pfp = File.open( "#{temp_vf}", 'r' )
39	while pfp.gets	50	while pfp.gets
40	break if $_ =~ /MAFFT v/	51	break if $_ =~ /MAFFT v/
41	end	52	end
Lines 114-148 Link Here
114	mafftopt += " " + $OPT_o + " "	125	mafftopt += " " + $OPT_o + " "
115	end	126	end
116		127
117	system "cat " + ARGV.to_s + " > /tmp/_if#{$$}"	128	system "cat " + ARGV.to_s + " > #{temp_if}"
118	ar = mafftopt.split(" ")	129	ar = mafftopt.split(" ")
119	nar = ar.length	130	nar = ar.length
120	for i in 0..(nar-1)	131	for i in 0..(nar-1)
121	if ar[i] == "--seed" then	132	if ar[i] == "--seed" then
122	system "cat #{ar[i+1]} >> /tmp/_if#{$$}"	133	system "cat #{ar[i+1]} >> #{temp_if}"
123	end	134	end
124	end	135	end
125		136
126	nseq = 0	137	nseq = 0
127	ifp = File.open( "/tmp/_if#{$$}", 'r' )	138	ifp = File.open( "#{temp_if}", 'r' )
128	while ifp.gets	139	while ifp.gets
129	nseq += 1 if $_ =~ /^>/	140	nseq += 1 if $_ =~ /^>/
130	end	141	end
131	ifp.close	142	ifp.close
132		143
133	STDERR.puts "Performing preliminary alignment .. "	144	if nseq >= 100 then
134	if nseq == 1 then	145	STDERR.puts "The number of input sequences must be <100."
135	system( "cp /tmp/_if#{$$}" + " /tmp/_pf#{$$}" )	146	exit
		147	elsif nseq == 1 then
		148	system( "cp #{temp_if}" + " #{temp_pf}" )
136	else	149	else
		150	STDERR.puts "Performing preliminary alignment .. "
137	if entiresearch == 1 then	151	if entiresearch == 1 then
138	# system( mafftpath + " --maxiterate 1000 --localpair /tmp/_if#{$$} > /tmp/_pf#{$$}" )	152	# system( mafftpath + " --maxiterate 1000 --localpair #{temp_if} > #{temp_pf}" )
139	system( mafftpath + " --maxiterate 0 --retree 2 /tmp/_if#{$$} > /tmp/_pf#{$$}" )	153	system( mafftpath + " --maxiterate 0 --retree 2 #{temp_if} > #{temp_pf}" )
140	else	154	else
141	system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} /tmp/_if#{$$} > /tmp/_pf#{$$}" )	155	system( mafftpath + " --maxiterate 1000 --localpair --core --coreext --corethr #{corethr.to_s} --corewin #{corewin.to_s} #{temp_if} > #{temp_pf}" )
142	end	156	end
143	end	157	end
144		158
145	pfp = File.open( "/tmp/_pf#{$$}", 'r' )	159	pfp = File.open( "#{temp_pf}", 'r' )
146	inname = []	160	inname = []
147	inseq = []	161	inseq = []
148	slen = []	162	slen = []
Lines 155-161 Link Here
155	end	169	end
156	pfp.close	170	pfp.close
157		171
158	pfp = File.open( "/tmp/_if#{$$}", 'r' )	172	pfp = File.open( "#{temp_if}", 'r' )
159	orname = []	173	orname = []
160	orseq = []	174	orseq = []
161	nin = 0	175	nin = 0
Lines 188-194 Link Here
188	#p act	202	#p act
189		203
190		204
191	afp = File.open( "/tmp/_af#{$$}", 'w' )	205	afp = File.open( "#{temp_af}", 'w' )
192		206
193	STDERR.puts "Searching .. \n"	207	STDERR.puts "Searching .. \n"
194	ids = []	208	ids = []
Lines 209-218 Link Here
209	end	223	end
210		224
211	if local == 0 then	225	if local == 0 then
212	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' > /tmp/_rid#{$$}"	226	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}"
213	system command	227	system command
214		228
215	ridp = File.open( "/tmp/_rid#{$$}", 'r' )	229	ridp = File.open( "#{temp_rid}", 'r' )
216	while ridp.gets	230	while ridp.gets
217	break if $_ =~ / RID = (.*)/	231	break if $_ =~ / RID = (.*)/
218	end	232	end
Lines 224-232 Link Here
224	while 1	238	while 1
225	STDERR.printf "."	239	STDERR.printf "."
226	sleep 10	240	sleep 10
227	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' > /tmp/_res#{$$}"	241	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}"
228	system command	242	system command
229	resp = File.open( "/tmp/_res#{$$}", 'r' )	243	resp = File.open( "#{temp_res}", 'r' )
230	# resp.gets	244	# resp.gets
231	# if $_ =~ /WAITING/ then	245	# if $_ =~ /WAITING/ then
232	# resp.close	246	# resp.close
Lines 247-263 Link Here
247	else	261	else
248	# puts "Not supported"	262	# puts "Not supported"
249	# exit	263	# exit
250	qfp = File.open( "/tmp/_q#{$$}", 'w' )	264	qfp = File.open( "#{temp_qf}", 'w' )
251	qfp.puts "> "	265	qfp.puts "> "
252	qfp.puts inseq[i]	266	qfp.puts inseq[i]
253	qfp.close	267	qfp.close
254	command = blastpath + " -p blastp -e #{eval} -b 1000 -m 7 -i /tmp/_q#{$$} -d #{localdb} > /tmp/_res#{$$}"	268	command = blastpath + " -p blastp -e #{eval} -b 1000 -m 7 -i #{temp_qf} -d #{localdb} > #{temp_res}"
255	system command	269	system command
256	resp = File.open( "/tmp/_res#{$$}", 'r' )	270	resp = File.open( "#{temp_res}", 'r' )
257	end	271	end
258	STDERR.puts " Done.\n\n"	272	STDERR.puts " Done.\n\n"
259		273
260	resp = File.open( "/tmp/_res#{$$}", 'r' )	274	resp = File.open( "#{temp_res}", 'r' )
261	while 1	275	while 1
262	while resp.gets	276	while resp.gets
263	break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ \|\| $_ =~ /(<Iteration_stat>)/	277	break if $_ =~ /<Hit_id>(.*)<\/Hit_id>/ \|\| $_ =~ /(<Iteration_stat>)/
Lines 310-326 Link Here
310	afp.close	324	afp.close
311		325
312	STDERR.puts "Performing alignment .. "	326	STDERR.puts "Performing alignment .. "
313	system( mafftpath + mafftopt + " /tmp/_af#{$$} > /tmp/_bf#{$$}" )	327	system( mafftpath + mafftopt + " #{temp_af} > #{temp_bf}" )
314	STDERR.puts "done."	328	STDERR.puts "done."
315		329
316	bfp = File.open( "/tmp/_bf#{$$}", 'r' )	330	bfp = File.open( "#{temp_bf}", 'r' )
317	outseq = []	331	outseq = []
318	outnam = []	332	outnam = []
319	readfasta( bfp, outnam, outseq )	333	readfasta( bfp, outnam, outseq )
320	bfp.close	334	bfp.close
321		335
322
323
324	outseq2 = []	336	outseq2 = []
325	outnam2 = []	337	outnam2 = []
326		338
Lines 356-360 Link Here
356	puts ">" + outnam2[i]	368	puts ">" + outnam2[i]
357	puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )	369	puts outseq2[i].gsub( /.{1,60}/, "\\0\n" )
358	end	370	end
359
360	system( "rm -rf /tmp/_if#{$$} /tmp/_vf#{$$} /tmp/_af#{$$} /tmp/_bf#{$$} /tmp/_pf#{$$} /tmp/_q#{$$} /tmp/_res#{$$} /tmp/_rid#{$$}" )