基于轉錄組學和機器學習算法的肺結核鐵死亡相關關鍵基因的研究

葉江娥 方雪暉 熊延軍 劉盛盛

結核病仍然是嚴重危害人類健康的慢性傳染性疾病。該病的治療以化學治療為主,但存在治療方案復雜、治療療程長、耐藥產生、藥物不良反應多等不足,亟需尋求新的治療方法[1-2]。因此,如何進一步從基因層面探索結核病的藥物治療靶點,對臨床治愈結核病及耐藥結核病具有重要現實意義。

近年來,有學者研究發現結核分枝桿菌(Mycobacteriumtuberculosis, MTB)在宿主體內會誘導巨噬細胞壞死,而這一過程可能是通過游離鐵、線粒體超氧化物和脂質過氧化增加、谷胱甘肽和谷胱甘肽過氧化物酶降低有關,上述參與巨噬細胞壞死的物質是鐵死亡進程的重要標志,且研究發現加入鐵抑制劑后MTB感染的巨噬細胞死亡得到了一定抑制,為肺結核調控機制提供了新的思路[3]。目前,國內外關于肺結核與鐵死亡相關性的機制與基于轉錄組學的分析相對較少,因此,本研究利用轉錄組學和機器學習算法,探索肺結核與鐵死亡相關差異基因,以及對機器學習篩選得到的關鍵基因的藥物調控網絡進行分析,并驗證關鍵基因的差異性,為今后進一步尋找治療靶點提供依據。

資料和方法

一、研究資料

