基于TCGA數(shù)據(jù)庫分析甲狀腺癌基因表達(dá)譜

趙國連，王冀邯，崔曉利

(1． 西安市胸科醫(yī)院 檢驗科，西安 710100 ；2.西北工業(yè)大學(xué) 醫(yī)學(xué)研究院，西安 710072)

甲狀腺癌(Thyroid cancer,THCA)是內(nèi)分泌系中最常見的惡性腫瘤，易受飲食、遺傳、環(huán)境等多種因素的影響[1]。近年來，中國的甲狀腺癌的發(fā)病率呈上升趨勢且女性高于男性[2]。基于甲狀腺癌術(shù)前診斷率低且晚期患者預(yù)后差的特點，探索其發(fā)病機制并尋找新型分子標(biāo)志物，對于早發(fā)現(xiàn)、早診斷、早治療具有重要意義[3]。近年來，隨著高通量測序技術(shù)及基因芯片技術(shù)的進(jìn)步，其在生命科學(xué)領(lǐng)域的應(yīng)用愈加廣泛。利用生物信息學(xué)方法在龐大的基因數(shù)據(jù)庫中篩選癌癥診斷的生物標(biāo)志物方法的有效性已經(jīng)被大量的臨床數(shù)據(jù)證實[4]。

目前已有學(xué)者[3]應(yīng)用基因表達(dá)綜合數(shù)據(jù)庫(The gene expression omnibus,GEO)對甲狀腺癌潛在的miRNAs 生物學(xué)標(biāo)志物及靶基因功能和信號通路進(jìn)行分析。Choi等通過腫瘤基因組圖譜(The Cancer Genome Atlas,TCGA)建立了一個12個基因預(yù)測模型(包括BCC8，CHI3L1，CLCNKA，F(xiàn)AM155B，GABRG1，LUM，MRO，MT1G，MT1H，SELV，SLC4A4和TMEM92)，用于預(yù)測甲狀腺乳頭狀瘤(Papillary thyroid carcinoma,PTC)中的淋巴結(jié)轉(zhuǎn)移[5]。此外，Lin等人使用腫瘤基因組圖譜(The Cancer Genome Atlas,TCGA)中與免疫相關(guān)的7個基因建立預(yù)后預(yù)測模型(包括AGTR1，CTGF，F(xiàn)AM3B，IL11，IL17C，PTH2R和SPAG11A)用于預(yù)測PTC預(yù)后情況[6]。因此，進(jìn)一步探索公共數(shù)據(jù)庫，將為尋找THCA發(fā)生發(fā)展的分子機制及挖掘疾病新型生物標(biāo)志物提供依據(jù)。本研究整合了TCGA中的THCA基因表達(dá)數(shù)據(jù)，應(yīng)用edgeR和limma兩種算法對診斷甲狀腺癌具有潛在應(yīng)用價值的基因標(biāo)志物做出預(yù)測，后續(xù)通過雙聚類分析及ROC分析進(jìn)一步驗證預(yù)測基因的可靠性。通過生物信息學(xué)分析鑒定出了11個THCA的差異表達(dá)基因(Differentially expressed genes,DEGs)及與疾病診斷相關(guān)的基因，以期為探索THCA發(fā)生發(fā)展的分子機制及挖掘疾病新型生物標(biāo)志物提供依據(jù)[2,7]。

1 資料與方法

1.1 數(shù)據(jù)來源及數(shù)據(jù)處理

