基于通路擾動相似性建模的抗類風濕關節炎藥對分子協同組合發現及機制分析＊

孫繼佳，王睿瑞，張磊，劉保成＊＊，呂愛平，3，4＊＊

（1.上海中醫藥大學中醫健康協同創新中心上海201203；2.上海中醫藥大學數學與物理教研室上海201203；3.上海中醫藥大學中西醫結合研究院上海201203；4.香港浸會大學中醫藥學院中國香港999077）

類風濕關節炎（Rheumatoid Arthritis，RA）是一種發生在滑膜關節及其他器官系統的慢性、全身性、炎癥性疾病，也是一種慢性進行性自身免疫性疾病[1]。近年來隨著分子生物學、免疫學、遺傳學的研究發展，顯示RA的發病機制與遺傳因素、細菌和病毒因素、T淋巴細胞、B淋巴細胞、細胞因子以及其他免疫細胞等的改變和影響有關[2]。RA的發病率隨著年齡的增長而增加，全世界每年約有0.3%-1.0%的人受其影響[3]，以35-50歲成年人多見，且女性發病率為男性的3倍左右[4]。目前，臨床上治療RA的藥物仍以西藥為主，但往往存在毒副作用較大[5-7]、治療費用較高等特點，因此進一步提高抗RA藥物的有效性、安全性以及降低治療成本，是當前亟需解決的問題。

中醫藥治療類風濕關節炎歷史悠久，積累了豐富的臨床應用經驗，RA屬于中醫“痹證”的范疇，具有痹癥共有的病因[8-9]?！端貑枴け宰C論》曰：“風寒濕三氣雜至，合而為痹也”，指出風寒濕是引起痹證的主要因素。RA中醫辨證分型主要有脾腎虧虛、風寒濕痹、痰濕阻絡、濕瘀痹阻等4種證型[10-11]，并認為肝腎不足、氣血虧虛是RA發生的內在因素，而風寒濕是RA發生的誘發因素。目前，中藥及復方在治療RA的過程中，許多有效成分之間存在的組合機制尚不完全清楚，有效成分之間的協同性很少能夠從分子調控網絡角度進行更為精準地進行定量描述；因此，需要建立更為適合的研究方法和技術，發現治療中起決定作用的主要成分，同時能夠篩選出具有協同治療作用的有效成分組合，以實現對中藥及復方成分的進一步“優化”，進而為中藥新組分開發提供理論和方法學上的支持。

藥對是中醫臨床遣藥組方常用的配伍形式，是歷代醫家逐漸積累經驗而來的精妙應用形式，藥對通過協同增效、相制減毒、相反相成等形式應用，介于中藥和方劑之間，是一個值得研究的群體[12-14]。那么，如何從分子機制角度深入探索藥對中不同成分之間的協同組合效應并進一步分析其作用機制，是亟需解決的實際問題。鑒于此，本研究從治療RA的典型藥對出發，綜合運用生物信息學、化學信息學、網絡藥理學等理論和技術，并創新性地提出一種基于通路譜擾動相似性的計算方法，闡明不同藥對中有效成分在抗RA中的協同組合特征；同時，進一步挖掘和分析藥對中有效成分作用的分子機制及功能，旨在為中藥藥對作用分子機理闡明及新組分開發提供有力的方法學上的支持。

1 材料與方法

1.1 類風濕關節炎（RA）差異基因篩選

在GEO（Gene Expression Omnibus，https://www.ncbi.nlm.nih.gov/geo/）[15]數 據庫中，輸 入“Rheumatoid Arthritis”查詢和收集類風濕關節炎相關的基因表達譜芯片，獲得GSE55235、GSE55457和GSE77298這3個與RA有關的表達譜數據；其中，GSE55235的平臺為GPL96[HG-U133A]Affymetrix Human Genome U133A Array，包含10個正常樣本與10個RA樣本；GSE55457的平臺為GPL96[HG-U133A]Affymetrix Human Genome U133A Array，包含10個正常樣本和13個RA樣本；GSE77298的平臺為GPL570[HG-U133_Plus_2]Affymetrix Human Genome U133 Plus 2.0 Array，包含7個健康樣本和16個RA樣本。

