RunningTBLASTN many times with python

Question

[trini1god]

Hello, I am trying to align a protein sequences to a nucleotide database so am using TBLASTN. Each protein sequence is contained in different files (like over 2000 files). So I opened the directory using glob and want to loop through each file and run the TBLASTN on them but it is giving me this error:

Command line argument error: Argument "query". File is not accessible:  `conSeqfile'

Below is my code:

import os
from os import path
import glob

datadir = '/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Muscle_files/Consensus'
os.chdir(datadir)
readfilenames = glob.glob(os.path.join(datadir, '*.fasta'))

db = '/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Mesculenta_305_v6.1.cds.fa'

#A directory for the output data
outdir = '/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Blast_output'


for conSeqfile in readfilenames:
     print conSeqfile      #printed here to make it is printing the file

    out_blast = path.join(outdir, "blastres_" + conSeqfile[4:] + ".txt")

    cmd = "tblastn -query conSeqfile -db db -out out_blast -outfmt 6 -evalue 0.00001"
    os.system(cmd)

I have been on it for hours and also sourced for similar problems online but still cant figure it out. Any help is appreciated. Thanks

[Andrzej Zielezinski]

Change the line:

cmd = "tblastn -query conSeqfile -db db -out out_blast -outfmt 6 -evalue 0.00001"

into this:

cmd = "tblastn -query {} -db {} -out {} -outfmt 6 -evalue 0.00001".format(conSeqfile, db, out_blast)

[trini1god]

Thanks @a.zielezinski, your answer helped me. I changed it and it worked

Answer 1

[st.ph.n]

You're on the right track, and your code can be simplified:

#!/usr/bin/env python

import os, glob, subprocess

os.chdir('/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Muscle_files/Consensus')

db = '/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Mesculenta_305_v6.1.cds.fa'

outfile = "/home/chisom/somtee/chisom/multifasta_files/multi_fasta/Blast_output/blastres_" + conSeqfile[4:] + ".txt"

for file in glob.glob("*.fasta"):
    print "Running tBLASTn for: " + file

    subprocess.call("tblastn -query " + file + " -db " + db + " -outfmt 6 -evalue 0.00001 -out " + outfile, shell=True)

Answer 2

[Joe]

Do you want to use python for this specifically? This would be less cumbersome as a simple shell script in my opinion. It looks like you're just trying to start the same command on several input fastas?

I'd just do (without any attempt to parallelise):

for file in /path/to/files/*.fasta ; do tblastn -query "$file" -db ~/somtee/chisom/multifasta_files/multi_fasta/Mesculenta_305_v6.1.cds.fa -out /path/to/outputdir/"${file%.*}"_out_blast.txt -outfmt 6 -evalue 0.00001 ; done

Answer 3

[Pierre Lindenbaum]

not python, but a simple Makefile. This can be parallelized using

make -j number_of_jobs

Answer 4

[Joe]

Invoke your python command via GNU parallel or explore the threading and multiprocessing modules within python if you want to handle it within the script.