機器學習預測內分泌干擾物水生生物毒性效應

王藝霖, 范俊韜 , 王書平, 黃國鮮, 閆振廣

1. 上海海洋大學海洋生態與環境學院,上海 201306

2. 中國環境科學研究院環境基準與風險評估國家重點實驗室,北京 100012

研究表明,含內分泌干擾物(endocrine disruptor chemicals, EDCs)類的化學品在農業、工業和日常生活中被廣泛使用[1],已在廢水、地表水、自來水中陸續檢出,表明其對水生生物乃至人類的影響正在逐漸擴大[2-5]。 EDCs 可以直接作用于內分泌系統,能夠以阻斷或模仿人類和動物體內自然激素的方式干擾激素行為,從而對心血管、代謝、免疫,尤其是生物的生殖系統造成影響,導致種群數量下降[6-9];大部分的EDCs 具有低劑量有效性、半衰期長和生物富集、生物放大等特點,因此會在環境中持久存在,造成較為長遠的影響[10-12]。 研究數據表明,我國多處水域均受到EDCs 污染,由此帶來的生態風險需要引起高度的重視[13-15]。

EDCs 生態風險的科學評估則依賴于繁殖毒性數據的獲取。 EDCs 的繁殖毒性數據主要來自與生物的生活史或部分生活史相關的實驗。 這些實驗周期長、成本高,難以在短期內積累足夠的EDCs 繁殖毒性數據,使得 EDCs 的生態風險評估非常困難[15-18],不利于以后科學開展生態風險評估和環境管理工作。 使用數學模型來預測毒性效應已成為國際生態毒理學研究熱點[19]。 數學建模工具可以在一定的框架下對現有的毒性實驗進行拓展,有利于深入了解劑量與反應關系之間的復雜性[15,20],從而保護生態系統,降低生態風險。 使用模型預測毒性效應數據相比實驗獲取也有一定的優勢,如擴充實驗數據、減少時間和物力消耗以及生物犧牲量[21-22],還可以對多種化學品的聯合作用進行分析[23]等。

定量構效關系(quantitative structure-activity relationship, QSAR)模型被廣泛應用于預測毒性效應。QSAR 是將一組化合物的某種性質或活性與這些化合物的化學成分或結構信息進行定量關聯的方法,可以用來預測化合物的毒性值、作用模式,篩選和排序化學品等[24-26],該方法通常與其他模型方法如機器學習耦合使用;其中機器學習在生態毒理學中得到了越來越多的應用,其一般原理是根據一定的規則將輸入變量與輸出變量之間的關系一般化,并用于預測未知的相似情況[27-28];機器學習方法可以更好地處理非線性問題,對于關系復雜或未知的輸入、輸出變量也有很好的適應性,且通常具有良好的精度,可以減少重復性試驗等[29-31]。 而 EDCs 繁殖毒性是慢性毒性的一種,急性毒性終點不適用于測量EDCs 的慢性繁殖毒性效應。 卵黃蛋白原(vitelloge-nin, VTG)、性腺指數(gonado-somatic index, GSI)、第二性征、血漿中的類固醇濃度和性腺組織病變被認為是用于評估EDCs 繁殖毒性終點的生物標志物,這些終點的變化需要長時間觀測,一般采用無觀察效應濃度(no observed effect concentration, NOEC)或最低可觀察效應濃度(lowest observed effect concentration, LOEC)指標表示[32],這就造成了EDCs 毒性數據較少,從而鮮見利用上述模型對EDCs 水生生物繁殖毒性進行預測[21]。

因此本文將首先對近年來應用機器學習方法預測化合物水生生物毒性效應的相關研究進展進行總結,并在搜集到的可靠數據的基礎上,利用QSAR建立用于預測EDCs 水生生物毒性效應的機器學習模型,從而為日后的化學品生態風險評估和檢測優先性等提供指導。

1 材料與方法(Materials and methods)

1.1 機器學習預測化學物質水生生物毒性文獻的檢索與評述

通過Web of Science 和中國知網數據庫對近年來國內外使用機器學習方法預測水生生物毒性文章進行檢索,采用的檢索詞如表1 所示。 對檢索到的文獻作如下分析:當前研究的主要目的;文獻中使用到的機器學習模型以及每種模型的使用頻率;對每項研究涉及的不同研究對象進行匯總,如化合物、歸屬于不同營養級的水生生物以及毒性終點等;另外還包括文獻內涉及到的研究手段與數據處理方法等。

