喉鱗狀細胞癌的轉錄組數據分析及生物標志物篩選

王娜娜，楊川

1延安大學附屬醫院耳鼻咽喉科，陜西延安716000

2河北醫科大學第二醫院耳鼻咽喉科，石家莊050000

3平昌縣人民醫院耳鼻咽喉科，四川巴中636040

喉癌是耳鼻咽喉-頭頸外科第二高發腫瘤，其中喉鱗狀細胞癌占全部病理分型的95%。喉鱗狀細胞癌的病因復雜，近年來其發病率不斷增加，患者的晚期生存率仍較低。隨著機器人輔助手術在頭頸部腫瘤治療中的開展及食管、氣管一期重建技術的不斷成熟，中國臨床腫瘤學會（Chinese Society of Clinical Oncology，CSCO）頭頸部腫瘤診療指南也強調手術治療對提高生存率的重要性，但仍面臨術后功能障礙、生活質量下降等嚴峻挑戰。另一方面，頭頸部腫瘤手術仍被認為是外科“皇冠”，手術水平直接關系到患者的生存預后。因此，進一步闡明喉鱗狀細胞癌的分子機制，尋找潛在的生物學靶點，對喉鱗狀細胞癌進行早期診斷和治療仍十分必要。信使RNA（messenger RNA，mRNA）、微小RNA（microRNA，miRNA）、長鏈非編碼RNA（long noncoding RNA，lncRNA）在喉鱗狀細胞癌的發生、發展中發揮重要作用。本研究通過腫瘤基因圖譜（The Cancer Genome Atlas，TCGA）轉錄組數據挖掘，分別闡明lncRNA、miRNA和mRNA在喉鱗狀細胞癌中的作用機制及其潛在聯系，現報道如下。

1 資料與方法

1.1 數據獲取及差異分析

在TCGA中下載喉鱗狀細胞癌轉錄組數據和相關臨床信息，其中腫瘤樣本111例，正常樣本12例。采用“edgR”包進行差異分析，設|logFC|﹥2，錯誤發現率（false discovery rate，FDR）﹤0.01，分別得到差異lncRNA、差異mRNA和差異miRNA，采用R軟件繪制可視化圖形并進行統計分析。

1.2 蛋白相互作用網絡分析

采用STRING11.0（https://string-db.org/cgi/input.pl）在線工具將喉鱗狀細胞癌Top500差異mRNA構建蛋白相互作用（protein-protein interaction，PPI）網絡，相互作用值=0.7；Cytoscape插件CentiScape2.2、cytoHubba0.1用于尋找關鍵基因，并進行可視化展示。

1.3 功能富集和通路富集分析

喉鱗狀細胞癌差異mRNA用Clusterprofiler包進行細胞組分（cellular component，CC）、生物學過程（biological process，BP）、分子功能（molecular function，MF）的基因文本（gene ontology，GO）富集及京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）信號通路富集。

1.4 MYHAS的共表達mRNA（cor-mRNA）

采用Pearson方法計算cor-mRNA，并取相關系數cor≥0.4。

1.5 構建競爭性內源RNA（competitive endogenous RNA，ceRNA）網絡和轉錄因子網絡

ceRNA網絡：采用本地版miRanda和Targetscan分別預測lncRNA結合的miRNA，并與差異miRNA取交集。再通過miRTarBase、Targetscan、miRBD分別預測miRNA結合的mRNA，并與差異mRNA取交集，建立lncRNA-miRNA-mRNA網絡聯系，并通過Cytoscape3.7.2軟件進行可視化展示。

轉錄因子網絡：將鈣信號通路相關的mRNA于DAVID6.8在線網站（https://david.ncifcrf.gov/）中進行轉錄因子富集，建立lncRNA相關mRNA的轉錄因子聯系，并通過Cytoscape3.7.2軟件進行可視化展示。

2 結果

2.1 一般資料

