Identification of Key Genes Related to Lung Squamous Cell Carcinoma Using Bioinformatics Analysis
文章思路
-典型建模类文章
-在本研究中，一方面
-数据：分析了GEO数据库的三个RNA芯片数据集和来自TCGA数据库的RNAseq数据；
-筛选：正常组织与肺鳞状细胞癌中的差异表达基因（DEG）
-分析：使用DEG进行PPI网络分析
-得到：获得关键基因。
另一方面，
-数据：使用TCGA数据库的DEG
-分析：进行WGCNA分析
-得到：获得与LUSC发生和发展有关的基因。
最终，将PPI网络分析与WGCNA分析得到的基因取交集得到关键基因作为对该疾病至关重要的中心基因。 此外，我们利用TCGA中DEG的患者临床信息进行单变量Cox回归分析，Lasso回归分析和多变量Cox回归分析，以鉴定与LUSC患者预后相关的关键基因。
下面我们来看看具体的解读吧！
01
02
肺癌是全球致命的恶性肿瘤之一，肺癌主要分为小细胞肺癌（SCLC）和非小细胞肺癌（NSCLC），其中NSCLC，占肺癌病例的约85％，其中NSCLC又可分为肺腺癌（LUAD）和肺鳞状细胞癌（LUSC）。
当前，许多治疗方法被用于治疗肺癌，早期NSCLC患者经常手术切除治疗为主，对于晚期病人，其肿瘤不能被手术切除，但靶向疗法或免疫疗法联合化疗是最好的方法之一。
然而，免疫疗法的治疗费用非常昂贵，最近有报道称这种疗法可导致严重的副作用，此外，LUSC占NSCLC病例的30％，在中老年男性中更常见，转移和复发率较高。非小细胞肺癌患者在晚期大多才被诊断，其5年生存率比早期患者低，因此，需要进一步研究预后标志物以个性化癌症治疗。
基因组芯片和高通量测序技术与生物信息学分析相结合成为发现疾病生物标记物强大工具。
03
1.数据下载
-GEO数据库：GSE2088，GSE6044，GSE19188，样本信息如表1所示
表1 GEO数据库三套数据集样本信息统计
-TCGA数据库：LUSC RNAseq数据以及相应的临床信息，包括502个肿瘤样本和49个正常样本。
2.数据预处理
使用R程序包sva去除不同数据集之间的批次效应，对TCGA数据库FPKM数据进行标准化，删除在所有样本平均表达值0.5的基因。
3.数据分析
筛选差异基因：R程序包limma ，阈值| log2FC | > 1.0和FDR<0.05
基因互作网络：STRING数据库，score>0.4, Cytoscape对结果进行可视化，R程序包WGCNA加权网络共表达分析。
生存分析：单因素Cox比例回归分析，Losso回归分析，多因素Cox比例回归分析，Kaplan–Meier曲线分析，诺莫图分析
分析过程中涉及R程序包：survival, caret, glmnet, survminer, survivalROC
04
1.筛选差异基因（DGEs）
———-
首先在三套GEO数据集中采用层次聚类去除离群值样本，共得到97个正常肺组织样本和84个肺鳞癌组织样本，然后对GSE2088、GSE6044和GSE19188三个数据集进行合并并去除批次效应。进一步在合并后的GEO数据集中进行差异分析，共得到485个上调基因119个下调基因（图1B），在从包含49个正常肺组织样本和499肺鳞癌样本的TCGA数据集中，共得到3348个上调基因和3387个下调基因。最后，分别对GEO和TCGA数据集中筛选的上下调基因取交集（图1C），得到337个显著上调和119个下调的基因进行下一步分析。
图1 正常组织与肿瘤组织中差异分析结果
2.差异基因PPI网络分析
———-
PPI网络是基于STRING数据库由Cytoscape构成的网络拓补图，由476个节点和4347个边缘组成，包括362个上调和114个下调基因。使用Cytoscape插件CytoHubba通过5种计算指标在PPI网络中筛选关键基因，在PPI网络中挑选每种指标得分前20个基因，最后提取出在5种计算得分均出现的基因作为hub基因，这些基因是：TOP2A、CCNA2、CDC20、AURKA、AURKB和FEN1，它们可能在LUSC进展中起重要作用（图2A）。然后使用Cytoscape插件MCODE进行模块分析，结果显示，上面5种指标计算的前20个基因大部分出现在模块1中，作者将模块1作为重要的模块。模块1包括54个节点和1380个边（图2B），这个模块中的基因都是上调基因。使用DAVID对模块中1的DEGs进行了富集分析。GO通路富集分析表明，模块1的基因参与了细胞分裂、有丝分裂核分裂、ATP和蛋白质的结合（图2C）等功能。
图2 PPI网络分析与差异基因富集分析
3.差异基因DEGs WGCNA 分析
———-
从TCGA RNA-seq数据集中提取差异基因表达谱，使用层次聚类法对样本进行聚类，排除两个离群样本，选择软阈值β=5建立基因调控网络。结果显示蓝色、黄色和绿松石色模块与样本的表型形状相关系数最大，分别为0.538、-0.542和-0.870（图3A）。此外，作者使用第一主成分分析和皮尔逊相关系数计算各模块之间的相关性，从这些结果表明，绿松石色、黄色模块中基因具有一致性，且两个模块与表型的相关系数均为负值。作者认为这两个模块所有基因的MM和GS值上四分位的是该模块的关键基因，蓝色模块的基因分布如图3B所示。值得注意的是，前面通过PPI网络筛选的关键基因大部分在蓝色模块中。最后，作者对这些基因进行富集分析，GO和KEGG分析结果表明，蓝色模块中基因参与与细胞周期、有丝分裂、核分裂、p53信号通路等途径，可能与癌细胞过度增殖有关，而且这些过程还参与了许多经典抗癌药物的代谢途径，因此，这个模块中的基因对药物开发具有指导意义。
图3 WGCNA分析与模块基因富集分析
4.LUSC中hub基因的免疫组化分析
———-
从PPI网络筛选出的关键基因和WGCNA分析中蓝色模块基因取成交集，主要包括CCNA2、AURKA、AURKB和FEN1，且GO分析表明蓝色模块与细胞周期相关，也与PPI子模块分析一致。这表明这四个基因可能在LUSC的发生发展中起着关键作用。所以作者将这四个基因CCNA2、AURKA、AURKB和FEN1定义为与LUSC相关的关键基因。为了进一步验证这些基因在LUSC的表达情况，使用免疫组化分析这些基因在正常组织和肿瘤组织中的表达情况，结果表示这些基因均在LUSC组织中高表达（图4）。
图4 hub基因的免疫组化分析
5.生存分析
———-
接下来，作者将TCGA数据集随机分为训练集数据和验证集数据，首先使用训练集携带的生存数据对前面筛选497个在TCGA样本和GEO样本均存在差异的基因进行单因素回归分析，共得到91个对预后有显著影响的基因（p<0.05），然后使用Losso回归分析并进行1000次的10倍较差验证，对基因的个数进行降维分析，最终选择λ=0.093，整个模型的误差最少，从91个基因中得到8个基因进行后续的分析，进一步，使用多变量Cox比例风险回归分析，共获得5个风险基因即MYEOV、LCE3E、PTGIS、OR2W3、RALGAPA2（图5A），并构建风险预后模型，即Risk score = (0.0137 × MYEOV) + (0.0152 × LCE3E) + (0.0223 × PTGIS) +(0.180 × OR2W3) + (0.0319 × RALGAPA2)，再根据风险得分公式对训练集和验证集中所有样本进行风险得分，并根据风险得分的中位数将患者分为高低风险两组，作者发现训练集和验证集中，低风险组患者的生存状况要优于高风险组，两组之间存在显著差异（图5B-C）。最后，作者使用在训练集和验证集使用ROC曲线对模型的预测对患者1年，3年和5年性能进行评估，结果显示训练集中风险预测模型得出的AUC分别为0.811、0.924和0.937），验证集数据中分别风险预测模型得出的AUC分别为0.692、0.722和0.651（图5D-F），并使用其他机器学习算法，包括决策树（DT）、Naïve贝叶斯算法（NB）和随机森林（RF）与多变量Cox模型的性能进行比较。可以得出结论多变量Cox模型在所有构建模型的方法中最有效。
图5 构建基因风险模型
此外，作者根据生存时间和基因风险模型得分绘制散点图。从图中可以看出随着风险得分的增加，在训练（图6A）和测试数据（图6B）中死亡的患者数量增加，患者存活时间逐渐减少。最后，在训练集和验证集中，与低风险组患者相比，基因MYOEV、PTGIS、OR2W3和RALGAPA2基因在高风险组患者中显著高表达，而基因LCE3E在两组患者中表达量并不很明显（图6C）。
图6.风险模型得分与基因表达量、临床信息的关系
6.模型验证
———-
美国癌症联合委员会（AJCC）定义了stage用来评价癌症患者的预后效果。另外，Carter等人用特异基因（包含25个基因）中鉴定出染色体不稳定性（CIN25），用这种方法可以评价各种癌症患者的预后情况。作者从TCGA数据库中收集LUSC患者的AJCC分期，然后计算患者的CIN25，并将其与AJCC分期和CIN25进行比较，计算了风险评分、AJCC分期和CIN25的C指数（表1）。结果显示，验证集风险评分C指数为0.642，显著高于AJCC分期（0.576，p<0.05）和CIN25（0.555，p<0.05）；训练集风险评分C指数为0.668，显著高于AJCC分期（0.527，p<0.05）和CIN25（0.545，p<0.05）。这说明，与其他预后方法相比，基因风险得分是患者预后状态的良好预测指标。
表1 基因风险评分、 AJCC分期和CIN25的C-index比较
最后，作者使用TIMER数据库分析风险模型中基因表达量与免疫细胞免疫浸润之间的相关性，结果显示OR2W3的表达情况与B细胞（Cor=0.183，p=6.36e-05）、CD4+T细胞（Cor=0.219，p=1.42e-06）、巨噬细胞（Cor=0.223，p=8.46e-07）和树突状细胞（Cor=0.208，p=4.92e-06；图7A）的免疫评分呈正相关。RALGAPA2的表达情况与B细胞（Cor=0.111，p=1.63e-02）、CD4+T细胞（Cor=0.418，p=1.52e-21）、巨噬细胞（Cor=0.315，p=2.00e-12）、中性粒细胞（Cor=0.311，p=3.71e-12）和树突状细胞（Cor=0.240，p=1.21e-07；图7B）的免疫评分呈正相关。PTGIS的表达情况与B细胞（Cor=0.314，p=2.93e-12）、CD8+T细胞（Cor=0.260，p=9.12e-09）、CD4+T细胞（Cor=0.414，p=3.81e-21）、巨噬细胞（Cor=0.464，p=7.53e-27）、中性粒细胞（Cor=0.329，p=1.85e-13）和树突状细胞（Cor=0.454，p=1.98e-25；图7C）的免疫评分呈正相关。LCE3E的表达情况与B细胞（Cor=-0.208，p=5.15e-06）、CD8+T细胞（Cor=-0.224，p=8.07e-07）、CD4+T细胞（Cor=-0.114，p=1.30e-02）、巨噬细胞（Cor=-0.11，p=1.67e-02）、中性粒细胞（Cor=-0.124，p=6.54e-03）和树突状细胞（Cor=-0.137，p=2.88e-03；图7D）的免疫评分呈负相关。MYEOV表达情况与CD4+T细胞（Cor=0.147，p=1.27e-03）、中性粒细胞（Cor=0.168，p=2.35e-04）和树突状细胞（Cor=0.195，p=1.96e-05；图7E）的免疫评分呈正相关。以上结果表明，风险评分公式中的5个基因与LUSC患者的免疫浸润过程密切相关，这可能是LUSC患者作为有效预后指标的原因之一。
图7 风险模型中基因表达量与免疫细胞免疫浸润的相关性
