GEO+TCGA双重挖掘！8图3表，3.9分肿瘤生信SCI全文复现

2023-06-25 09:06

大家好，我是四叶草~~今天带来的这篇文章题目是“Five hub genes contributing to the oncogenesis and trastuzumab-resistance in gastric cancer”是2022年10月发表在Gene期刊的生信文章，我们一起来看一下吧。

题目：Five hub genes contributing to the oncogenesis and trastuzumab-resistance in gastric cancer（五个在胃癌发生和曲妥珠单抗耐药中起关键作用的基因）

期刊：Gene

IF：3.913（2022年）

全文共8图3表（正文5图2表）。

1.期刊介绍：

2.文章概要

套路：肿瘤

基因：5个关键基因（筛选后获得）：GNGT1, KRT7, SOX9, TIMP1, KRT16

数据来源：

GEO数据库（GSE77346数据集）

TCGA数据库（TCGA- STAD RNA seq数据）

技术路线：该研究从 GEO数据库和 TCGA数据库中分别获取了GSE77346数据集以及胃癌（STAD）的RNA seq数据。使用 R 编程语言筛选出不同表达的基因 (DEGs)，分别采用GSEA、GO和 KEGG 富集分析来分析 DEGs 的功能和通路。接着，使用 STRING和 Cytoscape 软件建立蛋白质相互作用网络 (PPI)，并筛选出hub基因。采用ROC 分析方法来评估hub基因的诊断价值。最后通过Kaplan-Meier 分析和对数秩检验分析hub基因的表达和胃癌的总体生存期(OS)对应生存率的相关性。

3.图表简介

Fig. 1| TCGA-STAD和GEO数据集交集基因分析

Fig. 1A| TCGA-STAD数据集差异表达基因（DEGs）火山图

Fig. 1B| GSE77346数据集差异表达基因（DEGs）火山图

Fig. 1C| GSE77346数据集差异表达显著的基因热图

Fig. 1D| TCGA-STAD和GSE77346数据集差异表达基因（DEGs）韦恩图

Table 1| 差异表达基因的GO功能富集分析

Table 2| 差异表达基因的KEGG功能富集分析

Fig. 2| GSE77346数据集GSEA富集分析-经典可视化图

Fig. 3| 差异表达基因的GO-KEGG富集分析-柱状图

Fig. 4| PPI网络构建、子网络构建和hub基因筛选

Fig. 4A| STRING构建蛋白-蛋白互作网络

Fig. 4B| Cytoscape软件MCODE插件子网络分析

Fig. 4C| Cytoscape软件CytoHubba 关键基因（hub genes）筛选

Fig. 5| STAD中5个hub基因表达水平分组比较图和诊断ROC曲线

Fig. 5A-E| 5个hub基因mRNA表达水平的分组比较图

Fig. 5F-J| 5个hub基因的诊断ROC曲线

Fig. 5K| KRT16高低表达组的总生存期分析（生存曲线K-M）

4.分析工具：

仙桃学术（https://www.xiantao.love/）（新版）各种生信分析工具
GEO2R 对GSE23558数据集进行差异分析

5.复现流程

Fig. 1| TCGA-STAD和GEO数据集交集基因分析

Fig. 1A| TCGA-STAD数据集差异表达基因（DEGs）火山图

Fig. 1B| GSE77346数据集差异表达基因（DEGs）火山图

Fig. 1C| GSE77346数据集差异表达显著的基因的热图

Fig. 1D| TCGA-STAD和GSE77346数据集差异表达基因（DEGs）韦恩图

我们开始复现啦~~首先是Figure 1A复现。

我们先登录仙桃学术-生信工具页面。

进入仙桃学术（https://www.xiantao.love/）（新版）→选择“生信工具”

Fig. 1A复现：

第1步：TAGA-STAD数据集筛选分子

点击上方菜单栏【全部工具】→点击左侧导航栏【表达差异】→点击右侧导航栏[云]筛选分子→选择云端数据“TCGA-STAD RNA seq TPM格式表达谱”→选择【临床变量】为【（临床）status】→添加分组（左侧为参考组，右侧为实验组）→默认【分析参数】→点击【确认】，出现任务提交成功提醒→在【历史记录】中查询保存的【筛选分子】结果，更名并下载【筛选分子】数据。

