基于生物信息學的系統性硬化癥潛在miRNA-mRNA調控網絡的構建

單雨，卞華，2，郭嘉，陳爽，卞博

（南陽理工學院 1.河南省張仲景方藥與免疫調節重點實驗室；2.張仲景國醫國藥學院，河南 南陽 473004）

系統性硬化癥（systemic sclerosis，SSc）是以皮膚和內臟纖維化、炎癥和血管病變為特征的特發性全身性自身免疫性疾病，也是一種慢性炎性結締組織疾病，累及皮膚、肌肉骨骼系統和多個內臟器官［1-2］。我國SSc發病率和患病率居自身免疫性疾病的第三位，僅次于類風濕性關節炎和系統性紅斑狼瘡［3］。我國成年人SSc的年發病率較高，為（190～250）/100萬［4-5］。研究［6］顯示，SSc的整體發病率不高，但卻有較高的病死率。目前，SSc的發病機制還未明確，因此給臨床治療帶來了極大的困難。

微RNA（microRNA，miRNA）是一類內源性小分子非編碼RNA，由18～25個核苷酸組成［7］。miRNA通過直接降解mRNA或阻斷蛋白質的翻譯過程調控基因表達，從而參與細胞增殖、分化、凋亡等許多生物學過程［8-9］。miRNA通過靶向細胞外基質（extracellular matrix，ECM）蛋白、結締組織生長因子（connective tissue growth factor，CTGF）、轉化生長因子-β（transforming growth factor-β，TGF-β）通路、上皮-間質轉化（epithelial-mesenchymal transition，EMT）或肌成纖維細胞增殖來增強或抑制纖維化［10］。miRNA相對穩定，可在組織和體液中檢測到，因此有希望成為監測疾病進程和治療反應的生物標志物［11］。研究［12］發現miRNA-181b-5p、miRNA-223-3p以及miRNA-21-5p在SSc中的表達顯著增加。WUTTGE等［13］研究發現miR-20a-5p和miR-203a-3p可作為SSc相關的肺動脈高壓發展的預測因子。BORTOLOTT等［14］發現包括miRNA-409在內的特定miRNA可作為SSc診斷、進展評估和抗纖維化治療的生物標志物，這些都說明miRNA在SSc診斷和治療中具有極大潛力。本研究基于生物信息學分析構建SSc潛在的miRNA-mRNA調控網絡，希望能夠提供有前景的SSc診斷和治療靶點。

1 材料與方法

1.1 miRNA數據集獲取

搜索基因表達圖譜（Gene Expression Omnibus，GEO）數據庫高通量基因表達數據庫（https：//www.ncbi.nlm.nih.gov/geo/），關鍵詞選擇“systemic sclerosis”或“scleroderma”；研究類型選擇“array by non-coding RNA profiles”；種族選擇“home sapiens”，篩選SSc血漿相關miRNA數據集，見表1。選取GSE81293數據集進行后續分析。

表1 SSc相關數據集的信息Tab.1 Information about SSc related data sets

1.2 差異表達miRNA（differential expression miRNA，DEM）篩選、DEM靶基因獲取和靶基因上游轉錄因子的分析

GEO2R（https：//www.ncbi.nlm.nih.gov/geo/geo2r/）是基于R語言的網絡應用程序，可用于篩選SSc患者和正常健康者之間的DEM。采用Benjamini和Hochberg方法調整P值，糾正假陽性結果的發生，將調整后的P＜0.05和|log FC|＞1設置為閾值。利用TBtools工具包在樣本中篩選出DEM。利用 MiRNet（https：//www.mirnet.ca/）獲取DEM的潛在靶基因。利用Fun-Rich（https：//www.funrich.org/）分析靶基因的潛在上游轉錄因子，P＜0.05為差異有統計學意義。

1.3 靶基因基因本體（Gene Ontology，GO）和京都基因和基因組數據庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）分析

為研究DEM靶基因的生物學功能，利用DAVID在線注釋、可視化和集成發現數據庫（https：//www.david.ncifcrf.gov/）進行GO和KEGG路徑分析，P＜0.05為差異有統計學意義。

1.4 蛋白質-蛋白質相互作用（protein-protein interaction，PPI）網絡的構建以及核心（Hub）基因的篩選

利用STRING在線數據庫（https：//www.stringdb.org/）對DEM的靶基因進行PPI網絡構建和分析，并在Cytoscape軟件中進行可視化。利用Cytoscape的插件CytoHubba，通過幾種拓撲算法，根據預先加載的PPI網絡的特征對節點進行排序。采用最大Clique中心性（maximal clique centrality，MCC）方法，篩選出在PPI網絡中排名前10位的Hub基因。

