基于生物信息學分析LncRNA AC005479.2在甲狀腺乳頭狀癌中的表達及相關分析

梁新科，王鵬飛，霍勝男，李建英，侯慶田

(河北省邯鄲市邯鋼醫院甲狀腺科，河北 邯鄲 056000)

甲狀腺乳頭狀癌(papillary thyroid carcinoma，PTC)是最常見的甲狀腺癌類型，占甲狀腺惡性腫瘤的70%～85%，并且全球發病率逐年增加[1]。雖然甲狀腺乳頭狀癌普遍預后良好，但其腫瘤的復發和轉移亦嚴重影響患者的生存率。近年來，基因檢測得以廣泛普及，可以研究基因表達與疾病之間的關系，為腫瘤機制的研究提供了明確的方向。本研究通過生物信息學方法篩選分析PTC與健康人群的甲狀腺組織的差異表達基因，尋找其發揮作用的相關信號通路，探索PTC的發病機制，為PTC診斷及治療提供新的研究方向。

1 資料與方法

1.1一般資料 選取TCGA數據庫中2000年1月—2020年12月收錄并具有完整隨訪信息的甲狀腺乳頭狀癌樣本510例與癌旁組織樣本58例。納入標準：①腫瘤組織樣本來自甲狀腺乳頭狀癌患者，對照組為正常組織；②樣本須同時含有基因表達譜數據與臨床信息。

1.2TCGA數據的分析 使用Rstudio中的“Limma”軟件包對TCGA數據進行差異分析，尋找甲狀腺乳頭狀癌組織與正常組織之間是否存在差異長鏈非編碼RNA(long non-coding RNA, lncRNA)，并繪制火山圖和熱圖。差異lncRNA篩選標準為：|log2FC|>1，FDR<0.05。

1.3WGCNA構建及臨床相關性分析 使用“AnnotationDbi”R包將該數據中的探針ID轉換成基因符號以進行進一步分析。使用R軟件中的“FlsahClust”軟件包進行樣本層次聚類(implementation of optimal hierarchical clustering)，過濾變異最大的25%基因。利用WGCNA軟件包中的“Pick-Soft-Threshold”函數調節參數β的權重。使用WGCNA軟件包將數據中有相關性和相鄰關系的lncRNA計算成為拓撲重疊矩陣(topological overlap matrix，TOM) ，然后計算相應的相異度(1-TOM) 。使用1-TOM作為距離度量，使用動態樹切割法進行分層聚類以識別模塊，將高度相似的模塊通過聚類標記并合并。最小基因模塊>100，合并相似模塊的閾值=0.1用于搜索在心肌甲狀腺乳頭狀癌中起重要作用的模塊。獲得基因模塊后根據模塊化特征將甲狀腺乳頭狀癌患者的臨床特征和模塊內的基因聯合分析，探究模塊內基因的生物學意義。使用相關性熱圖探究基因與模塊的相關性，并進行聚類分析。此外，WGCNA算法還對在同一模塊中基因進行相互作用網絡的預測。Cytoscape是一種分析軟件，可以為生物學家提供生物網絡分析和二維可視化。

1.4模塊中差異lncRNA的篩選 Funrich是一種生物學分析軟件。本研究中，使用維恩圖，將blue模塊中的lncRNA與差異lncRNA取交集，得到blue模塊中差異lncRNA。

1.5統計學方法 本研究在Rstudio中使用“survival”包進行單因素和多因素Cox回歸分析，進行模型變量篩選。受試者工作特征(receiver operating characteristic curve， ROC)曲線用于評估基因的診斷價值。P<0.05為差異有統計學意義。

1.6功能富集分析 選取甲狀腺乳頭狀癌樣本，將其分為AC005479.2低表達組和AC005479.2高表達組，隨后導入GSEA中進行基因本體論(Gene Ontology，GO)分析和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Geomes，KEGG)分析。

2 結 果

2.1差異lncRNA的篩選 經TCGA數據庫篩選，共獲得差異表達lncRNA58個，其中包括高表達基因36個與低表達基因22個(P<0.05)，見圖1A。LncRNA表達情況見圖1B。

圖1 甲狀腺乳頭狀癌基因表達差異A.甲狀腺乳頭狀癌基因表達差異火山圖;B.甲狀腺乳頭狀癌LncRNA表達情況Figure 1 Differences in gene expression of thyroid papillary carcinoma