利用R軟件GEOquery包下載所有3個表達譜數據芯片，并對原始數據進行標準化、校正和基因名注釋等預處理；使用limma包對GSE55235、GSE55457和GSE77298分別進行差異表達基因分析（differentially expressed genes，DEGs），并根據|logFC|>1.0及校正adj.p.value<0.05篩選出每組芯片數據中的上調和下調差異表達基因。

1.2 類風濕關節炎（RA）治療靶點數據收集

我 們 依 次 從TTD（http://db.idrblab.net/ttd/）[16]、DrugBank（https://www.drugbank.ca/）[17]數據庫收集與類風濕關節炎相關的靶點數據；在TTD中輸入“Rheumatoid Arthritis”輸入查找找出與RA相關的基因；在DrugBank中也同樣通過輸入“Rheumatoid Arthritis”搜索與治療RA相關的已認證藥物及其作用的靶點。

再將基于GEO的差異基因分析篩選后的結果和通過數據庫查詢得到的靶點基因進行合并、去除重復；通過UniProt（https://www.uniprot.org/）[18]數據庫對所有收集到靶點信息進行確認，組成一個RA疾病相關靶點基因集合SRA。

1.3 抗類風濕關節炎（RA）典型藥對數據收集

根據課題組以往的文獻整理和分析，選擇了治療類風濕關節炎的兩種較為常見的藥對，即：肝腎不足型的典型藥對，獨活-桑寄生、濕熱痹阻型的典型藥對，黃柏-蒼術，以此進行下一步研究和分析。

從TCMID 2.0（http://www.megabionet.org/tcmid/）[19]、ETCM數 據 庫（http://www.nrc.ac.cn:9090/ETCM/index.php/Home/Index/）[20]分別收集到與這4種中藥相關的化學成分和靶點數據，所有的中藥成分化學信息以及Canonical SMILES通過PubChem（https://pubchem.ncbi.nlm.nih.gov/）獲得。

1.4 藥對成分潛在靶點預測與識別

藥對成分的靶點數據主要由SEA（http://sea.bkslab.org/）[21]、HitPick（http://mips.helmholtz-muenchen.de/proj/hitpick）[22]以及ETCM數據庫來確定。根據藥對成分所對應的Canonical SMILES，通過SEA和HitPick在線系統對抗RA藥對有效成分進行潛在靶點預測，其中，選取SEA數據庫中評分值Max Tc>0.6，HitPick在線預測系統Precision值>50%的預測結果作為藥對成分的潛在作用靶點；ETCM中根據潛在候選靶點（Candidate Target Genes）評分值大于0.9，篩選出相應成分的靶點，根據靶點基因集SRA，使用Venny 2.1（https://bioinfogp.cnb.csic.es/tools/venny/index.html）工具進行取交集，識別出藥對小分子潛在的抗RA靶點。

1.5 藥物小分子結構相似性計算

一般認為，如果兩個藥物小分子A和B的化學結構越相似，則具有協同組合治療的可能性就越大；因此，我們先通過DRAGON 7（http://www.talete.mi.it/index.htm）軟件計算出藥物小分子1024維的分子指紋，再利用Tanimoto相似性系數（Tanimoto Coefficientbased Similarity）計算藥物小分子之間的結構相似性即：

其中，FP(A)和FP(B)分別為藥物小分子A和B的分子指紋，SFPd(A,B)為分子結構相似度。

1.6 基于通路擾動相似性的中藥小分子協同組合機制發現

如果兩種藥物所作用的靶點屬于同一條信號通路，那么這兩種藥物分子協同組合的可能性就越高；根據現有方法，可以將每個藥物小分子所作用的靶點富集到所對應的信號通路上，構建藥物作用靶點的通路向量，以此來計算不同藥物分子作用通路的協同值。

