Blastp in linux

Here is python code for protein sequence blast in linux.

#!/usr/bin/python3

#now import modules to be used
import sys
import subprocess

#now get the filename and protein database name from the command line, the first argument
# str forces it to be a string, not a number or something else
file_name = str(sys.argv[1])
pro_db_name = str(sys.argv[2])

print(file_name, pro_db_name)

#now for neetness set up a subroutine to execute blast
#define function blastp_outfmt7 to blast in format 7
def blastp_outfmt7(seq, pro_db_name):
  cmdd = ("blastp -query " +  seq + " -db " + pro_db_name + " -outfmt 7")
  pss = subprocess.Popen(cmdd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
  output = str(pss.communicate()[0], 'ascii').split('\n')
  return output
#end def blastp_outfmt7

#define function blastp_outfmt6 to blast in format 6
def blastp_outfmt6(seq, pro_db_name):
    cmdd = ("blastp -query " +  seq + " -db " + pro_db_name + " -outfmt 6")
    pss = subprocess.Popen(cmdd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
    output = str(pss.communicate()[0], 'ascii').split('\n')
    return output
#end def blastp_outfmt6


#define function to find isoelectric point for each hit in database
def find_Isoelectric(seq, pro_db_name):
  cmdd =("seqret " + pro_db_name + ":" + seq  + 
    " -stdout -auto | pepstats -filter -auto -stdout | grep Isoe")
  pss = subprocess.Popen(cmdd, shell = True, stdout = subprocess.PIPE, stderr = subprocess.STDOUT)
  output = str(pss.communicate()[0], 'ascii').split('\n')
  return output
#end def find_Isoelectric


#main

# run blast
# in format 7
blastresult = blastp_outfmt7(file_name, pro_db_name) 

# in format 6
blastresult6 = blastp_outfmt6(file_name, pro_db_name)
 

#print results
print("blast of", file_name, "against Protein", pro_db_name, ":")

#print blast result
for i in range (len(blastresult)):
    print (blastresult[i])

#print isoelectric point
for i in range (len(blastresult6)):
    stuff = blastresult6[i].split( "\t")
    if(len(stuff)>1):
       gene_id = stuff[1]
       Isoelec = find_Isoelectric(gene_id, pro_db_name)
       #print (gene_id, Isoelec)
       print (gene_id, Isoelec[0])