基于原子特性知識增強的分子毒性預測方法

方舒言,劉 宇,2*,侯阿龍,秦歡歡,劉 嵩,4

(1.武漢科技大學 計算機科學與技術學院,湖北 武漢 430072;2.湖北省智能信息處理與實時工業系統重點實驗室,湖北 武漢 430072;3.武漢科技大學 醫學院,湖北 武漢 430072;4.湖北省職業危害識別與控制湖北省重點實驗室,湖北 武漢 430072)

0 引 言

新藥研發存在周期長、費用高和成功率低等特點。人工智能技術是近些年來的熱點技術之一,在很多領域都有非常廣泛的應用,多種人工智能方法已經成功應用于藥物的發現過程[1]。當前化學分子毒性預測方法主要從分子結構的序列表示中學習特征,如SMILES[2](Simplified Molecular Input Line Entry System)。SMILES通過對分子圖的生成樹實施深度優先的前序遍歷,為每個原子、鍵、樹遍歷決策和斷環生成符號,定義了分子的字符串表示,所得到的字符串對應于分子圖生成樹的序列化。因比其他表示結構的方法(包括圖)更緊湊,SMILES已被廣泛用于分子毒性預測。

目前,以Chemberta[3]為代表的分子毒性預測方法借鑒了自然語言處理技術的思想,直接使用大量無標注的SMILES作為語料庫來學習分子表征,并使用學習到的分子表示用于毒性預測。然而SMILES語法復雜且限制性強,常規字符集上的大多數序列不能明確定義分子。

由于上述限制,基于SMILES的深度學習模型只能將學習重點放在分子串的語法上,導致現有的毒性預測方法僅損失了分子中的原子特性等知識。針對上述問題,以Mol-BERT[4]為代表的方法選擇了使用摩根指紋標識符(Morgan Fingerprint Identifier,MFI)序列作為輸入來進行模型的訓練,并在MFI融入了原子特性的知識(如原子序數、電荷量、原子的度等)。摩根指紋(Morgan Fingerprints)是一種使用哈希算法生成的圓形分子指紋[5-6]。其中,原子的MFI是對原子及其屬性(如原子序數、電荷量等)進行哈希得到的一個定長整數,而子結構MFI則是對分子中某原子MFI及以該原子為中心一定半徑范圍內的其他原子的MFI進行聯合哈希得到的定長整數。根據設置的指紋半徑的不同,可以得到蘊含不同知識的MFI,比如以某原子為中心且半徑為0生成的MFI僅包含該原子的知識,以某原子為中心且半徑為1生成的MFI,則蘊含該原子及相鄰原子子結構知識。

根據上述對摩根指紋及MFI的說明可知,Mol-BERT雖在一定程度上融入了原子特性的知識,但是仍然存在局限。比如,Mol-BERT的輸入MFI序列中未顯式引入氫原子的知識,而且對于MFI序列的半徑信息Mol-BERT方法并沒有進行顯式的編碼。

針對目前基于分子序列的毒性預測方法均未考慮氫原子和摩根指紋半徑等領域知識的問題,該文提出了一種基于原子特性知識增強的分子毒性預測方法。一方面,改進的分子毒性預測方法考慮了分子中所有氫原子,輸入的MFI序列顯式引入了氫原子的MFI和以氫原子為中心的子結構的MFI。在考慮了分子中所有的氫原子之后,子結構MFI所代表的子結構粒度更小,模型的MFI字典大小變為引入氫原子之前的約一半。這是由于在同樣的半徑限制下,使得原本一些大而稀有的子結構因為氫原子的引入將被轉換成小而常見的子結構,如圖1(b)所示。該現象類似于自然語言處理領域的WordPiece[7],它將單詞拆成詞根和詞綴,一些生僻詞被拆成常見詞根和詞綴的組合,能夠減小自然語言處理模型的詞表大小。另一方面,受SongNet[8]的啟發,為了讓模型學習MFI序列時可以區分不同半徑的MFI,提出使用摩根指紋半徑的知識擴展嵌入層的方法,使模型學習到的原子MFI和子結構MFI有所區分。圖1(a)～(c)描述了氫原子和摩根指紋半徑知識增強MFI序列的過程及增強后MFI序列與原本MFI序列的區別。該文的主要貢獻包括3個方面：

(a)原本的MFI序列

(1)提出生成分子MFI序列時顯式地引入氫原子的方法,使得分子毒性預測模型可以學習到氫原子信息,同時減少了MFI字典的大小和模型的參數量;

