徐琪，謝樹仁，王嘉琪，才讓杰，黃佳雨，陳強(qiáng)，張婭，李永生

海南醫(yī)學(xué)院生物醫(yī)學(xué)信息與工程學(xué)院，海南 海口 571199

絕大部分傳染性疾病源于熱帶地區(qū)，近年來SARS、禽流感、埃博拉、寨卡熱的流行和新型冠狀病毒感染對人類健康都造成了重大影響[1]。登革熱(dengue fever)是由登革熱病毒感染、伊蚊傳播的一種急性傳染病，無法直接經(jīng)過呼吸道、消化道或接觸在宿主之間傳播。理解宿主-病原體間復(fù)雜的關(guān)系是開發(fā)新的診斷、預(yù)防和治療模型及建立動態(tài)防控體系的關(guān)鍵前提。隨著高通量測序技術(shù)的發(fā)展，基于多組學(xué)生物大數(shù)據(jù)已經(jīng)發(fā)現(xiàn)了一些與傳染病相關(guān)的關(guān)鍵基因及藥物靶點(diǎn)。例如，基于外周血的轉(zhuǎn)錄組測序，Nikolayeva等[2]識別了18個(gè)基因生物標(biāo)記能夠用于識別嚴(yán)重的登革熱患者。Robinson等[3]通過整合多套基因表達(dá)數(shù)據(jù)識別了20個(gè)基因標(biāo)記物用于預(yù)測登革熱患者的疾病進(jìn)展。然而對于登革熱感染關(guān)鍵分子標(biāo)記的識別仍然具有挑戰(zhàn)。

此外，登革熱感染經(jīng)常會伴隨著人類其他復(fù)雜疾病的產(chǎn)生，患者治愈后期會有多種并發(fā)癥。解析登革熱感染與人類其他復(fù)雜疾病的關(guān)系對于登革熱感染后期并發(fā)癥的早期預(yù)測具有十分重要的意義。同時(shí)，目前已經(jīng)開發(fā)了一些登革熱治療的方法，但是仍然需要新的治療方案的提出。如何基于多維組學(xué)數(shù)據(jù)整合分析識別新的候選藥物也是亟待解決的關(guān)鍵問題之一。本研究通過整合登革熱感染相關(guān)的多層次信息、挖掘潛在的分子標(biāo)記，并基于網(wǎng)絡(luò)分析技術(shù)識別登革熱感染與其他疾病間的關(guān)聯(lián)，優(yōu)選潛在的治療藥物，對于登革熱感染等疾病的預(yù)防和治療具有重要的理論意義。

1 資料與方法

1.1 微陣列芯片數(shù)據(jù)的獲取與預(yù)處理首先從GEO和ArrayExpress數(shù)據(jù)庫中共獲得了4套登革熱感染相關(guān)基因表達(dá)數(shù)據(jù)集，如表1所示。通過對芯片數(shù)據(jù)進(jìn)行對數(shù)化、探針過濾、補(bǔ)缺失值、標(biāo)準(zhǔn)化、探針注釋等預(yù)處理獲得基因表達(dá)矩陣。

表1 本研究使用數(shù)據(jù)集情況Table 1 Datasets used in this study

1.2 差異表達(dá)基因的篩選與整合根據(jù)樣本分組信息，基于Limma包對數(shù)據(jù)集進(jìn)行差異表達(dá)分析。首先，輸入基因表達(dá)矩陣，分組矩陣并構(gòu)建差異比較矩陣；然后將數(shù)據(jù)擬合到具有l(wèi)mFit函數(shù)的線性模型；使用contrasts.fit函數(shù)根據(jù)對比度模型計(jì)算差異；進(jìn)一步使用eBayes函數(shù)進(jìn)行貝葉斯測試；最后針對具有topTable功能的所有基因生成測試結(jié)果，使用Benjamini和Hochberg矯正方法校正P值。基于閾值FDR＜0.05，得到各個(gè)登革熱感染相關(guān)數(shù)據(jù)集的差異表達(dá)基因，并繪制對應(yīng)的火山圖與熱圖。