腫瘤患者的年齡、腫瘤T分期、病理分級、飲酒、吸煙情況對腫瘤患者總體生存率的影響不明顯。遠處轉移對患者總生存率的影響顯著，臨床轉移患者的5年生存率為0%，明顯低于無臨床轉移患者的50.33%，差異有統計學意義（χ=12.230，P=0.0005）；病理轉移患者的5年生存率為0%，明顯低于無病理轉移患者的44.83%，差異有統計學意義（χ=18.260，P=0.0001）。女性患者的5年生存率為30.00%，低于男性患者的50.11%，差異有統計學意義（χ=4.454，P=0.0348）。發生病理淋巴結轉移患者的5年生存率為48.08%，低于未發生病理淋巴結轉移患者的70.14%，差異有統計學意義（χ=4.539，P=0.0331），但臨床上發生淋巴結轉移與未發生淋巴結轉移患者的5年生存率比較，差異無統計學意義（χ=1.150，P=0.2835）。

2.2 lncRNA、miRNA、mRNA差異表達情況

lncRNA、miRNA和mRNA各自的主成分分析結果顯示，腫瘤組與對照組間表達譜數據相互獨立，組內重復性較好，數據可信度高（圖1A、1C和1E）。當|logFC|﹥2，FDR﹤0.01時，得到的差異lncRNA共687個，其中上調561個，下調126個（圖1B）；得到的差異miRNA共54個，其中上調28個，下調26個（圖1D）；得到的差異mRNA共1818個，其中上調950個，下調868個（圖1F）。

圖1 lncRNA、miRNA、mRNA的主成分分析、火山圖

2.3 差異lncRNA和miRNA表達的喉鱗狀細胞癌預后的單因素分析

根據RNA-Seq臨床數據計算了差異lncRNA表達的喉鱗狀細胞癌預后的單因素分析結果，取

P

﹤0.05，共有64個lncRNA滿足條件，根據是否能在NCBI中查到相關lncRNA信息，總生存率是否有統計學意義，篩選了11個lncRNA，其中

LINC00278

、

LINC00689

、

MYHAS

屬于低風險基因（HR﹤1），

MNX1-AS1

、

LINC02575

、

HOXB-AS4

、

LSAMP-AS1

、

LINC02086

、

IGFL2-AS1

、

LINC02253

、

CASC20

屬于高風險基因（HR﹥1）。同理，計算出了差異miRNA表達的喉鱗狀細胞癌預后的單因素分析結果，

P

﹤0.05的miRNA共7個，其中miRNA-383、miRNA-196a-2、miRNA-196a-1、miRNA-100、miRNA-4652是高風險基因（HR﹥1），miRNA-99a、miRNA-301a是低風險基因（HR﹤1）。（表1）

表1 差異lncRNA和miRNA表達的喉鱗狀細胞癌預后的單因素分析

2.4 差異mRNA功能分析和關鍵基因篩選

將1818個差異mRNA用Clusterprofiler進行GO分析和KEGG信號通路分析，取

P

值或

P.adjust

值﹤0.05。分子功能（MF）以通道活性和激活跨膜轉運活性為主，細胞組分（CC）以細胞外基質（extracelluar matrix，ECM）包含的膠原纖維和轉運復合體為主，生物學過程（BP）以細胞外基質形成和調節離子跨膜轉運為主（圖2A～2C）。KEGG信號通路中以細胞因子-細胞因子受體間相互作用和鈣信號通路為主要作用通路，其中也包括環磷酸腺苷（cyclic adenosine monophosphate，cAMP）信號通路、環磷酸鳥苷（cyclic guanosine monophosphate，cGMP）-cGMP依賴性蛋白激酶（cGMP-dependent protein kinase，PKG）信號通路和ECM受體相互作用等腫瘤發生、發展的經典信號通路（圖2D）。將差異倍數前500的mRNA進行String在線分析，構建PPI網絡，并通過Centiscape插件進行關鍵基因篩選，根據Degree進行篩選，綠色﹤10，青色≥10，紫色≥15，紅色≥20，其中

ACTN2

、

ACTN3

、

MYL1

、

MYL2

、

MYL3

、

MYH6

、

TTN

