楊浩藝,陳 微,姚澤歡，譚郁松，李 非

(1.國防科技大學(xué)計(jì)算機(jī)學(xué)院,湖南 長沙 410073；2.中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

1 引言

以侵襲性真菌感染為代表的感染性疾病是導(dǎo)致發(fā)病率高和死亡率高的重要原因，尤其對于艾滋病毒感染患者以及患有其他共病的免疫力低下的患者造成了重大威脅。這些真菌造成的死亡率與耐藥結(jié)核分岐桿菌造成的死亡率相當(dāng)，超過了瘧疾[1]。當(dāng)前已知的對人類有致病性的真菌多達(dá)300多種，根據(jù)抗真菌藥物的作用機(jī)制[2]，已用于臨床的抗真菌藥物大致可分為多烯類、三唑類、烯丙胺類、棘白菌素和其他抗真菌藥物。然而人類真菌感染的病理生理學(xué)研究仍然遠(yuǎn)落后于其他病原體引起的疾病[3]。此外，耐藥菌的出現(xiàn)和廣泛分布使得曾經(jīng)容易治愈的疾病變得再次致命[4]，以念珠菌為例，其對許多國家選擇的標(biāo)準(zhǔn)抗真菌藥物氟康唑以及新推出的棘孢菌素均具有耐藥性[5]。因此，安全有效的抗真菌藥物的研發(fā)顯得十分必要且迫切。

本文利用CMAP(Connectivity MAP)[6]和LINCS(Library of Integrated Network-based Cellular Signatures)[7]高通量轉(zhuǎn)錄組學(xué)數(shù)據(jù)，基于WTCS(WeighTed Connectivity Score)算法[8]，從已有抗真菌藥物出發(fā)，發(fā)現(xiàn)其他藥物潛在的抗真菌用途。本文通過將生物大數(shù)據(jù)應(yīng)用于快速藥物設(shè)計(jì)發(fā)現(xiàn)，基于生物醫(yī)藥數(shù)據(jù)特征構(gòu)建端到端的藥物預(yù)測策略，為加快新的抗真菌藥物的研發(fā)進(jìn)程提供計(jì)算方法。

2 數(shù)據(jù)驅(qū)動(dòng)的藥物重定位方法

根據(jù)Eroom定律[9]，新藥研發(fā)成本持續(xù)增長。面對從無到有的漫長傳統(tǒng)藥物研發(fā)過程，迅速積累的生物醫(yī)學(xué)高通量數(shù)據(jù)推動(dòng)了以生物信息學(xué)醫(yī)藥大數(shù)據(jù)為基礎(chǔ)的系統(tǒng)性藥物重定位[10]的發(fā)展，并引起研究人員廣泛關(guān)注。藥物重定位[11]針對已知的藥物識別和發(fā)現(xiàn)其新用途、新療效，是快速發(fā)現(xiàn)潛在藥物的不錯(cuò)選擇，既能夠高效地找到目標(biāo)藥物，也可提前預(yù)知藥物的副作用及用藥注意事項(xiàng)。以抗真菌藥物為例，在現(xiàn)有的、合成的或半合成的化合物庫中進(jìn)行抗真菌活性篩選，能夠提前確定這些篩選出的化合物的最低毒性[1]。

數(shù)據(jù)驅(qū)動(dòng)的藥物發(fā)現(xiàn)方法依賴高質(zhì)量的數(shù)據(jù)資源，GenBank數(shù)據(jù)庫[12]整合了來自所有可用公共來源的DNA序列；PharmGKB數(shù)據(jù)庫[13]為藥物研發(fā)提供了潛在的藥物-基因組關(guān)聯(lián)信息以及基因型-表型信息；CMAP和LINCS數(shù)據(jù)庫[6]提供了在不同細(xì)胞系中加入多種化合物所產(chǎn)生的基因表達(dá)譜；Drug Bank數(shù)據(jù)庫[14]作為一個(gè)獨(dú)特的生物信息學(xué)和化學(xué)信息學(xué)資源，結(jié)合了詳細(xì)的藥物/化學(xué)數(shù)據(jù)以及藥物靶標(biāo)/蛋白質(zhì)信息；蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)庫PDB(Protein Data Bank)[15]提供了目前最完整的蛋白質(zhì)三維結(jié)構(gòu)數(shù)據(jù)。

