Mercurial > repos > guerler > springsuite

comparison spring_package/Molecule.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 class Molecule:
2 def __init__(self, fileName = None):
3 self.calpha = dict()
4 self.biomol = dict()
5 self.rotmat = dict()
6 self.atoms = list()
7 if fileName is not None:
8 self.fromFile(fileName)
9
10 def fromFile(self, fileName):
11 biomolNumber = 0
12 biomolChains = list()
13 with open(fileName) as file:
14 for index, line in enumerate(file):
15 key = line[0:6].strip()
16 if key == "ATOM":
17 atom = line[12:16]
18 atomNumber = line[6:11]
19 chainName = line[21:22]
20 if chainName not in self.calpha:
21 self.calpha[chainName] = dict()
22 x = self.toFloat(line[30:38])
23 y = self.toFloat(line[38:46])
24 z = self.toFloat(line[46:54])
25 occupancy = self.toFloat(line[54:60], optional=True)
26 temperature = self.toFloat(line[54:60], optional=True)
27 residue = line[17:20]
28 residueNumber = self.toInt(line[22:26])
29 atomNumber = self.toInt(line[6:11])
30 atomName = line[12:16]
31 atomDict = dict(x=x, y=y, z=z,
32 residue=residue,
33 occupancy=occupancy,
34 temperature=temperature,
35 atomNumber=atomNumber,
36 atomName=atomName,
37 residueNumber=residueNumber,
38 chainName=chainName)
39 if atom.strip() == "CA":
40 self.calpha[chainName][residueNumber] = atomDict
41 self.atoms.append(atomDict)
42 biokey = "REMARK 350 BIOMOLECULE:"
43 if line[0:len(biokey)] == biokey:
44 biomolNumber = self.toInt(line[len(biokey):])
45 biokey = "REMARK 350 APPLY THE FOLLOWING TO CHAINS:"
46 nextLine = next(file)
47 while nextLine[:len(biokey)] != biokey:
48 nextLine = next(file)
49 biomolChains = nextLine[len(biokey):].split(",")
50 biomolChains = list(map(lambda x: x.strip(), biomolChains))
51 biokey = "REMARK 350 AND CHAINS:"
52 nextLine = next(file)
53 while nextLine[:len(biokey)] == biokey:
54 moreChains = nextLine[len(biokey):].split(",")
55 moreChains = list(map(lambda x: x.strip(), moreChains))
56 biomolChains = biomolChains + moreChains
57 nextLine = next(file)
58 biokey = "REMARK 350 BIOMT"
59 if nextLine[:len(biokey)] == biokey:
60 biomolMatId1, biomolMat1 = self.getFloats(nextLine)
61 nextLine = next(file)
62 biomolMatId2, biomolMat2 = self.getFloats(nextLine)
63 nextLine = next(file)
64 biomolMatId3, biomolMat3 = self.getFloats(nextLine)
65 if biomolMatId1 != biomolMatId2 or biomolMatId1 != biomolMatId3:
66 raise Exception("Invalid rotation matrix format [%s]." % biomolMatId1)
67 matrix = [biomolMat1, biomolMat2, biomolMat3]
68 biomolChains = [c for c in biomolChains if c]
69 if biomolNumber not in self.rotmat:
70 self.rotmat[biomolNumber] = list()
71 self.rotmat[biomolNumber].append(dict(chains=biomolChains, matrix=matrix))
72 removeChains = []
73 for chainName in self.calpha:
74 if len(self.calpha[chainName]) == 0:
75 removeChains.append(chainName)
76 for chainName in removeChains:
77 del self.calpha[chainName]
78 if not self.calpha:
79 raise Exception("Molecule has no atoms.")
80
81 def getFloats(self, nextLine):
82 matId = self.toInt(nextLine[20:23])
83 matLine = nextLine[23:].split()
84 matLine = list(map(lambda x: self.toFloat(x), matLine))
85 return matId, matLine
86
87 def createUnit(self, biomolNumber = 1):
88 molecule = Molecule()
89 chainCount = 0
90 for matrixDict in self.rotmat[biomolNumber]:
91 for chain in matrixDict["chains"]:
92 if chain in self.calpha:
93 chainCopy = dict()
94 for residue in self.calpha[chain]:
95 chainCopy[residue] = self.calpha[chain][residue].copy()
96 for atomNumber in chainCopy:
97 atom = chainCopy[atomNumber]
98 rotmat = matrixDict["matrix"]
99 self.applyMatrix(atom, rotmat)
100 if chain in molecule.calpha:
101 chainName = chainCount
102 else:
103 chainName = chain
104 molecule.calpha[chainName] = chainCopy
105 chainCount = chainCount + 1
106 return molecule
107
108 def applyMatrix(self, atom, rotmat):
109 newx = atom["x"] * rotmat[0][0] + atom["y"] * rotmat[0][1] + atom["z"] * rotmat[0][2] + rotmat[0][3]
110 newy = atom["x"] * rotmat[1][0] + atom["y"] * rotmat[1][1] + atom["z"] * rotmat[1][2] + rotmat[1][3]
111 newz = atom["x"] * rotmat[2][0] + atom["y"] * rotmat[2][1] + atom["z"] * rotmat[2][2] + rotmat[2][3]
112 atom["x"] = newx
113 atom["y"] = newy
114 atom["z"] = newz
115 return atom
116
117 def toFloat(self, x, optional=False):
118 try:
119 return float(x)
120 except:
121 if not optional:
122 raise Exception("Invalid float conversion [%s]." % x)
123 return 0.0
124
125 def toInt(self, x):
126 try:
127 return int(x)
128 except:
129 raise Exception("Invalid integer conversion [%s]." % x)
130
131 def saveChain(self, chainName, outputName):
132 print ("Writing PDB file to %s." % outputName)
133 f = open(outputName, "w")
134 for residueNumber in sorted(self.calpha[chainName].keys()):
135 ca = self.calpha[chainName][residueNumber]
136 if ca["residue"] is not None:
137 f.write(self.atomString(ca))
138 f.close()
139
140 def save(self, outputName, append=False, chainName=None):
141 print ("Writing atoms to PDB file to %s." % outputName)
142 fileFlag = "+a" if append else "w"
143 f = open(outputName, fileFlag)
144 for atom in self.atoms:
145 atom["chainName"] = chainName if chainName else atom["chainName"]
146 f.write(self.atomString(atom))
147 f.write("TER\n")
148 f.close()
149
150 def atomString(self, atom):
151 return "ATOM %5d %s %s %s%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n" % (atom["atomNumber"], atom["atomName"], atom["residue"], atom["chainName"], atom["residueNumber"], atom["x"], atom["y"], atom["z"], atom["occupancy"], atom["temperature"])
152