前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >OpenMM-组蛋白甲基转移酶分子动力学模拟-TIP3P

OpenMM-组蛋白甲基转移酶分子动力学模拟-TIP3P

作者头像
DrugScience
发布2021-04-23 14:39:10
1.4K0
发布2021-04-23 14:39:10
举报
文章被收录于专栏:DrugScience

组蛋白甲基转移酶模拟含水

介绍 步骤 step 1 使用pdbfixer来处理蛋白 step 2.1 溶剂化以及使用TIP3P进行模拟 step 2.2 溶剂化以及使用TIP4P-Ew进行模拟 step 3 导出文件

介绍

OpenMM包含一系列的水模型,例如:TIP3P, TIP4P-ew, TIP5P。 本例说明了OpenMM的使用不同的水模型,来模拟对于组蛋白甲基转移酶SET7/9(Uniprot:Q8WTS6)反应性十分重要的水通道。蛋白已经去除了配体。

步骤

step 1 使用pdbfixer来处理蛋白

获取1O9S的PDB文件,除去不是想要的链(将二体转化为单体,除去配体),加入缺失氨基酸(位于链中间),以及一些缺失的重原子。保留所有的水分子,因为我们对于水分子很感兴趣。

代码语言:javascript
复制
from pdbfixer import PDBFixer
from simtk.openmm.app import PDBFile

fixer = PDBFixer('1o9s.pdb')
# 看下fixer里面几条链
print(list(fixer.topology.chains()))
# 查看fixer里面链1有几个氨基酸,并转化为列表
residue_list = list(fixer.topology._chains[1].residues())
# 除去不相关的氨基酸
fixer.removeChains([1,2,3,4,5,7,9])
# 查看fixer的topo
print(fixer.topology)

# 寻找缺失残基
fixer.findMissingResidues()

# 在链的中间补充缺失缺失氨基酸,而不是末端
chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
missingResidues = dict()
for key in keys:
    chain = chains[key[0]]
    if not (key[1] == 0 or key[1] == len(list(chain.residues()))):
        missingResidues[key] = fixer.missingResidues[key]
fixer.missingResidues = missingResidues

fixer.findMissingAtoms()
fixer.addMissingAtoms()

# 输出
PDBFile.writeFile(fixer.topology, fixer.positions, open('1o9s_fixed.pdb', 'w'))

我们导入蛋白文件,并使用openmm app.Modeller为其加氢,增加溶剂环境. 我们将会使用以下水模型:TIP3P, TIP4P-ew, and TIP5P。 对于后两个水模型,我们将会执行一个额外的步骤,向晶体水中增加额外的水粒子。我们将会使用amber99sbildn力场-此力场通过XML文件创造。整个系统被创造,使用createSystem方法。接下来我们将会使用来源于Modeller对象的topology以及positions来设置LangevinIntegrator以及进行模拟。

在本案例中,我们将会使用CPU,混合精度。 此模拟将会进行能量最小化,平衡100步。 模拟10ps。

step 2.1 溶剂化以及使用TIP3P进行模拟

代码语言:javascript
复制
from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

# 导入准备好的蛋白文件
# 力场文件通过xml导入
pdb = app.PDBFile('1o9s_fixed.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml')
# 使用app.Modeller来进行加氢以及溶剂化
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
# 溶剂化,默认阳离子 'Na+' ,阴离子 'Cl-'
# 使用1ns的填充距离
modeller.addSolvent(forcefield, model='tip3p', padding=1.0*unit.nanometers)
# 输出文件
app.PDBFile.writeFile(modeller.topology, modeller.positions, open('1o9s_modeller_tip3p.pdb', 'w'))

# 准备系统以及integrator
system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
    2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

# 准备模拟
# CPU,混合精度
platform = mm.Platform.getPlatformByName('CPU')
properties = {'CpuThreads': 'mixed'}
simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)

# 最小化
print('Minimizing...')
simulation.minimizeEnergy()

# 平衡 100 steps 
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)

# 增加 reporters
simulation.reporters.append(app.DCDReporter('trajectory_tip3p.dcd', 100, enforcePeriodicBox=False))
simulation.reporters.append(app.pdbreporter.PDBReporter('trajectory_tip3p.pdb', 100,enforcePeriodicBox=False))
simulation.reporters.append(app.StateDataReporter('trajectory_tip3p.csv', 100, step=True,
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
    speed=True,totalSteps=5000, separator=','))


# 10ps的模拟
print('Running Production...')
simulation.step(5000)
print('Done!')

step 2.2 溶剂化以及使用TIP4P-Ew进行模拟

使用4点水模型TIP4P-Ew,需要我们使用modeller.addExtraParticles来模拟virtual sites 和上面十分相似,我就不加注释了

代码语言:javascript
复制
from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

# load in input PDB file and force field XML files
pdb = app.PDBFile('1o9s_fixed.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'tip4pew.xml')

# use app.Modeller to add hydrogens, extra particles, and solvent
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens()
modeller.addExtraParticles(forcefield)
modeller.addSolvent(forcefield, model='tip4pew', padding=1.0*unit.nanometers)
app.PDBFile.writeFile(modeller.topology, modeller.positions, open('1o9s_modeller_tip4pew.pdb', 'w'))

# prepare system and integrator
system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
    2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

# prepare simulation
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}
simulation = app.Simulation(modeller.topology, system, integrator, platform,
    properties)
simulation.context.setPositions(modeller.positions)

# minimize
print('Minimizing...')
simulation.minimizeEnergy()

# equilibrate for 100 steps
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)

# append reporters
simulation.reporters.append(app.DCDReporter('trajectory_tip4pew.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
    speed=True, totalSteps=25000000, separator='\t'))

# run 50 ns of production simulation
print('Running Production...')
simulation.step(25000000)
print('Done!')

step 3

你会得到一个文件trajectory_tip3p.dcd,一个文件trajectory_tip3p.pdb,trajectory_tip3p.csv,这些文件会包含模拟的信息

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-04-19,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 FindKey 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档