無論是傳統(tǒng)的藥物發(fā)現(xiàn)流程還是藥物重定位，其關(guān)鍵點(diǎn)在于確定化合物的作用模式MOA(Mode Of Action)及其非靶點(diǎn)效應(yīng)[16]?；诩?xì)胞轉(zhuǎn)錄反應(yīng)檢測藥物作用模式(MOA)進(jìn)而發(fā)現(xiàn)藥物的方法所需信息量最少，且可以快速應(yīng)用于新的化合物[11]。轉(zhuǎn)錄組數(shù)據(jù)能夠直觀反映在某一特定條件下基因表達(dá)、基因過表達(dá)或者基因沉默的情況，不同條件下轉(zhuǎn)錄組數(shù)據(jù)結(jié)果不盡相同[17]。針對藥物發(fā)現(xiàn)來說，一是可以通過比較各種化合物作用于細(xì)胞與正常條件下細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)差異，找到有效藥物；二是可以通過比較不同化合物在相同條件下作用于細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)，找到具有相同作用模式的化合物，進(jìn)而達(dá)到發(fā)現(xiàn)潛在治療藥物的目的。在此背景下，利用基因表達(dá)譜、轉(zhuǎn)錄組譜以及生成的關(guān)聯(lián)網(wǎng)絡(luò)進(jìn)行相關(guān)比對分析的方法以其快速高效、低成本的優(yōu)勢在臨床治療、藥物作用模式闡述、藥物重定位和系統(tǒng)生物學(xué)等多個(gè)研究領(lǐng)域得到了應(yīng)用和發(fā)展[12]。

Lamb等[6]創(chuàng)建了大型的藥物和基因標(biāo)簽公共數(shù)據(jù)庫CMAP和LINCS，使得通過多種模式匹配算法挖掘生物數(shù)據(jù)特征之間的關(guān)聯(lián)性成為可能[18]。他們采用L1000技術(shù)[19]基于大規(guī)模的統(tǒng)計(jì)分析辨識出人類細(xì)胞中978個(gè)基因作為全基因組的標(biāo)志基因，進(jìn)一步通過計(jì)算預(yù)測獲得全部基因的表達(dá)量[20],實(shí)現(xiàn)了低成本、高通量的實(shí)驗(yàn)數(shù)據(jù)獲取。到目前為止，LINCS計(jì)劃已獲得了77種典型細(xì)胞中的4 000多個(gè)沉默基因和7 000余種化學(xué)小分子刺激下的130余萬個(gè)全基因組表達(dá)譜[21]，為構(gòu)建不同藥物反應(yīng)之間的關(guān)聯(lián)關(guān)系奠定了牢固的數(shù)據(jù)基礎(chǔ)。

3 抗真菌藥物預(yù)測發(fā)現(xiàn)

本文針對抗真菌藥物進(jìn)行藥物預(yù)測發(fā)現(xiàn),以5種抗真菌分子化合物(氟胞嘧啶(flucytosine)、酮康唑(ketoconazole)、咪康唑(miconazole)、兩性霉素B(amphotericin B)和制霉菌素(nystatin))為基礎(chǔ)，基于CMAP數(shù)據(jù)庫對5種藥物與基因表達(dá)譜的聯(lián)通性分?jǐn)?shù)排序，針對每種藥物確定一組上調(diào)和下調(diào)基因，每組包括10個(gè)基因。運(yùn)用WTCS算法將查詢藥物的上、下調(diào)基因集與CMAP參考數(shù)據(jù)庫中的擾動(dòng)分子進(jìn)行富集分析，計(jì)算相似性分?jǐn)?shù)并對其進(jìn)行排序，得到每種目標(biāo)藥物的相似藥物列表。

考慮到選取的是不同種藥物作用機(jī)制下的代表性藥物，作用機(jī)制不同使得基因差異表達(dá)情況不同，綜合分析5種抗真菌的藥物特點(diǎn)及其相似藥物列表結(jié)果之后，本文選擇合并多個(gè)相似藥物列表，利用RankAggreg[22]中的交叉熵-蒙特卡洛算法，分別采用Spearman footrule distance和Kendall’s tau distance 2種距離函數(shù)對5種相似藥物列表進(jìn)行聚合，對比2種距離函數(shù)得到的聚合結(jié)果篩選出抗真菌藥物的預(yù)測藥物列表。