1.5 Hub基因表達的驗證

從GEO數據庫下載GSE146093數據集的表達矩陣及其GPL570注釋文件，根據注釋文件找到表達矩陣中相關的Hub基因，驗證并分析Hub基因的表達。顯著Hub基因需要滿足2個條件：（1）上調DEM的靶基因顯著下調或者下調的DEM的靶基因顯著上調；（2）P＜0.05為差異有統計學意義。

2 結果

2.1 DEM篩選結果

使用GEO在線工具對GSE81293數據集進行分析，該數據集包含15個SSc樣本和5個正常對照樣本，以調整后P＜0.05和|logFC|＞1為閾值。篩選得到26個DEM，其中上調有let-7e-5p、miR-455-3p，下調有miR-4459、miR-4507、miR-4417、miR-3687、miR-1246、miR-4668-5p、miR-4530、miR-4299、miR-4516、miR-663b、miR-762、miR-665、miR-4463、miR-4656、miR-3663-3p、miR-4505、miR-4440、miR-4685-5p、miR-4281、miR-149-5p、miR-4793-3p、miR-27a-5p、miR-1469、miR-4508。見表2、圖1。

圖1 DEM熱圖Fig.1 Heat map of DEM

2.2 DEM靶基因的篩選結果

通過miRNet數據庫共篩選出候選DEM的靶向基因5 450個，其中上調 1 407個，下調4 043個。每個DEM的靶基因數量見表3。

表3 DEM上調和下調潛在靶基因數量Tab.3 Number of potential target genes for up-and down-regulated DEM

2.3 靶基因轉錄因子的分析

利用FunRich軟件預測DEM靶基因的轉錄因子，結果顯示，上調DEM的轉錄因子包括NFYA、SP4、NRF1、ELK1、SP1、YY1；下調DEM的轉錄因子包括CACO、KLF7、E2F1、YY1、MYC、SP1。見圖2。

2.4 靶基因GO和KEGG分析

生物過程（biological process，BP）分析顯示，DEM靶基因主要富集在DNA轉錄及其調控、細胞分裂和細胞間黏附方面；細胞組分（cellular component，CC）分析顯示，靶基因主要富集于細胞膜、細胞溶質、核質以及蛋白質磺酰化等方面；分子功能（molecular function，MF）分析表明，靶基因在多聚腺嘌呤核糖核苷酸、DNA結合、轉錄因子結合、鈣黏蛋白結合參與細胞黏附等方面顯著富集。KEGG通路分析顯示，DEM靶基因顯著富集于P53信號通路、癌癥通路、FOXO信號通路、神經膠質瘤、RNA轉運、乙型肝炎、轉錄失調等。見圖3、4。

2.5 PPI和DEM-Hub基因網絡的構建

通過STRING數據庫分別建立了上調和下調DEM靶基因的PPI網絡。然后將上述的PPI網絡上傳到Cytoscape軟件，并使用Cytohubba插件篩選Hub基因。上調和下調的前10位Hub基因見圖5。

圖5 DEM靶基因的PPI網絡Fig.5 Protein-protein interaction network of DEM target genes

為了更好地研究DEM在SSc中的分子機制，利用Cytoscape軟件構建了DEM-Hub基因網絡。結果顯示，miR-455-3p與6個Hub基因（DDX55、RBM28、EBNA1BP2、IMP4、PWP2、RRP1B）相互作用；let-7e-5p與4個Hub基因（UTP15、RSL1D1、PDCD11、PES1）相互作用；miR-149-5p與6個Hub基因（DHX9、SRSF1、ILF3、HNRNPH1、MATR3、HNRNPA1）相互作用；miR-4440與PTBP1相互作用。

2.6 Hub基因表達的印證

根據DEM-Hub基因網絡，利用GSE146093數據集，印證前6個上調DEM靶基因（PWP2、UTP15、PDCD11、RSL1D1、RRP1B、EBNA1BP2）和下調DEM靶基因（HNRNPA1、EIF4A3、SRSFI、HNRNPH1、PTBP1、DHX9）的表達水平。對于上調的DEM，只有UTP15的表達呈現下降狀態，而PWP2、PDCD11、RRP1B、EBNA1BP2的表達則相反（圖6）。下調的DEM中，PTBP1的表達顯著增加，HNRNPA1、DHX9的表達反而下降（圖7）。因此認為let-7e-5p-UTP15和miR-149-5p-PTBP1是SSc潛在的調控途徑，見表4。

圖6 上調DEM中Hub基因的表達印證Fig.6 Validation of the expression of Hub gene in up-regulated DEM

