北鯤教程 | 利用Amber熱力學積分計算相對自由能變化


上周四,何博士為大家在北鯤云的直播間分享了Amber熱力學積分計算相對自由能變化(直播回放可在視頻號:北鯤云-直播回放中查看)。

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

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

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

本文所需的所有文件請在https://github.com/Xinheng-He/ti_toturial上下載。
北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖1
該應用場景解決將蛋白口袋內的小分子A變為小分子B所產生的相對自由能變。

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖2

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

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

   
   
save * my_case\ ben_ligand .mol2
   
   
save *\ my_case\ benfen_ligand .mol2

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

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖3

(保存pymol中的sele對象)

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

   
   
module add Anaconda3/2020.02
source /public/software/. local/easybuild/software/amber/aber20/amber.sh
ulimit -s unlimited
ulimit -l unlimited

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

   
   
antechamber -i ben_ligand.mol2 -fi mol2 -o ben_real.gaff2.mol2 -fo mol2 -rn BEN -at gaff2 -an  yes -dr  no -pf  yes -c bcc -nc  0  

生成frcmod力場參數文件。


parmchk2  -i  ben_real .gaff2 .mol2  -f  mol2  -o  BEN .gaff2 .frcmod   -s   gaff2   -a   yes

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

   
   
antechamber  -i  benfen_ligand .mol2  -fi  mol2  -o  benfen_real .gaff2 .mol2  -fo  mol2  -rn  FEN  -at  gaff2  -an  yes  -dr  no  -pf  yes  -c  bcc  -nc 0
parmchk2  -i  benfen_real .gaff2 .mol2  -f  mol2  -o  FEN .gaff2 .frcmod  -s  gaff2  -a  yes

根據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中的原子變成真實的樣子,這樣做是為了保證一樣原子的位置完全一致,減少不必要的變量。

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖4