(2)改進的分子毒性預測方法,額外考慮了各原子及相鄰原子構成的子結構的知識,并在模型的嵌入層強化了摩根指紋半徑知識;

(3)為了驗證改進方法的有效性,利用Pubchem[9]分子庫中的1 000萬分子預訓練了基于原子特性增強的分子毒性預測模型(Hydrogen and Radius enhanced Toxicity Prediction model,HRToxPred),并在主流的毒性預測數據集Tox21[10]上進行了驗證。

1 相關工作

在基于深度學習的分子毒性預測領域中,主流方法直接在分子的SMILES字符串上訓練模型來預測分子毒性。表1總結了基于分子序列的毒性預測模型的相關工作。FP2VEC[11]和Mol2vec[12]使用了自然語言處理領域的Word2vec[13]思想。其中,Mol2vec學習的是MFI的低維向量表示,然后輸入學習到的MFI向量表示至下游網絡中預測分子毒性,而FP2VEC學習的是比特分子指紋比特串的低維向量表示。然而,在自然語言處理領域中Word2vec的詞和向量是一對一的關系,所以無法解決多義詞的問題。在化學領域,同樣存在類似“一詞多義”的場景,例如相同的原子所處的環境不同其作用也不同,所以FP2VEC,Mol2vec存在著類似的問題。Chemberta和Mol-BERT使用了基于Transformer的掩碼語言模型。其中,Chemberta設計了一種SMILES的分詞方法,使用Roberta[14]的預訓練策略練了基于SMILES的分子毒性預測模型,然而SMILES的分詞方法需要化學領域的專家來設計,需要消耗大量人力。Mol-BERT使用分子的MFI序列作為語料庫,使用BERT[15]的MLM(Masked Language Model)策略預訓練了基于MFI序列的分子毒性預測模型。相比于Chemberta,Mol-BERT所使用的MFI序列不需要人工分詞。然而Mol-BERT使用的MFI序列未顯式引入氫原子且序列中原子指紋與子結構指紋混雜在一起,使模型學習到的分子表示并未區分原子指紋與子結構指紋。針對上述方法存在的問題,該文提出生成MFI序列時顯式引入氫原子并對半徑信息進行顯式編碼的知識增強方法。

表1 相關工作總結

2 方 法

2.1 摩根指紋標識符序列

為了便于解釋顯式引入氫原子后MFI的計算方法,現給出以下定義：

定義1 原子鄰域(Atom Neighborhood)：給定一個分子圖mg=(A,B),A為該分子中的所有原子的集合,B為該分子中所有的化學鍵的集合。對于一個原子a∈A,給定半徑r,以a為中心,r為半徑的原子鄰域定義為nbr(mg,a,r)={[a,at,b]|a∈A,at∈A,b∈B,a≠at,dist(a,at)≤r,b=},其中at為與a通過化學鍵b相連的原子,dist(·)表示分子圖中兩個原子之間的跳數。

將原子鄰域的半徑設為1,計算MFI時需考慮到兩種不同規模的原子鄰域nbr(mg,a,0)和nbr(mg,a,1)。如圖2所示,深色節點對應nbr(mg,a,0),而由深色節點a和淺色節點at組成的子結構對應nbr(mg,a,1)。由此,原子的MFI和子結構的MFI將由其對應的原子鄰域特性經過哈希得到。

圖2 原子鄰域、摩根指紋、摩根指紋標識符之間的關系

定義2 摩根指紋[5-6](Morgan Fingerprint)：摩根指紋是一個稀疏的比特串。一個分子可以通過散列分子中每個原子的鄰域信息得到一個索引,該索引位置的比特位被置為1。

定義3 摩根指紋標識符(Morgan Fingerprint Identifier,MFI)：摩根指紋標識符是一個固定長度的整數,用來表示原子鄰域的特征。通常情況下,一個原子鄰域的摩根指紋標識符可以對該原子鄰域的屬性(原子序數、電荷量等)進行哈希來獲得。記為mfinbr(mg,a,r)=Hash(Attrs(nbr(mg,a,r))),Hash(·)為哈希函數,Attrs(·)表示獲取原子鄰域的屬性。

原子鄰域、摩根指紋和摩根指紋標識符之間的關系如圖2所示。對于一個分子經過算法1的計算,可以獲得該分子的MFI字典。字典的鍵為MFI,字典的值為由該MFI對應原子鄰域中心原子的原子索引和原子鄰域半徑組成的元組[a.Index,r]。該分子的摩根指紋則是一個長的比特串,其中索引為MFI字典中鍵值的比特位被置為1,如圖2箭頭指示。

