利用Amber熱力學(xué)積分計(jì)算相對(duì)自由能變化

?

上周四，何博士為大家在北鯤云的直播間分享了Amber熱力學(xué)積分計(jì)算相對(duì)自由能變化）。

直播結(jié)束后有很多小伙伴來(lái)向我們要PPT資料，這里何博士也為大家準(zhǔn)備了文字版本的教程。將為大家介紹如何在北鯤云計(jì)算平臺(tái)上利用Amber熱力學(xué)積分計(jì)算相對(duì)自由能變化，體系包括小分子-蛋白（小分子改變），小分子-蛋白（蛋白突變），蛋白-蛋白相互作用。

本教程要求使用者一定程度了解Amber動(dòng)力學(xué)模擬程序。

Amber是美國(guó)加州大學(xué)Peter Kollman等開(kāi)發(fā)的一款著名的分子動(dòng)力學(xué)模擬軟件包。Amber主要適用于蛋白質(zhì)，小分子和多糖等生物分子體系的模擬。

本文所需的所有文件請(qǐng)?jiān)趆ttps://github.com/Xinheng-He/ti_toturial上下載。

該應(yīng)用場(chǎng)景解決將蛋白口袋內(nèi)的小分子A變?yōu)樾》肿覤所產(chǎn)生的相對(duì)自由能變。

?

將蛋白口袋內(nèi)的苯轉(zhuǎn)化為苯酚。

首先，使用pymol將分子打開(kāi)，并選中小分子，保存為mol2文件，如下圖所示，我們使用

save*my_caseben_ligand.mol2save*my_casebenfen_ligand.mol2

這兩個(gè)命令保存了變化前后的配體。

?

（保存pymol中的sele對(duì)象）

開(kāi)啟一個(gè)北鯤云管理節(jié)點(diǎn)加載環(huán)境。

moduleaddAnaconda3/2020.02
source/public/software/.local/easybuild/software/amber/aber20/amber.sh
ulimit-sunlimited
ulimit-lunlimited

對(duì)苯生成具有電荷的可用mol2文件，總電荷為0，殘基名為BEN。

antechamber-iben_ligand.mol2-fimol2-oben_real.gaff2.mol2-fomol2-rnBEN-atgaff2-anyes-drno-pfyes-cbcc-nc0

生成frcmod力場(chǎng)參數(shù)文件。

parmchk2-iben_real.gaff2.mol2-fmol2-oBEN.gaff2.frcmod-sgaff2-ayes

上述操作對(duì)苯酚再來(lái)一次。

antechamber-ibenfen_ligand.mol2-fimol2-obenfen_real.gaff2.mol2-fomol2-rnFEN-atgaff2-anyes-drno-pfyes-cbcc-nc0
parmchk2-ibenfen_real.gaff2.mol2-fmol2-oFEN.gaff2.frcmod-sgaff2-ayes

根據(jù)frcmod文件，生成兩個(gè)分子的文庫(kù)文件，該文件描述了分子內(nèi)部的原子類型和鍵連信息。

tleap-f-<<_EOF
source?leaprc.gaff2
loadamberparams?FEN.gaff2.frcmod
FEN?=?loadmol2?benfen_real.gaff2.mol2?
saveoff?FEN?FEN.lib
savepdb?FEN?FEN.pdb
quit
_EOF

tleap?-f?-?<<_EOF
source?leaprc.gaff2
loadamberparams?BEN.gaff2.frcmod
BEN?=?loadmol2?ben_real.gaff2.mol2?
saveoff?BEN?BEN.lib
savepdb?BEN?BEN.pdb
quit
_EOF

注意前后的力場(chǎng)要保持一致。

將兩個(gè)pdb文件（FEN.pdb和BEN.pdb）中同樣位置的全部重原子，保存成同樣的坐標(biāo)，注意名字要和lib中的一樣，放成一個(gè)lig.pdb，在下面的tleap過(guò)程中，tleap會(huì)自動(dòng)根據(jù)lib文件將complex中的原子變成真實(shí)的樣子，這樣做是為了保證一樣原子的位置完全一致，減少不必要的變量。

