KRAS突變對人肺腺癌轉錄組影響的生物信息學分析

，， ，， ，，，，

肺癌是我國常見惡性腫瘤之一，其發病率和病死率在惡性腫瘤中均居首位。非小細胞肺癌(non-small cell lung cancer,NSCLC)是肺癌最常見的類型，其中肺腺癌(lung adenocarcinoma,LA)是主要的NSCLC類型，超過40%的患者發現時已屬晚期[1]。表皮生長因子受體(epidermal growth factor receptor，EGFR)、間變淋巴瘤激酶(anaplastic lymphoma kinase,ALK)和鼠肉瘤病毒致癌基因(Kirsten rat sarcoma viral oncogene,KRAS)基因突變是LA中常見的突變類型。近年來，隨著分子靶向治療技術的不斷進步，以吉非替尼、克唑替尼為代表的靶向藥，為EGFR、ALK基因突變陽性的LA患者帶來了新的希望。盡管如此，依然有10%～30% KRAS突變型LA患者無法從現有治療方案中獲益[2]。到目前為止，KRAS基因突變所產生的耐藥機制仍然不明確。

隨著二代測序和大數據分析方法的不斷進步，通過全基因組轉錄組測序獲得的表達譜能鑒別腫瘤組織與正常組織的差異表達基因[3]。其中，以基因表達數據庫(gene expression omnibus，GEO)的數據量和類型最為豐富，提供了各種類型的轉錄組信息[4-5]。本研究擬采用GEO數據庫中的NSCLC基因芯片數據篩選出KRAS突變型NSCLC差異表達基因，并對差異表達基因進行生物信息學分析，預測和篩選出KRAS突變型LA的轉錄調控網絡，為進一步研究KRAS突變型LA的耐藥機制、早期診斷標志和潛在分子靶點提供理論依據。

1 資料和方法

1.1 資料

1.1.1 數據來源 以“KRAS”AND“lung cancer”和“KRAS”AND“lung”為檢索式在GEO公共數據庫進行檢索，查看所有“series”。具體數據篩選標準如下：①物種為人；②至少2個生物重復；③實驗設計思路清晰、數據質量好。在詳細查看所有數據的注釋文件，確定數據集中樣本類型及數量后進行分析。本研究使用的肺腺癌KRAS突變數據的具體描述如表1，其中，GSE是指一個實驗項目中的芯片實驗編號，GPL是芯片平臺。本研究所使用的數據集，主要為Affymetrix芯片數據。

表1 肺腺癌數據的基本信息

1.1.2 實驗設計 本研究所使用的數據集均為經GEO數據庫預處理后的轉錄組表達譜矩陣。根據數據來源實驗描述以及各樣本類型，將數據類型設定為對照與突變兩組，即：在未區分具體人種的情況下，正常的肺組織設為對照組，僅具有KRAS突變的LA組織設為突變組。進行差異分析并進行后續的分析研究。

1.2 方法