網絡連接最多，屬于關鍵基因（圖2E）。

圖2 差異mRNA的GO、KEGG信號通路富集和PPI網絡

2.5 MYHAS的cor-mRNA、關鍵基因篩選及功能富集

MYHAS的共表達mRNA（cor≥0.4）共253個，采用插件cytoHubba分別計算Closeness、Degree、EPC、MCC、MNC和Radiality 6種算法下前10的基因并取交集，得到關鍵基因：

TNN

、

MYL1

、

MYL3

、

ACTN2

、

MYH6

。GO富集顯示cor-mRNA主要參與調節金屬離子轉運、跨膜轉運復合體和肌動蛋白結合等（圖3A～3C）。信號通路分析顯示，cormRNA主要富集在鈣信號通路、cGMP-PKG信號通路和催產素信號通路等（圖3D）。其中鈣離子信號通路中富集的基因分別是PRKACA、CAMK2A、SLC8A3、CAMK2B、RYR1、PLN、CASQ2、PHKG1、HRC、MYLK3、TNNC1、CACNA1S、MYLK2、CASQ1、ATP2A1、TNNC2。

圖3 MYHAS 的cor-mRNA的GO和KEGG信號通路富集

2.6 ceRNA網絡和轉錄因子網絡

將單因素分析與生存分析中P﹤0.05的lncRNA取交集，得到21個lncRNA，并分別進行mi-Randa和Targetscan靶基因預測并取交集，得到2363個與lncRNA可以結合的miRNA，然后將預測的miRNA與差異miRNA的成熟體取交集，得到66個miRNA。對66個miRNA分別進行miRTarBase、Targetscan和MiRDB靶基因預測后與差異mRNA取交集，并構建ceRNA網絡（圖4A）。ceRNA網絡中有11個lncRNA（青色菱形），分別為LINC02576、LINC02086、AC020659.1、LINC00528、LINC00689、HOXB-AS4、MNX1-AS1、LINC00278、AC010624.1、AC016773.1、MYHAS；5個 miRNA（綠色三角形），分別為has-miRNA-206、has-miRNA-573、has-miRNA-3662、has-miRNA-133b、hasmiRNA-449a；12個mRNA（其中7個紅色6邊型為上調mRNA，5個藍色圓形為下調mRNA），分別為STC2、PAX3、NETO2、EIF5A2、ADPRHL1、SYNM、ACTA1、GPR156、CLIC5、BMP3、HNF4A、FOXL2。將cor-mRNA鈣信號通路中的16個基因進行轉錄因子富集，并繪制轉錄因子調控網絡，網絡中綠色橢圓形代表mRNA，藍色菱形代表轉錄因子，這些轉錄因 子 分 別 是 MYOD、RSRFC4、TAL1BETAITF2、YY1、P300、PAX3、PAX5、ZID和OLF1（圖4B）。

圖4 ceRNA網絡和轉錄因子網絡

2.7 DNA甲基化分析

為了進一步闡明關鍵基因在喉鱗狀細胞癌中的調控作用，進一步分析了TCGA和GEO中的喉鱗狀細胞癌DNA甲基化數據。設Δβ絕對值≥0.4，P﹤0.05，取二者交集，共得到75個甲基化位點；其中cg15700197和cg17892556分別為甲基化差異最大的位點，cg15700197在正常組織中高甲基化，在喉鱗狀細胞癌組織中低甲基化，cg17892556在正常組織中低甲基化，在喉鱗狀細胞癌組織中高甲基化。這75個差異甲基化位點主要位于1stExon、TSS1500和TSS200，其次為5'非翻譯區（untranslated region，UTR）和基因體區，3'非翻譯區和基因間隔區（intergenic region，IGR）分布較少。上述位點對應高甲基化基因56個，低甲基化基因15個，ZNF625是Δβ值最大的高甲基化基因，OR10J1是Δβ值最大的低甲基化基因。

3 討論