?（pdb文件的內(nèi)容）

使用pdb4amber，檢查蛋白是否有二硫鍵，或需要編輯的殘基。

pdb4amber pure_protein.pdb -o pure_check.pdb cat pure_check_sslink

沒(méi)有二硫鍵，之后使用pure_check.pdb。

接下來(lái)在tleap中加載配體和受體。

tleap-f-<<_EOF
source?leaprc.protein.ff14SB
source?leaprc.gaff2
source?leaprc.water.tip3p
loadAmberParams?frcmod.ionsjc_tip3p

loadoff?BEN.lib
loadoff?FEN.lib
loadamberparams?BEN.gaff2.frcmod
loadamberparams?FEN.gaff2.frcmod


ligands?=?loadpdb?lig.pdb
complex?=?loadpdb?pure_check.pdb
complex?=?combine?{ligands?complex}
check?complex

solvatebox?ligands?TIP3PBOX?15
addions?ligands?Na+?0
savepdb?ligands?ligands_vdw_bonded.pdb
saveamberparm?ligands?ligands_vdw_bonded.parm7?ligands_vdw_bonded.rst7

solvatebox?complex?TIP3PBOX?15
addions?complex?Cl-?0
savepdb?complex?complex_vdw_bonded.pdb
saveamberparm?complex?complex_vdw_bonded.parm7?complex_vdw_bonded.rst7?

quit
_EOF

注意根據(jù)實(shí)際電荷情況調(diào)整addions，如果ligand/complex帶負(fù)電，加Na+，反之加Cl-，離子類型可以自己選擇。

下載complex_vdw_bonded.pdb并檢查其中的配體分子區(qū)別，是否只是不一樣的配體部分不一樣，確定配體不一樣部分的原子。設(shè)置出發(fā)和結(jié)束原子。

?

紅框是有區(qū)別的原子，其他原子需要手動(dòng)編輯使其保持一致的坐標(biāo)，紅線是編輯后的原子。

修改initial中的in文件下述部分。

timask1=':BEN',timask2=':FEN',
scmask1=':BEN@H6',scmask2=':FEN@O1,H6'

使6個(gè)in文件符合實(shí)際更改的原子情況，只改變不同的原子，該步驟的目標(biāo)是優(yōu)化vdw變化過(guò)程中氫原子的位置。

使用sbatch run_v100.slurm，將任務(wù)提交到北鯤云的單卡V100集群上，該任務(wù)大概耗時(shí)1分鐘，注意該任務(wù)如果有報(bào)錯(cuò)，可能是以下問(wèn)題：

如果上述過(guò)程報(bào)關(guān)于ambmask的錯(cuò)（比如mask中不能用_符號(hào)），可以改寫ambmask來(lái)獲得正確可識(shí)別的mask。

ambmask-pcomplex_vdw_bonded.parm7-ccomplex_vdw_bonded.rst7-find:FEN

是ambmask的輸入方式，在parm7和rst7中尋找這些原子，并用pdb的形式輸出。

如果報(bào)錯(cuò)Error : Atom 12 does not have match in V1 ，說(shuō)明這個(gè)原子在兩個(gè)小分子配體中間的位置區(qū)別太大，TI不能識(shí)別這兩個(gè)分子作為一樣的背景，因此將這兩個(gè)原子（初始和結(jié)束配體的）加入坐標(biāo)中，就可以解決這個(gè)問(wèn)題。

運(yùn)行結(jié)束后，提取優(yōu)化過(guò)后的分子結(jié)構(gòu)。

pligands_vdw_bonded.rst7ligands_vdw_bonded.rst7.leap
cppress_lig.rst7ligands_vdw_bonded.rst7
cpcomplex_vdw_bonded.rst7complex_vdw_bonded.rst7.leap
cppress_com.rst7complex_vdw_bonded.rst7

