程 遠 張榮花 侯曉麗 韓向陽 盧鴻健 張 青 劉志勇

1華北理工大學臨床醫(yī)學院 河北唐山 063000；2華北理工大學基礎醫(yī)學院(河北省慢性疾病重點實驗室)；3華北理工大學醫(yī)學部

microRNAs(miRNAs)是一類大小為15～21個核苷酸的非編碼RNA小分子，通過與靶基因的3'utr區(qū)結(jié)合實現(xiàn)基因表達的轉(zhuǎn)錄后調(diào)控[1]，導致細胞生理功能失調(diào)，參與癌癥的發(fā)生、發(fā)展和轉(zhuǎn)移，是多種癌癥的促癌基因或抑癌基因。越來越多的研究表明miRNA可以作為潛在的癌癥診斷標志物[2]。

肺癌是世界上發(fā)病率最高的惡性腫瘤之一，據(jù)估計中國每年有50萬人死于肺癌[3]。根據(jù)特定病例的組織學形態(tài)，肺癌可分為小細胞肺癌和非小細胞肺癌(NSCLC)，NSCLC約占肺癌的80%～85%，包括肺腺癌、肺鱗狀細胞癌和大細胞癌[4]，其中肺鱗狀細胞癌占NSCLC的30%。目前，肺鱗狀細胞癌的診斷、預后及治療藥物迫切需要新的預測性生物標記物。本研究通過對腫瘤基因組圖譜(The cancer genome atlas，TCGA)數(shù)據(jù)庫中478例肺鱗狀細胞癌miRNAs表達數(shù)據(jù)進行分析，篩選與肺鱗狀細胞癌預后相關(guān)的miRNAs 及其靶基因預測，旨在為肺鱗狀細胞癌的診斷和治療提供理論依據(jù)。

1 資料與方法

1.1數(shù)據(jù)采集 本研究收集并處理TCGA數(shù)據(jù)庫中公布的肺鱗狀細胞癌miRNAs表達譜數(shù)據(jù)523例，其中正常樣本45例、肺鱗狀細胞癌樣本478例，同時收集肺鱗狀細胞癌患者的生存和淋巴轉(zhuǎn)移分期數(shù)據(jù)。

1.2肺鱗狀細胞癌 miRNAs差異分析 使用Deseq2R包、pheatmap程序包對收集的正常肺組織和肺鱗狀細胞癌組織樣本進行差異表達miRNAs分析，并使用ggplot2 R包繪制火山圖。正常肺組織與肺鱗狀細胞癌組織差異表達倍數(shù)(Fold Change)>2，即logFC>2且P<0.05，即為差異表達有統(tǒng)計學意義的miRNAs。

1.3差異表達的miRNAs與生存率相關(guān)性分析 利用Kaplan-Meier分析肺鱗狀細胞癌組織中差異表達的miRNAs與患者生存率的相關(guān)性，進一步采用GraphPad 5.0軟件制作Kaplan-Meier生存曲線。

1.4淋巴轉(zhuǎn)移分析 使用Excel處理肺鱗狀細胞癌淋巴轉(zhuǎn)移數(shù)據(jù)，并采用Graphpad 6.0軟件繪制散點圖。

1.5miRNAs靶基因分析 利用DIANA-miRPathv3.0生物信息網(wǎng)站預測與肺鱗狀細胞癌生存率相關(guān)的miRNAs靶基因并分析其功能。運用miRDB、miRwalk、TargetScan和microRNAs.org數(shù)據(jù)庫預測與肺鱗狀細胞癌生存率、淋巴轉(zhuǎn)移相關(guān)的miRNAs靶基因，取其交集，使用FunRich 3.1.3軟件繪制韋恩圖。根據(jù)miRDB數(shù)據(jù)庫對靶基因結(jié)合程度評分，篩選評分排名較高的靶基因。

1.6miR-30a-3p候選靶基因預后價值評估 使用GEPIA在線工具評估m(xù)iR-30a-3p候選基因的臨床預后意義并繪制Kaplan-Meier生存曲線，評估肺鱗狀細胞癌患者的總生存(OS)。通過log-rank檢驗計算危險比(HR)，評估m(xù)iR-30a-3p候選基因表達對肺鱗狀細胞癌患者生存的影響。P<0.05為差異有統(tǒng)計學意義。

1.7miR-30a-3p與候選靶基因相關(guān)性分析 Target Score在80分以上的靶基因與miR-30a-3p的相關(guān)性采用Pearson相關(guān)性分析，并利用GraphPad 5.0軟件繪圖。

2 結(jié)果