由于實(shí)驗(yàn)樣本的異質(zhì)性，使用不同的檢測平臺和數(shù)據(jù)處理方法均會導(dǎo)致結(jié)果不一致。因此，進(jìn)一步利用“RobustRankAggreg”[4]方法將來自不同數(shù)據(jù)集的分析結(jié)果整合，設(shè)閾值FDR＜0.05，fold change＜1.2，得到整合的差異表達(dá)基因。

1.3 基于隨機(jī)游走模型優(yōu)選生物分子標(biāo)記物隨機(jī)游走是給定源節(jié)點(diǎn)開始從當(dāng)前節(jié)點(diǎn)到隨機(jī)選擇的鄰居迭代walker轉(zhuǎn)換[5]。以RRA整合得到的生物分子標(biāo)記物中FDR＜0.01的基因作為種子節(jié)點(diǎn)，映射到HPRD數(shù)據(jù)庫的人類蛋白質(zhì)相互作用網(wǎng)絡(luò)進(jìn)行隨機(jī)游走分析。提取得分前5%的基因作為優(yōu)選的生物分子標(biāo)記物，利用Cytoscape從人類蛋白質(zhì)互作網(wǎng)絡(luò)中提取了439個(gè)優(yōu)選生物分子標(biāo)記物的PPI子網(wǎng)，并刪除了網(wǎng)絡(luò)中沒有互作的節(jié)點(diǎn)。最后將隨機(jī)游走得到的節(jié)點(diǎn)得分進(jìn)行排序后作為輸入，進(jìn)行基因集合功能富集分析(GSEA)。

1.4 登革熱病毒與人類蛋白之間的相互作用蛋白質(zhì)相互作用網(wǎng)絡(luò)的結(jié)構(gòu)和性質(zhì)是系統(tǒng)生物學(xué)中重要的研究內(nèi)容，通過網(wǎng)絡(luò)分析來評估相關(guān)疾病基因的特性是研究疾病最簡單方法之一，這有利于找到病毒潛在靶向的蛋白，并揭示有關(guān)其病因的重要線索。通過已有的登革熱公共數(shù)據(jù)庫(包括DenHunt[6]和DenvInt[7])和文獻(xiàn)[8-9]，收集到了登革熱病毒蛋白和人類蛋白質(zhì)之間的互作關(guān)系，得到8 526個(gè)PPIs，涉及11種登革熱病毒蛋白和2 863種人類蛋白。將隨機(jī)游走優(yōu)選后的節(jié)點(diǎn)映射到網(wǎng)絡(luò)中并使用“clusterprofiler”包對病毒蛋白靶向的人類蛋白進(jìn)行功能富集分析。

1.5 基于PPI的登革熱感染與復(fù)雜疾病的關(guān)聯(lián)為了揭示登革熱感染與人類復(fù)雜疾病之間的關(guān)系，通過查找文獻(xiàn)收集了多種復(fù)雜疾病相關(guān)的基因[10]，共涉及299種人類復(fù)雜疾病的3 173個(gè)基因，將這些疾病進(jìn)一步劃分為10個(gè)大類。隨后基于網(wǎng)絡(luò)之間的重疊評估了疾病相關(guān)基因和登革熱感染相關(guān)基因之間的相似度Svb值[11]，公式如下：

dvv和dbb分別代表登革熱感染相關(guān)蛋白網(wǎng)絡(luò)和人類復(fù)雜疾病蛋白網(wǎng)絡(luò)的平均最短路徑，dvb表示疾病b和登革熱v之間的成對平均最短距離，Svb＜0則表明登革熱v感染相關(guān)蛋白質(zhì)和與疾病b相關(guān)蛋白質(zhì)之間存在基于網(wǎng)絡(luò)的重疊模塊。

