超大規模藥物虛擬篩選的實現與應用

張寶花，李輝，劉倩，高美娜，黃荷，趙毅，于坤千+，金鐘

1.中國科學院計算機網絡信息中心，北京100190

2.中國科學院上海藥物研究所，上海201203

3.中國科學院大學，北京100049

惡性傳染病的間歇性暴發對人類社會產生了深遠而持久的影響。人類與病毒的較量，是科研與時間的賽跑，而快速藥物研發則是較量之關鍵。采用基于受體結構進行分子對接的高通量藥物虛擬篩選（high throughput virtual screening，HTVS）方法，有助于從海量小分子中篩選出與受體結構最強結合的部分，是先導化合物發現的重要步驟[1-2]。在選定了合適的分子對接軟件并進行了合理參數設置后，HTVS計算效率的提高主要依靠虛擬篩選規模和虛擬篩選計算的有效執行。

據化學家估計，滿足Lipinski類藥五規則的穩定小分子化合物數量達到1060左右[3]。研究表明，大規模篩選的對接得分隨虛擬篩選規模增加而提高[4-5]。增加虛擬篩選規模可從兩個方面來改善虛擬篩選的質量：通過識別具有更緊密結合親和力的配體，可降低劑量和減少脫靶效應；發現藥代動力學更有利的、細胞毒性較小的化合物[6]。增加虛擬篩選規模則主要有兩種途徑：一方面，考慮到化合物的可合成性和易購買性，可選擇以ZINC（＞5 億）和REAL（＞14 億）為代表的大型化合物庫[7]；另一方面，人們可通過碎片組合虛擬篩選庫、深度學習等方法生成海量的配體化合物[8-9]。

在可接受時間內完成海量化合物分子的虛擬篩選對高性能計算技術提出了巨大挑戰。幸運的是，目前的超級計算機已接近E級（Exascale，每秒可進行百億億次浮點運算）計算能力。美歐日均制定了E級計算計劃，我國也在研制下一代異構高性能計算系統，這均為超大規模虛擬篩選提供了可能。計算功能強大的超級計算機在藥物設計中的應用給藥物發現帶來了新機遇，使藥物分子設計在速度和成功率上均發生了突飛猛進的提高。然而，如何緊跟硬件的步伐，將超級計算機的計算能力發揮到最佳，還需要選擇合適的方法進行有效的大規模藥物虛擬篩選。

科研人員在如何將分子對接軟件應用于高通量藥物虛擬篩選方面開展了大量工作。例如使用AutoDock Vina 軟件來進行高通量虛擬篩選，可提前將配體小分子進行一定數量切割后存入不同的文件目錄。在計算資源充裕的條件下，不同的文件目錄可以同時獨立執行計算來最優化計算時間[10]。rDock與rxDock 軟件則通過將含有多個配體文件的sdf 文件，利用sdsplit 進行按需劃分，對每份文件在當前節點獨立發起分子對接操作或通過隊列調度系統調度到其他計算節點執行計算[11]。Vina MPI使用MPI包裝器，可支持同時啟動數千個并行AutoDock Vina 可執行文件[12]。POAP（parallelized Open Babel&AutoDock suite pipeline）將Open Babel 和AutoDock 軟件包連用，使用GNU Parallel將配體準備、受體準備、對接任務等并行化，通過Shell 命令行方式使用[13]。以上幾種方式，僅針對部分適合藥物發現的有機化合物（106～107個分子）[14]進行中小規模的藥物虛擬篩選計算，對多任務進行簡單分割，易引起負載不均衡，缺乏大規模藥物篩選的應對策略。李祖華等人利用網格計算技術構建了用于大規模藥物虛擬篩選的計算環境[15]，可有效整合散布于網絡上的各種異構計算資源，并能夠快速吸納新的網格結點和篩選工具。在整體架構上使用了支持大規模并發的主從模式，在網格監控上使用了支持多種異構軟硬件環境的適配器模式，滿足了真實環境大規模藥物虛擬篩選的要求，但存在數據轉移安全性問題，網絡通訊能力強弱對性能也有很大影響。Peng等人在“天河二號”上研發的D3DOCKxb軟件[16]則通過設計一種新型可擴展并行算法，利用CPU+MIC協同計算，采用分層、分段任務處理策略，分組休眠緩解IO壓力，“橫向并行、縱向分段”緩解通訊壓力，在8 000 節點上可獲得超過80%的加速，一天內篩選老藥化合物4 200 萬個。Gorgulla 等人在Nature上發表的超大規模藥物發現平臺[4]，其中的VFVS（virtual flow for virtual screening）使用Shell 腳本開發，文件列表形式進行多層級任務分發，數據壓縮形式遷移減少IO 問題，使用8 000 個CPU 核心，4 周時間篩選化合物約13 億。以上科研團隊在大規模藥物虛擬篩選的有效執行方面均取得了較好進展。但隨著異構超級計算機大量采用GPU 卡等加速部件，針對特定服務器架構設計的mD3DOCKxb 軟件和僅支持CPU 計算的VFVS 軟件已無法有效利用加速部件提升計算速度。此外，在超大規模虛擬篩選時的動態負載均衡和容錯管理方面亦存在不少優化空間。

