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

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

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

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

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

?

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

?

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

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

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

?

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

（保存pymol中的sele對象）

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

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力場參數(shù)文件。

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

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

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

?

根據(jù)frcmod文件，生成兩個(gè)分子的文庫文件，該文件描述了分子內(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 _EOFtleap?-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

?

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

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

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

使用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_tip3ploadoff?BEN.lib loadoff?FEN.lib loadamberparams?BEN.gaff2.frcmod loadamberparams?FEN.gaff2.frcmodligands?=?loadpdb?lig.pdb complex?=?loadpdb?pure_check.pdb complex?=?combine?{ligands?complex} check?complexsolvatebox?ligands?TIP3PBOX?15 addions?ligands?Na+?0 savepdb?ligands?ligands_vdw_bonded.pdb saveamberparm?ligands?ligands_vdw_bonded.parm7?ligands_vdw_bonded.rst7solvatebox?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ū)別的原子，其他原子需要手動編輯使其保持一致的坐標(biāo)，紅線是編輯后的原子。

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

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

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

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

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

ambmask?-p?complex_vdw_bonded.parm7?-c?complex_vdw_bonded.rst7?-find?:FEN

?

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

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

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

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.rst7cpptraj?-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 _EOFcpptraj?-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_tip3ploadOff?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.rst7decharge?=?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.rst7recharge?=?combine?{cphn?cphn?csolv} setbox?recharge?vdw savepdb?recharge?complex_recharge.pdb saveamberparm?recharge?complex_recharge.parm7?complex_recharge.rst7quit _EOF

?

生成好這些過程的文件后，使用setup_run.sh來產(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開始的那一串）。

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

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

python?run_mpi.py?ligands?500000?mpi

?

將會檢查所有l(wèi)igands下的文件，對于5分鐘內(nèi)沒更新，且info中運(yùn)行步驟在500000 以下的，會提交到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來分析結(jié)果，注意該文件在python2下運(yùn)行。

pip2?install?matplotlib pip2?install?scipy pip2?install?numpy pip2?install?pymbar==3.0.3

?

運(yùn)行：

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

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

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

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

去第162行。

?

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

將蛋白口袋內(nèi)的LEU轉(zhuǎn)化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，讀取野生型結(jié)構(gòu)。

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

?

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

?

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

?

同樣的方式讀取突變型結(jié)構(gòu)，使其保持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

?

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

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

?

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

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

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.rst7setBox?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運(yùn)算，此處的162為WT或者突變體的殘基數(shù)，@后面的內(nèi)容是python get_mutation.py L86Q得到的mapping結(jié)果，是在那個(gè)突變殘基上但也沒有變化的殘基，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 _EOFparmed?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文件轉(zhuǎn)成真正的ti文件放進(jìn)文件夾，同時(shí)使用slurm提交最小化，加熱和運(yùn)行步驟。

注意，這里直接使用cuda很容易斷，可以適當(dāng)自己修改之前的run_mpi.py來使用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 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 quitpython?dry_for_TI.py?box_check.pdb?wat1.pdb?omi_rbd.pdb,omi_ace2.pdb

同時(shí)生成一個(gè)小的水盒，用于跑蛋白部分的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 quitpython?dry_for_TI.py?ligands_recharge.pdb?rbd_wat.pdb?test_rbd.pdb

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

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

?

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

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.rst7setBox?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 _EOFparmed?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文件轉(zhuǎn)成真正的ti文件放進(jìn)文件夾，同時(shí)使用slurm提交最小化，加熱和運(yùn)行步驟（5ns）。

總結(jié)

以上是生活随笔為你收集整理的利用Amber热力学积分计算相对自由能变化的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。