第2步：数据整理

根据得到的筛选分子结果，筛选基因类型protein coding（mRNA），并保留gene name和log2FoldChange、padj几列数据，调整列顺序，得到的数据如下所示：

第3步：TCGA-STAD数据集差异表达基因（DEGs）火山图绘制

点击上方工具选择栏【全部工具】→点击左侧导航栏【表达差异】→点击【火山图】→上传第2步整理好的数据→点击【验证】→【主要参数】中【阈值】默认【logFC阈值】为1，【P值阈值】为0.05→默认其他主要参数→点击【确认】→下载【差异表格.csv】和火山图图片。

PS：原文标注了几个基因，应该是最终得到的hub基因，可以在完成全文的分析之后，再重新绘制火山图并标注hub基因。

可以设置【点】参数的点的大小和透明度，能够提高下面融合成一片的点之间的区分度。

Fig. 1B复现：

第1步：GSE77346数据集差异分析

1、进入仙桃学术的【数据集检索】模块页面

进入仙桃学术-数据集检模块（https://www.xiantao.love/gds）→数据集检索栏中输入数据集GSE77346→点击【检索】，可看到当前数据集的摘要及实验设计信息。

2、添加样本：找到对应数据集和平台，点击右下方的【选择样本】→点击第一行最前面的小方框，第一列都显示，表示选中了所有样本→点击【添加至样本库】。

我们在添加样本时发现，分组中样本数＜2，没有满足仙桃对GEO数据集进行差异分析的条件，所以我们尝试用GEO2R来进行差异分析。

GEO2R进行差异分析

1、在当前页面点击GEO2R

2、点击【Define groups】→输入分组名称并选择相应样本，再点击对应的分组，将样本添加到相应的分组中（先输入实验组后输入参考组）→点击下方的【Analyze】→待分析完成后显示结果图和数据表，点击下方页面的【Download full table】，下载差异分析表格→下载GEO2R分析结果的火山图。

PS：我们看到GEO2R差异分析之后会输出一些可视化图，其中有箱式图和火山图，我们看到箱式图比较整齐，说明数据是经过标准化的，需要再进行标准化设置，可以直接下载差异分析结果。我们可以直接下载GEO2R差异分析得到的火山图，也可以下载差异分析数据后进行整理，上传到仙桃生信工具的火山图模块进行差异分析。

Fig. 1D 复现：

第1步：数据准备：

根据TCGA-STDA RNA seq数据和数据集GSE77346差异分析结果，分别根据阈Padj< 0.05 & |log2(fold change)| >1筛选差异基因，得到两列差异基因列表：

TCGA-STDA RNA seq数据差异分析结果：

数据集GSE77346差异分析结果：

整理好的数据如下所示：

第2步：韦恩图绘制

点击【全部工具】→点击左侧的【基础绘图】→选择右侧的【韦恩图】→选择第1步中整理好的数据→点击【验证】→默认主要参数→下载【交集情况.xlsx】和韦恩图。

Fig.1C复现：

第1步：数据准备，下载表达矩阵

进入仙桃学术工具（https://www.xiantao.love/）→选择【数据集检索】→检索框中输入数据集GSE77346→点击【检索】→点击【数据下载】→跳转界面到GEO数据库→下载GSE77346_series_matrix.txt.gz。

下载的表达矩阵数据需要解压并提取，将文件拖动到打开的EXCEL文件中，并另存为新的EXCEL文件，数据如下所示：

第2步：数据整理，GEO2R差异分析top29显著差异基因筛选，探针替换及添加分组。

1、数据整理

删除1-62行，以及最后一行→整理根据GEO2R差异分析结果，用VLOOKUP函数匹配探针和Gene Symbol（GEO2R差异分析结果可以用来进行匹配），并最终用Gene Symbol替换探针→整理好格式之后另存为CSV格式文件。

2、GEO2R差异分析top29显著差异基因筛选