3.1 DPA藥物預(yù)測分析

基于Subramanian等[19]提出的WTCS計(jì)算相似性分?jǐn)?shù)方法，本文提出藥物預(yù)測分析DPA(Drug Prediction Analysis)方法，如圖1所示。首先，通過CMAP數(shù)據(jù)庫找到目標(biāo)藥物所對應(yīng)的上、下調(diào)基因集，經(jīng)過嘗試與計(jì)算，將每組上、下調(diào)基因集中的基因標(biāo)簽數(shù)量控制在10個(gè)；其次，通過計(jì)算目標(biāo)藥物與CMAP參考標(biāo)簽數(shù)據(jù)庫中每一個(gè)分子化合物擾動(dòng)的相似性分?jǐn)?shù)和差異顯著性值，將得到的藥物列表按照相似性分?jǐn)?shù)由高到低進(jìn)行排序，從而得到目標(biāo)藥物的相似藥物列表；最后，采用RankAggreg中的2種距離函數(shù)方法對5種抗真菌藥物的相似藥物列表進(jìn)行聚合，得到最可能的抗真菌藥物預(yù)測結(jié)果。為了驗(yàn)證該方法的正確性，本文利用Glaser等[23]公布的HDAC(Histone Deacetylase)抑制劑的基因表達(dá)標(biāo)簽，對所設(shè)計(jì)的藥物預(yù)測方法進(jìn)行正確性驗(yàn)證。

Figure 1 DPA drug discovery analysis process圖1 DPA藥物發(fā)現(xiàn)分析流程

3.2 WTCS相似性計(jì)算

經(jīng)典的基因集富集分析GSEA(Gene Set Enrichment Analysis)[21]方法以查詢目標(biāo)化合物分子的基因標(biāo)簽作為輸入，可以評估其與數(shù)據(jù)集中每個(gè)參考表達(dá)譜的相似性。給定所需要查詢計(jì)算的目標(biāo)化合物分子的基因標(biāo)簽(上調(diào)基因、下調(diào)基因)，將目標(biāo)化合物分子的基因標(biāo)簽與CMAP數(shù)據(jù)庫中的編目列表進(jìn)行比較分析，根據(jù)上調(diào)基因、下調(diào)基因在排序列表中的分布情況，可以將目標(biāo)化合物的基因標(biāo)簽與數(shù)據(jù)庫中的基因標(biāo)簽的關(guān)系分為正相關(guān)、負(fù)相關(guān)和無關(guān)3種。而正、負(fù)相關(guān)又可細(xì)分為強(qiáng)正(負(fù))相關(guān)和弱正(負(fù))相關(guān)。比對后可以得到目標(biāo)化合物分子基因標(biāo)簽與數(shù)據(jù)庫中化合物分子基因標(biāo)簽的聯(lián)通性分?jǐn)?shù)(Connectivity Score)，相似性分?jǐn)?shù)的取值在-1～1。

WTCS算法對傳統(tǒng)的GSEA富集分析方法進(jìn)行了改進(jìn)，通過計(jì)算不同化合物分子基因標(biāo)簽的富集分?jǐn)?shù)得到不同化合物的聯(lián)通性分?jǐn)?shù)，并通過聯(lián)通性分?jǐn)?shù)的高低找到與目標(biāo)藥物作用相似的藥物，進(jìn)而達(dá)到重定位藥物的目的。該算法原理相對簡單，易于操作實(shí)現(xiàn)，在尋找相似化合物分子計(jì)算中已有廣泛的應(yīng)用。

WTCS算法基于Kolmogorov-Smirnov富集統(tǒng)計(jì)ES(Enrichment Score)的非參數(shù)性相似性度量，對于輸入基因集(qup，qdown)，按照式(1)的計(jì)算方式得到與某一參考基因標(biāo)簽的相似性分?jǐn)?shù)Wq,r：

(1)

其中，ESup、ESdown分別是qup、qdown在參考基因標(biāo)簽下的富集分?jǐn)?shù)。

3.3 Rank aggregation加權(quán)聚合排序

