?

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

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

本教程要求使用者一定程度了解Amber動力學模擬程序。

Amber是美國加州大學Peter Kollman等開發的一款著名的分子動力學模擬軟件包。Amber主要適用于蛋白質，小分子和多糖等生物分子體系的模擬。

本文所需的所有文件請在https://github.com/Xinheng-He/ti_toturial上下載。

該應用場景解決將蛋白口袋內的小分子A變為小分子B所產生的相對自由能變。

?

將蛋白口袋內的苯轉化為苯酚。

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

save*my_caseben_ligand.mol2save*my_casebenfen_ligand.mol2

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

?

（保存pymol中的sele對象）

開啟一個北鯤云管理節點加載環境。

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

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

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

生成frcmod力場參數文件。

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

上述操作對苯酚再來一次。

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

根據frcmod文件，生成兩個分子的文庫文件，該文件描述了分子內部的原子類型和鍵連信息。

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

注意前后的力場要保持一致。

將兩個pdb文件（FEN.pdb和BEN.pdb）中同樣位置的全部重原子，保存成同樣的坐標，注意名字要和lib中的一樣，放成一個lig.pdb，在下面的tleap過程中，tleap會自動根據lib文件將complex中的原子變成真實的樣子，這樣做是為了保證一樣原子的位置完全一致，減少不必要的變量。

?（pdb文件的內容）

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

pdb4amber pure_protein.pdb -o pure_check.pdb cat pure_check_sslink

沒有二硫鍵，之后使用pure_check.pdb。

接下來在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

注意根據實際電荷情況調整addions，如果ligand/complex帶負電，加Na+，反之加Cl-，離子類型可以自己選擇。

下載complex_vdw_bonded.pdb并檢查其中的配體分子區別，是否只是不一樣的配體部分不一樣，確定配體不一樣部分的原子。設置出發和結束原子。

?

紅框是有區別的原子，其他原子需要手動編輯使其保持一致的坐標，紅線是編輯后的原子。

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

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

使6個in文件符合實際更改的原子情況，只改變不同的原子，該步驟的目標是優化vdw變化過程中氫原子的位置。

使用sbatch run_v100.slurm，將任務提交到北鯤云的單卡V100集群上，該任務大概耗時1分鐘，注意該任務如果有報錯，可能是以下問題：

如果上述過程報關于ambmask的錯（比如mask中不能用_符號），可以改寫ambmask來獲得正確可識別的mask。

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

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

如果報錯Error : Atom 12 does not have match in V1 ，說明這個原子在兩個小分子配體中間的位置區別太大，TI不能識別這兩個分子作為一樣的背景，因此將這兩個原子（初始和結束配體的）加入坐標中，就可以解決這個問題。

運行結束后，提取優化過后的分子結構。

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

生成好這些過程的文件后，使用setup_run.sh來產生三個步驟的輸入文件，在修改setup_run.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的都改成目標的ambmask，：前面的1或者2不要改。

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

該文件將直接生成所需的slurm文件，并提交到對應的機器上，默認使用g-v100-1，運行pmemd.cuda，大致運行時間1小時左右。

有時候pmemd.cuda會運行失敗，此時轉用cpu來運行，使用run_mpi.py，提交命令：

pythonrun_mpi.pyligands500000mpi

將會檢查所有ligands下的文件，對于5分鐘內沒更新，且info中運行步驟在500000 以下的，會提交到32核CPU機器上運行后續的模擬，直到結束。

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

運行結束后，使用alchemical-analysis/alchemical_analysis/alchemical_analysis.py來分析結果，注意該文件在python2下運行。

pip2installmatplotlib
pip2installscipy
pip2installnumpy
pip2installpymbar==3.0.3

運行：

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

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

如果見到pymbar的warning，只要注釋掉對應的assert就可以了。

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

去第162行。

該應用場景解決將蛋白口袋內的殘基A變為殘基B所產生的相對自由能變。

將蛋白口袋內的LEU轉化GLN。

首先，使用pymol將分子打開，并選中殘基，使用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，讀取野生型結構。

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

記錄盒子的范德華半徑，并將結構中的水提取出來備用。

pythondry_for_TI.pybox_check.pdbwat.pdbWT_receptor.pdb

同樣的方式讀取突變型結構，使其保持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

對比兩個去水后的文件，發現由于Amber重編號，突變的殘基變為84位，此時使用check_diff_online.py，來保證不是突變的殘基的位置都絕對一致，以允許進行TI過程。

pythoncheck_diff_online.pyL84QL86Q_dry.pdbWT_receptor.pdb84L86Q_check.pdbWT_check.pdb

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

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

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運算，此處的162為WT或者突變體的殘基數，@后面的內容是python get_mutation.py L86Q得到的mapping結果，是在那個突變殘基上但也沒有變化的殘基，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

來生成tmpl文件（run.tmpl文件需要上傳），并將tmpl文件轉成真正的ti文件放進文件夾，同時使用slurm提交最小化，加熱和運行步驟。

注意，這里直接使用cuda很容易斷，可以適當自己修改之前的run_mpi.py來使用cpu續跑中斷的模擬。

運行結束后的分析略。

本案例將實際運行一個蛋白-蛋白相互作用上的突變。我們計算新冠病毒受體結合區域（rbd）到人ACE2受體（ace2）復合物上發生Omicron的突變之一的返回突變A484E后的結合自由能變化。

生成連接了二硫鍵的大復合體水盒，記錄vdw盒子大小，額外加入0.15M/L NaCL （總水數量*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

同時生成一個小的水盒，用于跑蛋白部分的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區域之外沒有區別，注意輸入的順序，前面是突變后的蛋白，后面是原始的蛋白。

pythoncheck_diff_online.pyA152EA484E_rbd.pdbomi_rbd.pdb152A484E_check.pdbomi_check.pdb

設定正確的二硫鍵，分別加載兩種水盒，輸出protein和complex的拓撲學文件和坐標文件。

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

來生成tmpl文件（run.tmpl文件需要上傳），并將tmpl文件轉成真正的ti文件放進文件夾，同時使用slurm提交最小化，加熱和運行步驟（5ns）。

?
審核編輯黃昊宇