最终的top29 显著差异分子探针和Gene Symbol:

3、探针替换

4、添加分组

第3步：简易数值热图绘制

点击工具选择栏【全部工具】→点击左侧导航栏【表达差异】→点击右侧导航栏简易数值热图→上传第2步整理好的数据→点击【验证】→选择数据log2(value+1)转换和归一化设置（是否归一化根据数据情况选择）→设置【聚类（顺序）/分割】→设置【色阶】，添加色块描边→点击【确认】→下载热图图片。

Fig. 2| GSE77346数据集GSEA富集分析-经典可视化图

Fig. 2复现：

第1步，整理GSE77346数据集差异分析表格，保留所有分子的id和logFC两列数据，如下所示：

第2步：GSEA富集分析

进入仙桃学术生信工具→点击上方菜单栏【全部工具】→点击左侧导航栏【功能聚类】→点击右侧导航栏【GSEA分析】中的【[GSEA]富集分析】→上传第1步整理好的数据→点击【验证】→默认主要参数→点击【确认】→提示“任务提交成功”，点击【确定】→在【历史记录】中更名并下载GSEA分析结果。

下载的GSEA分析结果如下所示：

参考原文，我们用p <0.01, q value <0.25为阈值，筛选显著富集的通路。

第3步：GSEA经典可视化

进入仙桃学术生信工具→点击【全部工具】→点击左侧导航【功能聚类】→点击右侧导航栏【GSEA分析】中的【[GSEA]经典可视化】→选择第2步保存的结果→添加筛选的显著富集通路ID，我们选择了KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450→默认其他主要参数→点击【确认】→下载GSEA经典可视化图。

Table 1| 差异表达基因的GO功能富集分析

Table 2| 差异表达基因的KEGG功能富集分析

Table 1复现：

第1步：数据准备

根据韦恩图交集情况，提取交集基因，得到差异基因列表。

如下所示：

第2步：GO富集分析

进入仙桃学术生信工具→点击上方工具选择栏【全部工具】→点击左侧导航栏的【功能聚类】→点击右侧导航栏【GOKEGG】中的【[GOKEGG]分析】→上传第1步的数据→点击【验证】→【富集参数】的【条目】选择【全部GO】→默认其他参数→点击【确认】→点击【保存结果】并下载【GOKEGG.xlsx】。

根据原文，用padj＜0.05筛选显著富集的通路。提取相关参数，就可以完成Table 1。

Table 2是KEGG通路富集的分析结果，参考Table 1我们只要在GOKEGG分析中选择【条目】为【KEGG】即可。

Fig. 3| 差异表达基因的GO-KEGG富集分析-柱状图

Fig. 3 复现：

需要重新进行GOKEGG富集分析，【条目】选择【全部GO+KEGG】，过程参考Table 1复现。

接着进行GOKEGG富集分析可视化，绘制柱状图。

进入仙桃学术生信工具→点击上方工具选择栏【全部工具】→点击左侧导航栏【功能聚类】→点击右侧导航栏【GOKEGG】中的【[GOKEGG]柱状图】→选择GOKEGG富集分析保存结果→设置或默认【y轴映射】和【颜色映射】→默认其他参数→点击【确认】→下载图片。

PS：没有显示KEGG结果，是因为在 GOKEGG 富集分析模块中，最终的结果表格只保留了满足较宽的阈值（p<0.1 以及 qvalue<0.2）的结果，而不满足这一较宽值下的条目都会被过滤，如果整个类（BP、CC、MF、KEGG）都不满足这个阈值，那么最终的表格中就会缺少这个类。

Fig. 4|

Fig. 4A| STRING构建蛋白-蛋白互作网络

Fig. 4B| Cytoscape软件MCODE插件子网络分析

Fig. 4C| Cytoscape软件CytoHubba 关键基因（hub genes）筛选

Fig. 4A复现：

准备DEGs列表：