雖然，這種評估方法相對簡單，但沒有考慮到小分子對每個通路的擾動程度；因此，本文我們先利用clusterProfiler包[23]將藥物小分子作用的靶點進行KEGG通路富集分析，按p.value<0.05和q.value<0.05得到相關通路；根據注釋得到的相應通路集{TPk∣k=1,2,…,n}（n為通路個數）構建一個由包含所有注釋通路所構成的向量，即：通路向量TP=(tp1,tp2,…,tpn)，再根據以下流程（如圖1所示）構建協同組合評價模型：

圖1 基于通路擾動相似性方法的小分子協同機制發現方法示意圖

如果藥物小分子A作用的靶點集T(A)中的靶點屬于某一條通路TPk，那么通路譜TP對應位置上對這條通路的藥物擾動強度值為tpk(A)，否則為0，其中，tpk由下式計算得到，即：其中，l為通路TPk的靶點個數，m為靶點集T(A)的個數，dij為靶點之間在RA疾病相關PPI網絡上的最短路徑距離；由此得到藥物小分子A的靶點通路擾動向量TP(A)；同理，得到藥物小分子B作用的靶點通路擾動向量TP(B)。

根據所得到的向量TP(A)和TP(B)，就可以計算兩種藥物分子之間的通路譜擾動相似性評分STPd(A,B)，另外，為防止出現分母為0的情況，利用改進余弦距離公式進行計算，即：

1.7 GO功能注釋與KEGG Pathway富集分析

基因本體（Gene Ontology，GO）數據庫是GO組織（GO Consortium）在2000年構建的一個結構化的標準生物模型，涵蓋了基因的生物過程（Biological Process，BP）、分子功能（Molecular Function，MF）和細胞分組（Cellular Component，CC）。一個生物過程或通路通常是由一組基因共同參與，而不是由單個基因完成的，富集分析的主要依據是：如果一個生物學過程或通路在已知的研究中發生異常，則共同發揮功能的基因極可能被選擇出來作為與這一過程或通路相關的基因集合。

本文，我們進一步利用clusterProfiler工具包對所有藥對的潛在作用RA疾病的靶點集進行GO功能和KEGG Pathway富集分析，按照p.value<0.05和q.value<0.05進行篩選。

2 結果

2.1 類風濕關節炎疾病基因收集及PPI網絡構建

經過分析，GSE55235中篩選得到1071個差異基因，其中，上調基因599個、下調基因472個；GSE55457中篩選得到312個差異基因，其中，上調基因189個、下調基因123個；GSE77298中篩選得到432個差異基因，其中，上調基因237個、下調基因195個。另外，將只要同時出現在任意兩個芯片表達數據中的差異基因都列入RA疾病基因集，由此，共獲得與RA疾病相關的267個差異基因（圖2）。

從DrugBank數據庫中收集到與RA治療相關的73個已證實的藥物及其對應的193個靶點基因；又從TTD數據庫中收集得到139個與RA疾病相關的靶點基因；再將兩個數據庫所收集到的靶點集進行合并去除重復后，共獲得309個與RA相關的靶點基因。

最后，再將上述兩種方法（基于差異基因分析和靶點數據庫收集）獲得的所有基因進行合并后，總共得到555個與RA相關的靶點基因，即靶點基因集SRA；進一步再將SRA所有555個靶點基因導入STRING（https://string-db.org/）數據庫得到蛋白相互作用關系，其中，參數設置：Organism設為Homo Sapiens，Combine Score閾值設為0.7；同時，采用Gephi0.9.2軟件構建靶點蛋白的PPI網絡（圖2）。

圖2 基于差異基因分析和藥物靶點數據庫的抗RA靶點基因收集及PPI網絡構建

2.2 藥對數據收集及潛在靶點預測與識別

