Mercurial > repos > guerler > springsuite

comparison spring_cross.py @ 39:172398348efd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 26b4018c88041ee0ca7c2976e0a012015173d7b6-dirty"
author guerler
date Fri, 22 Jan 2021 15:50:27 +0000
parents 80a4b98121b6
children
comparison
equal deleted inserted replaced
38:80a4b98121b6 39:172398348efd
1 #! /usr/bin/env python3 1 #! /usr/bin/env python3
2 import argparse 2 import argparse
3 from os import system 3 from os import mkdir
4 from os.path import isdir
4 5
5 from spring_package.DBKit import createFile 6 from spring_package.DBKit import DBKit
6 from spring_package.Molecule import Molecule 7 from spring_package.Molecule import Molecule
8 from spring_package.Utilities import getName
7 9
8 10
9 def getId(line): 11 def hasInterface(mol, chainA, chainB, distance=10.0, contacts=5):
10 line = line.strip() 12 count = 0
11 if len(line) != 6 or line[4:5] != "_": 13 distance = distance ** 2
12 raise Exception("Invalid list entry (`PDB_CHAIN`): %s." % line) 14 for residueA in mol.calpha[chainA]:
13 return line[:4].upper() + line[4:6] 15 atomA = mol.calpha[chainA][residueA]
16 for residueB in mol.calpha[chainB]:
17 atomB = mol.calpha[chainB][residueB]
18 dist2 = ((atomA["x"] - atomB["x"]) ** 2 +
19 (atomA["y"] - atomB["y"]) ** 2 +
20 (atomA["z"] - atomB["z"]) ** 2)
21 if dist2 < distance:
22 count = count + 1
23 if count >= contacts:
24 return True
25 return False
14 26
15 27
16 def main(args): 28 def main(args):
17 logFile = open(args.log, "w") 29 logFile = open(args.log, "w")
18 system("mkdir -p %s" % args.temp) 30 if not isdir("temp"):
31 mkdir("temp")
19 pdbCount = 0 32 pdbCount = 0
20 partnerList = set() 33 partnerList = set()
21 entries = list() 34 entries = set()
22 with open(args.list) as file: 35 with open(args.index) as file:
23 for line in file: 36 for line in file:
24 entries.append(getId(line)) 37 entries.add(getName(line))
25 logFile.write("Found %s template entries.\n" % len(entries)) 38 logFile.write("Found %s template entries.\n" % len(entries))
26 for entryId in entries: 39 pdbDatabase = DBKit(args.index, args.database)
27 pdb = entryId[:4].lower() 40 for pdb in sorted(entries):
28 pdbChain = entryId[5:6] 41 print("Processing %s" % pdb)
29 pdbFile = "%s/temp.pdb" % args.temp 42 pdbFile = "temp/temp.pdb"
30 pdbDatabaseId = "%s.pdb" % pdb 43 pdbDatabaseId = "%s.pdb" % pdb
31 createFile(pdbDatabaseId, args.index, args.database, pdbFile) 44 pdbDatabase.createFile(pdbDatabaseId, pdbFile)
32 try: 45 try:
33 mol = Molecule(pdbFile) 46 mol = Molecule(pdbFile)
34 except Exception: 47 except Exception as e:
35 logFile.write("Warning: File '%s' not found.\n" % pdbDatabaseId) 48 logFile.write("Warning: Entry '%s' not found. %s.\n" % (pdbDatabaseId, str(e)))
36 continue 49 continue
37 pdbCount = pdbCount + 1 50 pdbCount = pdbCount + 1
38 logFile.write("Processing %s, chain %s.\n" % (pdb, pdbChain)) 51 for pdbChain in mol.calpha.keys():
39 logFile.write("Found %d biomolecule(s).\n" % len(mol.biomol.keys())) 52 logFile.write("Processing %s, chain %s.\n" % (pdb, pdbChain))
40 for biomolNumber in mol.biomol: 53 logFile.write("Found %d biomolecule(s).\n" % len(mol.biomol.keys()))
41 if biomolNumber == 0: 54 for biomolNumber in mol.biomol:
42 logFile.write("Processing biomolecule.\n")
43 bioMolecule = mol
44 else:
45 logFile.write("Processing biomolecule %d.\n" % biomolNumber) 55 logFile.write("Processing biomolecule %d.\n" % biomolNumber)
46 bioMolecule = mol.createUnit(biomolNumber) 56 bioMolecule = mol.createUnit(biomolNumber)
47 nChains = len(bioMolecule.calpha.keys()) 57 nChains = len(bioMolecule.calpha.keys())
48 print("Found %d chain(s)." % nChains) 58 if nChains > 1 and pdbChain in bioMolecule.calpha:
49 if nChains > 1 and pdbChain in bioMolecule.calpha: 59 for bioChain in bioMolecule.calpha:
50 for bioChain in bioMolecule.calpha: 60 if bioChain == pdbChain:
51 if bioChain == pdbChain: 61 continue
52 continue 62 if hasInterface(bioMolecule, pdbChain, bioChain):
53 partnerPdbChain = "%s_%s" % (pdb.upper(), bioChain[:1]) 63 corePdbChain = "%s_%s" % (pdb.upper(), pdbChain[:1])
54 partnerList.add("%s\t%s" % (entryId, partnerPdbChain)) 64 partnerPdbChain = "%s_%s" % (pdb.upper(), bioChain[:1])
55 else: 65 partnerList.add("%s\t%s" % (corePdbChain, partnerPdbChain))
56 logFile.write("Skipping: Chain not found or single chain [%s].\n" % pdbChain) 66 else:
57 logFile.flush() 67 logFile.write("Skipping: Chains have no interface [%s, %s].\n" % (pdbChain, bioChain))
68 else:
69 logFile.write("Skipping: Chain not found or single chain [%s].\n" % pdbChain)
70 logFile.flush()
58 with open(args.output, 'w') as output_file: 71 with open(args.output, 'w') as output_file:
59 for entry in sorted(partnerList): 72 for entry in sorted(partnerList):
60 output_file.write("%s\n" % entry) 73 output_file.write("%s\n" % entry)
61 74
62 75
63 if __name__ == "__main__": 76 if __name__ == "__main__":
64 parser = argparse.ArgumentParser(description='List filtering.') 77 parser = argparse.ArgumentParser(description='List filtering.')
65 parser.add_argument('-l', '--list', help='List of PDB chains [PDB_CHAIN]', required=True)
66 parser.add_argument('-i', '--index', help='PDB Database Index file (ffindex)', required=True) 78 parser.add_argument('-i', '--index', help='PDB Database Index file (ffindex)', required=True)
67 parser.add_argument('-d', '--database', help='PDB Database files (ffdata)', required=True) 79 parser.add_argument('-d', '--database', help='PDB Database files (ffdata)', required=True)
68 parser.add_argument('-o', '--output', help='Output file', required=True) 80 parser.add_argument('-o', '--output', help='Output file', required=True)
69 parser.add_argument('-t', '--temp', help='Temporary Directory', required=True)
70 parser.add_argument('-g', '--log', help='Log File', required=True) 81 parser.add_argument('-g', '--log', help='Log File', required=True)
71 args = parser.parse_args() 82 args = parser.parse_args()
72 main(args) 83 main(args)