#! /usr/bin/env ccp4-python # # Copyright (C) 2014 Ronan Keegan # # This code is distributed under the terms and conditions of the # CCP4 Program Suite Licence Agreement as a CCP4 Application. # A copy of the CCP4 licence can be obtained by writing to the # CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK. # # A MrBUMP wrapper for HHpred # Ronan Keegan - 26/2/14 # import os, sys import time import subprocess import shlex import shutil import string class Hit: def __init__(self): self.chainName="" self.pdbName="" self.chainID="" self.rank="" self.prob=0.0 self.evalue=0.0 self.pvalue=0.0 self.score=0.0 self.cols=0 self.alignments=dict([]) class HHpred: def __init__(self): self.hhblitsEXE="hhblits" self.hhsearchEXE="hhsearch" self.cpu=1 self.seqin="" self.seqout="" self.uniprotDB="" self.oa3mFile="" self.ina3mFile="" self.logfile="" self.iterations=1 self.termination=False self.workingDIR="" self.resultsList=[] self.resultsDict=dict([]) def runHHblits(self, seqin, uniprotDB, oa3mFile, cpu=1, iterations=1, logfile="hhblits.log", debug=False): """ Run HHblits """ self.cpu=cpu self.iterations=iterations self.seqin=seqin self.uniprotDB=uniprotDB self.oa3mFile=oa3mFile self.logfile=logfile self.termination=False #hhblits -cpu 4 -i seq.pir -d ${HHDB}/hhsuite_dbs/uniprot20/uniprot20_2013_03 -oa3m query.a3m -n 1 # Set the command line command_line=self.hhblitsEXE + " -cpu %d -i %s -d %s -oa3m %s -n %d" \ % (self.cpu, self.seqin, self.uniprotDB, self.oa3mFile, self.iterations) if debug == True: sys.stdout.write("HHblits command line:\n %s\n" % command_line) sys.stdout.write("\n") # Launch program if os.name == "nt": process_args = shlex.split(command_line, posix=False) p = subprocess.Popen(process_args, bufsize=0, shell="False", stdin = subprocess.PIPE, stdout = subprocess.PIPE, stderr = subprocess.STDOUT) else: process_args = shlex.split(command_line) p = subprocess.Popen(process_args, stdin = subprocess.PIPE, stdout = subprocess.PIPE, stderr = subprocess.STDOUT) (child_stdout, child_stdin) = (p.stdout, p.stdin) # Write the input child_stdin.close() # Open the log file for writing log=open(self.logfile, "w") # Watch the output for successful termination out=child_stdout.readline() while out: if debug == True: sys.stdout.write(out) if 'Writing A3M alignment' in out: self.termination=True log.write(out) log.flush() out=child_stdout.readline() child_stdout.close() log.close() def runHHsearch(self, seqout, pdbDB, ina3mFile, cpu=1, logfile="hhsearch.log", debug=False): """ Run HHsearch """ self.cpu=cpu self.seqout=seqout self.pdbDB=pdbDB self.ina3mFile=ina3mFile self.logfile=logfile self.termination=False CAPTURE=False #hhsearch -Ofas query.fasta -cpu 4 -i query.a3m -d ${HHDB}/hhsuite_dbs/pdb70/pdb70_08Feb14_hhm_db # Set the command line command_line=self.hhsearchEXE + " -Ofas %s -cpu %d -i %s -d %s" \ % (self.seqout, self.cpu, self.ina3mFile, self.pdbDB) if debug == True: sys.stdout.write("HHsearch command line:\n %s\n" % command_line) sys.stdout.write("\n") # Launch program if os.name == "nt": process_args = shlex.split(command_line, posix=False) p = subprocess.Popen(process_args, bufsize=0, shell="False", stdin = subprocess.PIPE, stdout = subprocess.PIPE, stderr = subprocess.STDOUT) else: process_args = shlex.split(command_line) p = subprocess.Popen(process_args, stdin = subprocess.PIPE, stdout = subprocess.PIPE, stderr = subprocess.STDOUT) (child_stdout, child_stdin) = (p.stdout, p.stdin) # Write the input child_stdin.close() # Open the log file for writing log=open(self.logfile, "w") # Watch the output for successful termination out=child_stdout.readline() while out: if debug == True: sys.stdout.write(out) if 'Done' in out: self.termination=True m=dict([]) if CAPTURE: hit=string.split(out) if len(hit) > 1: self.resultsList.append(hit[1]) self.resultsDict[hit[1]]=Hit() self.resultsDict[hit[1]].rank=hit[0] self.resultsDict[hit[1]].chainName=hit[1] self.resultsDict[hit[1]].pdbName=hit[1][0:4] if len(hit[1])>=6: self.resultsDict[hit[1]].chainID=hit[1][5:] else: self.resultsDict[hit[1]].chainID="A" self.resultsDict[hit[1]].prob=float(out[35:40]) self.resultsDict[hit[1]].score=float(out[57:63]) #self.resultsDict[hit[1]].prob=float(hit[len(hit)-9]) #self.resultsDict[hit[1]].score=float(hit[len(hit)-6]) #self.resultsDict[hit[1]].evalue=float(hit[len(hit)-8]) #self.resultsDict[hit[1]].pvalue=float(hit[len(hit)-7]) #self.resultsDict[hit[1]].cols=float(hit[len(hit)-4]) else: CAPTURE=False if "No Hit" in out and "E-value" in out: CAPTURE=True log.write(out) log.flush() out=child_stdout.readline() child_stdout.close() # Scan the output fasta file for the sequence alignments align=[] seqDict=dict([]) seq="" fastaFile=open(seqout, "r") line=fastaFile.readline() while line: if "# No" in line: if len(align) != 0: count=0 for i in align: if ">" in i: if seq != "": seqDict[title]=seq seq="" title=i if count==0: title=title+"T" if count==1: title=">Target_alignment" count=count+1 else: seq=seq+i if seq != "": seqDict[title]=seq ModelName=string.strip(title).replace(">","") self.resultsDict[ModelName].alignments=seqDict seq="" count=count+1 align=[] seqDict=dict([]) else: align.append(string.strip(line)) line=fastaFile.readline() if len(align) != 0: count=0 for i in align: if ">" in i: if seq != "": seqDict[title]=seq seq="" title=i if count==0: title=title+"T" if count==1: title=">Target_alignment" count=count+1 else: seq=seq+i if seq != "": seqDict[title]=seq ModelName=string.strip(title).replace(">","") self.resultsDict[ModelName].alignments=seqDict seq="" count=count+1 if __name__ == "__main__": seqin = "seq.pir" seqout = "results.fasta" uniprotDB = "/home/rmk65/Projects/hhpred/uniprot20_2013_03/uniprot20_2013_03" pdbDB = "/home/rmk65/Projects/hhpred/pdb70/pdb70_10Apr14_hhm_db" oa3mFile = "test.oa3m" hh=HHpred() hh.runHHblits(seqin=seqin, uniprotDB=uniprotDB, oa3mFile=oa3mFile) hh.runHHsearch(seqout=seqout, ina3mFile=oa3mFile, pdbDB=pdbDB) print hh.resultsList print hh.resultsDict["1wlt_A"].alignments[">1wlt_A"]