從TCMID、ETCM數據庫以及PubChem數據庫分別收集和整理獲得獨活、桑寄生、黃柏、蒼術、等4種中藥成分及其相關的化學成分數據；其中，獨活148個、桑寄生16個、黃柏62個、蒼術95個；說明：為了便于后續分析，我們還對獨活-桑寄生藥對164個成分進行一一編號DS1-DS164，黃柏-蒼術藥對157個成分也依次編號為HC1-HC157。

利用SEA、HitPick以及ETCM，經過靶點預測和閾值篩選得到相關成分的潛在作用靶點，進一步，通過靶點基因集SRA識別出抗RA有效成分及其潛在作用靶點；其中，獨活-桑寄生藥對有53個成分及其作用的52個潛在RA靶點，黃柏-蒼術藥對有32個成分及其作用的28個潛在RA靶點，合并得到67個潛在作用靶點信息（表1）；不同藥對的成分-靶點作用網絡（圖3）。

圖3 （B）黃柏-蒼術藥對成分-作用靶點網絡

圖3 （A）獨活-桑寄生藥對成分-作用靶點網絡

表1 兩種藥對所有67個潛在作用靶點信息

2.3 小分子協同組合評分結果

利用1.5所述，分別對兩種藥對中的小分子協同組合進行計算分析，根據協同分值Sd(A,B)值的大小來判斷兩種小分子是否存在可能協同作用及其強度；Sd(A,B)值越大，則對這對小分子對PPI網絡的擾動程度就越大、協同組合治療特征也越強；本文中，我們規定協同評分值Sd(A,B)≥0.15的小分子組合被認為是具有潛在的協同治療作用，兩種藥對中具有協同治療特性的小分子組合結果，如圖4所示，所有具有協同治療作用的中藥小分子信息（表2）。

表2 具有協同組合治療的2種藥對小分子信息

圖4 （B）黃柏-蒼術藥對協同小分子組合

圖4 （A）獨活-桑寄生藥對協同小分子組合

在獨活-桑寄生藥對中，我們發現存在15對具有通路協同治療特征的中藥小分子組合，其中，獨活中有4個小分子，分別為：DS51、DS75、DS120、DS137，桑寄生中有11個小分子，分別為：DS150、DS151、DS152、DS155、DS156、DS157、DS159、DS160、DS161、DS163、DS164。

在黃柏-蒼術藥對中，我們發現存在17對具有通路協同治療特征的中藥小分子組合，其中，黃柏中有14個小分子，分別為：HC9、HC16、HC23、HC24、HC33、HC39、HC60、HC9、HC16、HC23、HC24、HC33、HC39、HC60，蒼術中則有8個成分，分別是：HC66、HC70、HC71、HC126、HC66、HC70、HC71、HC126。

2.4 GO功能注釋與KEGG Pathway富集分析結果

利用clusterProfiler包進一步對兩種藥對的作用靶點集進行GO功能注釋和KEGG Pathway富集分析，其中，p.value<0.05和q.value<0.05，兩種藥對作用靶點的功能注釋與富集分析結果（圖5）。

圖5 （B）黃柏-蒼術藥對的作用靶點GO功能與KEGG通路富集分析結果

圖5 （A）獨活-桑寄生藥對的作用靶點GO功能與KEGG通路富集分析結果