圖7 下調DEM中Hub基因的表達印證Fig.7 Validation of the expression of Hub gene in down-regulated DEM

表4 Hub基因表達印證結果Tab.4 Validation of the expression of Hub gene

3 討論

隨著微陣列技術的發展，對各種疾病進展中的數千種基因的變化有了更深的研究。本研究提取GSE81293數據集數據，識別和篩選SSc與正常人血漿的DEM。結果顯示，共篩選出了2個上調的DEM（let-7e-5p、miR-455-3p）和24個下調的DEM（miR-4459、miR-4507、miR-4417、miR-3687、miR-1246、miR-4668-5p、miR-4530、miR-4299、miR-4516、miR-663b、miR-762、miR-665、miR-4463、miR-4656、miR-3663-3p、miR-4505、miR-4440、miR-4685-5p、miR-4281、miR-149-5p、miR-4793-3p、miR-27a-5p、miR-1469、miR-4508）。

研究［15］表明，miRNA表達可以被轉錄因子調控，本研究預測了SSc與正常人血漿的DEM的上游轉錄因子。其中，上調DEM的轉錄因子主要是SP1、SP4和NFYA，下調DEM的轉錄因子是SP1和KLF7，而且SP1在上行調控和下行調控中均占比最高。SP1是重要的轉錄因子，有研究［16-17］指出SP1可通過參與調控腫瘤細胞增殖相關轉錄因子FOXM1在卵巢癌中的表達，促進卵巢癌病情進展。KYPRIOTOU等［18］發現SP1和SP4轉錄因子的相互作用可以促進SSc成纖維細胞中Ⅰ型膠原的表達。NFYA在多種癌癥中的表達增加，與不良預后相關；并且在氧化固醇結合蛋白的基礎轉錄中起著重要的作用，可以通過P53/SRE BF2/NFYA信號通路來下調氧化固醇結合蛋白2的表達［19］。KLF7是神經系統發育和脂肪形成的關鍵因子，JIANG等［20］研究表明KLF7的表達可誘導內皮細胞的氧化應激、炎癥和細胞凋亡。

GO和KEGG分析顯示，DEM主要富集于DNA轉錄及其調控、細胞分裂和細胞間黏附方面以及P53信號通路、FOXO信號通路。P53可通過調控多種細胞免疫相關信號通路來參與機體的免疫應答［21］，通過上調Toll 樣受體增加促炎性細胞因子的產生［22］，P53功能失調與SSc 患者的成纖維細胞中的細胞凋亡途徑受損相關，可能導致SSc患者皮膚的慢性纖維化的發生［23］。FOXO是介導細胞存活的經典通路，PI3K/Akt/FOXO的相互作用對機體的免疫功能產生影響，在機體細胞的增殖、凋亡等方面發揮重要作用［24］。

通過構建DEM-Hub基因網絡，本研究發現大多數Hub基因與miR-455-3p、let-7e-5p、miR-149-5p聯系密切。在排名前12位的Hub基因中，只有UTP15、PTBP1的表達與GSE146093數據集一致，可能是樣本的來源不同所致。MOUILLESSEAUX等［25］發現UTP15的缺陷對細胞存活和血管系統的正常運行產生了影響，并且與P53信號通路聯系密切，而GUI等［26］研究結果顯示let-7e-5p可通過激活JAK1/STAT3通路抑制促炎性細胞因子表達。本研究驗證過程中發現，作為上調let-7e-5p的靶基因UTP15在SSc中下調表達，因此let-7e-5p-UTP15有可能作為SSc潛在的調控途徑。PTBP1 在人類幾乎所有的細胞類型中都有表達，并且介導多種細胞過程，包括神經元細胞的生長和分化以及免疫細胞的激活，其功能受多種分子調控，包括miRNA、長鏈非編碼RNA和RNA結合蛋白［27］。通過miR-4440-PTBP1途徑調控PTBP1的表達從而抑制免疫細胞的激活，可能成為控制SSc發展的一種有效方法。

綜上所述，基于GEO數據庫和生物信息學分析，本研究構建了SSc血漿中潛在的miRNA-mRNA調控網絡，發現了let-7e-5p-UTP15、miR-4440-PTBP1可作為SSc患者診斷和治療的潛在生物標志物。本研究的不足之處：（1）僅研究了SSc與正常對照之間的miRNA和mRNA，在SSc的不同階段結果可能會不同；（2）與生物標志物分析通常需要的樣本量比較，本研究采用的數據集的樣本量較小；（3）僅基于公共數據庫預測miRNA-mRNA的相互作用結果存在局限性，還需進一步通過體內和體外實驗驗證。