2.1肺鱗狀細胞癌組織中差異表達的miRNAs 45例正常肺組織和478例肺鱗狀細胞癌組織樣本中，共發(fā)現(xiàn)差異表達的miRNAs 585個，其中 484個上調(diào)，101個下調(diào)；利用 Origin 繪制火山圖，見圖1。按照FC>2或FC<-2且P<0.05標準篩選肺鱗狀細胞癌差異表達顯著的miRNAs共201個，其中173個上調(diào)，28個下調(diào)，見表1。

圖1 肺鱗狀細胞癌組織中miRNAs表達火山圖

表1 肺鱗狀細胞癌組織中差異表達的miRNAs

hsa-miR-526b-5p3.163.12E-07hsa-miR-487a-3p2.398.27E-08hsa-miR-338-3p-2.474.72E-44hsa-miR-31613.13.55E-07hsa-miR-744-3p2.391.90E-47hsa-miR-326-2.491.27E-33hsa-miR-149-5p3.097.73E-62hsa-miR-4745-5p2.384.18E-07hsa-miR-139-5p-2.522.80E-56hsa-miR-122-5p3.051.21E-03hsa-miR-4745-3p2.366.28E-09hsa-miR-139-3p-2.833.99E-52hsa-miR-4664-3p3.051.94E-19hsa-miR-5571-5p2.329.10E-08hsa-miR-30c-2-3p-2.843.23E-65hsa-miR-518f-5p3.021.04E-03hsa-miR-380-3p2.311.53E-05hsa-miR-30a-3p-2.841.39E-66hsa-miR-675-3p31.26E-17hsa-miR-33a-5p2.34.54E-34hsa-miR-133a-3p-2.861.05E-43hsa-miR-4665-5p31.37E-14hsa-miR-3200-3p2.32.49E-23hsa-miR-144-3p-2.932.08E-37hsa-miR-224-5p2.991.14E-43hsa-miR-4640-3p2.292.51E-10hsa-miR-30a-5p-2.932.67E-79hsa-miR-6512-5p2.989.47E-12hsa-miR-4652-3p2.291.35E-06hsa-miR-338-5p-3.026.14E-77hsa-miR-7112-3p2.965.54E-15hsa-miR-545-5p2.284.36E-12hsa-miR-133b-3.174.98E-33hsa-miR-888-5p2.953.05E-04hsa-miR-548aq-5p2.282.17E-06hsa-miR-451a-3.281.25E-51hsa-miR-4664-5p2.921.88E-15hsa-miR-135b-5p2.287.58E-26hsa-miR-4732-3p-3.353.96E-34hsa-miR-96-5p2.913.20E-72hsa-miR-6814-3p2.274.49E-09hsa-miR-490-3p-3.423.11E-15hsa-miR-514b-5p2.911.07E-06hsa-miR-3144-3p2.272.36E-05hsa-miR-144-5p-3.438.60E-61hsa-miR-183-5p2.912.64E-77hsa-miR-13052.277.99E-08hsa-miR-206-3.512.31E-10hsa-miR-873-5p2.892.12E-10hsa-miR-345-5p2.278.90E-50hsa-miR-486-5p-3.626.39E-61

2.2差異表達的miRNAs與患者生存率相關(guān)性分析 分析肺鱗狀細胞癌組織中差異表達的miRNAs與患者生存率的相關(guān)性，對表1的結(jié)果按FC的絕對值進行排序，選取FC絕對值較高的miRNAs繪制Kaplan-Meier生存曲線，見圖2。篩選出下調(diào)的miR-144-5p、miR-30a-3p及上調(diào)的miR-210-3p、miR-9-5p與患者生存率相關(guān)，具有統(tǒng)計學意義(P<0.05)，見圖2。

圖2 肺鱗狀細胞癌組織中差異表達的miRNAs的Kaplan-Meier生存曲線

2.3生存率相關(guān)的miRNAs與淋巴轉(zhuǎn)移相關(guān)性分析 通過TCGA數(shù)據(jù)集庫獲取肺鱗狀細胞癌患者淋巴轉(zhuǎn)移分期數(shù)據(jù)，利用Graphpad 5.0軟件繪制散點圖，統(tǒng)計學分析結(jié)果顯示，差異表達的miR-30a-3p與肺鱗狀細胞癌患者淋巴轉(zhuǎn)移相關(guān)(P<0.05)，而miR-144-5p、miR-210-3p、miR-9-5p與淋巴轉(zhuǎn)移無明顯相關(guān)(P>0.05)，見圖3。