本文提出了一種在超級計算機上進行大規模藥物虛擬篩選新的實現方式，即利用中央數據庫進行任務的多層動態分發。此外，對大規模虛擬篩選中的負載均衡、通訊、IO及容錯管理等問題也開展了研究。通過以上工作，有望利用高性能計算系統進行更大規模更快速的虛擬篩選，為應對暴發性惡性傳染病候選藥物的快速發現提供強大計算保障，大大地加速藥物研發。

文章將按照以下組織進行撰寫：系統設計部分將介紹系統的總體架構及在動態負載均衡、容錯管理、通訊和IO優化等方面的設計方案；系統實現部分將介紹系統實現中相關的技術方案，包括系統的運行流程和操作方式等；系統測試部分通過各項測試展示系統的實現效果；系統應用部分將介紹系統在新冠病毒藥物研發中的應用成效。最后部分為總結和展望。

1 系統設計

1.1 系統整體設計架構

藥物虛擬篩選是海量數據影響下的超高通量計算過程，由大批量松耦合的子任務構成。在超級計算機上進行此類計算模擬，需有效地利用各空閑資源，高效處理數十億化合物，確保高通量計算執行的速度和質量。為此，本文提出以中央數據庫進行集中任務分發的方案。考慮到單個數據庫進行任務分發將達到負載上限，如集中通訊量過大必然導致通訊速度減慢甚至異常。從效率和節點穩定性方面考慮須采用多數據庫分擔負載。為方便管理多個數據庫并及時補充待計算的任務，本文提出基于多層級數據庫的并行分發處理策略，設計了“樹”形的多層級任務分發框架。

如圖1所示，根服務節點控制整個計算任務的中央數據庫，計算任務從這里開始分配。T個“樹”的枝服務節點啟動枝數據庫，與根數據庫通訊獲取大任務并分割后，由圖2所示的任務執行單元來執行計算。

圖2 任務執行單元處理流程Fig.2 Processing flow of task execution unit

每個任務執行單元采取兩個層級的主從模式設計，第一層為枝數據庫與計算節點（圖中WG表示）的主從執行模式，數據庫負責管理和分發作業，WG 與數據庫通訊獲取和更新任務。在每個WG 里含第二層主從模式，主進程負責抓取和預處理任務并啟動多個子進程執行批量分子計算操作。子進程將計算任務全部執行完畢后，由主進程統計計算結果，對本組結果進行排序并將少量輸出結果寫入上級數據庫。之后去數據庫中抓取新任務，不斷輪詢“抓取任務-啟動子進程計算-結果處理”操作，直至數據庫中任務全部執行完畢，WG 退出，本數據庫作業執行完畢。然后由上級枝服務節點與根數據庫通訊獲取下一批次計算任務。

1.2 動態負載均衡策略

藥物虛擬篩選屬于分布式計算[17]，可把需要利用巨大的計算能力才能解決的藥物篩選問題進行“水平切分”，分布到不同的節點分別處理得到結果，并綜合起來得到最終結果。如何減少任務獲取、數據更新時的通訊負擔并保證計算任務在各種狀態下的負載均衡，是超大規模藥物虛擬篩選的主要痛點。在本文中，得益于多層級中央數據庫和“樹”形模式的任務分發框架，下級節點通過主動鏈接上級數據庫進行任務獲取和執行，計算過程中各節點都處于任務滿載狀態，并可根據節點服務能力動態調整負載。運行快的節點將從上級數據庫獲取更多的任務計算。反之，計算較慢的節點則會減緩與上級數據庫的通訊頻次。故障節點則不再與上級數據庫通訊獲取任務，由此可將慢節點和故障節點的影響降到最低。

1.3 高并發帶來的IO和通訊問題優化策略

超大規模虛擬篩選涉及到大量的數據通訊、文件讀寫等操作，對超級計算機的網絡通訊和IO 提出很大挑戰。