1.6 構(gòu)建藥物-靶點(diǎn)網(wǎng)絡(luò)優(yōu)選潛在抗病毒藥物利用DrugBank數(shù)據(jù)庫(v4.3)[12]、治療靶點(diǎn)數(shù)據(jù)庫(TTD)[13]、PharmGKB數(shù)據(jù)庫[14]、ChEMBL(v20)[15]、BindingDB[16]和IUPHAR/BPS藥理學(xué)指南[17]收集藥物-靶點(diǎn)相互作用信息。并從DrugBank[12]中提取每種藥物的化學(xué)結(jié)構(gòu)，采用SMILES格式。在這里，使用了符合以下三個(gè)標(biāo)準(zhǔn)的藥物-靶相互作用：(i)結(jié)合親和力，包括Ki、Kd、IC50、和EC50均小于等于10μmol/L；(ii)該目標(biāo)在UniProt數(shù)據(jù)庫[18]中標(biāo)記為“已審查”；(iii)人類靶標(biāo)由唯一的UniProt登錄號表示。本課題組最近的研究[19-21]提供了經(jīng)實(shí)驗(yàn)驗(yàn)證的藥物-靶點(diǎn)網(wǎng)絡(luò)的詳細(xì)信息。

進(jìn)一步基于互作網(wǎng)絡(luò)對候選藥物進(jìn)行了優(yōu)先級排序?；谂c病毒互作的宿主基因編碼蛋白質(zhì)V和藥物靶標(biāo)T，計(jì)算了V與每個(gè)候選藥物的靶標(biāo)蛋白質(zhì)T的網(wǎng)絡(luò)鄰近度。計(jì)算公式如下：

d(i,j)是人類蛋白質(zhì)相互作用網(wǎng)絡(luò)中蛋白質(zhì)i(病毒互作蛋白)和j(藥物靶標(biāo)蛋白)之間的最短距離?；陔S機(jī)檢驗(yàn)將網(wǎng)絡(luò)鄰近度結(jié)果進(jìn)一步轉(zhuǎn)換為Z-score：

2 結(jié)果

2.1 轉(zhuǎn)錄組分析識別登革熱感染的關(guān)鍵基因使用limma R包對四套數(shù)據(jù)集進(jìn)行差異分析，卡閾值FDR＜0.05篩選差異表達(dá)的基因，各個(gè)數(shù)據(jù)集差異表達(dá)基因數(shù)目如表2所示。并對四套數(shù)據(jù)集分別繪制了火山圖和熱圖(圖1)。

表2 不同數(shù)據(jù)集識別的登革熱感染上下調(diào)基因情況Table 2 Numbers of differentially expressed genes

圖1 登革熱感染四套數(shù)據(jù)集差異分析結(jié)果Figure 1 Differential expression of genes in dengue infections

通過熱圖觀察到基因表達(dá)模式具有異質(zhì)性，在GSE84331數(shù)據(jù)集中，轉(zhuǎn)錄組分析可以很明顯地將登革熱疾病樣本與正常樣本區(qū)分開，而在GSE18090和GSE51808的熱圖中，正常樣本和登革熱早期(DF)階段的疾病樣本之間并未檢測到明顯的差異，卻在DF到DHF發(fā)展的階段疾病轉(zhuǎn)錄組發(fā)生了急劇的變化(圖1)。在疾病恢復(fù)期，轉(zhuǎn)錄組水平與嚴(yán)重登革熱相比并未發(fā)生明顯改變。因此，轉(zhuǎn)錄組分析可區(qū)分早期登革熱患者和登革出血熱患者。

2.2 RRA整合分析識別穩(wěn)健的生物分子標(biāo)記隨后基于RRA方法整合了四個(gè)數(shù)據(jù)集的DEGs，并基于logFC的大小進(jìn)行排序，取FDR＜0.01的691個(gè)上調(diào)和1 071個(gè)下調(diào)DEGs繪制熱圖(圖2A)。CD1C陽性樹突狀細(xì)胞是真皮樹突狀細(xì)胞所有亞群中最豐富的，研究表明，人類皮膚樹突狀細(xì)胞可能是登革熱病毒感染的重要目標(biāo)，它們會將傳染性病毒從皮膚轉(zhuǎn)移到引流淋巴結(jié)，從而有效的激活全身免疫反應(yīng)，CD1C可能是引發(fā)早期適應(yīng)性抗登革病毒T細(xì)胞反應(yīng)的重要因素，這進(jìn)一步確定了CD1C細(xì)胞是登革熱的潛在靶標(biāo)[22]。隨后通過觀察在四個(gè)數(shù)據(jù)集里均差異表達(dá)的基因CASP71(圖2B)，以及在三個(gè)數(shù)據(jù)集中顯著差異的基因FEZ1、PITPNC1、SEMA4C、TXLNA(圖2C～E)的表達(dá)情況，發(fā)現(xiàn)差異基因在不同數(shù)據(jù)集中的差異程度都有所不同。

