Parse BLAST output against ChEMBL, to have the top hits for each query protein

format_blastp_output_for_chembl.py

import os
import sys
from collections import defaultdict
import FiftyHG_Chembl

#====================================================================#

def main():

    # find the blast output files:  
    blastoutput = defaultdict()
    mydir = os.getcwd() # the current directory
    subdirs = os.listdir(mydir) # there may be blast searches against ChEMBL for several query pathogen species
    for subdir in subdirs:
        species = subdir # e.g. schistosoma_mansoni
        subdir = os.path.join(mydir,subdir) 
        myfile = os.path.join(subdir,'blast.out')
        if os.path.exists(myfile):
            assert(species not in blastoutput)
            blastoutput[species] = myfile 

    # read the chembl fasta file of drug targets:
    chembl_file = 'none'
    # get a list of files in mydir:
    files = os.listdir(mydir)
    for file in files:
        #file = os.path.join(mydir,file)
        # check if it is the ChEMBL file (starting with chembl_ and ending with .fa):
        if str(file).endswith('.fa') and str(file).startswith('chembl_'):
            # check the file exists (and is not a broken symbolic link):
            if os.path.exists(file):
                chembl_file = file
    assert(chembl_file != 'none')
    (species_name,protein_seq) = FiftyHG_Chembl.read_chembl_file(chembl_file) 

    # format blastp output for each species:
    for species in blastoutput: 
        myfile = blastoutput[species]
        # parse and format this blast output file:
        output_file = '%s.txt' % species
        print('Making file',output_file)
        if not os.path.exists(output_file):  
            FiftyHG_Chembl.reformat_blast_output(myfile,species_name,protein_seq,output_file)                 

    print("FINISHED")
 
#====================================================================#

if __name__=="__main__":
    main()

#====================================================================#