直接存放和處理數十億級別的分子會導致文件系統卡頓甚至崩潰。本文采取“分而治之”的處理策略，將海量分子文件進行多層壓縮打包以減少存儲壓力。如圖3 所示，根數據庫存放多個頂級壓縮文件。執行計算時枝服務節點與根數據庫通訊，獲取頂級壓縮包，解壓得到多個二級壓縮包（內含多個小分子文件），并作為圖2 所示任務執行單元的處理對象。另外，從計算策略上進行優化，通過長時間間隔讀寫來減少IO 頻次。在計算節點上，每次分子對接計算執行完后，隨即進行能量提取操作，待同一批次全部分子執行完后，由主進程執行能量排序并根據用戶選擇挑選排名靠前的結果記錄到上級數據庫中。通過少量大文件長時間間隔讀寫減輕了系統IO壓力，使硬件系統在處理海量計算任務時流暢不卡頓。

圖3 海量小分子文件的多層級壓縮處理策略Fig.3 Multi-level compression strategy for massive small molecular files

通過文件壓縮和大文件遷移的方式，在緩解系統IO 壓力的情況下，還間接促進了通訊問題的改善。計算節點每次從數據庫獲取多個配體分子組成的壓縮包，進行二級分發后執行計算，如此可有效減少與上級數據庫的通訊次數。分子對接計算任務則完全在計算節點本地執行，計算節點間無任何通訊。

1.4 容錯管理策略

高通量任務運行時，不可避免地會出現節點異常、程序計算錯誤等問題。如何保障作業正確、有效、順利地執行，是超大規模藥物虛擬篩選的基本需求。

由前所述，通過采取中央數據庫動態管理和分發作業，相較于固定文件列表分發方式，中央數據庫服務在靈活性和易用性方面更有優勢，例如支持在計算過程中動態靈活增加數據、刪除數據、修改數據，容易統計任務執行狀態等，在訪問時可以對單條數據加鎖，且加鎖粒度更細，可將故障節點和慢節點的影響降到最小，這就要求中央數據庫服務必須持續運行。

并行分發系統配備的監控程序可定時檢查根數據庫服務和枝數據庫服務。如圖4，若檢測到服務節點異常掉線，則自動重啟服務，保證下級數據庫能夠鏈接獲取任務。同時，在每個節點中，對計算節點錯誤，如節點空掛無進程、加速卡故障等錯誤進行定時掃描，發現異常則刪除進程，釋放節點。對于失敗的計算作業，通過修改作業狀態可重新回到待計算任務數據庫中由其他正常節點獲取計算。

圖4 計算錯誤處置策略Fig.4 Calculation error handling strategy

2 系統實現

本工作采用Python語言實現了簡單易用的藥物虛擬篩選并行分發系統。MongoDB數據庫作為中央數據庫，是一款為現代Web 應用程序和互聯網基礎設施而設計的新型開源數據庫管理系統。作為一種典型的分布式NoSQL數據庫，它并發性能優異、支持水平擴展、模式自由靈活，且在大文件傳輸、故障自動轉移等方面具有優良的特點。藥物虛擬篩選涉及到的任務比較松散，非常適合使用MongoDB 數據庫來進行任務管理。

圖5展示了并行分發系統的啟動流程，用戶發起計算請求并通過隊列調度系統獲取到計算資源后，由主程序啟動一個根數據庫和T個枝數據庫。根服務節點啟動枝服務節點與根數據庫通訊獲取任務，并對任務分割后存于枝數據庫中。啟動計算節點，并開啟主進程與枝數據庫連接獲取任務，通過Python進程池調用節點上其他資源進行分子對接計算。在每個計算節點上，還配有一個心跳程序，每隔10 min連接枝數據庫對當前任務執行時間和狀態進行更新。主程序還會啟動監控程序來監控數據庫和計算節點的狀態，根據異常狀態來進行容錯。

圖5 并行分發系統啟動流程Fig.5 Parallel distribution system startup process

在進入虛擬篩選前，需對受體和配體分子進行處理，如受體補全殘基、格式轉換、加氫及加電荷等操作，并根據受體結合口袋位置和尺寸設置分子對接的配置文件。經過預處理的文件即可使用藥物虛擬篩選并行分發系統執行計算。圖6 展示了并行分發系統在超級計算機各節點上的映射關系。為管理方便，每個枝數據庫設置M個子庫負責管理和分發分組作業。在虛擬篩選執行時，為每個分組啟動N個計算節點，每個節點對應圖2中一個WG。在節點上，進程1作為主進程，負責執行數據抓取、后處理等任務，而其他進程則負責調用CPU 或競爭調用加速器執行具體的分子對接操作。