以“pulmonary tuberculosis”為關鍵詞從公共存儲庫NCBI GEO(http://www.ncbi.nlm.nih.gov/geo)中進行搜索,共得到647個肺結核轉錄組數據集,再根據測序類型(轉錄組學)、物種(HOMO sapiens)等條件對結果進行篩選,最終得到GSE153326和GSE67589兩套轉錄組數據集。下載GSE153326和GSE67589數據集,以GSE153326數據集為訓練組,以GSE67589數據集為驗證組,根據探針矩陣文件和注釋信息,分別對訓練組和驗證組數據集進行注釋,其中,前者由GPL21185注釋平臺生成,后者由GPL570注釋平臺生成,根據兩數據集中健康人群樣本和MTB陽性樣本分類特征獲得兩組臨床分組,即對照組1和對照組2(健康血液樣本)、實驗組1和實驗組2(MTB陽性樣本)。

二、研究方法及統計學處理

1.結核病差異表達基因(tuberculosis differentially expressed gene, TBDEG)的獲取:通過R軟件包“limma”分別從GSE153326和GSE67589數據集中獲取與健康人群比對后的TBDEG。被鑒定出來的TBDEG均需要滿足調整P<0.05和log2(Fold-change)>1。

2.鐵死亡相關結核病差異表達基因(TBFerDEG)的鑒定:FerrDb V2數據庫[4]主要收錄了與鐵死亡相關的驅動因子(driver)、抑制因子(suppressor)和標記因子(marker)三類調控蛋白與疾病關系的數據,是研究鐵蛋白調節劑和疾病關系的重要工具。參考FerrDb V2數據庫,使用R軟件包“limma”對GSE153326注釋平臺獲得的TBDEG進行TBFerDEG鑒定。借助corrplot軟件包對訓練組中416個TBFerDEG進行相關性分析,剔除差異無統計學意義的相關基因。

3.TBFerDEG功能富集分析:使用R語言的“enrichplot包”“ComplexHeatmap包”和“clusterProfiler包”分析通過轉錄組學篩選得到的TBFerDEG的基因本體(GO)功能富集與京都基因和基因組百科全書(KEGG)通路[5]。

4.機器學習算法篩選TBFerDEG關鍵基因:應用機器學習算法中的LASSO擬合廣義線性模型[6]和SVM支持向量機遞歸特征[7]篩選TBFerDEG中差異有統計學意義的基因,選二者的交集基因,獲得交叉和驗證的交集基因集合,并將其稱為TBFerDEG關鍵基因。

5.構建受試者工作特征曲線(receiver operating characteristic curve, ROC):將結核病相關鐵死亡關鍵基因置于R軟件包“pROC”中,繪制ROC曲線,以曲線下面積(area under curve,AUC)評估關鍵基因的診斷效能。以AUC>0.7為篩選標準,AUC越大,表明該基因作為診斷的準確度越高。

6.構建藥物調控網絡:由于1個靶點可能與多個藥物調控有關,將得到的TBFerDEG通過在線數據庫“https://dgidb.genome.wustl.edu”獲取藥物和鐵死亡關鍵基因的相互作用,再將藥物和基因作用關系的數據導入cytoscape v3.7.1軟件包,構建藥物調控網絡[8],推導出與其相關的調控藥物網絡。其中,在基因轉錄成mRNA時受到正向調控,促進表達為上調作用;受到抑制,表達量減少則為下調作用。

7.驗證組數據集的篩選、矯正和注釋:從GEO數據庫下載GSE67589數據集,再根據GSE67589上GPL570平臺文件獲得驗證組數據集的矩陣信息,以對篩選得到的TBFerDEG進行注釋,獲得統一規范的基因名。

8.驗證TBFerDEG關鍵基因:結合驗證組(GSE67589)獲得的TBDEG,檢驗訓練組(GSE153326)獲得的TBFerDEG關鍵基因的表達量在驗證組健康人群和MTB陽性人群樣本中是否也具有差異性,以驗證其在肺結核患者中的獨立性與普適性。

結 果

一、數據處理與鐵死亡基因表達分析

根據訓練組和驗證組數據集注釋結果和兩數據集中健康人群樣本和MTB陽性樣本分類特征臨床分組情況,獲得8份對照組1和30份對照組2健康血液樣本、52份實驗組1和27份實驗組2 MTB陽性血液樣本。根據臨床分組信息,通過R軟件包“limma”分別從訓練組和驗證組數據集中獲取與健康人群比對后的TBDEG,再根據鐵死亡相關基因數據庫(FerrDb V2, http://www.zhounan.org/ferrdb/current/)信息獲得416個訓練組TBFerDEG和66個驗證組TBDEG。

二、差異基因及其相關性分析

對訓練組中416個TBFerDEG進行相關性分析,剔除差異無統計學意義的相關基因,最終得到56個TBFerDEG(表1)。

表1 訓練組中56個TBFerDEG基因

三、GO和KEGG分析

通過“enrichplot包”“ComplexHeatmap包”和“clusterProfiler包”對轉錄組學篩選得到的訓練組56個TBFerDEG進行GO富集和KEGG通路分析,并以p.adjust進行升序排列,得到GO富集分析結果(表2)和前15位KEGG通路結果(表3)。

表2 訓練組中56個TBFerDEG的GO富集分析結果

四、LASSO與SVM篩選鐵死亡變量

通過LASSO回歸模型分析訓練組56個TBFerDEG,共獲得13個交叉驗證誤差最小的TBFerDEG關鍵基因(ALOX12、BID、PEX3、TRIM26、SRC、CHMP6、SREBF1、STK11、AR、OIP5-AS1、FABP4、TMSB4Y、EZH2)。通過SVM向量機算法,篩選得到5個TBFerDEG關鍵基因,即BID、AR、STK11、ALOX12和SRC。取二者篩選得到的變量交集基因,最終得到5個TBFerDEG關鍵基因,即BID、AR、STK11、ALOX12和SRC。

五、鐵死亡關鍵基因的ROC分析

對篩選得到的5個基因進行ROC分析,結果表明:ALOX12、BID、SRC、STK11和AR的AUC分別為0.840、0.807、0.880、0.734和0.858(圖1)。

圖1 肺結核與鐵死亡相關關鍵基因的ROC曲線

六、藥物調控網絡分析

對機器學習算法篩選得到的5個關鍵基因進行藥物相關基因篩選,將得到的相關調控藥物與靶點的關系導入cytoscape v 3.7.2軟件,結果顯示,ALOX12、SRC和AR作為靶點在調節藥物網絡中起下調作用,STK11則起上調作用,而BID未得到與之相關的藥物。起興奮劑(agonist)調節作用的藥物有5種,起拮抗劑(antagonist)調節作用的有5種,起抑制劑(inhibitor)作用的藥物有18種,起調節劑(modulator)調節作用的有1種,具體見表4。

表4 5個關鍵基因與藥物構建的調控網絡節點特征關系

七、TBFerDEG關鍵基因的驗證

將篩選出的訓練集5個TBFerDEG關鍵基因導入驗證組數據集(GSE67589)中進行注釋,并結合獲得的66個驗證組TBDEG,對5個TBFerDEG關鍵基因在驗證組中的對照組2與實驗組2進行驗證(圖2～6)。結果顯示,AR和SRC關鍵基因在對照組2與實驗組2中的差異均有統計學意義(P值均<0.05)。

注 橫坐標代表驗證組數據集GSE67589各組的臨床特征,縱坐標代表關鍵基因在驗證組各樣本中的表達量

討 論

長期抗結核藥物治療容易引起肝腎功能損傷、胃腸道反應,降低機體抵抗力,誘發其他感染,給患者治療和愈后帶來嚴峻挑戰。最新研究發現,MTB的存活可能與鐵離子之間存在一定聯系,即存在于細胞膜上的IrtAB蛋白可以幫助MTB將機體內游離鐵轉運至細菌中,使宿主體內游離鐵含量降低,引起發病[9]。這一機制的發現有助于探索鐵死亡相關

機制與MTB在宿主細胞內存活的關系,也有助于新機制和新的治療靶點的發現,因此,鐵死亡特征靶點有望成為今后相關探索的重點。

本研究通過對訓練組基因芯片數據與鐵死亡在線公共數據庫相結合進行篩選,共得到56個TBFerDEG。再對這些TBFerDEG開展GO功能富集和KEGG通路分析,以獲得其生物學過程和信號通路。GO功能富集分析表明,上述基因參與的生物學過程主要有細胞對化學應激的反應、自噬調節、對金屬離子的響應和凋亡信號通路的負調控等。研究發現,金屬離子作為重要的代謝輔助因子和生理信號在人體生理和病理進程中有著重要的意義,其濃度的調節也是致病細菌與其宿主之間相互作用的生理學核心[10]。當機體鐵離子、鎂離子進入細胞內后,其胞外濃度會發生一定變化,從而可引起相關毒力因素的加重[11],也可能加重MTB對宿主的危害,這提示TBFerDEG能夠在一定程度上參與到機體內金屬離子等的調節過程,如游離鐵等的增加是MTB感染宿主期間造成巨噬細胞凋亡的重要因素,而其中的相關過程可能與金屬調控過程有關[12]。KEGG信號通路分析表明,與上述差異基因相關的通路有內分泌抵抗和鐵死亡等,提示上述基因可能參與到細胞的鐵死亡進程中。鐵死亡是一種由游離鐵和有毒脂質過氧化物積累引起的調節性壞死,在宿主感染MTB后,誘導巨噬細胞壞死,這種壞死過程可能與游離鐵、線粒體超氧化物和脂質過氧化的增加有關[13]。正如有實驗研究發現,當給小鼠加入了鐵抑制劑后,小鼠肺部感染率和壞死率均明顯降低,這表明鐵死亡是MTB感染引起細胞壞死的主要機制,也是宿主導向治療結核病的靶點[3]。

LASSO和SVM-RFE是機器學習算法中兩類最常見模型。SVM-RFE模型廣泛用于模式識別、機器學習等領域,能夠對多個變量進行模式識別,以結構風險最小化為原則,并最小化經驗誤差,以此提高學習性能,構建預測模型;LASSO回歸分析是通過對一系列不同λ值進行擬合,并對其中的相似值進行10折交叉驗證,有助于將高維數據中的重要變量進行快速有效地提取[14-16]。本研究借助LASSO和SVM算法獲得了5個與鐵死亡有關的關鍵基因,分別是BID、AR、STK11、ALOX12和SRC。通過對5個基因的ROC分析,發現5個關鍵基因均具有較好的診斷性能(AUC值均>0.7),為今后結核病的診斷和后續研究提供了一定參考。繼而構建關鍵基因相關藥物調控網絡,結果發現有5種藥物對AR靶點起興奮作用、5種藥物對AR靶點起對抗作用、17種藥物對SRC靶點起抑制作用、1種藥物對STK11靶點起抑制作用。進一步將5個關鍵基因導入驗證集數據集中進行驗證,結果表明,僅AR和SRC基因在驗證集中的表達差異均有統計學意義,提示這兩個靶點可能與抗結核治療有關,可在今后的藥物調控研究中進一步分析。既往研究發現,調節巨噬細胞中的吞噬小體成熟和自噬調節SRC酪氨酸激酶活性的過程是強毒性MTB生存所必需的,若SRC的表達受到抑制,則可顯著抑制單核細胞THP-1誘導分化為巨噬細胞,影響其胞內的H37Rv及耐多藥(MDR)和廣泛耐藥(XDR)MTB菌株的存活率,還可控制豚鼠MTB感染發生[17];AR是雄激素受體超家族成員之一,可競爭性地與雌激素受體元件進行結合,從而誘導細胞凋亡[15],但當前尚未見AR與MTB存活相關的研究報道,提示應在今后的機制研究中進一步挖掘其潛在的治療作用[16-17]。

轉錄組學是基于臨床樣本開展的高通量測序,旨在借助于現代測序與生物信息學手段挖掘出更多與肺結核發病機制相關的潛在靶點與通路[18],但本研究是對已發表和上傳至公共數據庫的轉錄組學數據開展的二次分析,缺乏本地區臨床樣本代表性,也缺乏實驗驗證,只能為實際臨床研究提供參考,是本研究存在的短板與不足,今后應采集臨床樣本獨立開展轉錄組學數據檢測,以篩選和驗證相關差異基因。

本研究利用轉錄組學和機器學習算法,探索肺結核與鐵死亡相關差異基因和關鍵基因的藥物調控網絡進行分析,并驗證關鍵基因差異的統計學意義,結果提示,AR和SRC兩個特征靶點有望成為今后探索與MTB之間關系的重點。盡管當前對AR與MTB的機制研究較少,使得本研究無法對其直接進行分析,但今后可將其作為研究的潛在靶點之一,更好地服務于臨床。

利益沖突所有作者均聲明不存在利益沖突

作者貢獻葉江娥:數據分析、論文撰寫;方雪暉:數據分析;熊延軍:論文指導和校對;劉盛盛:課題指導、論文修改