2.2加權基因共表達網絡構建及臨床相關性分析 “pickSoftThreshold”函數的結果顯示，當權重參數β=12時，log(k)與log[p(k)]之間相關系數的平方>0.9。選取軟閾值β=12構建WGCNA，見圖2A。使用動態樹切割法識別基因模塊，設定模塊內最低基因個數為100，最終得到11個相應的模塊，見圖2B。結果顯示，模塊blue與甲狀腺乳頭狀癌顯著相關，見圖2C?；蚺c模塊的相關性分別展現在圖2D與圖2E中。預測出的Blue模塊中基因之間的相互作用見圖2F。

圖2 加權基因共表達網絡構建A.構建WGCNA;B.基因模塊;C.藍色模塊與甲狀腺乳頭狀癌顯著相關;D.基因與模塊的相關性;E.基因與模塊的相關性;F.藍色模塊中基因之間的相互作用Figure 2 Construction of weighted gene co-expression networks

2.3模塊中差異基因的篩選 通過維恩圖將得到甲狀腺乳頭狀癌組織與正常組織之間的58個差異lncRNA與blue模塊中174個模塊lncRNA取交集，得到藍色模塊內差異lncRNA共有34個，見圖3。

圖3 藍色模塊中的差異lncRNAFigure 3 Differential lncRNA in blue module

2.4預后模型的構建與驗證 對藍色模塊內34個差異lncRNA進行單因素和多因素Cox回歸分析(圖4A、B)，得到AC005479.2。隨后對AC005479.2繪制ROC曲線，曲線下面積(area under curve，AUC)值為0.838(圖4C)。

圖4 差異lncRNA的回歸分析A.單因素Cox回歸分析;B.多因素Cox回歸分析;C.ROC曲線Figure 4 Regression analysis of differential lncRNAs

2.5基因富集分析 通過GO和KEGG分析，揭示AC005479.2的潛在通路與功能，結果顯示AC005479.2通過脂肪細胞因子信號途徑、絲裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)信號通路、鈣黏蛋白結合參與細胞黏附、對凋亡信號通路的正調控等功能參與甲狀腺乳頭狀癌的病理過程，見圖5。

圖5 富集分析顯示AC005479.2的潛在通路與功能A.AC005479.2低表達組的GO分析;B.AC005479.2高表達組的GO分析;C.AC005479.2低表達組的KEGG分析;D.AC005479.2高表達組的KEGG分析Figure 5 Enrichment analysis showing the potential pathways and functions of AC005479.2

3 討 論

甲狀腺癌是最常見的內分泌腫瘤，占人類所有惡性腫瘤的0.5%～1%[2]。PTC被認為是復發率最高的甲狀腺癌類型，但同時也被認為是具有良好預后的惰性癌癥，占已分化甲狀腺癌的70%～90%[3]。對于大多數PTC患者，預后良好，然而，在某些情況下會出現侵襲性行為，30%～90%的患者表現出臨床或隱匿性頸部淋巴結轉移[4-5]。lncRNA已被證實參與PTC的細胞遷移和增殖過程[6]。因此，尋找腫瘤與正常組織之間的差異lncRNA有助于進一步了解甲狀腺乳頭狀癌的發病機制，為PTC的術前診斷提供生物標志物和治療靶點。

基于生物信息學分析，研究顯示，lncRNA AC005479.2在甲狀腺乳頭狀癌中表達與癌旁組織差異有統計學意義(P<0.05)，同時GO和KEGG結果顯示AC005479.2可能通過脂肪細胞因子信號途徑、MAPK信號通路、鈣黏蛋白結合參與細胞黏附、對凋亡信號通路的正調控等功能參與甲狀腺乳頭狀癌的病理過程。

近年來研究顯示，肥胖患病率的增加與包括乳腺癌、結腸癌和前列腺癌在內的各種人類癌癥的發病率、進展和預后較差有關[7]。最近，一些研究顯示甲狀腺癌的發病率隨著肥胖率的顯著上升而增加，并且越來越多的證據表明肥胖與甲狀腺癌風險增加之間存在關聯[8-9]。脂肪因子是脂肪細胞分泌的各種酶、激素、細胞因子、生長因子、蛋白質和其他生物活性物質，包括脂聯素、瘦素、抵抗素和白細胞介素等。據報道，體重指數越高的患者脂肪因子濃度越高[10]。由于營養過?；蜻\動不足所致肥胖時，脂肪因子將出現異常分泌，并參與多種病態的發生。脂肪因子參與以下病理和生理過程，如胰島素敏感性、炎癥和增殖，這些在腫瘤發生和發展過程中均發揮重要作用。因此，脂肪因子可能是將肥胖與甲狀腺癌聯系起來的目標之一。本研究顯示，AC005479.2可能通過脂肪因子參與了甲狀腺乳頭狀癌的發生發展過程。