圖6 并行分發系統在超級計算機節點上的映射關系Fig.6 Mapping of parallel distribution system on supercomputer nodes

本文并行分發系統屏蔽了復雜處理細節，用戶通過簡單設定配置文件相關參數即可開啟大規模并發計算模擬。通過編輯配置文件可使用服務器上已安裝的軟件執行分子對接計算，如選擇AutoDock GPU 程序利用GPU 加速計算或AutoDock Vina 程序利用CPU并行加速，解除對于軟件的高集成綁定，實現分子對接軟件與并行分發系統的松耦合。對于執行計算的分組數量、執行計算的進程等數量均可通過參數配置靈活調整以增強易用性。

3 系統測試

針對如何在超級計算系統上有效進行超大規模藥物虛擬篩選，本文設計了基于中央數據庫的“樹”形多層級任務并行分發系統。在動態負載均衡、高容錯技術、優化IO與通訊問題等方面進行的研究，為快速高效穩定地執行大規模藥物虛擬篩選提供了有力保證。

系統硬件配置：每節點配置1 顆國產X86 處理器，32核心，8*DDR4內存通道、128*PCI-E 3.0通道，4個國產加速器，16 GB HBM2顯存，帶寬1 TB/s。

計算參數設置：藥物虛擬篩選并行分發系統需要結合具體的分子對接程序執行計算。研究測試表明，雖開源軟件相較于商業軟件的算法不同，但構象排序和打分結果差異不大，并不存在一個對接軟件在各方面都優于其他軟件[18]。考慮到對接軟件的準確性和計算速度，本文選取了AutoDock GPU（OpenCL加速版本）[19]程序執行對接計算。

經測試，AutoDock GPU 計算并非全部在加速器上執行，加速卡在較少進程時出現了負載不滿的現象。節點內分子對接計算的進程數量會對計算效率產生影響，測試（圖7）發現每個節點啟動30個以上進程競爭使用加速器計算的執行時間最少，硬件能夠得到充分的利用。因此，在之后的計算中，設置每個WG中使用的子進程數量為30。

圖7 節點內計算進程數量對性能的影響Fig.7 Impact of the number of computing processes on performance intra node

計算受體采用由饒子和院士團隊公布的新冠病毒Mpro 蛋白及對接口袋位置[20]，使用AutoDock Tools 對受體進行加氫加電荷和格式轉換操作，對接口袋大小設置為22.5 ?，受體設置為剛性，使用AutoDock Tools 準備好.gpf 文件，進而使用autogrid程序生成需要的各個.map文件和.fld文件[21]。配體取自公開數據庫REAL和ZINC，使用Openbabel軟件進行了加氫加電荷和格式轉換處理，配體設置為柔性。

3.1 弱擴展性測試

對饒子和院士團隊報道的病毒RdRp 蛋白[22]與REAL 數據庫分子執行篩選計算，選取測試軟件為AutoDock GPU 程序，每個計算節點啟動30 個進程+4 加速器來執行計算，對接策略采取受體為剛性，配體為柔性的處理方案，每次對接產生10 個最優構象。以50節點對100萬分子執行半柔性對接進行測試，采用傳統劃分文件目錄后，順序調用加速器執行分子對接計算的處理方法，計算時間受任務劃分粒度、計算節點服務能力、小文件頻繁IO 等因素的影響，計算時間大于5 h。超過500節點的大規模模擬、計算不均衡等因素導致的影響則更為嚴重，甚至多次出現節點故障導致計算失效的情況。使用本文分發系統，計算時間可縮短到1.1 h。以50節點對100萬分子執行半柔性對接為基準測試計算效率。通過網絡傳輸效率的測試，本文采取了一個枝服務數據庫向500 個計算節點同時分發任務的方案。例如，在2 000節點計算時，采用了1個根數據庫和4個枝數據庫。由表1 所示，隨著并發規模的增大，計算效率在2 000節點時，效率超過95%。

表1 并行分發系統弱擴展性測試Table 1 Weak scalability test of parallel distribution system

對于更大并發規模的計算，可通過增加枝數據庫的數量來調度更多的計算節點。由于計算節點只與上級數據庫通訊，而彼此不進行通訊，每個分支計算互不影響，可實現計算規模的不斷擴展。

3.2 與VFVS計算性能對比測試

