:從分子對接到虛擬篩選與自由能計算指南)
看到“schrodinger 薛定諤”這個關鍵詞很多人的第一反應是那只既死又活的貓。但我作為一個常年跑計算模擬的人腦子里閃過的其實是另一個東西Schr?dinger Suite業(yè)內直接叫“薛定諤”一套做分子模擬、藥物設計、材料科學的高通量計算平臺。薛定諤方程是量子力學的基石而 Schr?dinger 這家公司把這套物理模型做成了工程化軟件讓人可以直接在圖形界面里做分子對接、虛擬篩選、自由能微擾計算算得上是藥物發(fā)現(xiàn)領域繞不開的“重型裝備”。這篇文章不聊量子力學科普就圍繞我在實際項目里用 Schr?dinger 套件做分子對接和虛擬篩選的經驗展開。我會從最初的軟件選型邏輯講起一步一步把結構準備、Glide 對接參數(shù)、MM-GBSA 重打分、FEP 邊界這些核心環(huán)節(jié)拆開再把實操中踩過的坑、總結的排查方法一并放出來。適合正在學計算化學、生物信息學或者剛拿到薛定諤 license 準備跑第一個真實課題的人參考。1. 先說清楚此薛定諤非彼薛定諤1.1 從量子力學的薛定諤到計算平臺Schr?dinger 這個名字源于物理學家 Erwin Schr?dinger他提出的薛定諤方程描述了微觀粒子的波函數(shù)演化?,F(xiàn)代分子模擬的底層邏輯本質上還是解這個方程——不過不是解析解而是通過各種近似方法分子力學、半經驗方法、從頭算、DFT在計算機上逼近真實體系的能量和運動狀態(tài)。Schr?dinger 這家公司成立于1990年核心就是用物理模型去預測分子的性質和行為把學術圈里散落的算法工具包裝成一套工業(yè)級工作流。很多人第一次接觸“薛定諤”是在藥企實習或者計算化學課程里打開 Maestro 界面看到一堆模塊菜單。整套套件里最常見的幾個模塊是Glide分子對接、Prime蛋白結構預測與 MM-GBSA 重打分、LigPrep配體準備、Desmond分子動力學模擬、FEP自由能微擾、Canvas化學信息學與構效關系分析。它們不是獨立的小工具而是共享同一套力場參數(shù)、同一套文件格式、同一個圖形界面這意味著從受體準備到結果分析整個流程都是串起來的不太需要在不同軟件之間玩“文件搬家”。1.2 誰在什么場景下會用到它我見過三類人會真正用到這套軟件。第一類是藥物研發(fā)團隊的計算化學家他們拿它做靶點的虛擬篩選、先導化合物優(yōu)化配合濕實驗驗證第二類是結構生物學出身的研究者手里有晶體結構或冷凍電鏡結構想快速看這個靶點能裝下什么形狀的分子第三類是跟計算沾邊的學生和跨界研究者比如做農藥設計、材料篩選、甚至酶催化機理的人也會借用對接和動力學模擬來補充實驗解釋。這些場景有一個共同點核心需求都是“預測”。結晶一個蛋白可能花掉幾個月合成一個化合物又需要時間和經費而計算可以通過打分函數(shù)先篩一遍把候選范圍從百萬級壓到幾十個。所以 Schr?dinger 解決的不是“算出真相”而是“用較少資源找到最可能正確的方向”。這點非常重要后面所有參數(shù)選擇、流程設計都是圍繞它來的。2. 用 Schr?dinger 搭一條分子對接流水線2.1 結構準備Protein Prep Wizard 和 LigPrep別跳過這一步我第一次用薛定諤跑對接的時候以為拿一個 PDB 結構直接 Glide 就能出結果結果對出來的打分爛得離譜。后來才明白PDB 文件里存的是實驗解析的原子坐標但這個坐標系里其實缺了很多化學上必須的信息氫原子基本都沒有組氨酸的質子化狀態(tài)不明確Asp、Glu 等殘基該帶幾個質子也沒定活性口袋里可能還混著結晶緩沖液分子、金屬離子、二聚體界面上的水。Protein Prep Wizard蛋白準備向導做的就是這些“補齊”工作。實際操作里我會依次做三件事。第一給蛋白加氫并分配正確的質子化狀態(tài)第二用 Prime 補缺失的殘基和側鏈同時處理掉不合理的原子位置沖突第三做一輪能量最小化讓整個結構在力場下達到一個穩(wěn)定構象。這里要特別提醒最小化的原子位移限制不要設太大默認 0.3 ? 就夠了否則會把活性位點壓變形后面的對接全廢。配體端的準備用 LigPrep。它會生成配體的 3D 構象計算特定 pH 下的質子化狀態(tài)還能產生合理的立體異構體和環(huán)構象。這一步很多人不重視直接拿 2D 結構的 smiles 去對接結果發(fā)現(xiàn) Glide 報錯或者產出一堆高能構象。常見做法是設定 pH 7.0 ± 2.0讓 Epik 模塊按生理條件計算所有可電離基團的質子化狀態(tài)再對每個輸入結構生成最多 32 個低能 3D 構象。聽起來簡單但這一步直接決定后續(xù)對接能不能找到“對的姿勢”。提示結構準備是整個流程的地基。模板結構里如果帶著共晶配體Prep Wizard 默認會保留它不要盲目刪掉。這個配體往往是網格生成的定位點也是判斷對接結果是否合理的錨定參照。2.2 Glide 分子對接參數(shù)怎么選Glide 是 Schr?dinger 的招牌對接工具它把配體放進受體結合口袋搜索所有可能的結合構象然后用 GlideScore 打分。具體用哪種精度取決于你要解決的問題。Glide 提供了三檔模式HTVS高通量虛擬篩選、SP標準精度、XP超高精度。HTVS 最粗糙也最快適合動輒百萬級的超大庫初篩SP 是日常最常用的精度兼顧速度和可靠性XP 會對配體-受體之間的形狀互補、疏水作用做更精細的懲罰和獎勵但極其耗時適合對幾十個命中做深度分析。網格生成這一步最容易被忽略。Glide 需要先定義一個受體網格Receptor Grid網格的中心通常放在共晶配體的質心上這樣能保證結合口袋被完整覆蓋。網格盒子的大小不需要貪大默認的 15 ? × 15 ? × 15 ? 一般夠用盒子里還可以指定某些殘基為柔性殘基讓側鏈在對接時微調但每加一個柔性殘基計算量會指數(shù)增加所以一般只對活性位點里最關鍵的幾個殘基開啟。如果做一個小規(guī)模的篩選我的習慣是先用 SP 跑一遍把打分前 20% 的配體拿出來再用 XP 重新對接這兩輪的結果綜合起來決定下一步實驗。很多人上來直接跑 XP一個配體跑十幾分鐘換來的是被噪聲淹沒的分數(shù)純屬浪費時間。2.3 打分函數(shù)怎么選、怎么看對接完成后你會在 Maestro 的結果表里看到一列 GlideScore。它表示的是配體結合到受體上時的一種“經驗性評分”不是真實的結合自由能但數(shù)值越負說明預測的結合越強。它由好幾項加在一起靜電相互作用、范德華力、氫鍵、疏水接觸、溶劑化效應、以及配體內部的應變能。理解“打分函數(shù)不是能量”這一點很重要。真實結合自由能要考慮熵、去溶劑化、蛋白構象變化這些在對接打分里要么是近似項要么根本沒有。所以分數(shù)差 0.3 分根本不算差距我自己一般以 1 分以上作為篩選閾值。另一點是分數(shù)高更負不代表活性好因為打分函數(shù)可能被某些“表面友好”的分子誤導比如分子太大、埋在口袋里強行填滿空隙導致范德華項虛高。我見過不少人直接按 GlideScore 排序取前 50 個去做活性測試結果命中率慘淡。我的做法是先按分數(shù)看排名再逐個目視檢查結合姿勢——看配體有沒有完全離開口袋、有沒有嚴重的原子碰撞、極性基團有沒有和水或缺電子區(qū)域形成合理相互作用。視覺檢查能過濾掉 30% 到 50% 的“分數(shù)好看但姿勢離譜”的候選。3. 進階從虛擬篩選到結合自由能計算3.1 虛擬篩選流程設計與富集率如果你面對的化合物庫有幾萬甚至上百萬個分子那就不是一個個跑 Glide 的事了而是一個系統(tǒng)流程設計問題。我習慣的虛擬篩選管線分四步。第一步是預過濾。所有配體先過一遍物理化學性質過濾器比如分子量、AlogP、氫鍵供體/受體數(shù)量、可旋轉鍵數(shù)參考類藥五規(guī)則。再跑一遍 ADMET 預測把有明顯毒性風險或透膜性差的分子提前踢掉。這一步能把庫縮到原來的三分之一甚至更少。第二步是構象和異構體準備。用 LigPrep 批量生成 3D 結構。注意控制立體異構體數(shù)量沒有手性中心的分子不要生成無謂的異構體否則后面對接時間直接翻倍。第三步是分級對接。所有化合物先跑 HTVS按分數(shù)保留前 20%30%剩下的跑 SPSP 結果里取前 1000 到 2000 個跑 XP。每一級都沒有必要把上一級的全部結果都拿過來篩的就是“低分基本不值得細看”的思路。第四步是重打分和聚類。用 Prime MM-GBSA 對 XP 命中的前 200 個分子重新打分再對化學結構聚類選每一類里分數(shù)最高的一兩個。這樣最終留下的幾十個化合物在結構和理化性質上都有代表性和差異避免 50 個命中全是同一個母核的不同尾巴。富集率是評價篩選流程優(yōu)劣的關鍵指標。操作上可以準備一組已知活性化合物和一組誘餌分子混入候選庫中跑一遍完整流程看看活性化合物是否比誘餌排在更前面。如果活性化合物大部分排在庫的頭部說明流程的富集能力靠譜如果活性化合物排名跟隨機差不多那即使表面分數(shù)好看也說明你的受體結構或準備步驟有問題。3.2 Prime MM-GBSA 的批量重打分GlideScore 是快速篩選的機槍但到了決賽圈你會需要精度更高一點的評價。Prime MM-GBSA 計算的是配體結合前后的能量差綜合了分子力學能量MM、連續(xù)溶劑化模型GBSA里的極性項和非極性項。它的核心優(yōu)勢在于部分考慮了溶劑化效應這是對接打分里處理得比較粗的部分。實際操作不復雜在 Glide 的結果列表里選中一批配體右鍵選擇 Prime MM-GBSA設置力場用 OPLS4、溶劑模型用 VSGB其他保持默認然后提交任務。它會自動對每個配體做一個小規(guī)模的構象采樣并給出一個 delta G bind 的估計值。通??梢园?Glide 的排名和 MM-GBSA 的排名做一個交叉驗證兩者都靠前的分子優(yōu)先做實驗。但這個手段也有明顯的邊界。它仍是一個“端態(tài)”計算忽略了配體結合過程中的熵效應與通路上的中間態(tài)因此它輸出的數(shù)值適合做排序不適合當真實自由能解讀。比如對比 A、B 兩個類似物說“B 比 A 預估強 5 kcal/mol”是合理的但說“A 的結合自由能就是 -48.6 kcal/mol”就沒有物理意義因為絕對數(shù)值受力場參數(shù)和模型簡化影響太大。在我自己的流程里MM-GBSA 永遠是“相對比較”工具絕不單獨用它拍板一個化合物能不能進合成清單。3.3 FEP 的思路和適用邊界如果項目推進到先導化合物優(yōu)化階段常見的問題是在某個核心骨架上我想把這個苯環(huán)換成吡啶把甲基換成乙基或者在這里加一個氟原子活性會變好還是變壞這時 MM-GBSA 的精度已經不夠了就需要引入自由能微擾FEP。FEP 的核心思路是用一系列中間態(tài)連接兩個結構高度相似的配體計算從一個配體“演化”到另一個配體過程中的自由能差值。因為兩個分子結構接近很多誤差會在求差的過程中相互抵消所以它預測的相對結合自由能精度可以做得非常高誤差通常在 1 kcal/mol 以內遠好于 MM-GBSA也更適合指導化學家往哪個方向修飾。但 FEP 不是拿來隨便跑的。它的第一個硬性要求是配體之間結構足夠相似最好是只在局部官能團上有差異如果你拿一個完全不同的骨架去做 FEP中間態(tài)很難收斂結果就失去了意義。第二個要求是計算資源充足每個配體需要搭一張合理的圖哪些配體之間建微擾邊每個邊都要跑分子動力學模擬動輒需要 GPU 集群跑幾天。第三個要求是你得有一個可靠的共晶結構或對接結構作為起點起點蛋白質構象如果不對后面一切都白搭。我自己項目里的分工是這樣的初篩用 Glide中篩用 MM-GBSA最后的十來個類似物決定具體合成順序時才上 FEP。這個三級遞進既保證了效率也把最貴的計算留給了最重要的決策點不會出現(xiàn)“GPU 集群跑了一禮拜結果發(fā)現(xiàn)模擬的體系根本不是活性構象”這種慘劇。4. 實操中我踩過的坑4.1 別拿原始 PDB 直接對接我剛開始做激酶項目時直接從 PDB 下載了一個分辨率不錯的晶體結構沒有用 Prep Wizard 處理只是刪掉了水分子和其他鏈就生成了網格開始對接。結果篩出來一個打分極高的化合物合成出來卻完全沒活性。事后回頭檢查才發(fā)現(xiàn)原始結構里 DFG 基序附近有一個關鍵殘基的側鏈密度不完整PDB 里存的坐標是扭曲的Prep Wizard 能自動檢測并修正而我跳過了這一步。正確操作是每次下載 PDB 后先看序列信息和配體信息再用 Protein Prep Wizard 按默認流程走一遍重點觀察 Log 窗口里有沒有“missing residues”或“alternate conformations”提示。如果有缺失長度超過五六個殘基的 loop看它離活性位點遠不遠遠的話直接保持缺口不要強行建模近的話要考慮換一個更完整的結構因為模型補出來的 loop 不確定性很大對接結果可信度低。4.2 質子化狀態(tài)和異構體的問題配體端最容易翻車的是異構體數(shù)量失控。有一次我想篩一個含有兩個手性中心的母核LigPrep 跑完提醒我生成了 4 個立體異構體各 32 個構象總共 128 個文件對接時間直接翻了 16 倍。后來我學會了在有明確實驗數(shù)據(jù)的前提下用 chiral 標記指定只有目標手性構型參與計算或者用 LigPrep 的“retain specified stereoisomers”功能把不需要的異構體刪掉。另一件常被忽略的事是水分子的處理。有些活性位點里會有一個水分子它同時在蛋白和配體之間搭橋就好像兩個陌生人之間的傳話人處于關鍵的位置。準備受體時如果無條件刪除所有水分子某些配體可能就少了一個重要的極性相互作用打分被低估。我的做法是保留所有太陽下看起來協(xié)調的水分子與蛋白極性殘基距離在 3.0 ? 左右且不與非極性殘基沖突的生成網格時再把它們當作受體的一部分參與打分。4.3 資源消耗與并行調優(yōu)Schr?dinger 的任務調度邏輯跟普通程序不太一樣很多新手會誤以為把所有 CPU 核心都填滿就能加速。實際上單個 Glide 對接任務本身對多核利用有限提速的關鍵是“同時跑多個任務”。我習慣在 Linux 服務器上把配體庫拆成幾個子文件每個子文件單獨提交一個 Glide 任務利用批處理腳本把它們并行調度起來。Desmond 和 FEP 則剛好相反它們強烈依賴 GPU 加速沒有 GPU 的話跑一個 100 ns 的分子動力學模擬可能要等一周而用單張中端 GPU 能縮短到一天以內。所以如果你準備長期跑薛定諤的自由能計算別省 GPU如果只是跑對接和虛擬篩選CPU 就夠用優(yōu)先保證內存不要爆掉就行。許可證連接也是一個隱蔽的坑。Schr?dinger 的 license 是網絡浮動式的如果服務器網絡不穩(wěn)定或者 license 服務沒啟動啟動 Maestro 或命令行工具時會一直卡在“could not checkout license”的報錯。排查時先看 license 服務器日志再確認本機環(huán)境變量 SCHROD_LICENSE_FILE 是否指向正確的端口和地址通常都能解決。5. 常見問題速查現(xiàn)象可能原因處理方法Prep Wizard 報錯缺少大量殘基PDB 文件本身分辨率差或截斷換用分辨率更高的結構或刪除無序區(qū)域盡力補全關鍵殘基Glide 網格生成失敗共晶配體坐標不完整或原子類型沖突檢查配體是否含金屬離子或非標準殘基編號必要時重新用 LigPrep 處理LigPrep 產出的構象太多立體異構體數(shù)量未限制指定手性構型控制每個結構的構象生成上限對接分數(shù)很好但結合姿勢奇怪配體沒進入口袋或原子碰撞嚴重用 Maestro 檢查受力分布關掉不合理的柔性殘基設置MM-GBSA 結果不收斂配體電荷設置不合理或結構構象過少重新檢查配體質子化狀態(tài)生成更多的初始構象FEP 自由能結果方差很大中間態(tài)不夠密或者模擬時間不夠增加兩個相鄰配體間的中間態(tài)數(shù)量延長平衡時長檢查是否選了結構差異過大的配體對還有一個值得單獨說的技巧批量任務報錯時優(yōu)先看 job 目錄下的 write.log 文件。Schr?dinger 的圖形界面經常把錯誤吞掉只顯示一個紅色的 “Job failed”但具體的失敗原因一定寫在日志里。學會從日志的第一條 error 開始處理比反復重跑要有效得多。實際項目里我第一次排查一個批量對接任務失敗的問題就是因為某個配體包含未被力場定義的鹵素原子類型日志里一行 “cannot type atom” 瞬間讓我定位到問題換了參數(shù)版本后整個隊列就順了。另外建議定期備份 Maestro 中構建好的自定義網格和結果視圖。薛定諤的項目文件結構里除了 .mae 主文件還有配套的 .maegz、.log、.grd 等文件。只拷貝主文件回家換個機器繼續(xù)看結果時經常發(fā)現(xiàn)結構缺失就是因為沒有把同一目錄下的附屬文件一起打包。把這些文件統(tǒng)一放進一個項目文件夾里算是一個做了就會感謝自己的好習慣。這套軟件不便宜學習曲線也不平緩但只要理順了“準備結構—生成網格—分級對接—重打分—自由能計算”這條主線它就能成為決策效率的強杠桿。我個人這幾年最深的體會是薛定諤的價值不在于幫你跑出更高的分數(shù)而在于逼著你想清楚每一步到底在算什么、誤差在哪里、結論能支撐多大范圍的判斷。想明白這層軟件選哪家、參數(shù)怎么調反而都是次要問題了。