File: Ensemble_phimap.py

package info (click to toggle)
pdb2pqr 1.8-1
  • links: PTS, VCS
  • area: main
  • in suites: wheezy
  • size: 8,576 kB
  • sloc: python: 24,897; sh: 12,005; cpp: 9,831; xml: 9,098; makefile: 355; ansic: 36
file content (322 lines) | stat: -rw-r--r-- 10,149 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
#!/usr/bin/env python
#
# This is a module in development for making conformationally averaged PBE maps
#
debug=False
import sys, os

print __file__
import os
try:
    file_name=__file__
    if file_name[:2]=='./':
        scriptpath=os.getcwd()
    else:
        scriptpath=os.path.join(os.getcwd(),os.path.split(file_name)[0])
        if scriptpath[-1] == "/":
            scriptpath=scriptpath[:-1]
except:
    scriptpath=os.path.split(sys.argv[0])[0]
    if scriptpath=='.':
        scriptpath=os.getcwd()
#
# Add to import path
#
pdb2pqr_path=os.path.split(scriptpath)[0]
sys.path.append(pdb2pqr_path)

import string
import math
import string
import getopt
import time
import shutil
from src.pdb import *
from src.utilities import *
from src.structures import *
from src.definitions import *
from src.forcefield import *  
from src.routines import *
from src.protein import *
from src.server import *
from StringIO import *
from src.hydrogens import *

