用於計算分子生物學的 Python 工具
用於計算分子生物學的 Python 工具
Mocapy++ 是一個用於訓練和使用貝氏網路的機器學習工具包。它已被用於開發生物分子結構的機率模型。此專案的目標是開發 Mocapy++ 的 Python 介面，並將其與生物Python整合。這將允許使用從資料庫中提取的資料來訓練機率模型。Mocapy++ 與生物Python的整合將為蛋白質結構預測、設計和模擬領域提供強有力的支援。
發現生物分子的結構是生物學中最大的問題之一。給定胺基酸或鹼基序列，其三維結構是什麼？生物分子結構預測的一種方法是建構機率模型。貝氏網路是一種由一組變數及其聯合機率分佈組成的機率模型，以有向無環圖表示。動態貝氏網路是一個代表變數序列的貝氏網路。這些序列可以是時間序列或符號序列，例如蛋白質序列。定向統計主要關注的是平面或三維空間中的單位向量的觀測值。樣本空間通常是一個圓或一個球體。必須有特殊的定向方法來考慮樣本空間的結構。圖形模型和定向統計的結合允許開發生物分子結構的機率模型。透過使用具有定向輸出的動態貝氏網路，可以建構序列和結構的聯合機率分佈。生物分子結構可以用幾何上自然的連續空間表示。Mocapy++ 是一個用於使用動態貝氏網路進行推斷和學習的開源工具包，它為定向統計提供了支援。Mocapy++ 非常適合建構生物分子結構的機率模型；它已被用於開發蛋白質和 RNA 原子級結構的模型。Mocapy++ 已在多篇具有高度影響力的出版物中使用，並將成為分子建模套件 Phaistos 的核心，該套件即將發布。此專案的目標是開發一個非常有用的 Mocapy++ Python 介面，並將該介面與生物Python專案整合。透過 Bio.PDB 模組，生物Python為生物分子結構資料庫的資料探勘提供了出色的功能。整合 Mocapy++ 和生物Python將允許使用從資料庫中提取的資料來訓練機率模型。將 Mocapy++ 與生物Python整合將為研究人員創建一個強大的工具包，以便快速實施和測試新想法、嘗試各種方法並完善他們的方法。它將為生物分子結構預測、設計和模擬領域提供強有力的支援。
Michele Silva (michele.silva@gmail.com)
指導老師
Thomas Hamelryck
Eric Talevich
'''了解 SEM 和定向統計'''
'''研究 Mocapy++ 的使用案例'''
'''使用 Mocapy++'''
'''設計 Mocapy++ 的 Python 介面'''
'''實作 Python 綁定'''
'''探索 Mocapy++ 的應用'''
已經有人嘗試使用 Swig 為 Mocapy++ 提供綁定。但是，如果需要效能，Swig 不是最佳選擇。Sage 專案旨在提供 Mathematica 或 Maple 的開源替代方案。Cython 是與 Sage 一起開發的（雖然它是一個獨立的專案），因此它基於 Sage 的需求。他們嘗試了 Swig，但由於效能問題而拒絕了它。根據Sage 程式設計指南，「這個想法是在 C++ 中為 SAGE 編寫需要快速的程式碼，然後將其封裝在 SWIG 中。這停滯不前，因為結果不夠快。首先，首先在 C++ 中編寫程式碼時存在開銷。其次，SWIG 在 Python 和執行實際工作的程式碼之間產生了多個程式碼層。」這是在 2004 年寫的，但似乎情況沒有太大變化。我考慮 Swig 的唯一原因是未來將 Mocapy++ 綁定包含在 BioJava 和 BioRuby 專案中。
Boost Python 是全面且廣泛被 Python 社群接受的。我會因為它廣泛的使用和測試而選擇它。我會因為它很難偵錯且建置系統複雜而拒絕它。我認為僅為了建立 Python 綁定而包含 boost 依賴項是不值得的，但由於 Mocapy++ 已經依賴 Boost，因此使用它成為一個更有吸引力的選項。以我個人的經驗而言，Boost Python 非常成熟，並且可以對其進行任何操作。就效能而言，Cython 仍然超越了它。看看Cython C++ 封裝基準測試並檢查 Cython 與 Boost Python 的計時。之前也有比較 Swig 和 Boost Python 的基準測試。
根據網路上多個基準測試，它比其他用於建立 C++ 程式碼的 Python 綁定的選項快得多。檢查 Cython 和 Boost.Python 之間的簡單基準測試。它也非常乾淨且簡單，但功能強大。Python 的關於移植擴充模組的文件提到了 cython：「如果您正在編寫一個新的擴充模組，您可能會考慮 Cython。」Cython 現在支援與 numpy 陣列進行有效的互動。它是一種年輕但不斷發展的語言，我絕對會因為它的精簡和速度而嘗試它。
由於 Boost 得到很好的支援，並且 Mocapy++ 已經依賴它，我們決定使用 Boost.Python 進行綁定。
如需更多資訊，請參閱Mocapy++生物Python - 想法箱。
原型的原始碼位於 gSoC11 分支：http://mocapy.svn.sourceforge.net/viewvc/mocapy/branches/gSoC11/bindings_prototype/
用於幾個 Mocapy++ 功能的綁定和幾個範例，以找出可能的實作和效能問題。
程序
Mocapy++ 的介面保持不變，因此測試看起來與 Mocapy/examples 中的測試相似。
在原型中，所有綁定都在單一模組中實作。對於實際實作，我們可以鏡像 src 套件結構，為每個套件（例如 discrete、inference 等）分離綁定。
可以實作執行範例所需的所有功能。在為 MDArray 的向量建立綁定時，無法使用 vector_indexing_suite。必須實作一些運算符（在 MDArray 中）才能將可索引的 C++ 容器匯出到 Python。
在 Python 中實作了兩個使用離散節點的 Mocapy++ 範例。公開 Mocapy 的資料結構和演算法沒有問題。Python 版本的效能非常接近原始 Mocapy++。
如需更多詳細資訊，請參閱 Mocapy++ 綁定原型報告。
'''資料結構'''
Mocapy 使用內部資料結構來表示陣列：MDArray。為了讓使用者更容易與 Mocapy 的 API 互動，決定提供一個接受 numpy 陣列的介面。因此，有必要實作 numpy 陣列和 MDArray 之間的轉換。
從 MDArray 到 Python 的轉換是透過使用 Boost.Python to_python_converter 完成的。我們實作了一個範本方法 convert_MDArray_to_numpy_array，它會將任何基本類型的 MDArray 轉換為對應的 numpy 陣列。為了執行轉換，原始陣列的形狀和內部資料會複製到新的 numpy 陣列中。
numpy 陣列是使用 Numpy 陣列 API 建立的。使用現有資料 (PyArray_SimpleNewFromData) 建立新的 PyArrayObject 並不會複製陣列資料，它只會儲存指向它的指標。因此，只有在沒有對物件的參考時才能釋放資料。這是透過使用膠囊來完成的。除了封裝資料外，膠囊還儲存一個在銷毀陣列時使用的解構函數。PyArrayObject 有一個名為「base」的欄位，它指向該膠囊。
從 Python 到 C++ 的轉換，即從 numpy 陣列建立 MDArray，稍微複雜一些。Boost.Python 將提供一塊記憶體，新的 C++ 物件必須在該記憶體中就地建構。有關更多詳細資訊，請參閱 如何編寫 boost.python 轉換器文章。
還實作了基本類型 (double、int…) 的 std::vector 和 Python 清單之間的轉換。對於自訂類型的 std::vector，例如 Node，未執行到 Python 清單的轉換。如果採用與基本類型相同的方式執行，則會引發類型錯誤：「TypeError: No to_python (by-value) converter found for C++ type」。使用 vector_indexing_suite 時，此問題已經解決。請參閱 封裝 std::vector<AbstractClass*>。使用 vector_indexing_suite 的唯一不便是建立諸如 vector_Node 之類的新類型，而不是使用標準的 Python 清單。
翻譯的程式碼位於 mocapy_data_structures 模組中。
’’’ 核心函數 ‘’’
mocapy Python 套件遵循 Mocapy 目前的原始碼樹狀結構。對於每個套件，都會建立一個包含綁定的共享程式庫。這樣可以加快編譯速度並簡化除錯。此外，如果只建立一個程式庫，將無法定義套件。
每個程式庫都稱為 libmocapy_
綁定程式碼可以在 Bindings 目錄中找到。
目前，正在開發針對剛建立的介面的測試。在 framework 套件下已經實作了一些測試：mocapy/framework/tests
'''資料結構'''
在實作其餘 Mocapy++ 功能的綁定時，遇到需要指向 mdarray 的指標和參考的方法的問題
DBN 表示九個連續的二面角，其中七個中心角來自單個核苷酸。每個切片 j（三個變數的列）對應於 RNA 片段中的一個二面角。每個切片中的變數為：一個角度識別符號 Dj、一個隱藏變數 Hj 和一個角度變數 Aj。角度識別符號用於追蹤切片表示哪個二面角，而角度節點則模擬實際的二面角值。隱藏節點會在序列中所有角度之間（而不僅僅是相鄰切片中的角度之間）引起依賴關係。Barnacle 的原始原始碼（其中包含以 Python 編寫的 Mocapy 嵌入版本）可以在 <http://sourceforge.net/projects/barnacle-rna> 找到。Barnacle 的修改版本已更改為與 Mocapy 綁定一起使用，可以在 <https://github.com/mchelem/biopython/tree/master/Bio/PDB/Barnacle> 找到。以下是一個使用範例： ``` python model = Barnacle("ACCU") model.sample() print("log likelihood = ", model.get_log_likelihood()) model.save_structure("structure01.pdb") ``` #### TorusDBN Wouter Boomsma、Kanti V. Mardia、Charles C. Taylor、Jesper Ferkinghoff-Borg、Anders Krogh 和 Thomas Hamelryck。局部蛋白質結構的生成機率模型。Proc Natl Acad Sci U S A. 2008 年 7 月 1 日；105(26)：8932–8937。 <http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2440424/> TorusDBN 旨在根據生物分子的胺基酸序列預測其 3D 結構。這是一個針對蛋白質局部序列-結構偏好的原子細節的連續機率模型。蛋白質的骨幹可以用二面角對 φ 和 ψ 的序列來表示，這在 [Ramachandran 圖](http://en.wikipedia.org/wiki/Ramachandran_plot) 中是廣為人知的。兩個角度（值範圍均為 −180° 到 180°）定義了環面上的一個點。因此，蛋白質的骨幹結構可以完全參數化為此類點的序列。 
圓形節點表示隨機變數。箭頭上的矩形框說明了它們之間條件機率分佈的性質。隱藏節點發射角度對、胺基酸資訊、二級結構標籤和順/反資訊。TorusDBN 模型最初作為 backboneDBN 套件的一部分實作，該套件可在 <http://sourceforge.net/projects/phaistos/> 免費取得。此專案中實作了新版本的 TorusDBN 模型，可以在 <https://github.com/mchelem/biopython/tree/master/Bio/PDB/TorusDBN> 找到。TorusDBNTrainer 可以用於使用給定的訓練集訓練模型： ``` python trainer = TorusDBNTrainer() trainer.train(training_set) # training_set 是一個檔案列表 model = trainer.get_model() ``` 然後可以使用該模型來取樣新序列： ``` python model.set_aa("ACDEFGHIK") model.sample() print(model.get_angles()) # 取樣的角度。 ``` 建立模型時，可以建立一個新的 DBN，指定隱藏節點的大小，或從檔案載入 DBN。 ``` python model = TorusDBNModel() model.create_dbn(hidden_node_size=10) model.save_dbn("test.dbn") ``` ``` python model = TorusDBNModel() model.load_dbn("test.dbn") model.set_aa("ACDEFGHIK") model.sample() print(model.get_angles()) # 取樣的角度。 ``` 也可以使用 find\_optimal\_model 方法來選擇隱藏節點的最佳大小： ``` python trainer = TorusDBNTrainer() hidden_node_size, IC = trainer.find_optimal_model(training_set) model = trainer.get_model() ``` IC 是 [貝葉斯資訊準則](http://en.wikipedia.org/wiki/Bayesian_information_criterion) (BIC) 或 [赤池資訊準則](http://en.wikipedia.org/wiki/Akaike_information_criterion) (AIC)（預設為 BIC。可以透過設定 use\_aic 標記來指定 AIC）。有關模型 API 的更多詳細資訊，請參閱測試檔案：<https://github.com/mchelem/biopython/blob/master/Tests/test_TorusDBNTrainer.py> 和 <https://github.com/mchelem/biopython/blob/master/Tests/test_TorusDBNModel.py>。 ### 效能 已進行了一些效能測量，比較了以 C++ 和 Python 實作的測試案例。測試在具有以下規格的電腦上執行：Core 2 Duo T7250 2.00GHz、記憶體雙通道 4.0GB (2x2048) 667 MHz DDR2 SDRAM、硬碟 200GB 7200RPM。沒有顯著的效能差異。對於兩種實作，消耗最多 CPU 時間的方法是相同的： 
效能分析測試使用 [Callgrind](http://valgrind.org/info/tools.html#callgrind) 進行，並使用 [Kcachegrind](http://kcachegrind.sourceforge.net/) 進行視覺化。以下是 Mocapy 提供的範例的平均執行時間（10 次執行）： | 測試名稱 | C++ (s) | Python (s) | |----------------------------|---------|------------| | hmm\_simple | 0.52 | 0.58 | | hmm\_discrete | 48.12 | 43.45 | | discrete\_hmm\_with\_prior | 55.95 | 50.09 | | hmm\_dirichlet | 340.72 | 353.98 | | hmm\_factorial | 0.01 | 0.12 | | hmm\_gauss\_1d | 53.97 | 63.39 | | hmm\_gauss | 16.02 | 16.96 | | hmm\_multinomial | 134.64 | 125.83 | | hmm\_poisson | 11.00 | 10.60 | | hmm\_vonmises | 7.22 | 7.36 | | hmm\_torus | 53.79 | 53.65 | | hmm\_kent | 61.35 | 61.06 | | hmm\_bippo | 40.66 | 41.81 | | infenginehmm | 0.01 | 0.12 | | infenginemm | 0.01 | 0.15 | #### TorusDBN 儘管 PDB 檔案已讀取、剖析並轉換為 mocapy 可以理解的格式，但最耗時的方法是在取樣過程中執行數學運算的方法（例如 Chebyshev 和 exp）。 
該模型已使用由約 950 個鏈組成的訓練集進行訓練，這些鏈的最大同源性為 20%，解析度低於 1.6 Å，R 因子低於 25%。讀取和訓練整個資料集大約需要 67 分鐘。產生的 DBN 可在 <https://github.com/mchelem/biopython/blob/master/Tests/TorusDBN/pisces_dataset.dbn> 取得，並且可以直接載入到模型中，如上面的 TorusDBN 部分所述。 ### 未來工作 這個夏天結束了，但工作仍在繼續... 我還有很多事情打算做： - 測試經過訓練的模型，以檢查它們在蛋白質結構預測中的有效性。- 嘗試減少動態配置，因為它是導致大量執行時間的原因。- 保證繫結中沒有記憶體洩漏。