圖2 穩(wěn)健差異表達(dá)基因的鑒定Figure 2 Robust differentially expressed genes

2.3 基于隨機(jī)游走優(yōu)選登革熱分子標(biāo)記蛋白質(zhì)-蛋白質(zhì)相互作用(PPI)網(wǎng)絡(luò)由無向圖表示(圖3A)。進(jìn)一步基于專門存儲經(jīng)過實(shí)驗(yàn)驗(yàn)證人類蛋白質(zhì)互作信息的HPRD數(shù)據(jù)庫中[23]的網(wǎng)絡(luò)進(jìn)行隨機(jī)游走分析，并取經(jīng)過RRA整合后FDR＜0.01的1 762個(gè)差異基因作為隨機(jī)游走模型的種子節(jié)點(diǎn)，通過迭代得到每個(gè)節(jié)點(diǎn)的在隨機(jī)游走模型中的得分，排序后提取前5%共439個(gè)節(jié)點(diǎn)，使用Cytoscape3.8.0將其映射到HPRD網(wǎng)絡(luò)來構(gòu)建優(yōu)選后的蛋白互作網(wǎng)絡(luò)(圖3A)。該網(wǎng)絡(luò)刪除了沒有互作的節(jié)點(diǎn)之后包括304個(gè)節(jié)點(diǎn)(node)，726條邊(edge)，同時(shí)對網(wǎng)絡(luò)的拓?fù)鋵傩赃M(jìn)行了分析如表3所示。

圖3 優(yōu)選登革熱分子標(biāo)記的功能富集Figure 3 Functional enrichment of prioritized biomarkers for dengue infection

表3 PPI網(wǎng)絡(luò)拓?fù)鋵傩訲able 3 Topological features of PPI network

TP53作為常見的與細(xì)胞周期控制密切相關(guān)的基因在該網(wǎng)絡(luò)中具有最高的連通度為36，已有研究證明，細(xì)胞周期停滯是p53激活后觀察到的主要生物學(xué)結(jié)果，可防止受損DNA的積累以及基因組不穩(wěn)定性[24]。其次CDK1基因編碼的蛋白質(zhì)作為Ser/Thr蛋白激酶家族的成員是高度保守的蛋白激酶復(fù)合物催化亞基，被稱為M期促進(jìn)因子(MPF)，它對真核細(xì)胞周期的G1/S和G2/M相變必不可少，這表明在登革熱疾病感染過程出現(xiàn)了細(xì)胞周期的改變[25]。

為了分析登革熱疾病中的分子生物學(xué)功能，將隨機(jī)游走的得分排序后作為GSEA的輸入，并對其進(jìn)行KEGG和GO功能富集分析得到了一些免疫反應(yīng)激活信號轉(zhuǎn)導(dǎo)通路，例如，體液和細(xì)胞適應(yīng)性免疫是人體抵御病毒感染的重要途經(jīng)，它在DENV感染期間被激活，眾多研究表明人體細(xì)胞被病毒感染后，免疫系統(tǒng)功能異常[26-28]，有研究通過對嚴(yán)重登革熱過度免疫激活進(jìn)行研究充分證實(shí)了免疫反應(yīng)在登革熱發(fā)病機(jī)制中的關(guān)鍵作用[29]。同時(shí)還發(fā)現(xiàn)差異基因主要富集在DNA復(fù)制、細(xì)胞周期蛋白體酶等通路中(圖3B～3C)。且已有研究表明，病毒感染細(xì)胞時(shí)，病毒編碼的蛋白或DNA會擾亂細(xì)胞周期通路，促進(jìn)細(xì)胞向G1/S期轉(zhuǎn)化或者使細(xì)胞靜息于G2/M期，這些現(xiàn)象表明在登革熱疾病感染過程出現(xiàn)了細(xì)胞周期的改變[30]。

