Skip to content

Added functionality to read CONECT entries from PDB. #55

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pbuslaev wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Added functionality to read CONECT entries from PDB. #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,4 @@ build
dist
*.egg-info
*.ipynb_checkpoints
.eggs/*
66 changes: 36 additions & 30 deletions src/pmx/ligand_alchemy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
# ================
# Helper functions
# ================
def reformatPDB(fname,num,randint=42,bStrict=False):
def reformatPDB(fname,num,randint=42,bStrict=False,bCONECT=False):
"""Higher level command to read, format and call a pdb writer.

Params
Expand All @@ -35,6 +35,8 @@ def reformatPDB(fname,num,randint=42,bStrict=False):
random number to mark the output pdb
bStrict : bool
limit atom names to 3 char (default False)
bCONECT : bool
read CONECT entries from pdb (default False)
Returns
-------
newname : str
Expand All @@ -46,7 +48,7 @@ def reformatPDB(fname,num,randint=42,bStrict=False):
"""

newname = "tempFormat_"+str(randint)+'_'+str(num)+".pdb"
m = Model().read(fname)
m = Model().read(fname,bCONECT=bCONECT)

# adjust atom names and remember the changes
atomNameID = {}
Expand All @@ -61,10 +63,11 @@ def reformatPDB(fname,num,randint=42,bStrict=False):
atomNameID[a.id] = a.name
a.name = newAtomName

writeFormatPDB(newname,m,bStrict=bStrict)
writeFormatPDB(newname,m,bStrict=bStrict,
_conect_entries=m._pdb_conect if bCONECT else None)
return(newname,atomNameID,sigmaHoleID)

def writeFormatPDB(fname,m,title="",nr=1,bStrict=False):
def writeFormatPDB(fname,m,title="",nr=1,bStrict=False,_conect_entries=None):
"""Writes formatted pdb of a ligand.

Params
Expand All @@ -76,9 +79,11 @@ def writeFormatPDB(fname,m,title="",nr=1,bStrict=False):
title : str
currently not used
nr : int
currently not used
currently not used
bStrict : bool
limit atom names to 3 char (default False)
_conect_entries : list
List of conect entries
Returns
-------
None
Expand All @@ -105,6 +110,9 @@ def writeFormatPDB(fname,m,title="",nr=1,bStrict=False):
print(foo,file=fp)
else:
print(atom,file=fp)
if _conect_entries is not None:
for line in _conect_entries:
fp.write(f"{line}\n")
fp.write('ENDMDL\n')
fp.close()

Expand Down Expand Up @@ -205,7 +213,7 @@ def assignFF(model, itp):
------
model : pmx Model
model
itp : itp TopolBase
itp : itp TopolBase
itp

Returns
Expand All @@ -229,7 +237,7 @@ def readPairsFile(fn):
Params
------
fn : str
filename
filename

Returns
-------
Expand Down Expand Up @@ -302,7 +310,7 @@ def superimposeStructures( fname1, fname2, m3, m2, plst, logfile ):
m3 : pmx Model
model m1 that will get the coordinates of m2
m2 : pmx Model
model m2 that will be fit on m1
model m2 that will be fit on m1
plst : array
pairs list
Returns
Expand Down Expand Up @@ -476,7 +484,7 @@ def _count_ringNonRing( self, n1,n2, mol1,mol2 ):
counter+=1
return(counter)


#*************************************************#
# MCS #
#*************************************************#
Expand Down Expand Up @@ -590,7 +598,7 @@ def mcs(self):
n1_foo = []
n2_foo = []
for n1,n2 in list(zip(n1_list,n2_list)):
n1,n2 = self._checkTop(n1,n2)
n1,n2 = self._checkTop(n1,n2)
n1_foo.append(n1)
n2_foo.append(n2)
n1_list = cp.copy(n1_foo)
Expand Down Expand Up @@ -646,7 +654,7 @@ def mcs(self):
n1,n2 = self._mapH( n1, n2 )

# one more final check for possible issues with the 1-2, 1-3 and 1-4 interactions
n1,n2 = self._checkTop(n1,n2)
n1,n2 = self._checkTop(n1,n2)

# remove sigma hole virtual particles
n1,n2 = self._removeSigmaHoles( n1,n2,self.sigmaHoleID1,self.sigmaHoleID2)
Expand Down Expand Up @@ -710,7 +718,7 @@ def _mapH( self, nfoo, nbar ):
nb1 = a1.GetNeighbors()
nb2 = a2.GetNeighbors()

# only allow mapping hydrogens in those cases,
# only allow mapping hydrogens in those cases,
# when n1 and n2 have an equal number of them bound
nH1 = self._countHydrogens(nb1)
nH2 = self._countHydrogens(nb2)
Expand Down Expand Up @@ -813,18 +821,18 @@ def _checkTop(self, n1, n2 ):
# 2) identify problematic mappings
# 3) fix the problems: discard the atom with fewer mapped neighbours

####### 1-2 #########
####### 1-2 #########
# 1a) 1-2 lists
dict12_mol1 = self._getList12(self.mol1,n1)
dict12_mol2 = self._getList12(self.mol2,n2)
# 2a) identify problems 1-2; and
# 2a) identify problems 1-2; and
# 3a) fix 1-2
rem12_mol2_start,rem12_mol2_end = self._findProblemsExclusions(n1,n2,dict12_mol1,dict12_mol2,case='1-2') # output: indeces of mol2
n2,n1 = self._fixProblemsExclusions(self.mol2,self.mol1,n2,n1,rem12_mol2_start,rem12_mol2_end)
rem12_mol1_start,rem12_mol1_end = self._findProblemsExclusions(n2,n1,dict12_mol2,dict12_mol1,case='1-2') # output: indeces of mol1
n1,n2 = self._fixProblemsExclusions(self.mol1,self.mol2,n1,n2,rem12_mol1_start,rem12_mol1_end)

####### 1-3 #########
####### 1-3 #########
# 1b) 1-3 lists
dict13_mol1 = self._getList13(self.mol1,n1)
dict13_mol2 = self._getList13(self.mol2,n2)
Expand All @@ -835,11 +843,11 @@ def _checkTop(self, n1, n2 ):
rem13_mol1_start,rem13_mol1_end = self._findProblemsExclusions(n2,n1,dict13_mol2,dict13_mol1,case='1-3') # output: indeces of mol1
n1,n2 = self._fixProblemsExclusions(self.mol1,self.mol2,n1,n2,rem13_mol1_start,rem13_mol1_end)

####### 1-4 #########
####### 1-4 #########
# 1b) 1-4 lists
dict14_mol1 = self._getList14(self.mol1,n1)
dict14_mol2 = self._getList14(self.mol2,n2)
# 2b) identify problems 1-4 and
# 2b) identify problems 1-4 and
# 3b) fix 1-4
rem14_mol2_start,rem14_mol2_end = self._findProblemsExclusions(n1,n2,dict14_mol1,dict14_mol2,case='1-4') # output: indeces of mol2
n2,n1 = self._fixProblemsExclusions(self.mol2,self.mol1,n2,n1,rem14_mol2_start,rem14_mol2_end)
Expand Down Expand Up @@ -985,7 +993,7 @@ def _fixProblemsExclusions(self,mol1,mol2,n1,n2,startList,endList):
a,b = self._getAttr(n1,n2,iEnd,bar)
if( (a!=None) and (b!=None) ):
endNeighb = endNeighb+1
elif( iEnd==bar ): # atom of interest
elif( iEnd==bar ): # atom of interest
a,b = self._getAttr(n1,n2,iEnd,foo)
if( (a!=None) and (b!=None) ):
endNeighb = endNeighb+1
Expand Down Expand Up @@ -1048,7 +1056,7 @@ def _disconnectedRecursive(self,mol1,mol2,ind1,ind2):
matchDict1 = {}
matchDict2 = {}
for id1,id2 in list(zip(ind1,ind2)):
# find id1
# find id1
key1 = ''
for i in range(0,len(n1_fragments)):
if id1 in n1_fragments[i]:
Expand Down Expand Up @@ -1710,7 +1718,7 @@ def _distance_based(self, mol1, mol2, id1=None, id2=None, calcOnly=False):
pairs1.append(keep1)
pairs2.append(keep2)
return(pairs1,pairs2)

# mcs
for ind1 in id1:
c1 = mol1.GetConformer()
Expand Down Expand Up @@ -1834,7 +1842,7 @@ def _selectOneMCS(self,n1_list,n2_list,mol1,mol2):
if bOK==True:
n1_largestB.append(nfoo)
n2_largestB.append(nbar)

# of the largest ones select the minRMSD
rmsdMin = 9999.999
for nfoo,nbar in list(zip(n1_largestB,n2_largestB)):
Expand Down Expand Up @@ -1984,8 +1992,8 @@ def _removeH(self, mol1,mol2,nfoo,nbar):
# when heavyAtom1 and heavyAtom2 have an equal number of them bound
heavyAtom1 = a1.GetNeighbors()[0]
heavyAtom2 = a2.GetNeighbors()[0]
nb1 = heavyAtom1.GetNeighbors()
nb2 = heavyAtom2.GetNeighbors()
nb1 = heavyAtom1.GetNeighbors()
nb2 = heavyAtom2.GetNeighbors()
nH1 = self._countHydrogens(nb1)
nH2 = self._countHydrogens(nb2)
if nH1!=nH2:
Expand Down Expand Up @@ -2061,15 +2069,15 @@ def __init__(self, **kwargs):
self.d = 0.05
self.logfile = None
self.bDist = True # this variable is always set to True

for key, val in kwargs.items():
setattr(self,key,val)


def _makeHybridTop( self ):
"""The main function of this class.
Generated hybrid topology

Parameters
----------
None
Expand Down Expand Up @@ -2144,7 +2152,7 @@ def _makeHybridTop( self ):
#####################################
# make exclusions
doLog(self.logfile,"Construct exclusions....")
self.newexclusions = []
self.newexclusions = []
self.bHasExclusions = False
self._make_exclusions()

Expand Down Expand Up @@ -2818,7 +2826,7 @@ def _sum_charge_of_states( self, itp ):

def _atoms_morphe( self, atoms ):
for atom in atoms:
if atom.atomtypeB is not None and (atom.q!=atom.qB or atom.m != atom.mB):
if atom.atomtypeB is not None and (atom.q!=atom.qB or atom.m != atom.mB):
return(True)
return(False)

Expand Down Expand Up @@ -3021,9 +3029,7 @@ def _write_FF_file(atoms,file):
for atype in atoms.keys():
at = atoms[atype]
fp.write(' %s %s %s %s %s %s\n' % (at.type,at.sigmaA,at.epsA,at.A,at.sigmaB,at.epsB) )

def _merge_FF_files( fnameOut, ffsIn=[] ):
atoms = _get_FF_atoms( ffsIn )
_write_FF_file( atoms, fnameOut )


28 changes: 17 additions & 11 deletions src/pmx/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -362,8 +362,10 @@ def make_residues(self):
r.chain = ch
r.chain_id = ch.id

def __readPDB(self, fname=None, pdbline=None):
def __readPDB(self, fname=None, pdbline=None, bCONECT=False):
"""Reads a PDB file"""
if bCONECT:
self._pdb_conect = []
if pdbline:
lines = pdbline.split('\n')
else:
Expand All @@ -374,6 +376,8 @@ def __readPDB(self, fname=None, pdbline=None):
self.atoms.append(a)
if line[:6] == 'CRYST1':
self.box = _p.box_from_cryst1(line)
if bCONECT and line[:6] == "CONECT":
self._pdb_conect.append(line.strip())
self.make_chains()
self.make_residues()
self.unity = 'A'
Expand All @@ -385,14 +389,14 @@ def __check_if_gap( self, atC, atN ):
if atN.name != 'N':
return(False)
d = atC - atN
if d > 1.7: # bond
if d > 1.7: # bond
return(True)
return(False)

def __compareWithoutLastChar(self, str1, str2):
def __compareWithoutLastChar(self, str1, str2):
if isinstance(str1,int) and isinstance(str2,int): # e.g. 52, 53
return(False)

if isinstance(str1,int): # e.g. 52, 52A
if str1==int(str2[0:-1]): # 52, 52A
return(True)
Expand All @@ -419,7 +423,7 @@ def _getChainIDs( self, lines ):
if( a.chain_id not in usedChainIDs ):
usedChainIDs.append( a.chain_id )
return( usedChainIDs )

# TODO: make readPDB and readPDBTER a single function. It seems like
# readPDBTER is more general PDB reader?
def __readPDBTER(self, fname=None, pdbline=None, bNoNewID=True, bPDBGAP=False, bPDBMASS=False):
Expand Down Expand Up @@ -467,7 +471,7 @@ def __readPDBTER(self, fname=None, pdbline=None, bNoNewID=True, bPDBGAP=False, b
bNewChain = True
if (prevAtomName == 'HH33') and ((prevResName=='NME') or (prevResName=='NAC') or (prevResName=='CT3')): # NME cap
bNewChain = True
# do not assign new chain IDs for waters, ions
# do not assign new chain IDs for waters, ions
if (a.resname=='WAT') or (a.resname=='SOL') or (a.resname=='TIP3') or (a.resname=='HOH') \
or (a.resname=='NA') or (a.resname=='CL') \
or (a.resname=='NaJ') or (a.resname=='Na') or (a.resname=='Cl') or (a.resname=='K') or (a.resname=='KJ') \
Expand Down Expand Up @@ -632,7 +636,7 @@ def assign_moltype(self):
else:
self.moltype = 'unknown'

def read(self, filename, bPDBTER=False, bNoNewID=True, bPDBGAP=False, bPDBMASS=False):
def read(self, filename, bPDBTER=False, bNoNewID=True, bPDBGAP=False, bPDBMASS=False, bCONECT=False):
"""PDB/GRO file reader.

Parameters
Expand All @@ -647,6 +651,8 @@ def read(self, filename, bPDBTER=False, bNoNewID=True, bPDBGAP=False, bPDBMASS=F
True. Default is True.
bPDBGAP : bool
whether search for gaps in the chain to assign new chain IDs.
bCONECT : bool
whether to read CONECT entries from PDB
"""
ext = filename.split('.')[-1]
if ext == 'pdb':
Expand All @@ -655,7 +661,7 @@ def read(self, filename, bPDBTER=False, bNoNewID=True, bPDBGAP=False, bPDBMASS=F
pdbline=None,
bNoNewID=bNoNewID, bPDBGAP=bPDBGAP, bPDBMASS=bPDBMASS)
else:
return self.__readPDB(fname=filename)
return self.__readPDB(fname=filename, bCONECT=bCONECT)
elif ext == 'gro':
return self.__readGRO(filename)
else:
Expand Down Expand Up @@ -808,7 +814,7 @@ def fetch_residue(self, idx, chain=None):
else:
valid_resids.append(r.id)

if chain is not None:
if chain is not None:
valid_chresids = []
for r in self.chdic[chain].residues:
if isinstance(r.id,str):
Expand Down Expand Up @@ -854,7 +860,7 @@ def fetch_residue(self, idx, chain=None):
else:
if r.id.replace(" ","")==idx:
return r


def fetch_residues(self, key, inv=False):
"""Gets residues using a list of residue names.
Expand Down Expand Up @@ -989,7 +995,7 @@ def assign_masses_to_model(model, topology=None):
topology : Topology
Topology object of the same molecule. When no topology provided, standard library masses are used.
'''

if topology!=None:
for ma, ta in zip(model.atoms, topology.atoms):
if ma.name != ta.name:
Expand Down
Loading