SPONGE常用教程：蛋白+配体模拟3

前序课程1
前序课程2

目录

1. 应用场景简述；- [ Done ]
2. DSDP:蛋白-配体对接；- [ Done ]
3. XPONGE：蛋白-配体建模，加溶剂；- [ Done ]
4. SPONGE：能量极小化-NVT-NPT-正式模拟；- [ Done ]
5. 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对齐到起始构象示意图（绿色：平均构象；黄色：初始构象；天蓝色：配体初始位置；紫色：氢键作用较强的残基）。

本教程为教学演示用，搭建水盒偏小，成品模拟时间较短，只进行了基础的数据分析，正式使用时建议适当调整。