cpptraj-pligands_vdw_bonded.parm7<<_EOF
trajin?ligands_vdw_bonded.rst7
strip?":1,2"
outtraj?ligands_solvated.pdb?onlyframes?1
unstrip
strip?":2-999999"
outtraj?ligands_BEN.pdb?onlyframes?1
unstrip
strip?":1,3-999999"
outtraj?ligands_FEN.pdb?onlyframes?1
_EOF

cpptraj?-p?complex_vdw_bonded.parm7?<<_EOF
trajin?complex_vdw_bonded.rst7
strip?":1,2"
outtraj?complex_solvated.pdb?onlyframes?1
unstrip
strip?":2-999999"
outtraj?complex_BEN.pdb?onlyframes?1
unstrip
strip?":1,3-999999"
outtraj?complex_FEN.pdb?onlyframes?1
_EOF

再次使用tleap生成decharge,vdw和recharge的文件,decharge是配體1，recharge是配體2，修改閱讀的pdb的名字即可。

tleap-f-<<_EOF
#?load?the?AMBER?force?fields
source?leaprc.protein.ff19SB
source?leaprc.gaff2
source?leaprc.water.tip3p
loadAmberParams?frcmod.ionsjc_tip3p

loadOff?BEN.lib
loadOff?FEN.lib
loadamberparams?BEN.gaff2.frcmod
loadamberparams?FEN.gaff2.frcmod

#?coordinates?for?solvated?ligands?as?created?previously?by?MD
lsolv?=?loadpdb?ligands_solvated.pdb
lbnz?=?loadpdb?ligands_BEN.pdb
lphn?=?loadpdb?ligands_FEN.pdb

#?coordinates?for?complex?as?created?previously?by?MD
csolv?=?loadpdb?complex_solvated.pdb
cbnz?=?loadpdb?complex_BEN.pdb
cphn?=?loadpdb?complex_FEN.pdb

#?decharge?transformation
decharge?=?combine?{lbnz?lbnz?lsolv}
setbox?decharge?vdw
savepdb?decharge?ligands_decharge.pdb
saveamberparm?decharge?ligands_decharge.parm7?ligands_decharge.rst7

decharge?=?combine?{cbnz?cbnz?csolv}
setbox?decharge?vdw
savepdb?decharge?complex_decharge.pdb
saveamberparm?decharge?complex_decharge.parm7?complex_decharge.rst7

#?recharge?transformation
recharge?=?combine?{lphn?lphn?lsolv}
setbox?recharge?vdw
savepdb?recharge?ligands_recharge.pdb
saveamberparm?recharge?ligands_recharge.parm7?ligands_recharge.rst7

recharge?=?combine?{cphn?cphn?csolv}
setbox?recharge?vdw
savepdb?recharge?complex_recharge.pdb
saveamberparm?recharge?complex_recharge.parm7?complex_recharge.rst7

quit
_EOF

生成好這些過(guò)程的文件后，使用setup_run.sh來(lái)產(chǎn)生三個(gè)步驟的輸入文件，在修改setup_run.sh時(shí)，注意以下部分。

decharge_crg=":2@H6"
vdw_crg=":1@H6|:2@O1,H6"
recharge_crg=":1@O1,H6"
decharge="ifsc=0,crgmask='$decharge_crg',"
vdw_bonded="ifsc=1,scmask1=':1@H6',scmask2=':2@O1,H6',crgmask='$vdw_crg'"
recharge="ifsc=0,crgmask='$recharge_crg',"

適配修改，注意H6的都改成初始的ambmask，O1,H6的都改成目標(biāo)的ambmask，：前面的1或者2不要改。

如有必要修改λ，改變prod.tmpl和setup.sh中的值（0.00922開(kāi)始的那一串）。

該文件將直接生成所需的slurm文件，并提交到對(duì)應(yīng)的機(jī)器上，默認(rèn)使用g-v100-1，運(yùn)行pmemd.cuda，大致運(yùn)行時(shí)間1小時(shí)左右。