VFVS 為美國哈佛大學最新發表的超大規模快速藥物虛擬篩選平臺。它使用Shell 腳本開發，采用文件列表進行任務分發，集成了包括QuickVina 和AutoDock Vina在內的多個CPU計算軟件，暫不支持GPU加速。利用該平臺，科研人員針對KEAPI受體，利用8 000CPU 核心花費約4 周時間，篩選了超13 億分子（包括3.3 億ZINC 數據庫分子及10 億REAL 數據庫分子），文中稱可利用Google 云平臺16 萬CPU 核心于15 h 篩選10 億化合物，是目前已知報道的篩選規模最大、速度最快的藥物虛擬篩選平臺。

本文針對新冠病毒Mpro 受體（參數設置同上），每個節點選擇30個CPU核心，使用QuickVina2軟件[23]對VFVS 和本文并發系統進行對比測試。首先測試單節點計算，本文選取了300 個分子，根據VFVS 軟件要求，分為30組（每組10個分子）寫入中央文件列表，使用30個進程執行QuickVina2軟件計算，每個進程串行執行。本文對300 個分子通過Python 進程池調度節點內30 個進程執行上述計算。測試發現：對于300 個分子結構全部相同的情況，VFVS 測試時間為513 s，本文程序為502 s，性能接近；若300 個分子結構不同，VFVS的計算時間為473 s，本文為295 s。同理，對于600和1 123個不同結構的分子，本文并發系統體現出了良好的速度優勢。

根據SlogP、MW、TPSA、HBA、HBD、RotB 等參數從REAL 數據庫中選擇了166 萬分子，分別使用VFVS 和本文并發程序進行測試，測試條件同上，測試規模為50 節點，VFVS 全部節點執行計算完成需19.5 h，本文為10.2 h。

分析計算時間出現差異的原因，VFVS在工作任務較少時，會發生進程的CPU利用率不高的狀況，由負載不均衡造成性能少許損耗。本文由于使用了進程池動態調度，CPU基本處于滿載狀態，有效克服了負載不均衡問題。

4 并行分發系統在抗新冠病毒藥物研發中的應用

新型冠狀病毒肺炎為感染新型冠狀病毒所致，該病毒在入侵細胞后首先利用細胞機器合成復制酶多肽，由主蛋白酶Mpro剪切后形成功能蛋白，進一步組裝成以RNA聚合酶RdRp為核心的復制轉錄機器，啟動自身遺傳物質的大量復制，干擾這兩種蛋白質的功能有望阻斷或抑制新冠病毒的復制過程。根據饒子和院士團隊解析Mpro和RdRp家族同源結構信息，本文提出了兩種不同活性位點的篩選方案，并在國產超級計算機上開展了對多個小分子數據庫的虛擬篩選。

采用本文的并行分發系統，結合AutoDock GPU程序（計算參數設置同上）篩選了超過20億個小分子化合物，化合物來自ZINC、REAL公開數據庫和自有數據庫等。如圖8，篩選出的活性分子有瑞德西韋類似物（圖中紫色部分對應的分子），還發現了若干新活性分子，目前已進入實驗活性測試階段（圖8 上圖為新冠病毒Mpro 與ZINC 數據庫對接后復合物之一，下圖為新冠病毒RdRp 蛋白與ZINC（紫色）及REAL數據庫對接后復合物之一（綠色））。

圖8 蛋白與活性分子復合物Fig.8 Protein and active ligands complex

5 結束語

本文開發了超大規模虛擬篩選并行分發系統，利用負載均衡優化和高容錯技術等實現了超大規模快速虛擬篩選計算，支持國產加速器，系統易用性好、穩定性強。系統采用“分而治之”的策略，中央數據庫及下級多層數據庫管理和分發作業，主從模式并發計算，可有效利用異構計算部件，實現近線性加速效果。通過多層級任務分發、文件打包遷移及計算和后處理流程式設計，實現節點間IO 和存儲的性能優化，可有效處理數十億級分子；自動掃描節點、快速重啟數據庫服務及相應計算程序，實現長時間穩定運行；屏蔽復雜處理細節，通過配置文件靈活調整程序輸入及運行參數等，有一定易用性。

未來工作主要包括以下三方面：一是進一步利用異構加速提升計算速度；二是拓寬程序應用領域和范圍，例如使用多個對接程序進行共對接等；三是進一步提升易用性，通過程序封裝并結合前端網絡應用提供更加方便易用的使用體驗。

致謝本文的部分計算得到“東方”超級計算系統的支持與幫助。