獨活-桑寄生藥對作用靶點的GO功能注釋結果顯示，所涉及的BP功能有526個，主要包括：有機陰離子運輸（GO:0015711）、類固醇激素反應（GO:0048545）、調節炎癥反應（GO:0050727）、營養水平反應（GO:0031667）、有機酸運輸（GO:0015849）等；所涉及的MF功能有74個，主要包括：跨膜轉運蛋白活性（GO:0022804）、核受體活性（GO:0004879）、轉錄因子活性（GO:0098531）、類固醇激素受體活性（GO:0003707）、一元羧酸結合性（GO:0033293）等；所涉及的CC功能有19個，主要包括：轉錄因子復合體（GO:0005667）、神經細胞體（GO:0043025）、RNA聚合酶II轉錄因子復合物（GO:0090575）、核轉錄因子復合體（GO:0044798）、膜筏（GO:0045121）等。KEGG通路富集分析顯示，獨活-桑寄生藥對作用靶點主要富集在14個相關信號通路，分別為：膽汁分泌通路（hsa04976）、乙肝通路（hsa05161）、ABC轉運通路（hsa02010）、花生四烯酸代謝通路（hsa00590）、PPAR信號通路（hsa03320）、利什曼 病 通 路（hsa05140）、PD-L1和PD-1通 路（hsa05235）、GnRH信號通路（hsa04912）、IL-17信號通路（hsa04657）、胰腺分泌通路（hsa04972）、甲狀旁腺激素通路（hsa04928）、Th17細胞分化信號通路（hsa04659）、葉酸生物合成通路（hsa00790）和抗葉酸通路（hsa01523）等。

續表

黃柏-蒼術藥對作用靶點的GO功能注釋結果顯示，所涉及的BP功能有373個，主要包括：細胞對外部刺激反應（GO:0071496）、化學突觸傳遞的調節（GO:0050804）、突觸信號轉導調控（GO:0099177）、有機陰離子轉運（GO:0015711）、多細胞生物體內平衡（GO:0048871）等；所涉及的MF功能有36個，主要包括：嘌呤 能 受 體 活 性（GO:0035586）、激 素 結 合（GO:0042562）、ATPase偶聯活性（GO:0042623）、跨膜轉運蛋白活性（GO:0022804）、ATP酶活性（GO:0016887）等；所涉及的CC功能有27個，主要包括：神經細胞體（GO:0043025）、膜筏（GO:0045121）、膜微區（GO:0098857）、膜區（GO:0098589）、樹突棘（GO:0043197）等。KEGG通路富集分析顯示，黃柏-蒼術藥對作用靶點主要富集在10個相關信號通路，分別為：神經活性配體-受體相互作用通路（hsa04080）、膽汁分泌通路（hsa04976）、鞘脂信號通路（hsa04071）、cAMP信號通路（hsa04024）、cGMP-PKG信 號 通路（hsa04022）、VEGF信 號 通 路（hsa04370）、腎 細 胞 癌 通 路（hsa05211）、胰腺癌通路（hsa05212）、氮素代謝通路（hsa00910）、近端小管碳酸氫鹽回收通路（hsa04964）等。

3 討論

近年來，許多學者從不同角度，尤其是基于對文獻的整理[24]、統計分析[25]、關聯規則挖掘[26]的方法來總結和歸納治療RA的核心藥對組合，發現川烏-桂枝、當歸-川芎、當歸-黃芪、黃柏-蒼術、獨活-桑寄生等藥對都具有較高使用頻率的用藥規律；此外，也有研究者利用網絡藥理學[27]、實驗手段[28,29]來分析和闡述抗RA藥對的藥理作用機制，進一步揭示了藥對配伍的特點。