Rank aggregation是一種對多個(gè)排序列表進(jìn)行整合得到一個(gè)綜合排序列表的算法[22]。在該算法中，Spearman footrule distance距離函數(shù)根據(jù)不同排序列表內(nèi)元素的排序位置進(jìn)行距離計(jì)算，該方法簡單且所需信息量少；而Kendall’s tau distance距離函數(shù)需要聯(lián)合不同排序列表中對應(yīng)的元素對距離進(jìn)行綜合計(jì)算，該方法復(fù)雜但最終得到的結(jié)果列表排序等級差異明顯。Rank aggregation算法在R語言中有現(xiàn)成的包RankAggreg可用，因此實(shí)驗(yàn)過程中只需要直接輸入需要聚合的排序列表，并調(diào)整迭代次數(shù)、距離函數(shù)等參數(shù)信息，使得結(jié)果收斂到最小。

4 結(jié)果

4.1 數(shù)據(jù)和方法

首先利用R語言下載CMAP數(shù)據(jù)集，該數(shù)據(jù)集中包含多種化合物作用下不同細(xì)胞系細(xì)胞的基因表達(dá)譜。這些化合物可以是蛋白、小分子化合物或者是復(fù)雜化合物。通過查詢氟胞嘧啶、酮康唑、咪康唑、兩性霉素 B和制霉菌素的ID號，可以確定在這些化合物作用下細(xì)胞的基因表達(dá)情況，根據(jù)不同基因表達(dá)結(jié)果的排序篩選得到5種化合物的基因標(biāo)簽集。本文的藥物相似性計(jì)算在整個(gè)CMAP數(shù)據(jù)集上針對所有化合物進(jìn)行相似性計(jì)算分析。GSEA和WTCS算法的計(jì)算流程由Bioconductor[24]中的piano包(使用各種統(tǒng)計(jì)方法從不同的基因統(tǒng)計(jì)水平和廣泛的基因集合進(jìn)行基因集分析)實(shí)現(xiàn)，通過分別對上調(diào)基因和下調(diào)基因進(jìn)行富集分析，可以得到所要查詢的基因標(biāo)簽與參考數(shù)據(jù)庫中的基因標(biāo)簽的聯(lián)通性分?jǐn)?shù)，以及差異顯著性值(P值)。藥物發(fā)現(xiàn)分析流程利用基于R語言實(shí)現(xiàn)的Bioconductor開源環(huán)境下的PharmacoGx[25]、piano、bioMaRt[26]和RankAggreg[22]等R語言包實(shí)現(xiàn)。

4.2 藥物相似性分析結(jié)果

在針對5種抗真菌藥物分別計(jì)算得到的排名前10的相似藥物列表結(jié)果(表1)中，酮康唑和咪康唑的相似藥物列表結(jié)果相似程度較高，酮康唑相似藥物列表中相似性分?jǐn)?shù)最高的可達(dá)0.82，而咪康唑相似藥物列表中相似性分?jǐn)?shù)最高可達(dá)0.88，這2種藥物的相似藥物列表在后續(xù)的實(shí)驗(yàn)中值得關(guān)注；與酮康唑和咪康唑的相似藥物列表相比，其余3種藥物的相似藥物結(jié)果相似程度相對較低，但前10位最相似藥物的相似性分?jǐn)?shù)均在0.6以上。從實(shí)驗(yàn)結(jié)果的準(zhǔn)確性上看，上述藥物計(jì)算結(jié)果中，酮康唑和咪康唑的相似藥物結(jié)果不僅在相似性分?jǐn)?shù)上表現(xiàn)較好，其結(jié)果的P值也幾乎全在0.01以內(nèi)；而兩性霉素B、制霉菌素和氟胞嘧啶的計(jì)算結(jié)果雖然在相似性分?jǐn)?shù)和差異顯著性值上表現(xiàn)相對較差，大部分計(jì)算結(jié)果的P值均大于0.01，但除去極個(gè)別藥物外，這些相似藥物結(jié)果的P值均小于0.05，仍然具有一定的參考價(jià)值。5種目標(biāo)藥物計(jì)算得到的P值結(jié)果如表2所示。從這些相似藥物結(jié)果也可以看出,由于不同的藥物作用機(jī)制不同，不同藥物的相似藥物以及其相似程度也存在明顯差異。

Table 1 List of inferred antifungal candidates表1 預(yù)測抗真菌候選藥物列表