class conf_avg:

    def __init__(self,options):
        """Initialize class and decide which kind of job to do"""
        self.options=options # Store options so we can access them anywhere
        potentials=[]
		# If directoryPath is specified then use that, otherwise use pdbfilename
        if options.directoryPath!='':
            listOfFiles=os.listdir(options.directoryPath)
            for currentPDB in listOfFiles:
                currentPDB=os.path.join(options.directoryPath,currentPDB) # Jens added this fix
                pots=self.process_one_pdb(currentPDB)
                potentials.append(pots)
        else:
            # Single file
            potentials.append(self.process_one_pdb(os.path.join(os.getcwd(),options.pdbfilename)))
        avgPotList=self.avg_pots()
        self.writeAvgPots(avgPotList)
        print "done"
        return
        
    #
    # ------
    #
    
    def process_one_pdb(self,pdbfilename):
        """Do everything for one input file"""
        print "Working on: %s" %pdbfilename
        pdbfile = getPDBFile(pdbfilename)
        
        if self.options.MD:
            #
            # Run an MD simulation for this PDB file and calculate potentials for all the snapshots
            #
            snapshots=self.run_MD(pdbfilename)
        else:
            snapshots=[pdbfilename]
        #
        # Get the potentials for everything
        #
        potentials=[]
        for pdbname in snapshots:
            pots=self.get_potentials(pdbfilename)
            potentials.append(pdbname)
        return potentials
        
    #
    # ------
    #
    
    def run_MD(self,pdbfilename):
        """Run an MD simulation and return a number of snapshots"""
        files=os.listdir(os.getcwd())
        addfiles=[pdbfilename]
        #for file in files:
            #addfiles.append(os.path.join(os.getcwd(),file))
        import Gromacs_class as Gclass
        G=Gclass.GROMACS(addfiles)
        #
        # Create Gromacs input file
        # 
        pdbfile=os.path.split(pdbfilename)[1]
        G.pdb2gmx(pdbfile,forcefield=1,ignore_Hs=True,auto_select_his=True)
        #
        # Set up the simulation box. The argument gives the distance between the box edges
        # and the protein in nm
        #
        G.center(1.5)
        G.solvate()
        #
        # Energy minimise and do the pre-MD
        #
        G.EM(2000,1000)
        params={}
        params['ref_t']='%d' %options.temperature
        G.PR_MD(250)
        # ----------------------------------
        # 500 ps timestep is 2 fs
        G.MD(self.options.MDtime*500,params)
        #
        # Get the snapshots
        #
        filenames=G.get_snapshots(self.options.numsnapshots)
        return filenames
        
        
    def get_potentials(self,currentPDB):
		"""Get the potentials by first running pdb2pqr and then apbs"""
		myProtein,apbs_inputfile=self.run_pdb2pqr(currentPDB)
		potentials=self.run_apbs(myProtein,apbs_inputfile)
		return potentials
        
        
    def run_pdb2pqr(self,currentPDB):
        """Run pdb2pqr, prepare input for apbs"""
        pdbfile = getPDBFile(currentPDB)
        pdblist, errlist = readPDB(pdbfile)
        #
        # Instantiate pdb2pqr
        #
        myDefinition = Definition()
        myProtein = Protein(pdblist, myDefinition)

        #
        # Setup everything
        #
        myRoutines = Routines(myProtein, verbose)
        myRoutines.updateResidueTypes()
        myRoutines.updateSSbridges()
        myRoutines.updateBonds()
        myRoutines.setTermini()
        myRoutines.updateInternalBonds()

        myforcefield=Forcefield(ff, myDefinition, None)
        myRoutines.applyNameScheme(myforcefield)

        myRoutines.findMissingHeavy()
        myRoutines.addHydrogens()
        myRoutines.debumpProtein()
        myProtein.reSerialize()
        #
        # Add and optimze hydrogens:
        # 
        from src.hydrogens import hydrogenRoutines
        myRoutines.updateInternalBonds()
        myRoutines.calculateDihedralAngles()
        myhydRoutines = hydrogenRoutines(myRoutines)
        #
        # Now optimize hydrogens
        #
        myhydRoutines.setOptimizeableHydrogens()
        myhydRoutines.initializeFullOptimization()
        myhydRoutines.optimizeHydrogens()
        myhydRoutines.cleanup()
        myRoutines.setStates()

        print "Created protein object (after processing myRoutines) -"
        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()

        #
        # Assign charges
        #
        for chain in myProtein.getChains():
            for residue in chain.get("residues"):
                for atom in residue.get("atoms"):
                    atomname = atom.get("name")
                    charge, radius = myforcefield.getParams1(residue, atomname)
                    atom.set("radius", radius)
                    atom.set("ffcharge", charge)
        #
        #
		method=""
		async=0
		split=0
		import pdb2pka.inputgen_pKa as IP
		igen = IP.inputGen(currentPDB)
		igen.maps=None
		igen.set_type('background')
		igen.pdie=8.0
		igen.sdie=80.0
		all_center,extent=igen.getCenter()
		igen.setfineCenter(all_center)
		print 'Center: %5.1fA %5.1fA %5.1fA' %(all_center[0],all_center[1],all_center[2])
		print 'Extent: %5.1fA %5.1fA %5.1fA'  %(extent[0],extent[1],extent[2])

		apbs_inputfile=igen.printInput()
		return myProtein, apbs_inputfile
        
    def run_apbs(self,myProtein,apbs_inputfile):
		"""runs apbs"""
		import pdb2pka.apbs 
		APBS=pdb2pka.apbs.runAPBS()
		potentials = APBS.runAPBS(myProtein, apbs_inputfile)
		# copies snapshots to separate directories so that they are not overwritten
		global run_no
		mydir="snapshot"+str(run_no)
		if not os.path.exists(mydir):
			os.mkdir(mydir)
		for i in range(0,4):
			myfile="potential"+str(i)+".dx"
			mypath=mydir+'/'+myfile
			shutil.copyfile(myfile,mypath)
		run_no+=1
		APBS.cleanup()
		return potentials

    def avg_pots(self):
		"""reads in potentials and averages them"""
		topdir=os.getcwd()
		for snapshot_no in range(0,run_no):
			os.chdir("snapshot"+str(snapshot_no))
			currPotFile=open("potential0.dx","r")
			print "starting snapshot no.", snapshot_no
			#discard first lines
			for j in range(0,10):
				currPotFile.readline()

			#get number of data items (in case a different grid than 65x65x65 is used)
			dataItems=currPotFile.readline().split()[9]
			dataLines=int(math.ceil(float(dataItems)/3))

			#initialize the list that will store current list of potentials read from the file
			currPotList=[]
			for j in range(0,dataLines):
				currLine=currPotFile.readline().split()
				currPotList.append(currLine)
				
			#flatten the currPotList
			flatCurrPotList=[item for sublist in currPotList for item in sublist]
			
			#initialize current sum list (only on the first go)
			if snapshot_no == 0:
				currSumPotList=[0]*len(flatCurrPotList)
			
			#add current potential list to sum of previous potential lists
			for k in range(0,len(flatCurrPotList)):
				currSumPotList[k]+=float(flatCurrPotList[k])

			os.chdir(topdir)
			currPotFile.close()

			print "done with snapshot no. ", snapshot_no
		
		#initialize averaged potentials list
		avgPotList=[0]*len(flatCurrPotList)
		#do the averaging
		for j in range(0,len(currSumPotList)):
			avgPotList[j]=float(currSumPotList[j])/(snapshot_no+1)

		return avgPotList

    def writeAvgPots(self,avgPotList):
		"""writes averaged potentials into a file"""
		avgPotsFile = open('averaged_pots', 'w')
		counter=0
		for each in avgPotList:
			avgPotsFile.write(str(each))
			avgPotsFile.write(" ")
			counter+=1
			if counter%3==0:
				avgPotsFile.write("\n")
		avgPotsFile.close()
		return


#
# ----
#

if __name__=='__main__':
    from optparse import OptionParser
    parser = OptionParser(usage='%prog [options]',version='%prog 1.0')
    parser.add_option('-p','--pdb',dest='pdbfilename',action='store',type='string',default='2lzt.pka.pdb',help='The PDB file. Default: %default')
    parser.add_option('-d','--dir',dest='directoryPath',action='store',type='string',default='',
                  help='Direcotry of the PDB files/snapshots. Default: %default')
    #
    # Flags
    #
    parser.add_option('--MD',dest='MD',action='store_true',default=False,help='Perform an MD simulation and use snapshots for calculating electrostatic potential')
    parser.add_option('-s','--MDsnapshots',dest='numsnapshots',action='store',type='int',default=100,help='Number of MD snapshots to use. Default: %default')
    parser.add_option('--MDtime',dest='MDtime',action='store',type='int',default=100,help='Time in picoseconds that MD should be run for. Default: %default')
    parser.add_option('-t','--temp',dest='temperature',action='store',type='float',default=310.15,help='Temperature for the MD run. Default: %default')
    #
    # We can think about adding flags for not solvating the structure etc here
    #

	# A global value keeping track of how many snapshot directories to create
    run_no=0


    (options, args) = parser.parse_args()

    verbose=True
    ff='parse'

    I=conf_avg(options)