1.2.1 數據處理及差異表達分析 在Bioconductor(http://www.bioconductor.org/)網站下載生物信息分析的R語言程序包[10]。利用Bioconductor中的Impute程序包對已獲取表達譜進行歸一化處理，使各個樣本的數據歸一化[11]。同時利用Bioconductor中的注釋包對數據進行注釋，將探針對應到基因上，如果多個探針對應一個基因，那么用所有探針表達值的平均值代表這個基因的表達值。得到3組數據集的表達譜用于后續分析。隨后利用Bionconductor中的Limma程序包對表達譜篩選KRAS突變LA樣本與正常樣本間的差異基因。首先使用P<0.05初步篩選差異基因。而后利用Bonferroni算法來校正假發現率(FDR)的原始P值并計算倍數變化(FC)[12]。本研究中，差異基因的篩選標準是FDR<0.05，Log2FC>1。3組數據集的差異結果保存用于后續分析。

1.2.2 篩選顯著性基因 通過上述分析得到的3組差異基因，隨后利用R語言中的RRA算法篩選3組差異基因中的顯著性基因。RRA方法主要是對不同數據集獲得的基因進行排序。如果一個基因在所有實驗中表達倍數越高，且P值越小，則其越顯著[13]。

1.2.3 構建顯著性基因的蛋白質互作網絡 利用STRINGv9.1(http://www.string-db.org/)開源數據庫構建蛋白質互作(protein-protein interaction,PPI)網絡[14]。為更好地理解重要基因的相互作用，使用STRING數據庫預測其編碼蛋白的PPI網絡，本研究使用的可靠性閾值為>0.4。Cytoscape軟件(http://cytoscape.org/)用于構建PPI網絡[15]。在PPI網絡中，每個節點代表基因、蛋白質或其他分子，節點之間的連線代表生物分子之間的相互作用。PPI網絡可用于識別由EGFR突變型LA中顯著基因編碼的蛋白質之間的相互作用和通路關系。核心節點，即具有重要生物功能的蛋白質通過計算蛋白質之間的連線數以及每個節點與特定節點有無直接連接來確定，并給出排名。本研究中利用R語言設計腳本篩選核心節點，結果輸出前10名基因。

1.2.4 富集分析 DAVID(The Database for Annotation，Visualization and Integrated Discovery，http://david.abcc.ncifcrf.gov/)[16]，是生物信息學常用的聚類分析數據庫，其中包含GO(Gene Ontology)[17]和KEGG pathway(The Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/pathway.html)[18]分析。本研究對出現在PPI網絡中的56個顯著性基因進行聚類，聚類結果以P<0.05為閾值進行篩選，最后利用R語言對結果進行可視化展示。

2 結果

2.1 差異表達分析及顯著性基因篩選 3個數據集GSE31210、GSE32863、GSE75037分別篩選出差異表達基因2 625個、902個、2 581個。其中上調表達的基因分別為1 115個、351個、1 191個，下調表達的基因分別為1 510個、551個、1 390個(圖1A、圖1B、圖1C)。利用RRA算法對3組差異結果綜合統計得到顯著性基因共120個，其中56個基因表達上調，64個基因表達下調，部分顯著性基因結果，見圖1D。

注：FDR<0.05，Log2FC>1圖1 各數據集差異表達分析及顯著性基因

2.2 顯著性基因蛋白質相互作用調控網絡構建 用STRING在線軟件對顯著性基因進行PPI網絡構建，其中有56個蛋白出現在PPI網絡中，顯示出與其他蛋白的相互作用關系，基因互作網絡見圖2A。而后對互作網絡進一步分析，依據基因連線數目篩選核心基因，統計結果圖形化展示見圖2B。

圖2 蛋白質互作網絡及核心基因

2.3 蛋白質互作網絡中出現的顯著性基因的富集分析 對120個顯著性基因進行富集分析，GO分析結果顯示，這些基因的富集情況分為3種類型，包括生物學進程、細胞成分、分子功能，圖3。其中，在生物學進程方面，這些基因大多與血管生成、免疫反應、細胞間信號轉導以及細胞粘附有關。在細胞成分方面，大多與藥物吸收、藥物結合以及耐藥有關。在分子功能方面，大多與細胞對腫瘤壞死因子的反應、T細胞遷移及細胞對白介素1的反應有關。

注：橫坐標代表GOID，縱坐標代表富集在各GO的基因數目圖3 GO富集分析結果

KEGG分析結果顯示這些顯著性基因主要富集在免疫相關通路、各種感染性疾病以及細胞周期、凋亡、轉錄調控等細胞增殖相關通路和藥物代謝相關通路上，圖4。

注：分析結果橫坐標代表富集在各通路的基因數目，縱坐標代表通路名稱，P<0.05圖4 KEGG富集分析結果

3 討論

KRAS突變型LA有難治、預后差、耐藥性強和容易復發等特點，是臨床治療最為棘手的肺癌類型。盡管目前針對KRAS突變LA有抑制KRAS基因(如法尼基蛋白轉移酶抑制劑、香葉基轉移酶、RAS轉換酶1等)、改變KRAS膜定位(如BKM120、GDC0941和XL147等)等治療方法[19]，但并沒有解決臨床的治療難題。本研究以KRAS突變相關的全轉錄組數據庫為研究對象，通過對KRAS突變型LA的轉錄組學數據進行生物信息學分析發現，KRAS突變產生的影響不僅限于改變細胞周期使腫瘤細胞無限增殖或使EGFR和ALK突變陽性患者出現靶向藥耐藥，實際上對整個基因組均具有不同程度的影響，且這種作用具有指向性。突變影響的基因表達包括對細胞增殖具有顯著影響的細胞周期以及腫瘤預后密切相關的免疫分子等，尤其以IL-6、MMP9、CDH5、TIMP1、TOP2A、TEK、CD36、CLDN5、LCN2、SPP1等基因最為顯著[20-29]，且這些基因與KRAS突變的關系已有研究證實。KRAS基因突變后對LA的預后可能與其表觀調控作用有關，如Caetano等[30]發現血清IL-6水平與LA存活率低及預后較差相關；Xu等[31]發現MMP9在LA的發生發展中發揮重要作用；Zhang等[29]發現SSP1上調PD-1使得LA細胞發生免疫逃逸等。以上均說明，采用轉錄組數據庫分析基因突變的方法學具有實用性，且可與DNA測序等相結合研究基因突變的病理生理學意義，指導基礎和臨床研究。

在基因本體方面，本研究發現，在KRAS突變影響的基因中，大量基因與腫瘤微環境有關，例如血管生成、細胞粘附、免疫反應等相關基因。包括TEK、IL-6、MMP9等在內的基因已經有明確的文獻證實其與腫瘤的轉移、預后不良有關[26,32]。同時，還有大量基因與藥物反應、藥物結合以及藥物代謝有關，這提示，KRAS突變所影響的基因中包含與藥物代謝相關的基因，證實了KRAS突變極易發生耐藥的特點，IL-6、CDH3、TOP2A為代表的KRAS突變影響與藥物代謝相關的分子，有望為未來研究KRAS突變LA的耐藥機制提供思路[33-34]。

在信號通路方面，本研究發現，KRAS突變影響的信號通路不僅限于KRAS突變后RAS通路的異常活化，還包括對原本具有抗腫瘤作用的免疫相關通路的抑制，例如對IL-17信號通路、白細胞內皮遷移等信號通路的抑制作用。這2個信號通路的抑制，可以改變腫瘤微環境，出現免疫耐受而導致預后不良。提示，KRAS突變對抗腫瘤免疫的抑制，可能是KRAS預后不良的重要機制之一。Akbay等[35]發現IL-17可以通過促炎反應促進肺癌的發展，同時這種促炎反應對PD-1的阻斷有抑制作用，導致肺癌細胞免疫耐受，這與我們的預測結果是一致的。另一方面，KRAS突變也可以使細胞周期相關通路異?；罨鰪娔[瘤細胞的增殖能力，同時抑制腫瘤細胞的凋亡信號通路，導致腫瘤細胞過度增殖。此外，KRAS突變還可以影響與人體藥物敏感性有關的細胞色素P450藥物代謝信號通路，發生耐藥。

由此可見，KRAS突變會產生全轉錄組范圍的影響，這種組學影響稱為“組學漣漪”。正是由于“組學漣漪”的出現，才使得KRAS突變成為LA的發生和發展中的“核心基因”。這種“組學漣漪”一方面通過促進細胞周期轉換，抑制凋亡而增強腫瘤細胞的生存能力，另一方面，抑制免疫系統的抗腫瘤作用，使腫瘤微環境有利于腫瘤細胞的生存。正是由于這一系列的組學改變，最終使得KRAS突變型LA出現耐藥和預后不良。

本研究表明，從全轉錄組范圍內進行研究，可以更加深入了解KRAS突變型LA的深層次病理機制，明晰KRAS突變在LA中的表觀調控機制，為KRAS突變型LA的臨床治療和醫學研究提供新的思路。