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的临床治疗和医学研究提供新的思路。