通過UCSC xean網(wǎng)站下載TCGA數(shù)據(jù)庫中的甲狀腺癌基因表達(dá)數(shù)據(jù)(https://gdc.xenahubs.net/download/TCGA-THCA.htseq_counts.tsv.gz)，該數(shù)據(jù)為Log2標(biāo)準(zhǔn)化后的數(shù)據(jù)。該數(shù)據(jù)集包含了510例腫瘤樣本和58例正常對照樣本。

在UCSC xean網(wǎng)站下載THCA對應(yīng)的ID/Gene Mapping (https://gdc.xenahubs.net/download/gencode.v22.annotation.gene.probeMap)，將基因ID 與基因名稱進(jìn)行匹配，當(dāng)有多個ID對應(yīng)同一個基因名稱時，求多個ID 的平均表達(dá)值。

1.2 統(tǒng)計學(xué)分析

分別運用R/Bioconductor中的edgeR包[8]和limma包[9]對預(yù)處理過后的THCA數(shù)據(jù)提取差異表達(dá)基因。選取腫瘤與正常對照組間表達(dá)差異倍數(shù)(Fold change，F(xiàn)C)大于2，P<0.05的基因作為差異表達(dá)基因(Differentially expressed genes,DEGs)，將兩種算下的DEGs取交集。運用R 中的pheatmap包對DEGs進(jìn)行雙聚類。運用Medcalc19.0.4統(tǒng)計軟件分析，檢驗所篩選的DEGs在鑒別腫瘤樣本和正常對照樣本的應(yīng)用效果，獲取敏感性、特異性、曲線下面積等指標(biāo)。

2 結(jié)果分析

2.1 THCA中差異表達(dá)基因的篩選

首先選取腫瘤與正常對照組間倍數(shù)改變大于2，P<0.05的基因。其中，利用edgeR包得到差異基因共2 768個(上調(diào)1 765個，下調(diào)1 003個)；利用limma包得到差異基因共2 699個(上調(diào)1 080個，下調(diào)1 619個)(見圖1)。將上述兩種算法的結(jié)果求交集并去除表達(dá)趨勢不一致的基因，最終得到差異基因共1 945個(上調(diào)1 033個，下調(diào)912個)。進(jìn)一步分析顯示，隨著組間差異倍數(shù)增大，差異基因主要表現(xiàn)為在腫瘤組織中上調(diào)(見圖2)。

圖1 腫瘤組與正常對照組間DEGs火山圖Fig.1 Volcanic diagram of DEGs between tumor group and normal control group

圖2 不同倍數(shù)改變的DEGs統(tǒng)計Fig.2 DEGs statistics with different multiples

2.2 候選基因?qū)HCA的診斷效能分析

分析顯示，隨著組間差異倍數(shù)的增大，腫瘤組織中DEGs絕大部分表現(xiàn)為上調(diào)的模式，我們進(jìn)一步篩選出組間差異倍數(shù)在32倍(log2(FC)=5)以上的DEGs進(jìn)行后續(xù)分析。該11個差異基因在兩種算法中的計算結(jié)果(見表1)。對11個DEGs和樣本進(jìn)行雙聚類分析，可以看出，基于組間的DEGs表達(dá)能夠較好的將腫瘤樣本和正常對照樣本進(jìn)行區(qū)分(見圖3)。

圖3 DEGs和樣本的雙聚類分析Fig.3 Biclustering analysis of DEGs and samples注：橫坐標(biāo)為樣本(紅色代表癌癥組，藍(lán)色代表正常組)，縱坐標(biāo)為差異表達(dá)基因.

表1 篩選出的DEGs匯總Table 1 Summary of screened DEGs

進(jìn)一步對篩選出的11個候選差異基因進(jìn)行顯示，基于基因表達(dá)值鑒別腫瘤組與對照組的敏感性和特異性均在70%以上，曲線下面積均大于0.8(見圖4及表2)。提示上述基因可以較好地鑒別THCA腫瘤組和正常組。

圖4 基于候選基因鑒別腫瘤樣本與正常對照組的ROC曲線Fig.4 ROC curves of tumor samples and normal control group based on candidate genes

表2 基于候選基因鑒別腫瘤樣本與正常對照組的應(yīng)用效果Table 2 Application effects of differentiating tumor samples from normal control group based on candidate genes

3 討 論

THCA是內(nèi)分泌系統(tǒng)常見的惡性腫瘤之一，尋找潛在的分子標(biāo)志物對于臨床與科研工作至關(guān)重要。TCGA作為全球最大的癌癥基因數(shù)據(jù)庫，其大量且規(guī)范的樣本及基因表達(dá)數(shù)據(jù)為研究探索THCA的發(fā)病機制及基因標(biāo)志物提供了平臺[10]。本文基于TCGA數(shù)據(jù)庫中的THCA基因表達(dá)數(shù)據(jù)，對edgeR算法和limma算法的處理結(jié)果取交集并選擇fold change>2、P<0.05且差異表達(dá)變化趨勢一致的基因為DEGs，最終得到了1 945個DEGs。且隨著差異倍數(shù)的不斷增大，腫瘤組織中DEGs主要表現(xiàn)為表達(dá)上調(diào)的改變模式。ROC結(jié)果顯示，11個差異顯著的DEGs在鑒別腫瘤與正常組具有較好的結(jié)果。預(yù)期由這11個表達(dá)差異的DEGs組合將為TCGA的診斷、預(yù)后及復(fù)發(fā)風(fēng)險評估有一定的應(yīng)用價值。

Jin Y等人發(fā)現(xiàn)GABRB2基因在甲狀腺腫瘤組織中過度表達(dá)，通過與正常組織為對照組的隊列研究中顯示GABRB2在PCT中過表達(dá)與淋巴結(jié)轉(zhuǎn)移相關(guān)，體外實驗表明GABRB2下調(diào)會顯著抑制三種PCT細(xì)胞系的集落形成，遷徙和侵襲[11]。說明其有作為分子診斷標(biāo)志物的潛力。HMGA2是一種非組蛋白的轉(zhuǎn)錄因子，可影響包括細(xì)胞周期過程、DNA損傷修復(fù)、細(xì)胞凋亡、衰老等生物學(xué)過程。Chiappetta G 等人通過免疫組織化學(xué)和定量RT-PCR分析，認(rèn)為HMGA2表達(dá)與人類甲狀腺腫瘤中的惡性表型相關(guān)[12]。Ivanamija通過對細(xì)針穿刺甲狀腺結(jié)節(jié)中HMGA2分析認(rèn)為其可以作為區(qū)分惡性和良性甲狀腺結(jié)節(jié)的輔助生物標(biāo)志物[13]。MUC21是一種從TA3-Ha細(xì)胞中鑒定出一種新型粘蛋白。它在甲狀腺癌中通過mRNA水平和抗體結(jié)合被發(fā)現(xiàn)，但在相鄰的正常上皮中卻沒有，這就進(jìn)一步說明這種粘蛋白有用作甲狀腺癌的組織或血清標(biāo)志物[14]。SYT12有相關(guān)研究證明，SYT12在甲狀腺癌中具有一定的預(yù)后意義，SYT12可用于PCT患者的病情進(jìn)展預(yù)測的過表達(dá)與癌癥的轉(zhuǎn)移有關(guān)。但SYT12子癌癥中的分子生物學(xué)作用仍不清楚[15]。一些研究表明ZCCHC12基因與某些疾病有關(guān)，但ZCCHC12在甲狀腺癌中的功能尚未確定。Wang O 的結(jié)論證明：ZCCHC12的表達(dá)在甲狀腺癌中顯著上調(diào)，該基因過表達(dá)與淋巴結(jié)轉(zhuǎn)移相關(guān)，說明該基因具有重要的生物學(xué)功能，并有作為甲狀腺癌癥中與轉(zhuǎn)移相關(guān)的癌基因的潛在價值[16]。

Li YDENG 等研究發(fā)現(xiàn)，LIPH在甲狀腺癌組織中的高表達(dá)與淋巴結(jié)轉(zhuǎn)移密切相關(guān)，其細(xì)胞功能實驗表明，LIPH與甲狀腺癌細(xì)胞系的惡性行為呈正相關(guān)，這可以作為甲狀腺癌診斷標(biāo)志物的有力證據(jù)[17]。Jarzab B在應(yīng)用基因芯片方法對23例甲狀腺癌患者基因表達(dá)譜分析中也明確RXRG的表達(dá)有顯著升高，但是該基因在甲狀腺癌發(fā)生發(fā)展中發(fā)揮具體作用的機制還未明確[18]。

除了以上7種預(yù)測基因在甲狀腺癌中的相關(guān)報道，目前尚未有對于PRR15、SLC22A31、SLIT1和SYTL54種基因在甲狀腺癌作用機制的報道，但是SYTL5和PRR15基因表達(dá)上調(diào)在其他癌癥中的有多次報道。Wright PK等人通過免疫組化顯示SYTL5在正常乳腺導(dǎo)管上皮細(xì)胞、原位導(dǎo)管癌和浸潤性乳腺癌細(xì)胞中表達(dá)[19]。Meunier D等人研究表明 PRR15在小鼠和人類胃腸道腫瘤中高表達(dá)，可能APC蛋白的缺失有關(guān)[20]。預(yù)測的11個基因中發(fā)現(xiàn)了4個以往沒有報道與甲狀腺癌相關(guān)的基因值得進(jìn)一步研究，但是這些基因用于甲狀腺癌診斷的可靠性還有待更加深入的機制研究。

綜上，本研究通過分析TCGA甲狀腺癌表達(dá)數(shù)據(jù)，鑒定出了與THCA發(fā)生發(fā)展相關(guān)的11種生物標(biāo)志物，鑒于此，在今后的臨床研究中可以以這些顯著表達(dá)差異的基因作為藥物治療的靶向治療點。本研究不足在于缺乏更深入的機制研究，首先轉(zhuǎn)錄組學(xué)的分析并不能完全代表機體總體變化，其次，由于缺乏體內(nèi)或體外試驗，該分子預(yù)測結(jié)果還需要進(jìn)一步的臨床樣本驗證。

4 結(jié)論

分析了TCGA中的甲狀腺癌表達(dá)譜數(shù)據(jù)，鑒定出了與疾病診斷顯著相關(guān)的11個差異表達(dá)基因，并通過雙聚類分析及ROC分析進(jìn)一步驗證顯示預(yù)測基因的可靠性，這將為探索甲狀腺腫瘤發(fā)生發(fā)展機制及尋找新型分子標(biāo)志物提供依據(jù)。