藥對小分子協同組合評價結果顯示，獨活-桑寄生藥對中的DS51（4-Methoxyacetophenone）與DS163（Quercitrin-7-olate）具有最高的協同性評分值，且DS51與桑寄生中的其他7個成分（DS160、DS157、DS159、DS152、DS161、DS156、DS164）等都存在一定的通路協同作用機制。王友慶[30]等人研究證實，槲皮素（DS150，Quercetin）可通過抑制NF-κB激活，減弱炎癥環境下軟骨細胞內MMP-13產生、基質降解和細胞凋亡，起到保護軟骨的作用，而通過計算結果顯示，槲皮素可能與DS120這個小分子化合物具有協同治療作用。鄧莉[31]等人研究發現，TNF-α能促進成纖維樣滑膜細胞IL-6及IL-1β的表達和分泌，并促進p38MAPK的磷酸化及NF-κB核轉位，而齊墩果酸（DS152，Oleanolic acid）呈濃度依賴性抑制TNF-α誘導的HFLS炎癥因子的表達（IL-6及IL-1β），并抑制TNF-α誘導p38MAPK的磷酸化及NF-κB核轉位，而計算結果也提示4-Methoxyacetophenone可能在藥對中與齊墩果酸具有組合治療效應。黃柏-蒼術藥對中的H24（28-Isofucostanol）與HC126化合物具有最高的協同評分，蒼術中的HC71（Icariside F2）與黃柏中的其他5個成分（HC33、HC16、HC9、HC39、HC60）均存在潛在的協同組合治療作用。研究顯示[32]，鹽酸藥根堿（HC9，Jatrorrhizine）能夠抑制體外滑膜細胞的增殖，遷移和分泌，并改善類風濕關節炎的大鼠模型；還有研究顯示[33]，Jatrorrhizine能抑制RANKL誘導的破骨細胞生成并防止磨損顆粒誘導的溶骨，已經具有發展成為治療RA的新型治療劑的巨大潛力。綠原酸（HC23，Chlorogenic acid）能夠通過下調核因子κB配體誘導的活化T細胞c1表達的受體激活劑來抑制破骨細胞分化和骨吸收，可能是破骨細胞相關疾病伴發炎性骨破壞的潛在治療選擇[34,35]；另有研究也證明[36]，綠原酸和木犀草素（Luteolin）通過調節NF-κB和JAK/STAT信號通路的活化，協同抑制白介素1β誘導的成纖維樣滑膜細胞的增殖。桑寄生中的兒茶素（HC66，Epicatechin）[37]具有抗氧化，抗炎，抗微生物，抗過敏，抗癌，抗血栓形成和保肝等治療作用。由此，可以發現，通過我們所提出的計算和分析方法，不僅可以找藥對中一些抗RA的主要活性成分，同時，也能夠發現藥對中的哪些其他成分與這些主要活性成分產生了“1+1>2”作用，從分子和定量角度揭示了中藥藥對配伍的特征。

從藥對所作用靶點富集的通路來看，獨活-桑寄生作用靶點所富集通路中，如已有研究報道IL-17能夠調節類風濕關節炎中SHP-2和IL-17RA/STAT-3所依賴的Cyr61、IL-23和GM-CSF的表達以及RANKL介導的成纖維細胞樣滑膜破骨細胞的生成，并可能揭示RA的發病機理的新機制[38]。已有報道還顯示[39]，可以通過調節成骨細胞中花生四烯酸的代謝以對類風濕關節炎的起到治療作用。黃柏-蒼術藥對中，例如，VEGF信號通路，目前已經明確報道與RA疾病多種發病分子機制有關[40,41]；還有研究顯示[42]，可以通過調節MMP14和VEGF介導的間充質細胞表達miRNA-150-5p外來體來治療類風濕關節炎。

此外，我們還發現獨活-桑寄生藥對主要涉及14個相關信號通路，黃柏-蒼術藥對則涉及10個相關信號通路，而兩種藥對中只有一個信號通路Bile secretion（hsa04976）是相同的，提示中醫在RA治療中針對不同證型采用不同的中藥配伍，因此，涉及到對不同信號通路的作用調控，從一定程度上反映了中醫“同病異治”的特征。

綜上所述，本文我們綜合了生物信息、化學信息學、網絡藥理學及計算建模等多種手段，通過GEO差異基因分析篩選RA關鍵基因，再結合已知數據庫中的RA的相關靶點信息，構建了一個較為可靠的RA疾病PPI網絡；選取了兩種不同證型的經典藥對：獨活-桑寄生和黃柏-蒼術，整理得到了兩種藥對的主要成分數據，并預測其潛在作用靶點；通過我們所提出的基于通絡擾動的相似性計算建模方法，挖掘出了不同藥對中具有協同組合治療的小分子組合；本文所提出的研究方法，從一個全新角度揭示了中藥藥對內在成分的作用機制和規律，為后續深入開展實驗研究提供技術上有力的支持。