定義4 原子的摩根指紋標識符：給定一個分子mg=(A,B)和一個原子a∈A,該原子的摩根指紋標識符可以通過哈希原子鄰域nbr(mg,a,0)的屬性得到。

算法1

輸入：分子圖mg=(A,B),半徑r

輸出：存儲該分子MFI的字典D

Initialization;

/*初始化一個容器存所有的mfi,偏移量為原子索引*/

Vectorinvars;

fora∈Ado

/*初始化一個容器存Attrs(nbr(mg,a,0))*/

Vectorcomponents;

components.pushback(a.AtomicNum);//原子序數

components.pushback(a.TotalDegree);//原子的度

ifa≠Hydrogen then

//非氫原子連接的氫原子數

components.pushback(a.TotalNumHs);

end

components.pushback(a.FormalCharge);//原子的電荷量

invars.pushback(boost：：hash(components));

end

/*將mfinbr(mg,a,0)的信息存入字典*/

fora∈Ado

D[invars[a.Index]].pushback(Tuple(a.Index,0));

end

/*初始化容器存mfinbr(mg,a,1)*/

Vectornbrs;

fora∈Ado

/*獲取從a出發的化學鍵*/

bonds = getChemicalBonds(a);

for ∈bonds do

ifa≠Hydrogen then

nbrs.pushback(Tuple(.type

invars[at.Index]));

end

end

invar = 0;

boost：：hashcom bine(invar, invars[a.Index]);

for nbr ∈ nbrs do

boost：：hashcombine(invar, nbr);

end

/*將mfinbr(mg,a,1)的信息存入字典*/

D[invar].pushback(Tuple(a.Index, 1));

end

該文考慮的原子鄰域nbr(mg,a,0)的屬性包含原子序數、電荷量、非氫原子連接的氫原子數、原子在分子圖中的度。算法1的第3～18行具體描述了該計算摩根原子MFI的過程。

定義5 子結構的摩根指紋標識符：給定一個分子mg=(A,B)和一個原子a∈A,以a為中心的子結構的摩根指紋標識符可以通過哈希原子鄰域nbr(mg,a,1)的屬性得到。

該文考慮的原子鄰域nbr(mg,a,1)的屬性包含該原子鄰域中化學鍵b的類型、該原子鄰域中原子a和at的摩根指紋。算法1的19～37行具體描述了計算過程。

定義6 摩根指紋標識符序列：摩根指紋標識符序列是由一個分子中所有的原子鄰域的摩根指紋標識符組成的,摩根指紋標識符的順序按照SMILES定義的原子順序排列。如果對于分子中某原子a有多個摩根指紋標識符,則按照r由小到大進行排序。

根據以上定義,算法1描述了計算最大半徑設置為1時的分子的MFI的過程。

2.2 數據預處理

數據預處理階段,因為同一個分子可以有多種不同的SMILES形式,所以首先要對分子SMILES進行標準化。SMILES之所以出現不同是因為可以隨意修改原子的“出場順序”以得到不同的SMILES。該文使用Rdkit[16-17]實現了標準化過程,并根據分子化合價態平衡原理補全SMILES中的氫原子。最后使用算法1計算并排序生產分子的MFI序列。

根據上述方法,該文選擇Pubchem分子庫中1 000萬分子的SMILES,并使用算法1分子的MFI并排序形成MFI序列。整個過程中有關分子的處理都使用Rdkit分子處理工具實現。統計1 000萬分子的是否顯式引入氫原子的MFI序列,得到不同MFI序列的MFI字典信息,如表2所示。

表2 是否引入氫的MFI字典大小比較

顯式引入氫原子MFI之后,分子中的某些相對較大的子結構標識符可以分解為一些相對較小的子結構的組合,例如圖1中以該羥基(-OH)上氧原子為中心半徑為1的子結構會被拆分成以羥基中氫原子為中心半徑為1的羥基結構和以氧原子為中心半徑為1的結構。這類似于自然語言處理領域把單詞拆成詞根與詞綴的方法。因此,引入顯式表示氫原子之后,生成的MFI字典大小更小且MFI對應的原子鄰域粒度更細。該文嘗試擴展算法1的最大半徑限制,當最大半徑達到2時,即便引入了氫原子一定程度上減小了MFI字典的規模,然而最大半徑限制為2時MFI字典的大小超過了100萬,因此該文僅考慮最大半徑限制為1的情況。