Table 2 P-value list of inferred antifungal candidates表2 預(yù)測抗真菌候選藥物P值列表

4.3 抗真菌藥物預(yù)測結(jié)果

在5種藥物的相似藥物結(jié)果的基礎(chǔ)上，本文利用RankAggreg中的交叉熵-蒙特卡洛算法，分別采用Spearman footrule distance和Kendall’s tau distance 2種距離函數(shù)對5種藥物的相似藥物列表結(jié)果進(jìn)行聚合排序。在4.2節(jié)得到的相似藥物排序結(jié)果中，根據(jù)化合物分子與目標(biāo)藥物的正負(fù)相關(guān)性最終得到的相似性分?jǐn)?shù)存在0～1,0,-1～0這3種情況。由于本實(shí)驗(yàn)中只考慮正相關(guān)的情況，對于每種藥物的相似藥物列表只選擇相似性分?jǐn)?shù)大于或等于0的化合物分子進(jìn)行聚合排序，在實(shí)驗(yàn)過程中不斷調(diào)整迭代次數(shù)、分位數(shù)和生成樣本數(shù)量等參數(shù)以獲得更加精確的收斂結(jié)果。由于計(jì)算方法的差異，采用Spearman footrule distance距離函數(shù)和Kendall’s tau distance距離函數(shù)最終收斂得到的結(jié)果數(shù)值水平存在一定差異(如圖2所示)，但利用2種距離函數(shù)聚合得到的藥物列表排序結(jié)果相似。通過排除掉5種目標(biāo)查詢藥物，最終選擇了10種化合物分子作為抗真菌藥物預(yù)測藥物結(jié)果(如表3所示)，分別是伊利替康(irinotecan)、氯硝柳胺(niclosamide)、舍他康唑(sertaconazole)、普尼拉明(prenylamine)、銀杏內(nèi)酯A(ginkgolide A)、莫索尼啶(moxonidine)、槲皮素(quercetin)、舒洛地爾(suloctidil)、卡米達(dá)佐(calmidazolium)和STOCK1N-35874(詳細(xì)計(jì)算分析結(jié)果參見https://github.com/yeaouh/antifungal)。

Table 3 Prediction results of antifungal candidates 表3 抗真菌候選藥物預(yù)測結(jié)果

Figure 2 Convergence results of RankAggreg圖2 RankAggreg收斂結(jié)果

5 結(jié)束語

本文基于WTCS、GSEA算法以及所設(shè)計(jì)的藥物發(fā)現(xiàn)分析方法對5種具有代表性的抗真菌藥物進(jìn)行了相似藥物預(yù)測，并通過綜合排序得到了抗真菌藥物最相似藥物列表，最高相似度為88%，得到的抗真菌相似藥物列表具有一定的參考價(jià)值。在最終計(jì)算得到的抗真菌藥物預(yù)測結(jié)果中，經(jīng)過文獻(xiàn)查證，已經(jīng)證明舍他康唑[27]、舒洛地爾[28]能夠應(yīng)用于臨床抗真菌藥物治療。另外，根據(jù)文獻(xiàn)顯示，一些抗氧化(如槲皮素)、抗癌(如伊利替康、槲皮素)等藥物的作用機(jī)理可能同樣適用于真菌[1]，從而達(dá)到治療真菌感染的目的，但這些猜測顯然需要實(shí)驗(yàn)研究進(jìn)行確認(rèn)。

利用高通量組學(xué)數(shù)據(jù)對現(xiàn)有藥物進(jìn)行篩選分析能夠大大縮短藥物發(fā)現(xiàn)流程。通過計(jì)算得到的理論結(jié)果有待進(jìn)一步的細(xì)胞、動(dòng)物及臨床實(shí)驗(yàn)驗(yàn)證。下一步將基于更大規(guī)模的數(shù)據(jù)集進(jìn)行探索，尋找使得預(yù)測結(jié)果更加可靠的方法，以及探索深度學(xué)習(xí)技術(shù)在組學(xué)數(shù)據(jù)分析處理、藥物設(shè)計(jì)發(fā)現(xiàn)中的應(yīng)用，實(shí)現(xiàn)更為精準(zhǔn)的抗真菌藥物預(yù)測發(fā)現(xiàn)。