(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并檢查其中的配體分子區別,是否只是不一樣的配體部分不一樣,確定配體不一樣部分的原子。設置出發和結束原子。

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖5

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

修改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  -p  complex_vdw_bonded .parm7  -c  complex_vdw_bonded .rst7  -find  :FEN

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

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

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

   
   
p  ligands_vdw_bonded .rst7  ligands_vdw_bonded .rst7 .leap
cp  press_lig .rst7  ligands_vdw_bonded .rst7
cp  complex_vdw_bonded .rst7  complex_vdw_bonded .rst7 .leap
cp  press_com .rst7  complex_vdw_bonded .rst7

cpptraj  -p  ligands_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,提交命令:

   
   
python  run_mpi .py  ligands 500000  mpi

將會檢查所有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下運行。

   
   
pip2  install matplotlib
pip2  install scipy
pip2  install numpy
pip2  install pymbar== 3.0 .3

運行:

   
   
mkdir -p ana_recharge &&  cd ana_recharge
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py -a AMBER -d . -p ../recharge/[01]*/ti00[1-9] -q out -o . -t 300 -v -r 5 -u kcal -f 50 -g -w
mkdir -p ../ana_decharge &&  cd ../ana_decharge
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py -a AMBER -d . -p ../decharge/[01]*/ti00[1-9] -q out -o . -t 300 -v -r 5 -u kcal -f 50 -g -w
mkdir -p ../ana_vdw &&  cd ../ana_vdw
../../alchemical-analysis/alchemical_analysis/alchemical_analysis.py -a AMBER -d . -p ../vdw_bonded/[01]*/ti00[1-9] -q out -o . -t 300 -v -r 5 -u kcal -f 50 -g -w

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

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖6

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

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

第162行。

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖7

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖8
該應用場景解決將蛋白口袋內的殘基A變為殘基B所產生的相對自由能變。
北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖9
將蛋白口袋內的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
source leaprc.protein.ff14SB
source leaprc.gaff2
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
zn = loadpdb WT.pdb
check zn
solvateBox zn TIP3PBOX 10
addIons zn Cl- 0
savepdb zn box_check.pdb
quit

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

北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖10

   
   
python  dry_for_TI .py  box_check .pdb  wat .pdb  WT_receptor .pdb

同樣的方式讀取突變型結構,使其保持Amber的原子順序。

   
   
tleap
source leaprc.protein.ff14SB
source leaprc.gaff2
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
zn = loadpdb L86Q.pdb
check zn
solvateBox zn TIP3PBOX 10
addIons zn Cl- 0
savepdb zn L86Q_leap.pdb
quit
python dry_for_TI.py L86Q_leap.pdb wat1.pdb L86Q_dry.pdb

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

   
   
python  check_diff_online .py  L84Q  L86Q_dry .pdb  WT_receptor .pdb 84  L86Q_check .pdb  WT_check .pdb

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

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

   
   
tleap
source leaprc.protein.ff14SB
source leaprc.gaff2
loadOff FEN.lib
loadamberparams FEN.gaff2.frcmod
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
ligand = loadmol2 FEN.gaff2.mol2 
m1 = loadpdb L86Q_check.pdb
m2 = loadpdb WT_check.pdb
w = loadpdb wat1.pdb
protein = combine {m1 m2 w}
complex = combine {m1 m2 ligand w}
set default nocenter on
setBox protein vdw {39.415 43.577 52.292}
savepdb protein protein.pdb
saveamberparm protein protein.parm7 protein.rst7

setBox complex vdw {39.415 43.577 52.292}
savepdb complex complex.pdb
saveamberparm complex complex.parm7 complex.rst7
quit

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

   
   
parmed  protein .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

正確獲得這些文件后,使用:

   
   
python  auto_gene_inp_run .py 84 162  CA, C, O, N, H, HA, CB  L86Q

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

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

運行結束后的分析略。
北鯤教程 | 利用Amber熱力學積分計算相對自由能變化的圖11
本案例將實際運行一個蛋白-蛋白相互作用上的突變。我們計算新冠病毒受體結合區域(rbd)到人ACE2受體(ace2)復合物上發生Omicron的突變之一的返回突變A484E后的結合自由能變化。

生成連接了二硫鍵的大復合體水盒,記錄vdw盒子大小,額外加入0.15M/L NaCL (總水數量*0.002772)。

   
   
tleap
source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
zn = loadpdb omi_SS.pdb
bond zn.333.SG zn.358.SG
bond zn.376.SG zn.429.SG
bond zn.388.SG zn.522.SG
bond zn.477.SG zn.485.SG
bond zn.637.SG zn.645.SG
bond zn.848.SG zn.865.SG
bond zn.1034.SG zn.1046.SG
check zn
solvateBox zn TIP3PBOX 10
addIons zn Na+ 0
addIons zn Na+ 80
addIons zn Cl- 0
savepdb zn box_check.pdb
quit

python dry_for_TI.py box_check.pdb wat1.pdb omi_rbd.pdb,omi_ace2.pdb

同時生成一個小的水盒,用于跑蛋白部分的TI。

   
   
tleap
source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
m1 = loadpdb omi_rbd.pdb
bond m1.4.SG m1.29.SG
bond m1.47.SG m1.100.SG
bond m1.59.SG m1.193.SG
bond m1.148.SG m1.156.SG
check m1
solvateBox m1 TIP3PBOX 10
addIons m1 Na+ 0
addIons m1 Na+ 28
addIons m1 Cl- 0
savepdb m1 ligands_recharge.pdb
quit

python dry_for_TI.py ligands_recharge.pdb rbd_wat.pdb test_rbd.pdb

檢查輸入的是否在TI區域之外沒有區別,注意輸入的順序,前面是突變后的蛋白,后面是原始的蛋白。

   
   
python  check_diff_online .py  A152E   A484E_rbd .pdb  omi_rbd .pdb 152  A484E_check .pdb  omi_check .pdb

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

   
   
tleap
source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadamberparams frcmod.ionsjc_tip3p
ligand = loadpdb omi_ace2.pdb
bond ligand.308.SG ligand.316.SG
bond ligand.519.SG ligand.536.SG
bond ligand.705.SG ligand.717.SG
m1 = loadpdb omi_check.pdb 
bond m1.4.SG m1.29.SG
bond m1.47.SG m1.100.SG
bond m1.59.SG m1.193.SG
bond m1.148.SG m1.156.SG
m2 = loadpdb A484E_check.pdb
bond m2.4.SG m2.29.SG
bond m2.47.SG m2.100.SG
bond m2.59.SG m2.193.SG
bond m2.148.SG m2.156.SG
w1 = loadpdb rbd_wat.pdb
w2 = loadpdb wat1.pdb
protein = combine {m1 m2 w1}
complex = combine {m1 m2 ligand w2}
set default nocenter on
setBox protein vdw {43.215 53.421 59.922}
savepdb protein protein.pdb
saveamberparm protein protein.parm7 protein.rst7

setBox complex vdw {64.171 75.490 114.587}
savepdb complex complex.pdb
saveamberparm complex complex.parm7 complex.rst7
quit

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

   
   
parmed  protein .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

正確獲得這些文件后,使用:

   
   
python  auto_gene_inp_run .py 152 193  CA, C, O, N, H, HA, CB  A484E

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

登錄后免費查看全文
立即登錄
App下載
技術鄰APP
工程師必備
  • 項目客服
  • 培訓客服
  • 平臺客服

TOP