序
HTMD是一个HTMD基于Python的分子可编程环境,用于准备、处理、模拟、可视化和分析分子系统。用户使用HTMD可以在几行代码内完成非常复杂的协,且HTMD是开源的,所以你可以添加你自己的应用程序
由于其基于Python,其可视化及分析作图是其非常大的优势,可以很好地利用Python的作图功能进行动力学数据的分析作图。本文主要介绍蛋白构象变化的分析,以及计算RMSD、自由能RMSF等。
Conformational analysis
#!/usr/bin/python3
# coding: utf-8
# In[1]:导入模块,开启视图与可视化
from htmd.ui import *
config(viewer='ngl')
%pylab inline
# In[2]:过滤轨迹
fsims = simlist(glob('./filtered/*/'), './filtered/filtered.pdb')
# In[3]:计算metrics,选择10到54的残基,设置时间步长
metr = Metric(fsims)
metr.projection(MetricDihedral(protsel='protein and resid 10 to 54', sincos=True))
data = metr.project()
data.fstep = 0.1
# In[4]:绘制轨迹长度
data.plotTrajSizes()
# In[5]:降维
tica = TICA(data, 20)
dataTica = tica.project(3)
# In[6]:聚类
dataTica.cluster(MiniBatchKMeans(n_clusters=200), mergesmall=5)
# In[7]:计算MaximumLikelihoodMSM并绘图
model = Model(dataTica)
model.plotTimescales(lags=list(range(1,1000,50)))
# In[8]:建立马尔科夫模型
model.markovModel(600, 8)
eqDist = model.eqDistribution()
print(eqDist)
# In[9]:可视化
#we can now visualize representatives for each of the equilibrium species
from htmd.config import config
config(viewer='ngl')
model.numsamples=1
model.viewStates(protein=True)
# In[10]:统计分析
# In[11]:过滤数据并可视化
# we can visualize which residues are different between states
filtered = Molecule('./filtered/filtered.pdb')
filtered.view(sel='protein',style='NewCartoon',hold=True)
filtered.view(sel='resid 16 17 24 25 32 33 44 45',style='Licorice')
# In[12]:获取一段轨迹
np.where(model.macro_ofmicro == 6)
# In[13]:输出获取轨迹信息
_,rel = model.sampleStates(10, 10, statetype='micro')
print(rel)
# In[14]:获取数据文件路径
# In[15]:提取2000帧轨迹
simid = 232
parent = None
input = []
trajectory = ['./filtered/9x9/9x9-GERARD_VERYLONG_CXCL12_confAna-0-1-RND2283_9.filtered.xtc']
molfile = ('./filtered/filtered.pdb')
numframes = [2000]
# In[16]:计算metrics
# The first selection corresponds to beta-sheet 2 carbons alpha, the second one to beta-sheet 3 CA.
# We specify metric='contacts' to create contact maps instead of proper distances,
# this means: create an interatom matrix and put 1 if the distance is below cutoff; 0 otherwise.
metr = Metric(fsims)
metr.projection(MetricDistance('resid 38 to 42 and noh', 'resid 22 to 28 and noh', metric='contacts'))
data3 = metr.project()
data3.fstep = 0.1
# In[17]:降维
# tICA projection (dimensionality reduction along the slow process)
tica3 = TICA(data3, 20)
dataTica3 = tica3.project(3)
# In[18]:聚类
# Clustering
dataTica3.cluster(MiniBatchKMeans(n_clusters=200), mergesmall=5)
# In[19]:#绘制timescales
model3 = Model(dataTica3)
model3.plotTimescales(lags=list(range(1,1000,50)))
# In[20]:建立马尔科夫模型
model3.markovModel(600, 4)
eqDist = model3.eqDistribution()
print(eqDist)
# In[21]:可视化
# Visualize states
model3.viewStates(protein=True, numsamples=1)
# In[22]:
# Map back the trajectory/ies that originated the macro. Substitute 1 for the macro that showed the pocket opening.
# This function is giving you the microclusters that are inside a given macrocluster
np.where(model3.macro_ofmicro ==1)
# In[23]:
# substitute 48 for the micro number from the previous step
# This function gives you trajectory-frame pairs that visited a given micro
_, rel = model3.sampleStates(48, 5, statetype='micro')
print(rel)
# In[24]:获取数据文件路径
print(model3.data.simlist[277])
# In[25]:提取2000帧轨迹
simid = 277
parent = None
input = []
trajectory = ['./filtered/10x23/10x23-GERARD_VERYLONG_CXCL12_confAna-0-1-RND9861_9.filtered.xtc']
molfile = ('./filtered/filtered.pdb')
numframes = [2000]
# In[26]:获取文件集合
simus = simlist(glob('./filtered/10x23/'), './filtered/filtered.pdb')
# In[27]:计算与参考轨迹的RMSD
refmol = Molecule('./filtered/filtered.pdb')
metr = Metric(simus)
metr.projection(MetricRmsd(refmol, 'resid 38 to 42 or resid 22 to 28 and noh', trajalnstr='protein'))
rmsd = metr.project()
# In[28]:绘制口袋变化的RMSD图
# Do you see the pocket opening at 50ns?
plt.plot(rmsd.dat[0])
plt.xlabel('Simulation length (frames; 0.1ns)', fontsize=10)
plt.ylabel('RMSD (Angstroms)', fontsize=10)
# In[29]:可视化轨迹,建立video
# You can also visualize the trajectory from your browser
refmol.read('./filtered/10x23/10x23-GERARD_VERYLONG_CXCL12_confAna-0-1-RND9861_9.filtered.xtc')
refmol.align('protein')
refmol.view()
参考资料
1.https://software.acellera.com/docs/latest/htmd/userguide/introduction.html
2.https://software.acellera.com/docs/latest/htmd/userguide/building.html
3.https://software.acellera.com/docs/latest/htmd/userguide/running.html
4.https://software.acellera.com/docs/latest/htmd/userguide/analysing.html
5.http://gainstrong.net/works/hudong/
领取专属 10元无门槛券
私享最新 技术干货