有時(shí)候pmemd.cuda會(huì)運(yùn)行失敗，此時(shí)轉(zhuǎn)用cpu來(lái)運(yùn)行，使用run_mpi.py，提交命令：

pythonrun_mpi.pyligands500000mpi

將會(huì)檢查所有l(wèi)igands下的文件，對(duì)于5分鐘內(nèi)沒(méi)更新，且info中運(yùn)行步驟在500000 以下的，會(huì)提交到32核CPU機(jī)器上運(yùn)行后續(xù)的模擬，直到結(jié)束。

這個(gè)運(yùn)行步驟非常緩慢，使用32核CPU算完1納秒的步驟可能需要7-8天，可以嘗試先運(yùn)行一段CPU，再運(yùn)行一段GPU，改為python run_mpi.py ligands 500000 cuda即可。

運(yùn)行結(jié)束后，使用alchemical-analysis/alchemical_analysis/alchemical_analysis.py來(lái)分析結(jié)果，注意該文件在python2下運(yùn)行。

pip2installmatplotlib
pip2installscipy
pip2installnumpy
pip2installpymbar==3.0.3

運(yùn)行：

mkdir-pana_recharge&&cdana_recharge
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py-aAMBER-d.-p../recharge/[01]*/ti00[1-9]-qout-o.-t300-v-r5-ukcal-f50-g-w
mkdir-p../ana_decharge&&cd../ana_decharge
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py-aAMBER-d.-p../decharge/[01]*/ti00[1-9]-qout-o.-t300-v-r5-ukcal-f50-g-w
mkdir-p../ana_vdw&&cd../ana_vdw
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py-aAMBER-d.-p../vdw_bonded/[01]*/ti00[1-9]-qout-o.-t300-v-r5-ukcal-f50-g-w

此時(shí)只能輸出三種變化的結(jié)果，將其TI 一列加和得到最終的結(jié)果。

如果見(jiàn)到pymbar的warning，只要注釋掉對(duì)應(yīng)的assert就可以了。

vim/home/cloudam/.local/lib/python2.7/site-packages/pymbar/timeseries.py

去第162行。

該應(yīng)用場(chǎng)景解決將蛋白口袋內(nèi)的殘基A變?yōu)闅埢鵅所產(chǎn)生的相對(duì)自由能變。

將蛋白口袋內(nèi)的LEU轉(zhuǎn)化GLN。

首先，使用pymol將分子打開(kāi)，并選中殘基，使用wizard-mutagenesis-protein完成突變，或者使用命令行完成突變：

load*.pdb
cmd.wizard("mutagenesis")
cmd.do("refresh_wizard")
cmd.get_wizard().set_mode("GLN")
cmd.get_wizard().do_select("86/")
cmd.get_wizard().apply()
cmd.set_wizard("done")
save*out_name.pdb,enabled

將突變后的蛋白文件保存為pdb文件。

將突變前的蛋白（WT.pdb）和突變后的蛋白（L86Q.pdb）的pdb文件上傳。使用tleap，讀取野生型結(jié)構(gòu)。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff2
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
zn=loadpdbWT.pdb
checkzn
solvateBoxznTIP3PBOX10
addIonsznCl-0
savepdbznbox_check.pdb
quit

記錄盒子的范德華半徑，并將結(jié)構(gòu)中的水提取出來(lái)備用。

pythondry_for_TI.pybox_check.pdbwat.pdbWT_receptor.pdb

同樣的方式讀取突變型結(jié)構(gòu)，使其保持Amber的原子順序。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff2
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
zn=loadpdbL86Q.pdb
checkzn
solvateBoxznTIP3PBOX10
addIonsznCl-0
savepdbznL86Q_leap.pdb
quit
pythondry_for_TI.pyL86Q_leap.pdbwat1.pdbL86Q_dry.pdb