进入STRING数据库（https://string-db.org/）→左侧栏选择【Multiple proteins】→右侧上传DEGs列表→物种选择【homo sapiens】→点击【search】跳转界面→点击【CONTINUE】→跳转界面→点击【Settings】设置【minimum required interaction score】，这里默认了0.40，如果修改阈值，需要点击【UPDATE】更新设置→点击【Exports】，选择SVG格式的PPI网络图片和TSV格式的可以导入Cytoscape软件进一步分析及可视化的文件。

Fig. 4B复现：

文件导入：打开Cytoscape软件→点击【File】→点击【Import】→点击【Network from file】导入STRING数据库下载的TSV格式的互作网络数据。

MCODE插件子网络分析：

在Cytoscape上方【APP】中选择【MCODE】，默认设置Degree Cutoff=2， Max. Depth=100， K-Core=2，Node Score Cutoff=0.2→点击【Analyze Current Network】→可以看到每个子网络的得分，分别点击每个子网络→点击右侧的三个横线符号→选择【Create Cluster Network】。

设置Layout和Style，大家可以根据自己的不要设置不同的布局和样式。

子网络图片导出：点击【File】→【Export】→【Network to image】，保存子网络图片。

我们导出了分数最高的子网络，大家也可以按照相同的方法导出其他的子网络。

Fig. 4C复现：

1、CytoHubba插件鉴定hub基因

点击Cytoscape软件左侧的【cytoHubba】→选择目标PPI网络→点击【Calculate】→选择TOP【10】→默认MCC拓扑算法→点击【Submit】

2、设置Layout

3、保存hub基因的Excel表和hub基因网络图片：点击右侧的【Save Current Rank】即可保存hub基因的Excel表→点击【File】→点击【Export】→点击【Network to image】保存图片。

Fig. 5| STAD中5个hub基因表达水平分组比较图和诊断ROC曲线

Fig. 5A-E| 5个hub基因mRNA表达水平的分组比较图

Fig. 5F-J| 5个hub基因的诊断ROC曲线

Fig. 5K| KRT16高低表达组的总生存期分析（生存曲线K-M）

通过STRING构建PPI网络和Cytoscape的cytoHubba插件，我们得到了10个hub基因，接下来应该对这10个hub基因先进行基因表达水平分组比较图比较它们在STAD和正常组织中的mRNA表达水平的差异，接着通过诊断ROC曲线判断它们对STAD诊断的准确性，最后用K-M曲线进行生存分析。

最后选定的5个基因，是在STAD组织中的mRNA表达水平的显著高于正常组织，并且对STAD诊断的准确性高的基因。

Fig. 5A-E：

基因表达水平分组比较图复现：

点击工具选择栏【全部工具】→点击左侧导航栏【表达差异】→点击右侧【 [云]疾病vs非疾病】→选择云端数据集（TCGA-STAD RNA seq TPM格式数据集）→输入hub基因，以AURKB为例→默认主要参数→点击【确认】→下载基因表达水平的分组比较图。

这里建议大家，在做和文章类似的分析时，先对10个hub基因做基因表达水平的分组比较图，根据诊断ROC和K-M分析之后，确定选定基因的数量，例如原文是最终鉴定了5个hub基因，这可能是因为，不是每个hub基因的mRNA表达水平都在正常组和疾病组之间存在显著差异，也不是每个hub基因用于疾病的诊断准确率高，也不是每个hub基因的表达量对生存结局有显著的影响，所以需要综合这几个分析的结果，鉴定出最终的hub基因。

Fig. 5F-J：

诊断ROC曲线复现：点击工具选择栏【全部工具】→点击左侧导航栏【临床意义】→点击右侧【[云]诊断ROC】→选择云端数据集（TCGA-STAD RNA seq TPM格式数据集）→输入hub基因，以AURKB为例→默认主要参数→点击【确认】→下载诊断ROC图片。

Fig. 5K：

K-M生存曲线复现：点击工具选择栏【全部工具】→点击左侧导航栏【临床意义】→点击右侧【[云]生存曲线（KM图）】→选择云端数据集（TCGA-STAD RNA seq TPM格式数据集）→输入hub基因，以AURKB为例→【预后类型】选择【OS】，【时间】选择【月】→选择分组方法（我们选择了P值最小分组）→点击【确认】→下载生存曲线（KM图）图片。