特發性肺纖維化相關的微陣列鑒定基因和microRNA分析

宋 鵬，張 猛，姜淑娟，楊棟梁

(1.山東第一醫科大學附屬省立醫院呼吸與危重癥醫學科，山東 濟南 250000；2.山東省濟南市人民醫院全科醫學科， 山東 濟南 271100；3.滄州醫學高等?？茖W校公共課教學部，河北 滄州 061000)

特發性肺纖維化(idiopathic pulmonary fibrosis，IPF)是特發性間質性肺炎最常見的形式，其特征在于臨床癥狀為咳嗽和呼吸困難，氣體交換受限以及進行性肺瘢痕形成[1]。目前己批準2種用于治療IPF的有效藥物(尼達尼布和吡非尼酮)[2-3]。但是，由于藥物的可及性及療效的限制，大多數IPF患者最終死于呼吸衰竭。IPF的預后仍然充滿挑戰，因此需要更加全面地了解其疾病發病機制。微小RNA(microRNA，miRNA)是一類非編碼RNA，可以通過與mRNA結合來抑制靶基因的表達[4]。 由于其在細胞凋亡、增殖和分化等細胞過程中的關鍵作用而引起了廣泛關注。在過去的十年中越來越多的研究尋找IPF潛在的生物標志物，多種miRNA可以作為IPF預測的潛在生物標記，包括miR-92a、miR-210、miR-29、miR-326、miR-98 和miR-let-7d[5-7]。但是單個miRNA的研究不能滿足IPF快速發展的要求，構建miRNA-mRNA聯合網絡可能有助于探索適用于臨床的生物標記物和治療靶標的進一步開發。本研究旨在通過基因表達綜合(gene expression omnibus，GEO)數據庫鑒定miRNA表達譜中的差異表達miRNA(differentially expressed miRNAs，DEMis)和差異表達mRNA(differentially expressed mRNAs，DEMs)，以探索IPF的生物學過程。使用全面的生物信息學分析檢查DEMs和DEMis之間的相關性，并結合從DEMis和DEMs數據中檢索到的信息，通過GO及KEGG途徑富集分析篩選出IPF的潛在生物標志物。

1 資 料 與 方 法