對(duì)比兩個(gè)去水后的文件，發(fā)現(xiàn)由于Amber重編號(hào)，突變的殘基變?yōu)?4位，此時(shí)使用check_diff_online.py，來(lái)保證不是突變的殘基的位置都絕對(duì)一致，以允許進(jìn)行TI過(guò)程。

pythoncheck_diff_online.pyL84QL86Q_dry.pdbWT_receptor.pdb84L86Q_check.pdbWT_check.pdb

print的是空字典，說(shuō)明所有的原子都是匹配的。

再次用tleap讀取受體和配體（之前第一部分保存的mol2文件），并讀取水盒子，其中l(wèi)igand是剛才保存的配體（不變部分），m1和m2分別是突變前后的部分（注意突變只改變側(cè)鏈，不改變主鏈），注意這里使用了剛才記錄的范德華半徑。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff2
loadOffFEN.lib
loadamberparamsFEN.gaff2.frcmod
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
ligand=loadmol2FEN.gaff2.mol2
m1=loadpdbL86Q_check.pdb
m2=loadpdbWT_check.pdb
w=loadpdbwat1.pdb
protein=combine{m1m2w}
complex=combine{m1m2ligandw}
setdefaultnocenteron
setBoxproteinvdw{39.41543.57752.292}
savepdbproteinprotein.pdb
saveamberparmproteinprotein.parm7protein.rst7

setBoxcomplexvdw{39.41543.57752.292}
savepdbcomplexcomplex.pdb
saveamberparmcomplexcomplex.parm7complex.rst7
quit

使用parmed處理protein.parm7和complex.parm7，以保證正確的所改變的位置提供給TI運(yùn)算，此處的162為WT或者突變體的殘基數(shù)，@后面的內(nèi)容是python get_mutation.py L86Q得到的mapping結(jié)果，是在那個(gè)突變殘基上但也沒(méi)有變化的殘基，84是突變位置，246是84+162。

parmedprotein.parm7<<_EOF
loadRestrt?protein.rst7
setOverwrite?True
tiMerge?:1-162?:163-324?:84&!@CA,C,O,N,H,HA,CB?:246&!@CA,C,O,N,H,HA,CB
outparm?merged_L86Q_protein.parm7?merged_L86Q_protein.rst7
quit
_EOF

parmed?complex.parm7?<<_EOF
loadRestrt?complex.rst7
setOverwrite?True
tiMerge?:1-162?:163-324?:84&!@CA,C,O,N,H,HA,CB?:246&!@CA,C,O,N,H,HA,CB
outparm?merged_L86Q_complex.parm7?merged_L86Q_complex.rst7
quit
_EOF

正確獲得這些文件后，使用：

pythonauto_gene_inp_run.py84162CA,C,O,N,H,HA,CBL86Q

來(lái)生成tmpl文件（run.tmpl文件需要上傳），并將tmpl文件轉(zhuǎn)成真正的ti文件放進(jìn)文件夾，同時(shí)使用slurm提交最小化，加熱和運(yùn)行步驟。

注意，這里直接使用cuda很容易斷，可以適當(dāng)自己修改之前的run_mpi.py來(lái)使用cpu續(xù)跑中斷的模擬。

運(yùn)行結(jié)束后的分析略。

本案例將實(shí)際運(yùn)行一個(gè)蛋白-蛋白相互作用上的突變。我們計(jì)算新冠病毒受體結(jié)合區(qū)域（rbd）到人ACE2受體（ace2）復(fù)合物上發(fā)生Omicron的突變之一的返回突變A484E后的結(jié)合自由能變化。

生成連接了二硫鍵的大復(fù)合體水盒，記錄vdw盒子大小，額外加入0.15M/L NaCL （總水?dāng)?shù)量*0.002772）。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
zn=loadpdbomi_SS.pdb
bondzn.333.SGzn.358.SG
bondzn.376.SGzn.429.SG
bondzn.388.SGzn.522.SG
bondzn.477.SGzn.485.SG
bondzn.637.SGzn.645.SG
bondzn.848.SGzn.865.SG
bondzn.1034.SGzn.1046.SG
checkzn
solvateBoxznTIP3PBOX10
addIonsznNa+0
addIonsznNa+80
addIonsznCl-0
savepdbznbox_check.pdb
quit

