User:Tlinnet
Author of following PyMOL scripts
PROPKA for PyMOL 12200
Displacement map for conformational changes 16400
F?orster distance calculator 7600
Rotation matrix calculator 6300
Color by displacement 4200
Surface Cysteine pKa predictor 4000
Tutorial
Test
from __future__ import print_function
from pymol import cmd
import os
import sys
import math
from time import localtime, strftime
# Thx for inspiration from Simple scriptin PymMOl http://www.pymolwiki.org/index.php/Simple_Scripting
# Made by Ma.Sc student. Troels Linnet, 2011-08. troels.linnet@bbz.uni-leipzig.de
# Based solely on the work by:
# Maik H. Jacob, Dan Amir, Vladimir Ratner, Eugene Gussakowsky, and Elisha Haas
# Predicting Reactivities of Protein Surface Cysteines as Part of a Strategy for Selective Multiple Labeling. (Biochemistry 2005, 44, 13664-13672)
# Example of pymol script: Directory "predict_reactivity" has script file cyspka.py and cysteine residue pdb file: cys.pdb
# import cyspka
# fetch 4AKE, async=0
# create 4AKE-A, /4AKE//A and not resn HOH
# delete 4AKE
# hide everything
# show cartoon, 4AKE-A
# cyspka 4AKE-A, A, 18
def cyspka(molecule, chain, residue, SeeProgress='yes', pH=7.2, MoveSGatom='no', SGatom=str((0, 0, 0))):
# If SeeProgress='yes', computation time will take 10-20% extra, but nice to follow.
cmd.refresh()
RotationRange = 360
RotationDegree = 1
# For error checking, the energies can be printed out
printMC = 'no'
printSC = 'no'
# Parameters
DieElecSpheDist = 7.0
DieElecWaterDist = 1.4
DieElecWater = 78.5
DieElecCore = 4.0
BornPenaltyB = 1.0
AvogadroR = 8.31446216
Temp = 298
DeltapKMCSC = 0
pK1 = 9.25
pK2 = 8.0
NotPopuDist = 2.4
PopEnergyPenalty = 10000000
# Side chain discrete charges
DieElecSC = 40.0
SCchargeASP = -1
SCchargeGLU = -1
SCchargeOXT = -1
SCchargeARG = +1
SCchargeHIS = +1
SCchargeLYS = +1
SCchargeMET1 = +1
# Main chain partial charges
NrMainchainNeighBours = 5
DieElecMC = 22.0
MCchargeC = +0.55
MCchargeO = -0.55
MCchargeN = -0.35
MCchargeH = +0.35
MCchargeProCA = +0.1
MCchargeProCD = +0.1
MCchargeProN = -0.2
# Loading an Cys residue, give it a logic name, and aligning it. The oxygen atom can not be aligned in many cases, and are skipped.
# We use only this molecule, to find the initial position of the SG atom, and to rotate the SG atom around the CA-CB bond. The molecule atom positions are not used for electric potential calculatons.
Cysmolecule = str(molecule) + str(residue) + "Cys"
cmd.fragment("cys")
cmd.set_name('cys', Cysmolecule)
# We use pair_fir, since align and super gets unstable with so few atoms
pairfitCys(Cysmolecule, molecule, chain, residue)
# Give nice representations quickly
cmd.show("sticks", Cysmolecule)
cmd.select(str(molecule) + str(residue) + "Res", "/" + molecule + "//" + chain + "/" + residue)
print("/" + molecule + "//" + chain + "/" + residue)
cmd.show("sticks", str(molecule) + str(residue) + "Res")
cmd.disable(str(molecule) + str(residue) + "Res")
# Find out what is the residuename we are investigating for
Respdbstr = cmd.get_pdbstr(str(molecule) + str(residue) + "Res")
Ressplit = Respdbstr.split()
residueName = Ressplit[3]
print("")
print("# Hello, PyMOLers. It should take around 1 minute per residue.")
print("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH))
# Determine the range of neighbour residues possible.
Maxresidues = cmd.count_atoms("/" + molecule + "//" + chain + " and name CA")
for i in range(NrMainchainNeighBours + 1):
if int(residue) - i >= 1:
Minresidue = int(residue) - i
else:
break
for i in range(NrMainchainNeighBours + 1):
if int(residue) + i <= Maxresidues:
Maxresidue = int(residue) + i
else:
break
# Get the position and the vector for the CA->CB bond.
dihedN = "/" + Cysmolecule + "//" + "/" + "/N"
dihedCA = "/" + Cysmolecule + "//" + "/" + "/CA"
dihedCB = "/" + Cysmolecule + "//" + "/" + "/CB"
dihedSG = "/" + Cysmolecule + "//" + "/" + "/SG"
dihedralPosCA = cmd.get_atom_coords(dihedCA)
dihedralPosSG = cmd.get_atom_coords(dihedSG)
dihedralVector = AtomVector(dihedCA, dihedCB)
# To compare with article, we can move the SGatom to a starting position. The rotation is still determined around the CA-CB bond.
if MoveSGatom == 'yes':
SGatom = [float(SGatom[1:-1].split(",")[0]), float(SGatom[1:-1].split(",")[1]), float(SGatom[1:-1].split(",")[2])]
Translate = [(SGatom[0] - dihedralPosSG[0]), (SGatom[1] - dihedralPosSG[1]), (SGatom[2] - dihedralPosSG[2])]
cmd.translate(Translate, dihedSG, state=0, camera=0)
dihedralPosSG = cmd.get_atom_coords(dihedSG)
# Create a pymol molecule, that in the end will hold and show all SG atoms. Gives the representation of the rotameric states.
SGName = str(molecule) + str(residue) + "SG"
cmd.create(SGName, "None")
# Create a pymol molecule, that in the end will hold and show all Amide protons. Gives a nice representation, and easy to delete.
AmideName = str(molecule) + str(residue) + "NH"
cmd.create(AmideName, "None")
# Check if there are any nearby SG atoms, which could make a SG-SG dimer formation. The
breakDimer = "no"
breakDimer = CheckDimer(dihedSG, molecule, chain, residue)
# Create a list for appending the calculated energies.
ListofEnergies = []
ListofRotamerDiscarded = []
# print "Angle before rotation", cmd.get_dihedral(dihedN,dihedCA,dihedCB,dihedSG)
# Enter into the loop of rotameric states
for i in range(int(math.floor(RotationRange / RotationDegree))):
Angle = i * RotationDegree
# Create pymol molecule/SG atom for which we will calculate for.
SGNameAngle = str(residue) + "SG" + str(Angle)
cmd.create(SGNameAngle, dihedSG)
# Calculate new coordinates for rotation around CA->CB bond. Then translate the created SG atom.
SGNewPos = fRotateAroundLine(dihedralPosSG, dihedralPosCA, dihedralVector, Angle)
Translate = [(SGNewPos[0] - dihedralPosSG[0]), (SGNewPos[1] - dihedralPosSG[1]), (SGNewPos[2] - dihedralPosSG[2])]
cmd.translate(Translate, SGNameAngle, state=0, camera=0)
# If one wants to "see it happen" while its making the states. But it will take extra computation time.
if SeeProgress == 'yes':
cmd.refresh()
# Calculate number of neighbours within 2.4 Angstrom. Amide hydrogens are not considered, and are actually not build yet.
nameselect = "(((/" + molecule + "//" + chain + " and not /" + molecule + "//" + chain + "/" + residue + ") or /" + molecule + "//" + chain + "/" + residue + "/N+CA+C+O) within " + str(NotPopuDist) + " of /" + SGNameAngle + "//" + "/" + "/SG) and not resn HOH"
# print nameselect
cmd.select("NotPop", nameselect)
NotPopNr = cmd.count_atoms("NotPop")
# print Angle, NotPopNr, cmd.get_dihedral(dihedN,dihedCA,dihedCB,SGNameAngle)
# If no neighbours, then proceed calculating
if NotPopNr == 0:
SumAllWMC = 0.0
# Now calculate the electric potential due to the side chains.
SumWSC = fSumWSC(molecule, SGNameAngle, chain, residue, DieElecSC, SCchargeASP, SCchargeGLU, SCchargeOXT, SCchargeARG, SCchargeHIS, SCchargeLYS, SCchargeMET1, printSC)
# Now we calculate for the flanking 5 peptide groups on each side of the Cysteine CA atom.
# For the first residue, only calculate for the tailing C,O atom in the peptide bond. No test for Proline.
SumWMCFirst = fSumWMCFirst(molecule, SGNameAngle, chain, residue, Minresidue, DieElecMC, MCchargeC, MCchargeO, printMC)
# For the residue itself, we dont test for PRO, since it should be a Cysteine.
SumWMCresidue = fSumWMCresidue(molecule, SGNameAngle, chain, residue, int(residue), DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, AmideName, printMC)
# For the last residue, we test for Proline. We only calculate for the N,H atom, or if Proline, N,CA and CD atom.
SumWMCLast = fSumWMCLast(molecule, SGNameAngle, chain, residue, Maxresidue, DieElecMC, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC)
# Then loop over rest of the residues in the chain.
for j in (list(range(Minresidue + 1, int(residue))) + list(range(int(residue) + 1, Maxresidue))):
MCNeighbour = j
# print "Looking at neighbour", j
SumWMC = fSumWMC(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC)
SumAllWMC = SumAllWMC + SumWMC
# print "Rotation: %s Neighbour: %s " % (Angle, j)
# Since the SG atom is negative, we multiply with -1.
SumMCSC = -1 * (SumWSC + SumWMCFirst + SumWMCresidue + SumWMCLast + SumAllWMC)
# Makes the neighbour count. Everything in 'molecule" within 7 ang of aligned SG atom. Not counting 'residue'. Adding 5 for 'residue' N,CA,C,O,CB
ListNeighbourCount = fNeighbourCount(molecule, SGNameAngle, chain, residue, DieElecSpheDist)
# Calculate the weighted electric potential and alter the b factor for coloring. Then add the rotated SG into bucket of SG atoms.
SG_MCSC_Weight = fBoltzSingleState(SumMCSC, AvogadroR, Temp) * SumMCSC
cmd.alter(SGNameAngle, 'b="%s"' % SG_MCSC_Weight)
cmd.alter(SGNameAngle, 'name="S%s"' % Angle)
cmd.create(SGName, SGName + " + " + SGNameAngle)
# Then save the calculated values
ListofEnergies.append([Angle, SumMCSC, ListNeighbourCount, NotPopNr, SG_MCSC_Weight, cmd.get_atom_coords(SGNameAngle)])
cmd.delete(SGNameAngle)
else:
SumMCSCPenalty = PopEnergyPenalty
ListNeighbourCount = fNeighbourCount(molecule, SGNameAngle, chain, residue, DieElecSpheDist)
ListofRotamerDiscarded.append([Angle, SumMCSCPenalty, ListNeighbourCount, NotPopNr, 0, cmd.get_atom_coords(SGNameAngle)])
cmd.delete(SGNameAngle)
# Now show all the SG atoms as the available rotameric states.
cmd.show("nb_spheres", SGName)
cmd.delete("NotPop")
cmd.spectrum("b", selection=SGName)
AvailRotStates = len(ListofEnergies)
# print "Available Rotational States: ", AvailRotStates
# Do the calculations according to eq 5.
# Find the partition function
BoltzPartition = 0.0
for i in range(len(ListofEnergies)):
Boltz = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp)
BoltzPartition = BoltzPartition + Boltz
# Find the summed function
BoltzSumNi = 0.0
for i in range(len(ListofEnergies)):
BoltzNi = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp) * ListofEnergies[i][1]
BoltzSumNi = BoltzSumNi + BoltzNi
# Check if there was any possible rotamers
nostates = "no"
if len(ListofEnergies) == 0:
print("####################################################")
print("########### WARNING: No states available ###########")
print("########### Did you mutate a Glycine? ###########")
print("####################################################")
BoltzSumNi = 0
BoltzPartition = 0
BoltzMCSC = 0
DeltapKMCSC = 99
NeighbourCount = 0
nostates = "yes"
else:
# Final calculation
BoltzMCSC = (BoltzSumNi) / (BoltzPartition)
DeltapKMCSC = fDeltapK(BoltzMCSC, AvogadroR, Temp)
# Find average number of neighbours
NCSum = 0.0
NCWeightedSum = 0.0
for i in range(len(ListofEnergies)):
NCi = ListofEnergies[i][2]
NCSum = NCSum + NCi
NCWeightedi = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp) * ListofEnergies[i][2] / BoltzPartition
NCWeightedSum = NCWeightedSum + NCWeightedi
# print "Weighted neighbour", int(round(NCWeightedSum))
#NeighbourCount = int(round(NCSum/len(ListofEnergies)))
NeighbourCount = round(NCWeightedSum, 1)
# If we found dimers
if breakDimer == "yes":
print("####################################################")
print("########### WARNING: Dimer formation? ###########")
print("####################################################")
BoltzSumNi = 0
BoltzPartition = 0
BoltzMCSC = 0
DeltapKMCSC = 99
NeighbourCount = 0
# Calculate the BornPenalty based on the neighbour count. It's a wrapper script for equation 13, 12, 11.
EnergyBornPenalty = fEnergyBornPenalty(DieElecSpheDist, DieElecWaterDist, NeighbourCount, DieElecWater, DieElecCore, BornPenaltyB)
DeltapKB = fDeltapK(EnergyBornPenalty, AvogadroR, Temp)
# Do the calculations according to eq 3 and 9.
pKm1 = fpKm1(DeltapKMCSC, pK1)
pKm2 = fpKm2(DeltapKMCSC, DeltapKB, pK2)
FracCysm1 = fFracCys(pKm1, pH)
FracCysm2 = fFracCys(pKm2, pH)
# Lets make a result file, and write out the angle, the SumMCSC, and the number of neighbours for this state.
Currentdir = os.getcwd()
Newdir = os.path.join(os.getcwd(), "Results")
if not os.path.exists(Newdir):
os.makedirs(Newdir)
filename = os.path.join(".", "Results", "Result_" + molecule + "_" + chain + "_" + residue + ".txt")
filenamelog = os.path.join(".", "Results", "Result_log.log")
logfile = open(filenamelog, "a")
outfile = open(filename, "w")
timeforlog = strftime("%Y %b %d %a %H:%M:%S", localtime())
logfile.write("# " + timeforlog + "\n")
logfile.write("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH) + "\n")
logfile.write("# BoltzSumNi: BoltzPartition: BoltzMCSC" + "\n")
logfile.write(("# %.4f %.4f %.4f" + '\n') % (BoltzSumNi, BoltzPartition, BoltzMCSC))
logfile.write("# Res NC States pKmcsc pK1 pKB pK2 pKm1 pKm2 f(C-)m1 f(C-)m2" + "\n")
logfile.write(("; %s %s %s %s %.4f %s %.4f %s %.4f %.4f %.6f %.6f" + '\n') % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
if nostates == "yes":
logfile.write("##### ERROR; No states available ###" + "\n")
if breakDimer == "yes":
logfile.write("##### ERROR; Dimer formation ###" + "\n")
logfile.write('\n')
outfile.write("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH) + "\n")
outfile.write("# BoltzSumNi: BoltzPartition: BoltzMCSC" + "\n")
outfile.write(("# %.4f %.4f %.4f" + '\n') % (BoltzSumNi, BoltzPartition, BoltzMCSC))
outfile.write("# Res NC States pKmcsc pK1 pKB pK2 pKm1 pKm2 f(C-)m1 f(C-)m2" + "\n")
outfile.write(("; %s %s %s %s %.4f %s %.4f %s %.4f %.4f %.6f %.6f" + '\n') % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
if nostates == "yes":
outfile.write("##### ERROR; No states available ###" + "\n")
if breakDimer == "yes":
outfile.write("##### ERROR; Dimer formation ###" + "\n")
outfile.write('\n')
outfile.write("#Ang SumMCSC NC rNC MCSC_Weight SG[X,Y,Z]" + "\n")
for i in range(len(ListofEnergies)):
outfile.write("%4.1d %10.3f %2.1d %1.1d %10.3f [%8.3f, %8.3f, %8.3f]" % (ListofEnergies[i][0], ListofEnergies[i][1], ListofEnergies[i][2], ListofEnergies[i][3], ListofEnergies[i][4], ListofEnergies[i][5][0], ListofEnergies[i][5][1], ListofEnergies[i][5][2]) + '\n')
for i in range(len(ListofRotamerDiscarded)):
outfile.write("%4.1d %10.3f %2.1d %1.1d %10.3f [%8.3f, %8.3f, %8.3f]" % (ListofRotamerDiscarded[i][0], ListofRotamerDiscarded[i][1], ListofRotamerDiscarded[i][2], ListofRotamerDiscarded[i][3], ListofRotamerDiscarded[i][4], ListofRotamerDiscarded[i][5][0], ListofRotamerDiscarded[i][5][1], ListofRotamerDiscarded[i][5][2]) + '\n')
outfile.close()
# Now, we are done. Just print out. The ; is for a grep command to select these lines in the output.
print("# residue: %s %s. Average NeighbourCount NC= %s " % (residueName, residue, NeighbourCount))
print("# From residue %s to residue %s" % (Minresidue, Maxresidue))
print("# BoltzSumNi: BoltzPartition: BoltzMCSC")
print("# %.4f %.4f %.4f" % (BoltzSumNi, BoltzPartition, BoltzMCSC))
print("# Result written in file: %s" % (filename))
print("# Res NC States pKmcsc pK1 pKB pK2 pKm1 pKm2 f(C-)m1 f(C-)m2")
print("; %s %s %s %s %.4f %s %.4f %s %.4f %.4f %.6f %.6f" % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
if nostates == "yes":
print("##### ERROR; No states available ###")
if breakDimer == "yes":
print("##### ERROR; Dimer formation ###")
cmd.extend("cyspka", cyspka)
def loopcyspka(molecule, chain, residue, SeeProgress='no', pH=7.2, MoveSGatom='no', SGatom=str((0, 0, 0))):
residue = residue.split('.')
residueList = []
for i in residue:
if '-' in i:
tmp = i.split('-')
residueList.extend(list(range(int(tmp[0]), int(tmp[-1]) + 1)))
if '-' not in i:
residueList.append(int(i))
print("Looping over residues")
print(residueList)
for i in residueList:
cyspka(molecule, chain, str(i), SeeProgress, pH, MoveSGatom, SGatom)
cmd.extend("loopcyspka", loopcyspka)
def fNeighbourCount(molecule, Cysmolecule, chain, residue, DieElecSpheDist):
nameselect = "(((/" + molecule + "//" + chain + " and not /" + molecule + "//" + chain + "/" + residue + ") or /" + molecule + "//" + chain + "/" + residue + "/N+CA+C+O) within " + str(DieElecSpheDist) + " of /" + Cysmolecule + "//" + "/" + "/SG) and not resn HOH"
# print nameselect
cmd.select(residue + "NC", nameselect)
# Adding 1 for CB
Neighbours = cmd.count_atoms(residue + "NC") + 1
cmd.delete(residue + "NC")
return Neighbours
def fNeighbourWater(DieElecSpheDist, DieElecWaterDist, NeighbourCount):
Waters = 0.74 * math.pow(DieElecSpheDist, 3) / math.pow(DieElecWaterDist, 3) - NeighbourCount
return Waters
def fDieElecEF(NeighbourWater, DieElecWater, NeighbourCount, DieElecCore):
DieElecEF = (NeighbourWater * DieElecWater + NeighbourCount * DieElecCore) / (NeighbourWater + NeighbourCount)
return DieElecEF
def fBornPenalty(BornPenaltyB, DieElecEF, DieElecWater):
BornPenalty = (1.39 * math.pow(10, 6)) / (2 * BornPenaltyB) * (1.0 / DieElecEF - 1.0 / DieElecWater)
return BornPenalty
def fEnergyBornPenalty(DieElecSpheDist, DieElecWaterDist, NeighbourCount, DieElecWater, DieElecCore, BornPenaltyB):
NeighbourWater = fNeighbourWater(DieElecSpheDist, DieElecWaterDist, NeighbourCount)
DieElecEF = fDieElecEF(NeighbourWater, DieElecWater, NeighbourCount, DieElecCore)
BornPenalty = fBornPenalty(BornPenaltyB, DieElecEF, DieElecWater)
return BornPenalty
def fDeltapK(Energy, AvogadroR, Temp):
DeltapK = -1 * math.log10(math.exp(-Energy / (AvogadroR * Temp)))
return DeltapK
def fRotateAroundLine(OriPoint, ThroughLinePoint, LineVector, AngleDeg):
# See http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/. Section 6.1
AngleRad = math.radians(AngleDeg)
x = OriPoint[0]
y = OriPoint[1]
z = OriPoint[2]
a = ThroughLinePoint[0]
b = ThroughLinePoint[1]
c = ThroughLinePoint[2]
u = LineVector[0]
v = LineVector[1]
w = LineVector[2]
L = math.pow(u, 2) + math.pow(v, 2) + math.pow(w, 2)
xPos = ((a * (math.pow(v, 2) + math.pow(w, 2)) - u * (b * v + c * w - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * x * math.cos(AngleRad) + math.sqrt(L) * (-c * v + b * w - w * y + v * z) * math.sin(AngleRad)) / L
yPos = ((b * (math.pow(u, 2) + math.pow(w, 2)) - v * (a * u + c * w - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * y * math.cos(AngleRad) + math.sqrt(L) * (c * u - a * w + w * x - u * z) * math.sin(AngleRad)) / L
zPos = ((c * (math.pow(u, 2) + math.pow(v, 2)) - w * (a * u + b * v - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * z * math.cos(AngleRad) + math.sqrt(L) * (-b * u + a * v - v * x + u * y) * math.sin(AngleRad)) / L
NewPos = [xPos, yPos, zPos]
return NewPos
def fWSC(charge, DieElecSC, DistR):
# print charge, DistR
WSC = 1.39 * math.pow(10, 6) * charge / (DieElecSC * DistR)
return WSC
def fSumWSC(molecule, SGNameAngle, chain, residue, DieElecSC, SCchargeASP, SCchargeGLU, SCchargeOXT, SCchargeARG, SCchargeHIS, SCchargeLYS, SCchargeMET1, printSC):
SumWSC = 0.0
SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
# Sidechain ASP
nameselect = "/" + molecule + " and resn ASP and name CG and not resi " + residue
cmd.select("SC", nameselect)
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distASP', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeASP, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC ASP ", str(SClist[i]), " ", SCchargeASP, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distASP')
# Sidechain GLU
nameselect = "/" + molecule + " and resn GLU and name CD and not resi " + residue
cmd.select("SC", nameselect)
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distGLU', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeGLU, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC GLU ", str(SClist[i]), " ", SCchargeGLU, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distGLU')
# print "GLU", cmd.count_atoms("SC"), SumWSC
# Sidechain OXT
nameselect = "/" + molecule + " and byres name OXT and not resi " + residue
cmd.select("SC", nameselect)
cmd.select("SC", "SC and name C")
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distOXT', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeOXT, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC OXT ", str(SClist[i]), " ", SCchargeOXT, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distOXT')
# print "OXT", cmd.count_atoms("SC"), SumWSC
# Sidechain ARG
nameselect = "/" + molecule + " and resn ARG and name CZ and not resi " + residue
cmd.select("SC", nameselect)
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distARG', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeARG, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC ARG ", str(SClist[i]), " ", SCchargeARG, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distARG')
# print "ARG", cmd.count_atoms("SC"), SumWSC
# Sidechain HIS
nameselect = "/" + molecule + " and resn HIS and name CD2 and not resi " + residue
cmd.select("SC", nameselect)
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distHIS', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeHIS, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC HIS ", str(SClist[i]), " ", SCchargeHIS, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distHIS')
# print "HIS", cmd.count_atoms("SC"), SumWSC
# Sidechain LYS
nameselect = "/" + molecule + " and resn LYS and name NZ and not resi " + residue
cmd.select("SC", nameselect)
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distLYS', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeLYS, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC LYS ", str(SClist[i]), " ", SCchargeLYS, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distLYS')
# print "LYS", cmd.count_atoms("SC"), SumWSC
# Sidechain MET1
nameselect = "/" + molecule + " and resn MET and res 1 and not resi " + residue
cmd.select("SC", nameselect)
cmd.select("SC", "SC and name N")
SClist = cmd.identify("SC")
for i in range(len(SClist)):
ResDist = cmd.dist(residue + 'distMET1', SGnameselect, molecule + " and id " + str(SClist[i]))
WSC = fWSC(SCchargeMET1, DieElecSC, ResDist)
SumWSC = SumWSC + WSC
if printSC == 'yes':
print("SC MET1 ", str(SClist[i]), " ", SCchargeMET1, " ", DieElecSC, " ", ResDist, " ", WSC)
cmd.delete(residue + 'distMET1')
# print "MET1", cmd.count_atoms("SC"), SumWSC
cmd.delete("SC")
return SumWSC
def fWMC(charge, DieElecMC, DistR):
WMC = 1.39 * math.pow(10, 6) * charge / (DieElecMC * DistR)
return WMC
def fSumWMCFirst(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, printMC):
# print "First", MCNeighbour
SumWMCFirst = 0.0
SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
cmd.select("MC", NBnameselect)
MCpdbstr = cmd.get_pdbstr("MC")
MCsplit = MCpdbstr.split()
residueName = MCsplit[3]
# print NBnameselect, residueName
# Mainchain C
Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
ResDist = cmd.dist(residue + 'distFirstC', SGnameselect, Cnameselect)
WMC = fWMC(MCchargeC, DieElecMC, ResDist)
SumWMCFirst = SumWMCFirst + WMC
if printMC == 'yes':
print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain O
Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
ResDist = cmd.dist(residue + 'distFirstO', SGnameselect, Onameselect)
WMC = fWMC(MCchargeO, DieElecMC, ResDist)
SumWMCFirst = SumWMCFirst + WMC
if printMC == 'yes':
print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
cmd.delete(residue + 'distFirstC')
cmd.delete(residue + 'distFirstO')
cmd.delete("MC")
return SumWMCFirst
def fSumWMCresidue(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, AmideName, printMC):
# print "residue", MCNeighbour
SumWMCresidue = 0.0
SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
cmd.select("MC", NBnameselect)
MCpdbstr = cmd.get_pdbstr("MC")
MCsplit = MCpdbstr.split()
residueName = MCsplit[3]
# print NBnameselect, residueName
AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
cmd.h_add(HbuildSelect)
cmd.create(AmideName, AmideName + " + " + AmideProt)
cmd.remove(AmideProt)
# Mainchain AmideH
ResDist = cmd.dist(residue + 'distResH', SGnameselect, Hnameselect)
WMC = fWMC(MCchargeH, DieElecMC, ResDist)
SumWMCresidue = SumWMCresidue + WMC
if printMC == 'yes':
print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain C
Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
ResDist = cmd.dist(residue + 'distResC', SGnameselect, Cnameselect)
WMC = fWMC(MCchargeC, DieElecMC, ResDist)
SumWMCresidue = SumWMCresidue + WMC
if printMC == 'yes':
print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain O
Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
ResDist = cmd.dist(residue + 'distResO', SGnameselect, Onameselect)
WMC = fWMC(MCchargeO, DieElecMC, ResDist)
SumWMCresidue = SumWMCresidue + WMC
if printMC == 'yes':
print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain N
Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
ResDist = cmd.dist(residue + 'distResN', SGnameselect, Nnameselect)
WMC = fWMC(MCchargeN, DieElecMC, ResDist)
SumWMCresidue = SumWMCresidue + WMC
if printMC == 'yes':
print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
cmd.delete(residue + 'distResH')
cmd.delete(residue + 'distResC')
cmd.delete(residue + 'distResO')
cmd.delete(residue + 'distResN')
cmd.show("nb_spheres", AmideName)
cmd.delete("MC")
return SumWMCresidue
def fSumWMCLast(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC):
# print "Last", MCNeighbour
SumWMCLast = 0.0
SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
cmd.select("MC", NBnameselect)
MCpdbstr = cmd.get_pdbstr("MC")
MCsplit = MCpdbstr.split()
residueName = MCsplit[3]
# print NBnameselect, residueName
if residueName == "PRO":
# Proline CA
CAnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CA"
ResDist = cmd.dist(residue + 'distLastProCA', SGnameselect, CAnameselect)
WMC = fWMC(MCchargeProCA, DieElecMC, ResDist)
SumWMCLast = SumWMCLast + WMC
if printMC == 'yes':
print("MC ProCA ", MCNeighbour, " ", MCchargeProCA, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline CD
CDnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CD"
ResDist = cmd.dist(residue + 'distLastProCD', SGnameselect, CDnameselect)
WMC = fWMC(MCchargeProCD, DieElecMC, ResDist)
SumWMCLast = SumWMCLast + WMC
if printMC == 'yes':
print("MC ProCD ", MCNeighbour, " ", MCchargeProCD, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline N
Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
ResDist = cmd.dist(residue + 'distLastProN', SGnameselect, Nnameselect)
WMC = fWMC(MCchargeProN, DieElecMC, ResDist)
SumWMCLast = SumWMCLast + WMC
if printMC == 'yes':
print("MC ProN ", MCNeighbour, " ", MCchargeProN, " ", DieElecMC, " ", ResDist, " ", WMC)
else:
AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
cmd.h_add(HbuildSelect)
cmd.create(AmideName, AmideName + " + " + AmideProt)
cmd.remove(AmideProt)
# Mainchain AmideH
ResDist = cmd.dist(residue + 'distLastH', SGnameselect, Hnameselect)
WMC = fWMC(MCchargeH, DieElecMC, ResDist)
SumWMCLast = SumWMCLast + WMC
if printMC == 'yes':
print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain N
Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
ResDist = cmd.dist(residue + 'distLastN', SGnameselect, Nnameselect)
WMC = fWMC(MCchargeN, DieElecMC, ResDist)
SumWMCLast = SumWMCLast + WMC
if printMC == 'yes':
print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
cmd.delete(residue + 'distLastProCA')
cmd.delete(residue + 'distLastProCD')
cmd.delete(residue + 'distLastProN')
cmd.delete(residue + 'distLastH')
cmd.delete(residue + 'distLastN')
cmd.show("nb_spheres", AmideName)
cmd.delete("MC")
return SumWMCLast
def fSumWMC(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC):
# print "chain", MCNeighbour
SumWMC = 0.0
SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
cmd.select("MC", NBnameselect)
MCpdbstr = cmd.get_pdbstr("MC")
MCsplit = MCpdbstr.split()
residueName = MCsplit[3]
# print NBnameselect, residueName
if residueName == "PRO":
# Proline CA
CAnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CA"
ResDist = cmd.dist(residue + 'distProCA', SGnameselect, CAnameselect)
WMC = fWMC(MCchargeProCA, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC ProCA ", MCNeighbour, " ", MCchargeProCA, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline CD
CDnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CD"
ResDist = cmd.dist(residue + 'distProCD', SGnameselect, CDnameselect)
WMC = fWMC(MCchargeProCD, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC ProCD ", MCNeighbour, " ", MCchargeProCD, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline N
Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
ResDist = cmd.dist(residue + 'distProN', SGnameselect, Nnameselect)
WMC = fWMC(MCchargeProN, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC ProN ", MCNeighbour, " ", MCchargeProN, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline C
Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
ResDist = cmd.dist(residue + 'distProC', SGnameselect, Cnameselect)
WMC = fWMC(MCchargeC, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC ProC ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
# Proline O
Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
ResDist = cmd.dist(residue + 'distProO', SGnameselect, Onameselect)
WMC = fWMC(MCchargeO, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC ProO ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
else:
AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
cmd.h_add(HbuildSelect)
cmd.create(AmideName, AmideName + " + " + AmideProt)
cmd.remove(AmideProt)
# Mainchain AmideH
ResDist = cmd.dist(residue + 'distH', SGnameselect, Hnameselect)
WMC = fWMC(MCchargeH, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain C
Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
ResDist = cmd.dist(residue + 'distC', SGnameselect, Cnameselect)
WMC = fWMC(MCchargeC, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain O
Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
ResDist = cmd.dist(residue + 'distO', SGnameselect, Onameselect)
WMC = fWMC(MCchargeO, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
# Mainchain N
Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
ResDist = cmd.dist(residue + 'distN', SGnameselect, Nnameselect)
WMC = fWMC(MCchargeN, DieElecMC, ResDist)
SumWMC = SumWMC + WMC
if printMC == 'yes':
print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
cmd.delete(residue + 'distProCA')
cmd.delete(residue + 'distProCD')
cmd.delete(residue + 'distProN')
cmd.delete(residue + 'distProC')
cmd.delete(residue + 'distProO')
cmd.delete(residue + 'distH')
cmd.delete(residue + 'distC')
cmd.delete(residue + 'distO')
cmd.delete(residue + 'distN')
cmd.show("nb_spheres", AmideName)
cmd.delete("MC")
return SumWMC
def fBoltzSingleState(SumMCSC, AvogadroR, Temp):
BoltzSingleState = math.exp(-SumMCSC / (AvogadroR * Temp))
return BoltzSingleState
def fpKm1(DeltapKMCSC, pK1):
pKm1 = DeltapKMCSC + pK1
return pKm1
def fpKm2(DeltapKMCSC, DeltapKB, pK2):
pKm2 = DeltapKMCSC + DeltapKB + pK2
return pKm2
def fFracCys(pKm, pH):
FracCys = 1.0 / (math.pow(10, (pKm - pH)) + 1)
return FracCys
def AtomVector(AtomStart, AtomEnd):
PosStart = cmd.get_atom_coords(AtomStart)
PosEnd = cmd.get_atom_coords(AtomEnd)
VectorDiff = [(PosEnd[0] - PosStart[0]), (PosEnd[1] - PosStart[1]), (PosEnd[2] - PosStart[2])]
return VectorDiff
def pairfitCys(Cysmolecule, molecule, chain, residue):
RN = "/" + Cysmolecule + "//" + "/" + "/N"
PN = "/" + molecule + "//" + chain + "/" + residue + "/N"
RCA = "/" + Cysmolecule + "//" + "/" + "/CA"
PCA = "/" + molecule + "//" + chain + "/" + residue + "/CA"
RC = "/" + Cysmolecule + "//" + "/" + "/C"
PC = "/" + molecule + "//" + chain + "/" + residue + "/C"
RCB = "/" + Cysmolecule + "//" + "/" + "/CB"
PCB = "/" + molecule + "//" + chain + "/" + residue + "/CB"
cmd.select("CBatom", PCB)
CBatomNr = cmd.count_atoms("CBatom")
# If PRO or GLY, then only fit N, CA, C atoms
if CBatomNr == 0:
cmd.pair_fit(RN, PN, RCA, PCA, RC, PC)
else:
# cmd.pair_fit(RN,PN,RCA,PCA,RC,PC,RCB,PCB)
cmd.pair_fit(RN, PN, RCA, PCA, RC, PC)
cmd.delete("CBatom")
def CheckDimer(dihedSG, molecule, chain, residue):
breakDimer = "no"
nameselect = "(/" + molecule + "//" + chain + " and name SG and not /" + molecule + "//" + chain + "/" + residue + ") within 5 of " + dihedSG
cmd.select(str(molecule) + str(residue) + "Dimer", nameselect)
DimerSG = cmd.count_atoms(str(molecule) + str(residue) + "Dimer")
if DimerSG > 0:
print("####################################################")
print("########### WARNING: SG in near detected ###########")
print("########### Is this a dimer? ###########")
print("####################################################")
cmd.select(str(molecule) + str(residue) + "Dimer", "byres " + str(molecule) + str(residue) + "Dimer")
cmd.show("sticks", str(molecule) + str(residue) + "Dimer")
breakDimer = "yes"
else:
cmd.delete(str(molecule) + str(residue) + "Dimer")
return breakDimer