喉鱗狀細胞癌是頭頸部常見的腫瘤之一，病因較多，影響因素復雜。本研究針對喉鱗狀細胞癌患者的一般情況進行分析，結果顯示，患者年齡、T分期、病理分級、飲酒、吸煙情況對總生存率的影響不明顯，而遠處轉移、淋巴結轉移會降低患者的總生存率。關于腫瘤分期（T分期）和病理分級對腫瘤預后無影響，可能與目前分期系統的局限性有關；因此針對喉鱗狀細胞癌的分期，根據腫瘤萌芽活性（tumor budding activity，BA）和細胞巢大?。╟ell nest size，CNS）進行分期對預后似乎更具指導意義。本研究結果中淋巴結浸潤對喉鱗狀細胞癌預后的影響與Ho等研究結果一致，尤其是病理組織上的淋巴結轉移對指導臨床治療意義重大，當然這也可能與從臨床到病理診斷的過程增加了“開始治療時間（time to treatment initiation，TTI）”，從而降低了預后有關。

近年來有關lncRNA在喉鱗狀細胞癌中作用的研究較多，尤其是miRNA“海綿”靶基因是喉癌發生、發展的重要機制之一。作為ceRNA，lncRNA通過靶向結合miRNA間接調控靶基因，是較好的生物標志物和潛在的治療靶點。生存分析顯示，miRNA-210、miRNA-455和miRNA-4326的表達量與喉鱗狀細胞癌預后相關（數據未展示），其中有關miRNA-210在肺癌、膠質瘤、口腔癌、骨肉瘤等惡性腫瘤中的作用已廣泛報道，并且miRNA-210在組織、血清和外泌體中均可存在。miRNA-455也可通過自噬、p21活化蛋白激酶2（p21-activated kinase 2，PAK2）軸和 Janus激酶 1（Janus kinase 1，JAK1）軸等作用方式抑制腫瘤生長，并且miRNA-455上述功能的發揮有賴于不同lncRNA的作用，這些lncRNA作為ceRNA競爭性結合miRNA-455；當然，也有較多miRNA-455單獨發揮作用的報道。同時，Xu等亦報道miRNA-4326通過抑制APC2基因表達促進肺癌細胞增殖。但遺憾的是，目前尚未見上述3個miRNA在喉鱗狀細胞癌中被報道，有待于進一步驗證。

在具有抑癌作用的lncRNA中，LINC00689、LINC00278已分別在膠質瘤和食管鱗狀細胞癌中報道，尤其是LINC00278編碼的多肽YY1BM，可抑制YY1與雄激素受體結合從而降低真核生物延伸因子2激酶（eukaryotic elongation factor 2 kinase，eEF2K）的表達，促進腫瘤細胞凋亡。喉結是第二性征標志，雄激素受體高表達，由此可推斷LINC00278在喉鱗狀細胞癌中是一個較好的潛在靶點。在具有促癌作用的lncRNA中，IGFL2-AS1、LSAMP-AS1和MNX1-AS1已分別在胃癌、前列腺癌和食管鱗狀細胞癌中報道，為下一步喉鱗狀細胞癌的研究提供了較好的參考意義。