MAPK通路代表一種普遍存在的信號轉導通路，在細胞活力、凋亡和侵襲的調控中發揮重要作用。MAPK信號通路能夠被受體酪氨酸激酶、G蛋白偶聯受體和部分細胞因子受體激活，將細胞外的信號刺激傳遞給細胞，引起細胞增殖、凋亡、分化等一系列反應，進而促進腫瘤的發生發展。先前的研究表明，MAPK信號通路在多種癌癥中被激活，包括胃癌[11]、肺癌[12]、卵巢癌[13]和肝癌[14]。這一途徑已被證明在腫瘤發生發展中發揮關鍵作用，因此，成為包括甲狀腺癌在內的惡性腫瘤的重要治療靶點[15]。MAPK屬于絲氨酸/蘇氨酸蛋白激酶，MAPK信號級聯是腫瘤細胞增殖、分化、凋亡和耐藥性的關鍵途徑。此外，MAPK激活可以磷酸化核轉錄因子和蛋白激酶，調節相關基因的轉錄并參與各種生理過程[16]。MAPK信號通路由4個子通路組成：extracellular signal-regulated kinase(ERK)、c-Jun N-terminal kinase(JNK)、p38 和big MAPKinase(BMK)。研究表明，ERK信號通路的激活促進細胞生長、分化、遷移和存活[17]。JNK和p38的激活涉及細胞分化、凋亡和細胞存活的變化[18]。研究顯示，MAPK通路中，BRAF的突變及RET/ PTC重排能夠促進甲狀腺濾泡細胞向乳頭狀甲狀腺癌的轉化，這種改變被認為是甲狀腺癌發生和進展的早期表現[19]。據Nucera等[20]報道，根據甲狀腺癌的遺傳學變異，有28%～86%的患者伴隨BRAF基因突變。其中BRAF V600E突變占甲狀腺乳頭狀癌的45%[19]。Zhang等[21]研究顯示BRAF基因激活的非編碼RNA(BRAF-activated non-coding RNA，BANCR)通過激活ERK-MAPK信號通路進而對乳頭狀甲狀腺癌細胞增殖和遷移產生影響。

MAPK還可通過調節上皮-間質轉化(epithelial mesenchymal transformation EMT)信號通路發揮作用，上皮-間質轉化被認為是癌細胞從原發部位向遠處組織或淋巴結轉移過程中的一個重要步驟，通過該過程上皮細胞逐漸獲得間充質特征并增強遷移和侵襲的特性。甲狀腺乳頭狀癌的多灶性病變顯示出淋巴結轉移傾向[22]，PTC的淋巴結轉移導致預后不良的結果，研究顯示EMT參與了PTC細胞的局部轉移過程。EMT的分子特征包括將極化上皮表型標志物(如E-cadherin)轉換為間充質標志物(如N-cadherin、纖連蛋白等)。轉錄因子(如Snail、Twist、ZEB1和ZEB2等)的激活觸發EMT的發生。據報道，MAPK/ERK的激活會增加Snail、Slug和N-鈣黏蛋白的水平，并降低E-鈣黏蛋白的表達水平，從而導致EMT[23-25]。EMT通路已被廣泛證明在腫瘤進展和轉移中起關鍵作用，可作為甲狀腺癌的臨床治療靶點[15,26]。本研究顯示AC005479.2可能通過MAPK通路，改變細胞黏附狀態，形成EMT環境，促進PTC的淋巴結轉移，導致PTC的進展。

綜上所述，本研究基于TCGA甲狀腺乳頭狀癌及癌旁組織的基因數據，顯示AC005479.2在甲狀腺乳頭狀癌和癌旁組織中的表達差異有統計學方法意義，而且通過GO和KEGG分析其可能通過MAPK通路、脂肪因子通路，調節EMT及凋亡的發生，從而參與甲狀腺乳頭狀癌的發生發展過程，為驗證AC005479.2的作用提供理論依據。