ACEMD之利用HTMD进行构象分析

序

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/