pythondry_for_TI.pybox_check.pdbwat1.pdbomi_rbd.pdb,omi_ace2.pdb

同時(shí)生成一個(gè)小的水盒，用于跑蛋白部分的TI。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
m1=loadpdbomi_rbd.pdb
bondm1.4.SGm1.29.SG
bondm1.47.SGm1.100.SG
bondm1.59.SGm1.193.SG
bondm1.148.SGm1.156.SG
checkm1
solvateBoxm1TIP3PBOX10
addIonsm1Na+0
addIonsm1Na+28
addIonsm1Cl-0
savepdbm1ligands_recharge.pdb
quit

pythondry_for_TI.pyligands_recharge.pdbrbd_wat.pdbtest_rbd.pdb

檢查輸入的是否在TI區(qū)域之外沒(méi)有區(qū)別，注意輸入的順序，前面是突變后的蛋白，后面是原始的蛋白。

pythoncheck_diff_online.pyA152EA484E_rbd.pdbomi_rbd.pdb152A484E_check.pdbomi_check.pdb

設(shè)定正確的二硫鍵，分別加載兩種水盒，輸出protein和complex的拓?fù)鋵W(xué)文件和坐標(biāo)文件。

tleap
sourceleaprc.protein.ff14SB
sourceleaprc.gaff
sourceleaprc.water.tip3p
loadamberparamsfrcmod.ionsjc_tip3p
ligand=loadpdbomi_ace2.pdb
bondligand.308.SGligand.316.SG
bondligand.519.SGligand.536.SG
bondligand.705.SGligand.717.SG
m1=loadpdbomi_check.pdb
bondm1.4.SGm1.29.SG
bondm1.47.SGm1.100.SG
bondm1.59.SGm1.193.SG
bondm1.148.SGm1.156.SG
m2=loadpdbA484E_check.pdb
bondm2.4.SGm2.29.SG
bondm2.47.SGm2.100.SG
bondm2.59.SGm2.193.SG
bondm2.148.SGm2.156.SG
w1=loadpdbrbd_wat.pdb
w2=loadpdbwat1.pdb
protein=combine{m1m2w1}
complex=combine{m1m2ligandw2}
setdefaultnocenteron
setBoxproteinvdw{43.21553.42159.922}
savepdbproteinprotein.pdb
saveamberparmproteinprotein.parm7protein.rst7

setBoxcomplexvdw{64.17175.490114.587}
savepdbcomplexcomplex.pdb
saveamberparmcomplexcomplex.parm7complex.rst7
quit

使用parmed處理protein.parm7和complex.parm7，以保證正確的位置。

parmedprotein.parm7<<_EOF
loadRestrt?protein.rst7
setOverwrite?True
tiMerge?:1-193?:194-386?:152&!@CA,C,O,N,H,HA,CB?:345&!@CA,C,O,N,H,HA,CB
outparm?merged_A484E_protein.parm7?merged_A484E_protein.rst7
quit
_EOF

parmed?complex.parm7?<<_EOF
loadRestrt?complex.rst7
setOverwrite?True
tiMerge?:1-193?:194-386?:152&!@CA,C,O,N,H,HA,CB?:345&!@CA,C,O,N,H,HA,CB
outparm?merged_A484E_complex.parm7?merged_A484E_complex.rst7
quit
_EOF

正確獲得這些文件后，使用：

pythonauto_gene_inp_run.py152193CA,C,O,N,H,HA,CBA484E

來(lái)生成tmpl文件（run.tmpl文件需要上傳），并將tmpl文件轉(zhuǎn)成真正的ti文件放進(jìn)文件夾，同時(shí)使用slurm提交最小化，加熱和運(yùn)行步驟（5ns）。

?
審核編輯黃昊宇