圖3 差異表達的miRNAs對淋巴轉(zhuǎn)移的影響注： N0：無淋巴結(jié)轉(zhuǎn)移。N1：支氣管周圍或同側(cè)肺門淋巴結(jié)轉(zhuǎn)移。N2：同側(cè)縱隔或隆突下淋巴結(jié)轉(zhuǎn)移。

2.4miRNAs靶基因分析 通過DIANA-miRPath v3.0以及TarBase v7.0數(shù)據(jù)庫預測miR-144-5p、miR-30a-3p、miR-210-3p和miR-9-5p的靶基因，進一步通過GO分析和KEGG評估預測靶基因功能和可能參與的信號通路，并繪制熱圖，見圖4。結(jié)果顯示，miR-30a-3p預測靶基因參與多種癌癥相關(guān)通路，包括細胞外基質(zhì)受體交互作用、轉(zhuǎn)化生長因子(TGF)信號通路、病毒致癌作用、河馬信號通路、肺小細胞肺癌、賴氨酸退化等。進一步通過miRDB、miRwalk、TargetScan和microRNAs.org數(shù)據(jù)庫預測miR-30a-3p靶基因并取其交集，并使用FunRich3.1.3軟件繪制韋恩圖，見圖5。篩選出53個候選靶基因，根據(jù)miRDB數(shù)據(jù)庫的評分將靶基因降序排列，見表2。

圖4 預測靶基因的KEGG通路和GO富集的熱圖注： a為KEGG通路熱圖，行代表MIRNAS，列代表PATHWAY；b為GO富集熱圖，行代表miRNAs，列代表GO項。單元格的顏色富集的P值決定。KEGG，京都基因和基因組百科全書；GO，基因本體論。

圖5 數(shù)據(jù)庫預測miR-30a-3p靶基因交集

表2 miR-30a-3p的候選靶基因

2.5miR-30a-3p候選靶基因預后價值評估 運用GEPIA數(shù)據(jù)庫分析miR-30a-3p候選靶基因與肺鱗狀細胞癌OS的相關(guān)性，結(jié)果顯示高表達USP38、SLMAP與更差的OS顯著相關(guān)(HR= 1.4、P=0.031，HR= 1.3、P=0.033)，見圖6。而其他候選靶基因生存分析無統(tǒng)計學意義，還需要進行更多的研究。

圖6 候選靶基因的Kaplan-Meier生存曲線

2.6miR-30a-3p與候選靶基因相關(guān)性分析 Pearson相關(guān)性分析顯示，Target Score在80分以上的靶基因中有3個與miR-30a-3p表達呈顯著負相關(guān)(P<0.01)：RCN2(r=0.1006)、C8orf44 (r=0.09696)、NAA25(r=0.1038)，見圖7。

圖7 MiR-30a-3p與候選靶基因相關(guān)性

3 討論

目前，癌癥相關(guān)死亡最常見的原因是肺癌，占全球癌癥相關(guān)死亡的27%以上。肺鱗狀細胞癌是一種常見的肺癌病理類型，發(fā)病率有上升的趨勢[5]。miRNA參與真核細胞重要的轉(zhuǎn)錄后調(diào)控機制[6]，越來越多的證據(jù)表明，miRNAs以其在腫瘤特異性表達中的獨特模式，成為尋找生物標記物的重要靶點[7]。篩選與肺鱗狀細胞癌預后相關(guān)的miRNAs可為臨床治療方案的優(yōu)化提供理論支持。

TCGA數(shù)據(jù)庫是目前最大規(guī)模測序結(jié)果數(shù)據(jù)庫，包含廣泛腫瘤類型及亞型、多維度的基因組變的圖譜[8]。TCGA 數(shù)據(jù)庫通過其海量數(shù)據(jù)為腫瘤研究者提供極大便利，其中蘊含的大量關(guān)于腫瘤規(guī)律仍等待挖掘[9]。

本研究收集 TCGA 數(shù)據(jù)庫的45例正常肺組織和478例肺鱗狀細胞癌組織樣本，發(fā)現(xiàn)差異表達的miRNAs 585個。按照logFC>2、P<0.05的標準篩選出肺鱗狀細胞癌差異表達顯著的miRNAs共201個，其中173個表達上調(diào)、28個表達下調(diào)。利用Kaplan-Meier繪制肺鱗狀細胞癌組織中l(wèi)ogFC值較高的差異表達miRNAs生存曲線，結(jié)果顯示差異表達的miR-144-5p、miR-30a-3p、miR-210-3p和miR-9-5p與肺鱗狀細胞癌患者的預后生存率相關(guān)。有研究證實miR-144-5p低表達與腎細胞癌患者的不良預后顯著相關(guān)[10]；miR-144-5p/ATF2軸去調(diào)控在NSCLC細胞放射敏感性中發(fā)揮了重要的作用，為NSCLC的治療提供了新的潛在靶點[11]。miR-210-3p在NSCLC組織中上調(diào)，且具有良好的敏感性和可接受的特異性，可以區(qū)分癌癥組織形成非癌性組織，可能成為NSCLC診斷的潛在生物標志物[12]；miR-210-3p也可促進口腔鱗狀細胞癌細胞增殖[13]。miR-9-5p通過下調(diào)TGFBR2的表達促進NSCLC細胞的增殖、轉(zhuǎn)移和侵襲[14]。