2.3 模型架構

HRToxPred包含Embedding模塊、Transformer模塊和下游任務模塊,如圖3所示。圖中M為輸入的MFI,F為MFI嵌入(Token嵌入),P為位置編碼(Position嵌入),R為半徑嵌入(Radius嵌入)。因為生成的每個MFI序列代表一個分子,而每個分子都是單獨的個體,無需考慮兩個分子之間的關系,所以文中的預訓練去掉了BERT模型預訓練策略中“下一句預測”(Next Sentence Prediction,NSP)預訓練任務,輸入樣本變為單個MFI序列,因此BERT模型原本用來區分兩個句子的Segment編碼失去了意義。該文修改了BERT的Embedding模塊,將Segment嵌入改為摩根指紋半徑嵌入來區分MFI的類別,最終Embedding模塊包含半徑嵌入、位置嵌入和MFI編碼嵌入。

圖3 模型架構示意圖

經過上述步驟輸入數據就處理完成,接下來進行嵌入。對于輸入特征,首先要獲得MFI編碼、半徑編碼、位置編碼的嵌入,通過把所有嵌入向量相加的方式,公式如下：

Ht=Emt+Ept+Ert

(1)

式中,t為序列中MFI的索引,Emt,Ept,Ert分別為MFI嵌入、位置嵌入和半徑嵌入,即通過半徑嵌入Ert區分了原子MFI和子結構MFI。處理完輸入數據和嵌入之后,HRToxPred使用自注意力機制進行學習。

qt=WQHt

(2)

kt=WKHt

(3)

vt=WVHt

(4)

ct=self-att(qt,kt,vt)

(5)

式中,WQ,WK,WV為參數矩陣,self-att(·)為注意力函數,Ct為該注意力頭輸出。根據上述方法將得到的所有注意力頭的Ct向量進行拼接,然后與參數矩陣W求積得到輸出c,再將c輸入到激活函數中得到對應的輸出o。在下游任務中只需將o輸入到下游任務的分類器網絡中即可。式1～7描述了模型的計算過程。

c=Concat(ct)W

(6)

o=max(0,ZW1+b1)W2+b2

(7)

3 實驗與分析

3.1 基準模型

選取了當前被廣泛使用的基于分子序列的毒性預測模型Chemberta和Mol-BERT作為對比模型。其中,Chemberta基于Roberta模型實現。預訓練時,使用1 000萬分子的SMILES作為語料庫,隨機遮蔽了每個輸入字符串中15%的Token。在學習恢復Mask Token時,該模型形成了一個化學空間的分子結構表示。Mol-BERT使用未顯式引入氫原子MFI和摩根指紋半徑特征的MFI序列作為輸入。Mol-BERT同樣使用了1 000萬分子的MFI序列,并利用了BERT的架構及預訓練策略進行了預訓練,最后將預訓練好的模型進行下游的毒性預測任務。

3.2 預訓練

HRToxPred與主流的自然語言處理預訓練語言模型的方式不同,去掉了NSP預訓練任務,僅留MLM任務。具體來說,該文使用bert-base的模型設置,并將分子MFI序列中隨機15%的標識符屏蔽為[MASK]。此外,預訓練使用了動態掩碼,即同一個MFI序列會產生不同的掩碼序列,這樣可以重復使用預訓練樣本,而且可以防止過擬合。

3.3 微 調

選擇主流毒性預測數據集Tox21的12個分類任務進行微調。數據集詳情見表3。

表3 Tox21各任務細節

3.4 評估標準

由于該文選取的數據集正負樣本不平衡,而精度、準確率、召回率等指標受正負樣本不平衡的影響,不能客觀反映出模型的性能,且在分子毒性預測中區分正負樣本的概率閾值并不確定,所以選擇AUC[18](Area Under Curve)作為評測指標可避免以上因素帶來的影響。

AUC被定義為ROC(Receiver Operating Characteristic)曲線下的面積。ROC曲線,又稱接受者操作特征曲線。二分類任務分類閾值的設定不同會得出不同的真陽率(True Positive Rate,TPR)和假陽率(False Positive Rate,FPR),將同一模型每個閾值的(FPR,TPR)坐標都畫在ROC空間里,就成為特定模型的ROC曲線。AUC計算公式如下：

(8)

式中,P為正樣本數量,N為負樣本數量,pi為正樣本預測得分,nj為負樣本預測得分。

3.5 實驗設置

HRToxPred包含Embedding模塊、Transformer模塊和下游任務模塊,模型的實現框架使用UER-py[19]。