1.1數據集的獲取和分析 IPF相關miRNA和mRNA表達譜中的微陣列數據可從GEO數據庫(http://www.ncbi.nlmNih.gov/geo)中檢索和下載。使用“idiopathic pulmonary fibrosis”作為關鍵字進行查詢。搜索僅限于以下特定領域：研究類型，按陣列進行的表達譜分析和種類(智人)。下載miRNA表達微陣列數據集GSE27430和mRNA表達微陣列數據集GSE135065。使用R studio中的limma軟件包獲得IPF患者與健康對照組之間差異表達的miRNA和mRNA。貝葉斯方法糾正批量效應。如果有多個探針映射到同一基因中，則平均表達值被用來等于該基因的表達值。采用t檢驗過濾差異表達的基因。DEMis和DEMs通過P<0.05和log FC>平均(logFC)+ 2×SD(logFC)進行篩選。為了顯示DEMis和DEMs在不同樣品中的差異表達，通過在R studio中應用圖和Pheatmap包繪制火山圖和熱圖。

1.2轉錄因子 靶標的GO富集分析將DEMis上傳到FunRich(3.1.3)，這是用于對轉錄因子途徑的富集靶標(基因/mRNA)進行GO功能富集分析的常用工具。GO富集分析還用于開發miRNA、基因/mRNA和轉錄因子之間的相互作用網絡分析。

1.3DEMis靶向mRNA的預測 使用miRDB(http://www.mirdb.org)、miRTarBase(http://mirtarbase.mbc.nctu.edu.tw)和TargetScan(http://www.targetscan.org)數據庫來預測DEMis的靶向mRNA。之后，通過匹配之前選擇的DEMs進一步過濾DEMis的預測mRNA，然后得到DEMis-DEMs的配對組合。

1.4miRNA-mRNA調控網絡的構建 將上面選擇的所有DEMis-DEMs通過配對來構建miRNA-mRNA網絡，同時使用Cytoscape軟件(3.7.2版)對其進行可視化。

1.5網絡中mRNA的GO和KEGG富集分析 在R Studio中使用了richplot和ggplot2軟件包對網絡中的mRNA進行GO和KEGG途徑分析。

1.6統計學方法 應用R語言V3.6.0 統計軟件分析數據。計量資料比較采用t檢驗,為了進一步控制錯誤率，使用了Benjamini和Hochberg方法來計算調整后的P值。P<0.05為差異有統計學意義。

2 結 果

2.1IPF中的DEMis篩選結果 根據前所述排除標準，miRNA的log FC的臨界值為1.5。確定了21個上調或下調的miRNAs，包括14個上調的miRNA以及7個下調的miRNA，見表1。通過log10(P值)和log(FC)相關性的火山圖直觀地說明了IPF和健康對照之間miRNA表達差異的分布(圖1)。并繪制熱圖以直觀顯示IPF和健康對照組之間的差異(圖2)。

圖1 DEMis的火山圖。紅點:上調的基因；綠點:下調的基因

圖2 DEMis的熱圖。左側12個樣品來自對照組，右側13個樣品來自IPF組;紅色：高表達; 藍色：低表達

2.2差異miRNA的轉錄因子富集分析 32個轉錄因子中總共映射748個基因。通過轉錄因子靶標的富集，篩選了與miRNA密切相關的前10個轉錄因子，分別為特異性蛋白1(specific protein 1，SP1)、早期生長應答因子1(early growth response 1，EGR1)、NK6同源框蛋白1重組蛋白(recombinant NK6 homeobox protein 1，NKX6-1)、結構域 2 類轉錄因子 1(poudomain class 2 transcriptionfactor 1，POU2F1)、富含 AT 的交互式含域蛋白 3A(AT-rich interactive domain-containing protein 3，ARID3A)、心肌細胞特異性增強因子2A重組蛋白(recombinant myocyte specific enhancer factor 2A，MEF2A)、同源異型盒基因D8(homeobox gene D8，HOXD8)、同源異型盒基因A9(homeobox gene A9，HOXA9)、肝細胞核因子4α(hepatocyte nuclear factor 4-alpha，HNF4A)、兔抗人受體相關孤兒受體α(rabbit anti-human retinoic acid receptor related orphan receptor A，RORA)(圖3)，表明這些轉錄因子與DEMis有調節關系。

表1 IPF樣本中21種差異表達的miRNA(DEMis)Table 1 21 differentially expressed miRNAs(DEMis) in IPF samples

2.3IPF中的DEMs篩選結果 本研究探討了9例IPF患者和9例健康對照的mRNA表達水平。 mRNA的log(FC)閾值為0.47；確認了56個上調的mRNA和127個下調的mRNA。上調和下調的前20個基因(P<0.05為差異有統計學意義)見表2?；鹕綀D和熱圖見圖4～5。

2.4在IPF中重建miRNA-mRNA網絡 構建了一個帶有3個mRNA和3個miRNA的miRNA-mRNA調控網絡(圖6)，以進一步展示DEMis和DEMs之間的相互作用，有助于理解miRNA在IPF中的作用。網絡中的度數、緊密度和中間度參數由Cytoscape軟件(版本3.7.2)中的cytoHubba插件計算。調控網絡顯示miR-338-3p、miR-199a-5p和miR-299-5p可能作為關鍵miRNA在IPF發生和發展中發揮關鍵作用?？缒? 超家族成員 1(transmembrane 6 superfamily member 1，TM6SF1)是miR-388-3p的靶標mRNA；軸突生長誘向因子G1(Netrin G1，NTNG1)是miR-199a-5p的靶標mRNA；鈉通道電壓門控Ⅱ型α(sodium channel voltage-gated type Ⅱ alpha，SCN2A)是miR-299-5p靶標mRNA。這3條途徑可能與IPF的發展有關。

圖3 DEMis的轉錄因子富集

圖4 DEMs的火山圖。紅色:上調的基因;綠點:下調的基因

圖5 DEM的熱圖。紅色:高表達;綠色:低表達

表2 前40個差異表達的mRNATable 2 The first 40 differentially expressed mRNAs

2.5調節網絡中mRNA的功能富集分析 對調節網絡中mRNA的GO分析表明，CC terms在電壓門控鈉通道復合物、突觸膜的錨固成分、軸突旁節區、突觸前、突觸膜、谷氨酸能突觸、突觸膜的固有成分、突觸前膜、突觸前膜的固有成分、郎飛氏結等部位顯著富集。分子功能途徑包括細胞間黏附介體活性、鈉通道活性、電壓門控鈉通道活性及細胞內黏附介體活性(圖7)，未做出BP途徑富集。mRNAs與富集途徑之間的關系見圖8。KEGG最重要的3條路徑(圖9)包括味覺轉導、細胞黏附分子及軸突導引。mRNA與KEGG富集途徑的關系見圖10。

圖6 miRNA-mRNA網絡。紅色：上調的基因；綠色：下調的基因；mRNA用橢圓形表示，miRNA用三角形表示

圖7 IPF中DEMs富集的GO生物過程

圖8 mRNAs與富集途徑之間的關系

圖9 網絡中mRNA的KEGG途徑富集分析

圖10 KEGG中富集的mRNA之間的關系途徑

3 討 論

IPF是一種慢性纖維化性肺疾病，其特征是成纖維細胞的增殖和活化增加，包括成纖維細胞的積累、膠原蛋白的合成、細胞外基質蛋白和糖蛋白的沉積。目前該病無法治愈，具有復雜的臨床表現，已批準吡非尼酮和尼達尼布用于減緩IPF的進展。miRNA是一類非編碼的小RNA，其長度約為22個核苷酸，是基因調控中的重要調控因子。先前的研究表明，肺纖維化的發病機制與多種因素有關，包括DEMis和miRNA控制的差異基因表達[8]。因此，篩選和鑒定在IPF中差異表達的miRNA和基因，以及研究DEMis和DEMs之間的相關性，有助于闡明IPF發病的分子機制，并為臨床新藥的靶點研發提供指導。

GEO是基于R編程語言的分析工具，用于研究DEMs。本研究分析來自IPF和正常組織樣品(GSE27430和GSE135065)的miRNA表達微陣列數據，并鑒定了21種差異表達的miRNA。這些miRNA中有7個下調，14個上調。在21種已鑒定的DEMis和183種常見DEMs調控的前20條信號通路中，至少有3條與IPF的進展有關，表明21個DEMis可能與IPF的進展有關。為了更好地了解DEMis靶標mRNA的機制，本研究篩選了可能的轉錄因子。其中SP1是最常見的轉錄因子，SP1參與調控與細胞外基質生成和降解密切相關基因(MMP14、PAI等)的表達，而細胞外基質過度沉積是肺纖維化的重要特征[9]。第2個主要的轉錄因子是EGR1，EGR1是肺纖維化常見的激活轉錄因子，可與SP1競爭結合MMP14啟動子的GC序列，促進其轉錄表達[10]。而NKX6-1、POU2F1、ARID3A、MEF2A、HOXD8、HOXA9、HNF4A、RORA等轉錄因子在肺纖維化的發病機制中亦有不同的作用。由此可見，參與IPF的轉錄因子及其信號傳導途徑非常復雜，各種轉錄因子間會形成繁雜的信號網絡系統，參與調控IPF相關基因的表達。

在肺上皮細胞和成纖維細胞中，miRNA是約22個堿基的小型非編碼RNA，其可以通過靶向TGF-β信號傳導事件，α平滑肌肌動蛋白和膠原蛋白，Ⅰ型基因表達影響纖維化活性，由此引發各種信號傳導途徑和生物進展，包括Smad介導和非氨介導的途徑，分化、增殖、遷移、上皮間充質轉換和細胞外基質[11-12]。迄今為止，已在人類中鑒定出約2 000種miRNA。有研究顯示，在博來霉素誘導的肺纖維化小鼠模型中miR-199a-5p過度表達。然后在IPF患者的肺樣品中，也通過在成纖維細胞灶內原位雜交，證實了其過表達[13]。表明miR-199a-5p可能與這種疾病的病理生理學有關，特別是在肺成纖維細胞的活化及其分化為成肌纖維細胞方面。實際上，miR-199a-5p在肺成纖維細胞中的過表達誘導了這些細胞的增殖、遷移和侵襲能力的增加，以及分化成肌纖維母細胞的能力。有研究分析了IPF肺中miRNA的表達，并通過定量RT-PCR驗證了miR-299-5p的表達增加，并確定其在胚胎肺中類似地表達[14]。目前未檢索到miR-388-3p在IPF發展中的作用，值得進一步探索。NTNG1和SCN2A是miRNA-mRNA網絡中的關鍵基因。NTNG1是一種軸突導向分子，包含層粘連蛋白型EGF樣結構域，并通過GPI錨點與質膜結合[15]。SCN2A基因位于染色體2q24.3，包含26個外顯子，編碼2 005個氨基酸。SCN2A基因編碼電壓門控鈉離子通道α2亞基，主要分布在興奮性神經元軸突起始段、未髓鞘化軸突部位，在大腦中分布廣泛。到目前為止，與NTNG1相關的證據集中在精神和神經疾病上，與SCN2A相關的證據多集中在兒童癲癇[16]，而兩者與IPF之間的關聯尚不清楚，是一個值得探索的新方向。

生物學過程和信號傳導途徑的富集分析表明，本研究中183個差異表達的基因與一系列生物學過程顯著相關，例如細胞黏附和生物學黏附。ECM受體相互作用和細胞黏附是DEGs豐富的重要途徑，已被證實與IPF的調節密切相關[17-18]。這些發現表明，183個差異表達基因在IPF基因表達譜數據集中呈顯著差異表達，并且可能通過參與細胞黏附，生物黏附，ECM-受體相互作用等過程而參與IPF的進展。

本研究有一些局限性。首先應進一步解釋IPF的潛在機制，并建立基因，miRNA，轉錄因子和mRNA組成的多重網絡。此外，由于缺乏實驗驗證，需要進一步探索這些結論。