miR-30a-3p在癌細胞中下調(diào)[15]；下調(diào)的miR-30a-3p/5p通過抑制Wnt2和Fzd2，激活Wnt信號通路，促進食管鱗狀細胞癌細胞增殖[16]。以上研究均提示miR-144-5p、miR-30a-3p、miR-210-3p和miR-9-5p是惡性腫瘤的調(diào)控因子，與本研究TCGA數(shù)據(jù)庫篩選的分析結(jié)果一致。目前尚未見miR-30a-3p對肺鱗狀細胞癌調(diào)控作用的報道。肺癌淋巴結(jié)轉(zhuǎn)移程度是影響患者預后的重要因素，分析差異表達的miRNAs與肺鱗狀細胞癌淋巴轉(zhuǎn)移的關(guān)系，發(fā)現(xiàn)miR-30a-3p與肺鱗狀細胞癌淋巴轉(zhuǎn)移分期相關(guān)(P<0.05)，淋巴轉(zhuǎn)移的程度越高、miR-30a-3p表達越低。利用DIANA-miRPath v3.0生物信息網(wǎng)站預測并分析miR-144-5p、miR-30a-3p、-miR-210-3p和-miR-9-5p的靶基因，對每個miRNAs靶基因的功能進行富集分析，結(jié)果顯示miR-30a-3p預測靶基因參與癌癥相關(guān)通路最多，包括細胞外基質(zhì)受體交互作用、TGF信號通路、病毒致癌作用、河馬信號通路、肺小細胞肺癌、賴氨酸退化等。因此，推測miR-30a-3p在肺鱗狀細胞癌發(fā)生發(fā)展過程中發(fā)揮重要調(diào)控作用，可能成為判斷預后的新指標。另外，本研究使用miRDB、miRwalk、TargetScan和microRNAs.org數(shù)據(jù)庫預測miR-30a-3p靶基因并取其交集，篩選出53個候選靶基因。為了進一步分析miR-30a-3p與候選靶基因的靶向關(guān)系，實驗中運用GEPIA在線工具評估53個候選基因的臨床預后意義，并分析其相關(guān)性。結(jié)果顯示高表達的USP38和SLMAP與更差的OS顯著相關(guān)(P<0.05)，但其表達量與miR-30a-3p無顯著負相關(guān)；RCN2 、C8orf、NAA25和與miR-30a-3p表達呈顯著負相關(guān)(P<0.01)，但與OS無相關(guān)性。提示預后意義顯著的靶基因與miR-30a-3p無顯著相關(guān)，與miR-30a-3p顯著負相關(guān)的靶基因預后意義不顯著。考慮可能與幾個因素有關(guān)：收集臨床樣本的數(shù)量不夠充分，統(tǒng)計結(jié)果存在抽樣誤差；靶基因的功能及其與miR-30a-3p的靶定關(guān)系需要實驗進行驗證；另外，miRNAs發(fā)揮作用有一種途徑是直接抑制蛋白的表達而不影響靶基因mRNAs的水平?？傊?，生物信息學的預測結(jié)果對后續(xù)的實驗研究有一定的參考價值，但仍需大量實驗對預測結(jié)果加以驗證。

綜上所述，通過TCGA數(shù)據(jù)庫提供的肺鱗狀細胞癌 miRNAs組數(shù)據(jù)，發(fā)現(xiàn)4種miRNAs與肺鱗狀細胞癌生存率相關(guān)，其中miR-30a-3p與預后相關(guān)性較高；其機制可能是通過負調(diào)控靶基因USP38、SLMAP、RCN2 、C8orf44和NAA25參與肺鱗狀細胞癌的發(fā)生發(fā)展過程，有望成為新的預后標志物。后續(xù)我們將通過實驗驗證miR-30a-3p及其靶基因的功能，為臨床肺鱗狀細胞癌診斷和治療提供新的靶標。