Embedding模塊包含Token嵌入、Position嵌入和半徑嵌入。HRToxPred設置的輸入序列長度與BERT保持一致均為512個Token。

Transformer模塊選擇bert-base結構,bert-base包含12個Transformer層和12個自注意力頭,隱藏層維度為768。

下游任務模塊的微調網絡使用UER-py框架中默認的網絡,即在預訓練完成的模型后接2層全連接層。對于每個數據集,按照80/10/10的比例隨機劃分為訓練集/驗證集/測試集,并進行了32個epoch的微調,具體參數見表4。

表4 模型參數設置

3.6 消融實驗分析

為了進一步驗證HRToxPred中引入的摩根指紋半徑特征及考慮氫原子的作用,進行了消融實驗,即用Pubchem中的1 000萬分子預訓練如下四個模型：

?Base：輸入MFI序列不顯式引入氫原子及其子結構MFI,嵌入層沒有額外的摩根指紋半徑嵌入;

?AddHS：輸入MFI序列不顯式引入氫原子及其子結構MFI,嵌入層沒有額外的摩根指紋半徑嵌入;

?AddRadius：輸入MFI序列不顯式引入氫原子及其子結構MFI,嵌入層有額外的摩根指紋半徑嵌入;

?HRToxPred：MFI序列顯式引入氫原子及其子結構MFI,嵌入層有額外的摩根指紋半徑嵌入。

然后將上述四個模型在Tox21上微調。如圖4所示,使用顯式引入氫原子MFI的模型預訓練收斂得更快,主要原因是HRToxPred考慮了氫原子及其子結構MFI造成的詞表減小。另外,從表5中Tox21不同通道上的上的AUC分數對比來看,僅使用考慮了氫原子MFI序列作為輸入(AddHs)和僅引入摩根指紋半徑嵌入(AddRadius)都可以提升模型的分類能力,而兩者的組合對模型分類能力的提升更顯著。HRToxPred相比于Base模型,在sr_p53,nr_ahr,nr_er,nr_ar,sr_are,nr_ar_lbd,nr_er_lbd,sr_hse等任務上取得了顯著的提升。具體而言,sr_p53任務提升了4.5%,nr_ahr任務提升了1.7%,nr_er任務提升了3.1%,nr_ar任務提升了7.5%,sr_are任務提升了4.8%,nr_ar_lbd任務提升了2.6%,nr_er_lbd任務提升了1.5%,sr_hse提升了3.3%。而在nr_aromatase,nr_ppar_gamma,sr_atad5,sr_mmp等任務上的表現不如預期。

表5 消融實驗結果比較(ROC-AUC↑)

圖4 Base和HRToxPred預訓練Loss曲線

3.7 對比實驗分析

為了檢驗HRToxPred的有效性,將HRToxPred與Chemberta和Mol-BERT模型進行了比較。表6報告了各個模型Tox21數據集12個通道上的AUC分數。HRToxPred在sr_p53,nr_ahr,nr_er,nr_ar,sr_are,nr_ar_lbd,nr_aromatase,nr_er_lbd,sr_hse等任務上取得了最佳的表現。同時,如圖5所示,HRToxPred在Tox21的雷達圖中擁有最大的面積。

表6 對比實驗結果比較(ROC-AUC↑)

圖5 Tox21各任務AUC分數雷達圖

4 結束語

該文提出了引入摩根指紋半徑以及顯式引入氫原子的知識增強方法,并預訓練了一個原子知識增強的分子毒性預測模型HRToxPred。該模型與現有的基于序列的分子毒性預測模型有兩個方面的優勢：一方面是顯式地引入了氫原子的知識,使得MFI字典規模更小,進而減少了模型參數量并提升了毒性預測的性能;另一方面是額外引入了摩根指紋半徑的知識,增強了MFI序列的順序特征。從對比實驗和消融實驗的結果可見,這兩方面能使模型性能有了顯著提升。

雖然該方法在一些毒性預測任務數據集上展現了良好的性能,但仍有一些局限性。從消融實驗的結果來看,HRToxPred在Tox21數據集中nr_aromatase,nr_ppar_gamma,sr_atad5,sr_mmp等任務上的表現不如預期,后續工作會結合深度學習理論和結構化學的知識來詳細分析引入氫原子及半徑對不同毒性預測任務的影響。此外,根據化學的領域知識可知,分子中的官能團對分子的特性有較大影響[20],因此后續工作會考慮引入相關領域知識。