目录
- 应用场景简述;- [ Done ]
- DSDP:蛋白-配体对接;- [ Done ]
- XPONGE:蛋白-配体建模,加溶剂;- [ Done ]
- SPONGE:能量极小化-NVT-NPT-正式模拟;- [ Done ]
- XPONGE:数据简单后处理。
5. XPONGE:数据简单后处理
经过1ns的SPONGE分子动力学模拟,得到了轨迹文件"mdcrd.dat"和盒子边长信息"mdbox.txt"。
cp ../1ae5_pro_lig_sol.pdb .
#复制pdb文件到文件夹内进行分析
import MDAnalysis
import Xponge.analysis.md_analysis as xmda
from MDAnalysis import Universe
u = Universe("1ae5_pro_lig_sol.pdb","mdcrd.dat",box="mdbox.txt",format=xmda.SpongeTrajectoryReader)
#calculate RMSD
ca = u.select_atoms('name CA')
import MDAnalysis.analysis.rms as rms
R = rms.RMSD(ca, ca, select='all', ref_frame=0)
R.run()
R.rmsd.shape
import numpy as np
np.savetxt('rmsd.txt',R.rmsd)
#calculate Rg
import numpy as np
from MDAnalysis.analysis.base import (AnalysisBase, AnalysisFromFunction, analysis_class)
def radgyr(atomgroup, masses, total_mass=None):
# coordinates change for each frame
coordinates = atomgroup.positions
center_of_mass = atomgroup.center_of_mass()
# get squared distance from center
ri_sq = (coordinates-center_of_mass)**2
# sum the unweighted positions
sq = np.sum(ri_sq, axis=1)
sq_x = np.sum(ri_sq[:,[1,2]], axis=1) # sum over y and z
sq_y = np.sum(ri_sq[:,[0,2]], axis=1) # sum over x and z
sq_z = np.sum(ri_sq[:,[0,1]], axis=1) # sum over x and y
# make into array
sq_rs = np.array([sq, sq_x, sq_y, sq_z])
# weight positions
rog_sq = np.sum(masses*sq_rs, axis=1)/total_mass
# square root and return
return np.sqrt(rog_sq)
protein =u.select_atoms('protein')
rog = AnalysisFromFunction(radgyr, u.trajectory, protein, protein.masses, total_mass=np.sum(protein.masses))
rog.run()
np.savetxt('Rg.txt',rog.results['timeseries'])
#calculate Hbonds
u.add_TopologyAttr('charges')
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
hbonds = HBA(universe=u, between=['resname HUP', 'protein'])
hbonds.hydrogens_sel ="type H"
hbonds.acceptors_sel = "type O" or "type N"
hbonds.donors_sel = "type H"
hbonds.d_h_cutoff = 2.0
hbonds.run()
hbonds.results.hbonds
import numpy as np
np.savetxt('hbond.txt',hbonds.results.hbonds)
out=np.zeros((hbonds.results.hbonds.shape[0],3),dtype=float)
count=np.zeros(u.atoms.n_residues)
for i in range(0, hbonds.results.hbonds.shape[0]-1):
tempDresid=u.atoms.resids[int(hbonds.results.hbonds[i,1])]
tempAresid=u.atoms.resids[int(hbonds.results.hbonds[i,3])]
out[i,0]=hbonds.results.hbonds[i,0]
out[i,1]=tempDresid
out[i,2]=tempAresid
count[tempDresid-1]=count[tempDresid-1]+1
count[tempAresid-1]=count[tempAresid-1]+1
np.savetxt('res_hbond.txt',count)
np.savetxt('t_hbond.txt',hbonds.count_by_time())
#calculate RMSF
from MDAnalysis.analysis import align
average = align.AverageStructure(u, u, select='name CA', ref_frame=0,in_memory=True).run()
ref = average.universe
aligner = align.AlignTraj(u, ref, select='name CA', in_memory=True).run()
chain = ca.select_atoms('all')
R = rms.RMSF(chain).run()
import numpy as np
np.savetxt('rmsf.txt',R.rmsf)
#after RMSF, had better to exit, 'alignTraj' change the coord
#save average structure
from MDAnalysis.analysis import align
average = align.AverageStructure(u, u, select='protein', ref_frame=0,in_memory=True).run()
average = average.universe
average = average.select_atoms('all')
average.write('average.pdb')
对所得文本数据进行分析,得到下图:
数据存储格式如下:
RMSD(图左上,rmsd.txt): time_frame;time_frame;RMSD(unit:A)
Rg(图左中,Rg.txt): Rg(unit:A); Rg_x; Rg_y; Rg_z
Hbond(图左下,t_hbond.txt):hbond_counts
以上横坐标为时间t(/ps);
RMSF(图右上,rmsf.txt):RMSF(unit:A)
Hbond(图右中,res_hbond.txt):hbond_counts
以上横坐标为蛋白残基序号。
右下为平均构象average.pdb对齐到起始构象示意图(绿色:平均构象;黄色:初始构象;天蓝色:配体初始位置;紫色:氢键作用较强的残基)。
本教程为教学演示用,搭建水盒偏小,成品模拟时间较短,只进行了基础的数据分析,正式使用时建议适当调整。
标签:教程,sq,SPONGE,配体,np,import,txt,hbonds,average From: https://www.cnblogs.com/bgalang/p/18338363