AmberTools将被用来生成输入文件,OpenMM 将被用来运行模拟,模拟平台为在线百度AI Stuio, 并使用GPU加速。
First thing’s first, 到PDB 蛋白数据库下载一需要模拟的靶点晶体,备用。
1. H++ web server ,输入蛋白编码,预处理蛋白,下载 后缀为 crd和top 两个文件备用 (Note: If you use it regularly, please register)
(如果第一步生成失败了,请阅读一下原因,如果是因为氨基酸序列缺失,backbone boken,or missing residue,调到第5步,使用openmm-setup修复下载,再回来走完整合流程。)
2. 在线(AI Studio)使用conda 安装一个独立的python3.9 环境,取名可为AmberTools (可能需要重启命令终端) ,如有询问y/n,输入y
(请自行用手机号提前注册一个百度AI studio平台的账号,新建一个"项目",点击"启动环境",这样,你在云端就有了一台性能不错的虚拟机)
conda create --name AmberTools python=3.9
- 激活python环境, 并安装AmberTools21,如有询问y/n,输入y (!注意,如果你是用本地计算机需要用conda activate 但是在百度的平台我们需要用sourcre activate)
source activate AmberTools
conda install -c conda-forge ambertools=21 compilers
运行AmberTools 所包含的一些列工具处理小分子,合并小分子和蛋白,生成新的crd和top文件(相对于第1步)
- Windows本地电脑安装Anaconda,然后启动Anaconda命令行,输入(与第3步类似),并安装一个用于准备OpenMM动力学输入文件的图形界面工具,如有询问y/n,输入y
conda create --name openmm python=3.9 # 这一句命令后,可能需要关掉再打开
conda activate openmm
(openmm)conda install -c conda-forge openmm-setup
5.启动openmm-setup图形界面,你的默认浏览器将会自动弹出成为一个动力学准备软件界面, 输入openmm-setup,回车即可
(openmm)openmm-setup
- 回到百度AI Stuido, 在公开项目中搜索”OpenMM-蛋白-小分子动力学模拟“项目,喜欢,点赞后,fork到你自己的账号中,点击,”运行一下“。上传第5步生成的输入文件,开始模拟。
- 下载模拟生成的output.pdb或者output.dcd到本地计算机,使用VMD进行分析,分析的教程请自行搜索,再次不在赘述。
步骤3:(在继续之前,请验证你已经成功完成了ambertools的安装)
conda activate ambertools
注意,第1步下载的两个文件的名称分别为 0.15_80_10_pH6.5_1ODX.top,和0.15_80_10_pH6.5_1ODX.crd, 数字代表了生成这些文件的H++web server条件。这些名字是什么不重要,你可以随意修改。下面一条命令是利用这两个文件来重新生成一个蛋白pdb,我们后面的模拟就用这个新的蛋白。新蛋白的名字为0.15_80_10_pH6.5_1ODX.pdb
ambpdb -p 0.15_80_10_pH6.5_1ODX.top -c 0.15_80_10_pH6.5_1ODX.crd > 0.15_80_10_pH6.5_1ODX.pdb
接下来,我们需要用到最最开始我们下载的那个PDB结构了,虽然这个结构中的蛋白部分我们不再需要(因为我们刚刚新生成了一个),但是小分子还是要从这个结构中提取。并对提取的这个小分子进行加氢,或者还原,并对其格式进行正规化矫正。从这里开始,我们下面将一共生成6个与小分子有关的结构,即ligand1到ligand6,注意命名方式。
awk '$4=="0E8"' 1ODX.pdb > ligand1.pdb # 该命令意思为提取pdb文件的第四列,小分子在这里
reduce ligand1.pdb > ligand2.pdb # 还原,也就是向小分子加氢
#如果加氢失败,请使用open babel 加氢。
#conda install -c conda-forge openbabel
#obabel -ipdb ligand1.pdb -opdb -O ligand2.pdb -h
pdb4amber -i ligand2.pdb -o ligand3.pdb # 格式美化,或者格式矫正
(选择性步骤:用txt编辑器打开确认你的clean版本的小分子pdb文件已经加氢,并且其格式不同于格式标准化之前的版本)
接下来我们需要生成小分子的mol2格式文件
(注意!如果小分子含有氧原子如酮,需要特别注意加氢是否正确,特别是如果使用pymol代替reduce加氢的情况下)
antechamber -fi pdb -i ligand3.pdb -fo mol2 -o ligand4.mol2 -c bcc -pf y
#This step may take a little bit of time, upon finish, check "sqm.out", the last line should be
# -------------calculation complete---------------------
antechamber -i ligand4.mol2 -fi mol2 -o ligand5.prepi -fo prepi -pf y
parmchk2 -f prepi -i ligand5.prepi -o ligand6.frcmod
(ATTENTION: You DON’T have to generate prepi file, you could move on with ONLY mol2 file, but then you need modify “http://tleap.in” file accordingly.)
……………………………………………………………………………………………………………………………………
合并小分子和蛋白,并对复合物pdb格式进行标准化
cat 0.15_80_10_pH6.5_1ODX.pdb 0E8_clean_H.pdb > 1ODX_H.pdb
pdb4amber -i 1ODX_H.pdb -o 1ODX_clean_H.pdb
新建一个文件,名称与格式为 tleap.in, 写入下面的内容,保存
source leaprc.protein.ff14SB #Source leaprc file for ff14SB protein force field
source leaprc.gaff #Source leaprc file for gaff
source leaprc.water.tip3p #Source leaprc file for TIP3P water model
loadamberprep ligand5.prepi #Load the prepi file for the ligand
loadamberparams ligand6.frcmod #Load the additional frcmod file for ligand
mol = loadpdb 1ODX_clean_H.pdb #Load PDB file for protein-ligand complex
solvatebox mol TIP3PBOX 8 #Solvate the complex with a cubic water box
addions mol Cl- 0 #Add Cl- ions to neutralize the system
saveamberparm mol 1ODX.prmtop 1ODX.inpcrd #Save AMBER topology and coordinate files
quit #Quit tleap program
(As we have mentioned earlier, if you do not have a “prepi” file, that is fine, what you need do is delete the “loadamberprep” line and add a new line for loading the mol2 file directly.)
在命令终端使用AmberTools的tleap 运行这个文件
tleap -s -f 1ODX_tleap.in > 1ODX_tleap.out
这样我们就得到了新的两个文件1ODX.prmtop 1ODX.inpcrd。这两个文件将是第5步的输入文件。
步骤3是最容易出错也是最重要的一个步骤,请参考
Tutorial for the LEaP Program
选择 Amber 文件,然后上传上面得到的两个文件,OpenMM 会要求用户输入模拟的步数,即时长,模拟的溶剂类型,以及很多其他参数,正常情况下是都有一个默认值已经填好的,依次点击下一步,平台 platform 选择 CUDA,如果你可用的平台没有GPU则更改CUDA为CPU即可。最后保存, ”save all files" ,得到一个压缩文件。里面包含了我们openmm 模拟需要的文件和脚本。如果你不知道如何修改界面上的参数,请全部保持默认即可!
AmberInpcrdFile('SYS_gaff2.crd')
你需要把上一步生成的你的prmtop文件和inpcrd文件(或者crd)替换到以上两个括号里。
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
prmtop = AmberPrmtopFile('SYS_gaff2.prmtop')
inpcrd = AmberInpcrdFile('SYS_gaff2.crd')
system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
if inpcrd.boxVectors is not None:
simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('amber_output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
potentialEnergy=True, temperature=True))
simulation.step(10000)
你可以使用上面这段代码,也可以使用下面这一段,作用是一样的,模拟的参数文件稍有区别。 如果使用下面这一段,我们需要替换的两个文件是 “complex.prmtop" 和 “complex.inpcrd”
from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
# Input Files
prmtop = AmberPrmtopFile('complex.prmtop')
inpcrd = AmberInpcrdFile('complex.inpcrd')
# System Configuration
nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001
hydrogenMass = 1.5*amu
# Integration Options
dt = 0.004*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25
# Simulation Options
steps = 1000000
equilibrationSteps = 1000
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}
dcdReporter = DCDReporter('trajectory.dcd', 10000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 10000)
# Prepare the Simulation
print('Building system...')
topology = prmtop.topology
positions = inpcrd.positions
system = prmtop.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance, hydrogenMass=hydrogenMass)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(positions)
if inpcrd.boxVectors is not None:
simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
# Minimize and Equilibrate
print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)
# Simulate
print('Simulating...')
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.step(steps)
最理想的情况是,你自己根据自己的需要或者理解用Openmm生成适合自己蛋白体系的代码区块然后替换上面的例子。如果你会使用本地Jupyter notebook, 理论上也是可以打开以上脚本运行的,但墙裂不建议这么做,因为那样速度将非常非常慢,还容易出错。自己有Linux系统的,可以考虑把步骤2和3挪到本地,这样会更容易一些。如:
https://github.com/quantaosun/Ambertools-OpenMM-MD