表1 用于檢索使用機器學習預測內分泌干擾物水生生物毒性效應文獻的關鍵詞Table 1 Key words for searching papers that applied machine learning to predict the toxicity effects of endocrine disruptor chemicals on aquatic organisms

1.2 基于機器學習和QSAR 的EDCs 毒性預測的模型構建

1.2.1 數據獲取與預處理

參考文獻中描述的毒性數據篩選方法[33],在美國環境保護局(US EPA) ECOTOX 數據庫檢索了以NOEC、LOEC 等作為毒性終點,與黑頭軟口鰷(Pimephales promelas)繁殖毒性相關的數據。 若搜集所得數據集內的相同化學品在相同毒性終點上存在不同的數據點,則取幾何平均值;篩選后得到了83種不同化學品對黑頭軟口鰷的繁殖毒性數據,考慮到數據量的因素,未對化學品繼續篩選[34]。

分子描述符是一組將分子的不同屬性(如物理化學、拓撲和結構等)進行量化表示的數值[35-36]。 為了獲得分子描述符,首先需要收集不同化學物質對應的簡化分子線性輸入規范(simplified molecular input line entry specification, SMILES);SMILES 數據收集自PubChem 網站(https://pubchem.ncbi.nlm.nih.gov/);使用了 PaDEL-descriptor 軟件[37]的 python 接口用于計算分子描述符,該軟件可以根據SMILES為每種化合物計算出共1 875 種分子描述符。

在獲得的描述符數據集中,并不是所有的描述符對于模型構建都是必要的。 具體篩選方法如下。

(1) 一些化合物的某些分子描述符的計算值可能為空值或無窮值(體現在excel 或csv 文件中即為無數據和Inf/Infinity),這些數值無法被輸入至機器學習模型中用于訓練,由于數據集中化合物的數量較少,因此刪除了具有非法值的描述符[38]。

(2) 常數項或半常數項(該系列的80%及以上數值都相等)的描述符通常對模型的貢獻較小,因此采取方差過濾法并選取0.01 作為過濾界限[39-40]。

(3) 一些分子描述符之間具有線性相關性,若成對的描述符之間的Pearson 相關系數＞0.99,則只留下其中一個[34]。

(4) 經過上述篩選,大多數冗雜特征被去除,但仍需要選擇最優子集。 這個選擇過程被認為是比較困難的,因為沒有合適的規則作為指導,通常以個人經驗與其他算法相結合的方式進行[41-42]。 本文使用了遞歸特征消除(recursive feature elimination,RFE)[43],RFE 可以結合具有判斷變量重要性的機器學習算法,重復建模為特征的重要性進行排序并逐漸刪除指定個數特征,直到剩余規定數量的特征為止。 為了消除數據之間由于數量級差異帶來的影響,首先對所有描述符作了標準化,公式如下所示:

式中:Xi為第n個描述符的第i個數值,μn為第n個描述符的平均值,Sn為第n個描述符的標準差;然后使用結合隨機森林的RFE 法選擇最終特征子集。

了解化合物的可能毒性范圍有利于開展初步生態風險評估工作[44]。 根據中華人民共和國國家標準《化學品水生環境危害分類指導第3 部分:水生毒性》(GB/T 36700.3—2018),對于慢性毒性不大于100 μg·L-1的物質,認為其毒性較高,反之則認為其毒性較低;在此標準的指導下,選取了100 μg·L-1作為分類界限,NOEC 小于等于該值的化合物為類別“1”,大于該值的為類別“0”。 數據集被以4∶1 的比例劃分為訓練集和測試集,測試集用于模型的效果評價,不用于模型的訓練。

1.2.2 機器學習模型的構建

采用的支持向量機(support vector machine,SVM)模型與線性神經網絡(linear neural network,LNN)模型,分別由 scikit-learn[45]和 Keras 搭建。SVM 模型可以執行線性和非線性的分類與回歸任務,且被認為非常適用于中小型數據集[46],其中應用到的核函數為高斯徑向基(Gaussian radial basis function, RBF),該核函數常被應用于SVM 的構建中。LNN 模型中,每個神經元都代表一個多元線性函數,如下式所示。

式中:Y為該神經元的輸出值,X1~Xn為輸入特征,W1~Wn為權重,b為偏置值,采用了單隱藏層結構[47];Sigmoid 函數為激活函數,可以將輸出的數值范圍變為0 ~1,即“預測為正類”的概率值;二元交叉熵作為損失函數。

1.2.3 模型評估標準

在二元分類中,模型的預測性能根據真陽性(true positives, TP)、真陰性(true negatives, TN)、假陽性(false positives, FP)、假陰性(false negatives, FN)的數量以及敏感性(sensitivity, SE)、特異性(specificity,SP)和預測準確度(accuracy, Acc)來判定[44];此外還應用了受試者工作特征(receiver operating characteristic, ROC)曲線與曲線下面積(area under curve,AUC)來評價模型的分類性能;ROC 曲線的x軸為假陽性率(false positive rate),y軸為真陽性率(true positive rate);AUC 取值為 0.5 ~1.0,當 AUC=1.0 時表示這是一個完美的分類器,而AUC=0.5 時說明該分類器沒有分類能力[48-49]。 所涉及到的評價參數的含義和計算式如表2 所示。

表2 二元分類模型能力判定標準Table 2 Assessment standard of binary classification models

1.2.4 應用領域

經濟合作與發展組織關于QSAR 模型的指導文件[50]中指出,“一個(Q)SAR 模型需要定義其應用域(application domain, AD)”,即根據模型訓練集中化學物質的結構或物理化學等信息確定模型的預測能力限制范圍,對超出該范圍的化學物質(與訓練集中物質的相似性不足)的預測結果被認為可靠程度較低。 由于相似性有很多不同的表達方式(一般通過理化性質來定義),因此AD 的評估也可以是多樣化的,如杠桿方法[51]和基于Euclidean 距離的AD 分析法[52-54]。 其中Euclidean 方法將化學分子表示為多維向量中的一點(維數等于每種描述符中的變量數量),并以Euclidean 距離計算任意2個分子之間的相 似 性。 Ambit Discovery 軟 件 (http://ambit.sourceforge.net/download_ambitdiscovery.html)可以直接構建基于Euclidean 距離的AD 分析,并顯示處于AD 之外的化合物,因此 AD 分析將使用該軟件進行。

2 結果(Results)

2.1 機器學習預測化學物質毒性進展

根據檢索詞共篩選出英文文獻61 篇,中文文獻2 篇,發文數量與年份增長之間的關系如圖1 所示。由圖1 可知,結合機器學習方法來預測化合物對水生生物毒性的文章數量從2009年開始增多并且呈現明顯的上升趨勢,說明這種策略正得到越來越多的認可。 這一方面是由于機器學習方法所具備的優勢,另一方面也和計算機技術的發展為機器學習的應用提供了更優秀的條件有關[55]。

圖1 近年來使用機器學習或建模方法預測化學品水生生物毒性的文章數量和趨勢Fig.1 The number and trend of papers that used machine learning or modeling methods to predict the toxicity of chemicals on aquatic organisms in recent years

每種算法的使用次數與應用方式(用于預測離散、連續型數據,或者變量篩選)如圖2 所示。 其中,使用次數最多的是SVM,共25 次,且在回歸與分類問題上的使用較為均衡,一定程度上體現了其廣泛適用性[56-58];線性回歸的使用次數僅次于SVM,并與神經網絡一起更多地被應用于回歸問題;遺傳算法幾乎僅被用于輔助作用,即作為一種選擇描述符子集的手段,而不用于預測化合物的毒性效應;決策樹、隨機森林和k最近鄰等算法被較多地應用于分類問題[59-60]。

圖2 被用于預測化學品水生生物毒性的算法及其應用的頻率與目的Fig.2 Algorithms used to predict the toxicity of chemicals on aquatic organisms and their frequency and purpose of application

文獻中涉及的水生生物、化合物和毒性終點如圖3 所示。 涉及的水生生物包括脊椎生物、無脊椎生物和藻類,其中脊椎生物即魚類,如黑頭軟口鰷(Pimephales promelas)、斑馬魚(Brachydanio rerio)和虹鱒(Oncorhynchus mykiss)等;無脊椎生物中較多的是浮游生物,如梨形四膜蟲(Tetrahymena pyriformis)、大型溞(Daphnia magna)等。 所探究的化合物種類也較多:按照結構信息,有取代苯類化合物、芳香族化合物和酚類化合物等;根據作用,包含農藥(如生物殺滅劑、除草劑等)、個人護理產品(如抗抑郁藥、降壓藥和麻醉藥等)和工業化學品等。 根據危害方式,大多數文獻所研究的毒性終點為急性毒性,如半抑制生長濃度[61]、半致死濃度[62]和半數效應濃度[63]等,這可能與其實驗周期短、數據量較多、誤差較低以及當前管控優先度較高等因素有關。 而在慢性毒性當中,以 NOEC 作為毒性終點的研究較少[34,64],且模型的性能也相對較差,如Sheffield 和Judson 等[34]的研究中為該終點構建了回歸模型,評估回歸模型常用的標準之一是由實際值與預測值所計算出的決定系數(R2),在其研究中所構建的部分模型的R2為0.6 左右,盡管在QSAR 領域中R2＞0.5時模型即被認為具有預測性能[65],但相較于大多數其他學者的研究而言則處于較低水平[66-68]。

圖3 各文獻中使用到的水生物種與毒性終點注:IC50 表示半抑制濃度;IGC50 表示半抑制生長濃度;LC50 表示半數致死濃度;EC50 表示半數效應濃度;NOEC 表示無觀測效應濃度。Fig.3 Aquatic creatures and toxicity endpoints applied in papersNote: IC50 stands for 50% inhibitory concentration; IGC50 stands for 50% impairment growth concentration; LC50 stands for lethal concentration 50%;EC50 stands for concentration for 50% of maximal effect; NOEC stands for no observed effect concentration.

2.2 模型性能評估

2.2.1 描述符選擇及AD 評估

經過RFE 方法篩選,最終選擇了ATSC0m、ATSC7p、MATS3i 和 TpiPC 作為輸入變量。 其中 ATSC0m、ATSC7p 和 MATS3i 是 2D 自相關描述符,ATSC0m 和ATSC7p 分別為原子質量加權和原子極化率加權的 Broto-Moreau 中心自相關描述符,MATS3i 是電離勢加權的 Moran 中心自相關描述符,分別表征了原子質量、極化率與電離勢的影響;TpiPC 則與步進計數的常規鍵序 ID 號相關[69-70]。使用Ambit Discovery 構建的AD 部分表征如圖4所示,軟件計算結果顯示訓練集與測試集中均無化合物落在AD 之外,這說明選取的訓練集具有良好的代表性。

圖4 基于Euclidean 距離的應用域表征Fig.4 Application domain based on Euclidean distance

分子描述符的數值變化對毒性帶來的影響如圖5 所示,圖 5 中(a)、(b)、(c)和(d)分別為 ATSC0m、ATSC7p、MATS3i 和 TpiPC。 藍色柱狀條代表標準化后的每個化合物的分子描述符的數值;橙色柱狀條代表毒性,存在與否表示該化合物是否具有較高毒性。 可以看出,對于描述符ATSC0m 和TpiPC,隨著數值的增大,橙色柱狀條開始變得相對密集,即化合物傾向于具有高毒性;ATSC7p 則與之相反,隨著其數值增大,更多的化合物毒性較低;MATS3i 顯示出了不同的趨勢,其增大與減小時化合物毒性均較低,而在均值附近時較多的化合物具有較高毒性。

圖5 分子描述符數值大小與毒性之間的關系注:橫坐標表示不同的化合物,縱坐標表示標準化后的化合物毒性值。Fig.5 Relationship between molecular descriptors and toxicityNote: Abscissa represents different chemicals, and ordinate represents the toxicity of chemicals after standardization.

2.2.2 性能評估

數據集中化合物名稱、CAS 號和模型的預測結果如表3 所示,其中模型Ⅰ為SVM,模型Ⅱ為LNN。

續表3

續表3

訓練集和測試集的評估如表4 所示。 其中,SVM 在訓練集和測試集上的預測準確率分別為0.91 和0.88 左右,均達到了較好的水平,說明預測能力可以接受;模型對測試集的預測結果中,對高毒性與低毒性化合物的召回率,即SE 與SP 分別為1.00 與0.67,相比訓練集中的0.93 與0.88 來說不夠均衡,這可能是由于測試集數據量較少導致的,但是SE 較高可以減少實際有毒化合物漏檢的可能性;訓練集與測試集的預測準確率差距不大,說明模型沒有發生過擬合。 SVM 與 LNN 構建模型得到的ROC 曲線分別如圖6 和圖7 所示,其中SVM 的訓練集與測試集的AUC 分別為0.93 和0.88,遠大于下限0.5,因此這是一個較好的分類器。

圖6 由SVM 構建模型得到的訓練集與測試集受試者工作特征(receiver operating characteristic, ROC)曲線注:AUC 表示曲線下面積。Fig.6 Receiver operating characteristic (ROC) curve for training set and test set based on SVMNote: AUC stands for area under curve.

圖7 由線性神經網絡(linear neural network, LNN)構建模型得到的訓練集與測試集ROC 曲線Fig.7 ROC curve for training set and test set based on linear neural network (LNN)

表4 最終模型預測性能表征Table 4 Statistical results of developed models

LNN 在訓練集和測試集上的預測準確率均為0.82 左右,未出現過擬合現象,SE 分別為0.88 與1.00,SP 分別為0.73 與0.50;該模型的預測結果同樣有不均衡的SE 與SP 分布,可能進一步說明該問題的出現與數據集有關;訓練集與測試集的AUC分別為0.87 與0.88,說明分類性能良好。

2.2.3 模型對比

(1) SVM 比 LNN 穩定。 如圖 8 所示,保持超參數等條件不變,SVM 可以通過訓練得到恒定最優解;而對于LNN,若訓練次數不斷增加,結果也在逐漸發生變化,如圖9 所示,訓練集預測準確率(Acc)上升,測試集預測準確率(val_acc)不變,但測試集損失函數(val_loss)卻與訓練集損失函數(loss)呈現相反趨勢,說明模型傾向于朝過擬合發展,這可能與數據集較小有關。 SVM 的預測結果也略優于LNN,一定程度上說明SVM 較LNN 更適合于小數據集。

圖8 經過10 次相互獨立的訓練后SVM 的預測準確率Fig.8 The prediction accuracy of SVM after trained for ten times separately

圖9 LNN 的訓練過程中結果持續變化Fig.9 The result of LNN kept changing while training

(2) SVM 的訓練難度相對較低。 如上所述,隨著訓練的進行,SVM 可以得到恒定最優解,而LNN不能;另外,在相同的訓練次數內,LNN 的預測準確率也會呈現不同的變化趨勢或規律,結束訓練時得到的結果也可能不同,如圖10 所示。

圖10 相互獨立的LNN 訓練過程中出現不同結果Fig.10 Separate training process of LNN led to different results

(3) SVM 的訓練耗時相較于LNN 更短:SVM得到本實驗中最優解的訓練時間遠＜1 s,對LNN 每訓練1 000 輪則需要20 s 左右(具體耗時與進行訓練所使用的設備以及模型的超參數有關,此處僅針對本實驗條件作討論)。

3 討論(Discussion)

本文對機器學習模型方法在水生毒性預測領域的應用研究進行了概括與總結,并使用 SVM 與LNN 結合QSAR,使用較少被其他研究者采用的EDCs 繁殖毒性的NOEC 作為終點,在黑頭軟口鰷數據集上構建了預測毒性高低的二分類模型;SVM在該領域中的使用頻率最高;對急性毒性的研究多于慢性毒性;描述符子集的篩選是非常重要的步驟,結合了隨機森林的RFE 方法較好地篩選出了合適的描述符子集,篩選結果說明化合物對黑頭軟口鰷的繁殖毒性可能與分子質量、極化率、電離勢和相鄰原子成鍵強度有關;根據準確率與ROC 曲線等分類模型評定標準可知,本文中所構建的模型均具有可接受的預測能力,其中SVM 的預測能力和訓練表現等相較于LNN 更優,驗證了SVM 更適用于中小數據集。 本實驗中所使用的方法和構建的模型可為日后的AD 內未知化合物的檢測優先性起到指導作用,并且為水生生物毒性領域中對EDCs 的繁殖毒性的研究提供了一定的支撐。