基于目前報道的有關lncRNA-miRNA-mRNA的相互作用關系，本研究利用差異lncRNA、差異miRNA和差異mRNA構建了ceRNA網絡，通過該網絡發現，網絡中miRNA既無生存意義，又無單因素意義；網絡中mRNA亦不屬于Centiscape中的關鍵基因；同時，與miRTarBase疾病搜索中喉鱗狀細胞癌的miRNA無交集。因此認為lncRNA并不能完全主導miRNA和mRNA發揮作用，lncRNA與miRNA在喉鱗狀細胞癌發病過程中各自存在獨立的作用機制。為了進一步說明lncRNA的功能，本研究選擇了MYHAS作為代表，通過分析MYHAS共表達基因cor-mRNA，發現cor-mRNA得到的5個關鍵基因全部與差異mRNA得到的關鍵基因吻合；并且cor-mRNA的通路分析中主要信號通路鈣信號通路、cGMP-PKG信號通路和cAMP信號通路亦是差異mRNA的主要信號通路。有關cAMP、cGMP-PKG和鈣信號通路在不同腫瘤中的調控作用已有很多報道，有研究報道了cGMP-PKG在頭頸部鱗狀細胞癌中的作用，認為cGMP-PKG在頭頸部鱗狀細胞癌中是較好的治療靶點，并且全球第一個鳥苷酸環化酶激動劑利那洛肽已于2019年初在中國上市。此外，Park和Juhnn亦報道cAMP信號激活可通過增加去乙酰化酶8的表達促進順鉑對肺癌細胞的凋亡作用。而關于鈣信號通路在腫瘤中的作用，Raynal等研究表明，鈣信號通路中一系列激酶相互作用后，引起鈣/鈣調素依賴蛋白激酶Ⅱ（calcium/calmodulin-dependent protein kinaseⅡ，CAMKⅡ）激活，使甲基化結合域蛋白甲基化CPG結合蛋白2（methyl-CPG binding protein 2，MeCP2）出核，導致基因組發生去甲基化，使已經甲基化的抑癌基因重新激活，以發揮抑癌作用。MYHAS屬于抑癌基因，其通過相關基因cormRNA在喉鱗狀細胞癌組織中參與鈣信號通路的失活，從而引起下游靶基因高甲基化。基于上述結果，本研究又分析了數據庫中喉鱗狀細胞癌的DNA甲基化數據，共篩選出了71個基因；本實驗室已經對ZNF667基因進行了試驗驗證，其在喉鱗狀細胞癌組織和細胞中均呈高甲基化、低表達，發揮抑制喉鱗狀細胞癌細胞增殖、侵襲的作用。由此認為MYHAS可以較好地代表喉鱗狀細胞癌中的差異lncRNA，在喉鱗狀細胞癌中是一個較好的生物標志物和治療靶點，具有潛在研究價值。

同時，本研究針對cor-mRNA中參與鈣信號通路的16個基因構建了轉錄因子調控網絡，擬進一步說明非編碼RNA在喉鱗狀細胞癌信號通路中的上下游調控關系。這些轉錄因子與MYHAS相互作用，直接在轉錄水平調控靶基因表達。網絡中PAX5是B淋巴細胞特異性轉錄因子，在B細胞淋巴瘤的發生、發展中具有重要作用，但也見于在其他惡性腫瘤中的報道，并且Brazao等報道了PAX5對不同階段B細胞中lncRNA的調控作用。而轉錄因子PAX3的功能較為復雜，這是因為PAX3具有7個可變剪切體，不同剪切體在腫瘤中表現的功能不同。P300被認為是轉錄共激活因子，是一種組蛋白乙酰轉移酶，通常被招募到轉錄增強子中，并通過乙酰化染色質來調節基因表達。Liang等研究報道，P300可在食管鱗狀細胞癌中調控LINC00460的表達，這亦為喉鱗狀細胞癌的研究提供了可行性。轉錄因子YY1是研究較多的腫瘤相關轉錄因子之一，其可調控lncRNA的表達，引起直腸癌、惡性膠質瘤、乳腺癌和肝癌等腫瘤的發生、發展，其中有lncRNANPCCAT1通過上調YY1促進鼻咽癌發生發展的報道，這對研究YY1在喉鱗狀細胞癌中的作用具有重要參考意義。MYOD屬于肌源性轉錄因子，在橫紋肌肉瘤中與同類轉錄因子MYF5共同發揮重要作用；此外Dodd等亦報道了MYOD通過直接結合miRNA-182的啟動子增加miRNA-182的表達，從而發揮促腫瘤作用。因此，在喉鱗狀細胞癌中轉錄因子MYOD與lncRNA的作用還是值得期待。

綜上所述，公共數據分析顯示喉鱗狀細胞癌中存在眾多差異lncRNA和miRNA，其中lncRNA MYHAS通過抑制鈣信號通路引起下游靶基因甲基化可能是喉鱗狀細胞癌發生、發展的主要機制之一，lncRNA MYHAS是一個較好的生物標志物和治療靶點，具有潛在研究價值。