Mercurial > repos > guerler > springsuite

comparison spring_model.py @ 17:c790d25086dc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit b0ede77caf410ab69043d33a44e190054024d340-dirty"
author guerler
date Wed, 28 Oct 2020 05:11:56 +0000
parents
children
comparison
equal deleted inserted replaced
16:16eb2acaaa20 17:c790d25086dc
1 #! /usr/bin/env python3
2 import argparse
3 import os
4
5 from spring_package.Alignment import Alignment
6 from spring_package.Energy import Energy
7 from spring_package.Molecule import Molecule
8
9 def buildModel(resultFile, templateFile, chainName, outputName):
10 template = Molecule(templateFile)
11 if chainName not in template.calpha:
12 raise Exception("Chain not found in template [%s]" % chainName)
13 chain = template.calpha[chainName]
14 alignment = Alignment(resultFile)
15 alignment.createModel(chain)
16 template.saveChain(chainName, outputName)
17 os.system("./build/pulchra %s" % outputName)
18
19 def TMalign(fileA, fileB):
20 baseA = os.path.basename(fileA)
21 baseB = os.path.basename(fileB)
22 baseA = os.path.splitext(baseA)[0]
23 baseB = os.path.splitext(baseB)[0]
24 tmName = "temp/tmalign.%s.%s" % (baseA, baseB)
25 os.system("build/TMalign %s %s -m %s.mat > %s.out" % (fileA, fileB, tmName, tmName))
26 rotmat = list()
27 with open("%s.mat" % tmName) as file:
28 line = next(file)
29 line = next(file)
30 for i in range(3):
31 line = next(file)
32 rotmatLine = line[1:].split()
33 rotmatLine = list(map(lambda x: float(x), rotmatLine))
34 rotmatLine = [rotmatLine[1], rotmatLine[2], rotmatLine[3], rotmatLine[0]]
35 rotmat.append(rotmatLine)
36 with open("%s.out" % tmName) as file:
37 for i in range(14):
38 line = next(file)
39 try:
40 tmscore = float(line[9:17])
41 line = next(file)
42 tmscore = max(tmscore, float(line[9:17]))
43 except:
44 raise Exception("TMalign::Failed to retrieve TMscore.")
45 molecule = Molecule(fileA)
46 for atom in molecule.atoms:
47 molecule.applyMatrix(atom, rotmat)
48 return tmscore, molecule
49
50 def main(args):
51 os.system("mkdir -p temp")
52 os.system("rm -f temp/*.*")
53 buildModel(args.a_result, args.a_template, args.a_chain, "temp/monomerA.pdb")
54 buildModel(args.b_result, args.b_template, args.b_chain, "temp/monomerB.pdb")
55 interfaceEnergy = Energy()
56 templateMolecule = Molecule(args.template)
57 modelCount = 0
58 for biomolNumber in range(len(templateMolecule.rotmat.keys())):
59 os.system("rm -f temp/template*.pdb")
60 if biomolNumber == 0:
61 bioMolecule = templateMolecule
62 else:
63 bioMolecule = templateMolecule.createUnit(biomolNumber)
64 if len(bioMolecule.calpha.keys()) > 1 and args.template_core in bioMolecule.calpha:
65 for chainName in bioMolecule.calpha.keys():
66 bioMolecule.saveChain(chainName, "temp/template%s.pdb" % chainName)
67 coreTMscore, coreMolecule = TMalign("temp/monomerA.rebuilt.pdb", "temp/template%s.pdb" % args.template_core)
68 maxScore = -9999
69 maxMolecule = None
70 for chainName in bioMolecule.calpha.keys():
71 if chainName != args.template_core and len(bioMolecule.calpha[chainName]) > 0:
72 print("Evaluating chain %s..." % chainName)
73 try:
74 partnerTMscore, partnerMolecule = TMalign("temp/monomerB.rebuilt.pdb", "temp/template%s.pdb" % chainName)
75 except:
76 continue
77 minTM = min(coreTMscore, partnerTMscore)
78 print("min-TMscore: %5.5f" % minTM)
79 energy = interfaceEnergy.get(coreMolecule, partnerMolecule)
80 print("Interaction: %5.5f" % energy)
81 springScore = minTM * args.wtm - energy * args.wenergy
82 print("SpringScore: %5.5f" % springScore)
83 if springScore > maxScore:
84 maxMolecule = partnerMolecule
85 modelName = "temp/model%s_%s.pdb" % (modelCount, chainName)
86 coreMolecule.save(modelName, chainName="A")
87 maxMolecule.save(modelName, chainName="B", append=True)
88 modelCount = modelCount + 1
89
90 if __name__ == "__main__":
91 parser = argparse.ArgumentParser(description='Create a 3D model from HH-search results.')
92 parser.add_argument('-ar', '--a_result', help='First HHR target file result', required=True)
93 parser.add_argument('-ac', '--a_chain', help='First template chain name', required=True)
94 parser.add_argument('-at', '--a_template', help='First template PDB', required=True)
95 parser.add_argument('-br', '--b_result', help='Second HHR target file result', required=True)
96 parser.add_argument('-bc', '--b_chain', help='Second structure chain name', required=True)
97 parser.add_argument('-bt', '--b_template', help='Second template PDB', required=True)
98 parser.add_argument('-ts', '--template', help='Structure template', required=True)
99 parser.add_argument('-tc', '--template_core', help='Core template chain name', required=True)
100 parser.add_argument('-o', '--output', help='Output PDB file', required=False)
101 parser.add_argument('-wt', '--wtm', help='Weight TM-score', type=float, default=12.4, required=False)
102 parser.add_argument('-we', '--wenergy', help='Weight Energy term', type=float, default=-0.2, required=False)
103 args = parser.parse_args()
104 main(args)
105