2.4 登革熱相關(guān)基因與病毒蛋白高頻率互作接下來利用已有的登革熱公共數(shù)據(jù)庫(包括DenHunt[6]和DenvInt[7])和文獻(xiàn)[8-9]收集了病毒與人類蛋白的互作信息，將隨機(jī)游走優(yōu)選后的304個(gè)節(jié)點(diǎn)映射到網(wǎng)絡(luò)中，生成病毒與人類蛋白互作網(wǎng)絡(luò)(圖4A)。由于節(jié)點(diǎn)的拓?fù)湫再|(zhì)可以通過其在網(wǎng)絡(luò)中的位置來顯示，因此分別比較受病毒蛋白靶向的人類蛋白和HPRD網(wǎng)絡(luò)中其他蛋白的介數(shù)、緊密度、連通度三個(gè)拓?fù)鋵傩?圖4B～4D)。發(fā)現(xiàn)病毒靶向的人類蛋白介數(shù)、緊密度、連通度均明顯高于其他蛋白，表明病毒與人類蛋白互作網(wǎng)絡(luò)具有很強(qiáng)的模塊性，且位于網(wǎng)絡(luò)的中心位置，進(jìn)一步證實(shí)病毒可能會針對位于網(wǎng)絡(luò)中心或在信息傳播中起重要作用的蛋白質(zhì)進(jìn)行復(fù)制。隨后利用“clusterprofiler”包對病毒靶向的人類蛋白進(jìn)行功能富集分析，發(fā)現(xiàn)大部分顯著富集的功能都與病毒感染相關(guān)，而且差異蛋白主要富集在蛋白酶體蛋白質(zhì)分解代謝過程、肽鏈內(nèi)切酶復(fù)合體、單鏈DNA結(jié)合等通路。已有研究表明，在登革熱感染過程中，自噬有助于病毒復(fù)制，而病毒可以通過蛋白酶體降解逐漸降低自噬受體的水平，說明蛋白酶體分解代謝通路與登革熱病毒的發(fā)展密切相關(guān)[31]。

圖4 登革熱相關(guān)基因與病毒蛋白互作Figure 4 PPI interactions among dengue proteins and human proteins

2.5 登革熱感染與多種復(fù)雜疾病關(guān)聯(lián)隨后從文獻(xiàn)中收集了涉及299種人類復(fù)雜疾病的3 173個(gè)基因，來探索登革熱感染與人類復(fù)雜疾病之間的關(guān)系。通過計(jì)算疾病相關(guān)基因和登革熱感染相關(guān)基因之間基于網(wǎng)絡(luò)的的相似度Svb值，發(fā)現(xiàn)了登革熱感染可能與癌癥、消化系統(tǒng)腫瘤、哮喘、類風(fēng)濕性關(guān)節(jié)炎，以及一些心血管疾病相關(guān)(圖5A)。特別是已有研究證實(shí)類風(fēng)濕性關(guān)節(jié)炎[32]、消化系統(tǒng)腫瘤[33-34]、腦血管疾病[35]與登革熱感染相關(guān)。通過對這三類疾病進(jìn)行分子網(wǎng)絡(luò)可視化發(fā)現(xiàn)，MIF、TP5、HNRNPD等基因不僅靶向登革熱病毒還與復(fù)雜疾病有關(guān)。其中巨噬細(xì)胞遷移抑制因子MIF(圖5B)是一種多效性炎癥細(xì)胞因子，在先天性和適應(yīng)性免疫反應(yīng)的調(diào)節(jié)中非常重要，它與自身免疫性疾病的發(fā)病機(jī)制有關(guān)[36]，并且該基因的中和會導(dǎo)致DENV感染巨噬細(xì)胞的遷移能力提高，與DENV感染的嚴(yán)重程度有關(guān)[37]。有研究表明，p53(圖5C)介導(dǎo)的干擾素刺激基因在體外和體內(nèi)都能抑制病毒的復(fù)制[38]，并且TP53的突變可作為患者的預(yù)后預(yù)測指標(biāo)[39]。HNRNPD(圖5D)作為一種通用的黃病毒宿主因子與AUF1蛋白結(jié)合可以促使登革熱病毒復(fù)制[40]。

圖5 登革熱與人類復(fù)雜疾病之間的關(guān)系Figure 5 Relationships of dengue infection and human diseases

2.6 登革熱感染潛在藥物的挖掘根據(jù)ZdVT對優(yōu)選藥物進(jìn)行排序，共獲得39個(gè)與登革熱疾病相關(guān)的藥物。經(jīng)篩選后獲得了23種至少靶向兩個(gè)登革熱蛋白的藥物(圖6A)。已有研究表明，DENV可能導(dǎo)致人體多種疾病，包括癌癥、哮喘、潰瘍性結(jié)腸炎等[41]，其中哮喘可以作為預(yù)測登革熱發(fā)展的標(biāo)志[42]，而藥物關(guān)聯(lián)性分析獲得的oxtriphylline已被FDA批準(zhǔn)用于治療哮喘[43](圖6B)。此外，登革熱病毒感染常伴有關(guān)節(jié)炎發(fā)生[44]，Tenoxicam、Acemetacin、Mefenamic acid、Etoricoxib、Loxoprofen作為常見的非甾體抗炎藥(圖6C)，常用于緩解骨關(guān)節(jié)炎的癥狀[45-49]。同時(shí)，Omega-3 fatty acids(圖6D)可通過降低甘油三酯血癥來降低動脈粥樣硬化性心血管疾病的發(fā)病風(fēng)險(xiǎn)[50]。Romidepsin(圖6E)作為一種組蛋白去乙?；敢种苿?，對多種實(shí)體瘤具有抗腫瘤作用[51]。此研究可以為DENV的治療及潛在藥物的發(fā)掘提供參考。

圖6 潛在藥物的挖掘Figure 6 Potential drugs for dengue infection

3 討論

登革熱感染作為一種傳染性極強(qiáng)、發(fā)病率極高的蚊媒傳染性疾病，目前還沒有發(fā)現(xiàn)有效的抗病毒療法可以預(yù)防和治療，因此進(jìn)一步研究登革熱的分子生物學(xué)機(jī)制，探索更加有效的預(yù)防和治療方法是熱帶疾病研究領(lǐng)域亟待解決的關(guān)鍵。本研究旨在探討登革熱感染相關(guān)的生物分子標(biāo)記、潛在的藥物組合療法的挖掘，為登革熱感染發(fā)病機(jī)制的探索與診療提供了新的思路。

本文通過RRA方法將從四個(gè)數(shù)據(jù)集中獲取的差異基因信息整合到一起，再基于隨機(jī)游走模型優(yōu)選出關(guān)鍵生物分子標(biāo)記物，并建立全面的登革熱感染相關(guān)功能調(diào)控網(wǎng)絡(luò)，經(jīng)過功能富集分析發(fā)現(xiàn)大部分差異基因都富集在DNA復(fù)制、細(xì)胞周期相關(guān)的通路上。這表明在登革熱感染能影響細(xì)胞周期相關(guān)通路，感染登革熱病毒時(shí)會影響細(xì)胞活動，此外這些基因還富集在免疫相關(guān)通路上，而患者的免疫反應(yīng)與后期疾病發(fā)展的嚴(yán)重程度相關(guān)。

另外，對于登革熱(DF)和登革出血熱(DHF)兩組也做了差異分析，發(fā)現(xiàn)與DF和DHF相關(guān)基因主要富集在T細(xì)胞受體信號通路、適應(yīng)性免疫應(yīng)答、T細(xì)胞激活、調(diào)節(jié)淋巴細(xì)胞活化等免疫相關(guān)通路，進(jìn)一步表明人體細(xì)胞被病毒感染后，會出現(xiàn)免疫功能失調(diào)，導(dǎo)致免疫系統(tǒng)功能異常。最后通過分析登革熱與299種復(fù)雜疾病的關(guān)系進(jìn)一步尋找與登革熱相關(guān)的潛在藥物，在探究登革熱的預(yù)防和治療道路上取得新的進(jìn)展。

綜上所述，本研究通過建立全面的登革熱感染相關(guān)功能調(diào)控網(wǎng)絡(luò)，優(yōu)選登革熱感染相關(guān)的基因分子標(biāo)記，解析了與病毒感染相關(guān)的功能調(diào)控通路。通過探究登革熱與其他復(fù)雜疾病的關(guān)聯(lián)，和潛在藥物的挖掘，為登革熱感染預(yù)防和診療提供了新的思路。