大家好呀！今天小编和大家分享的是2021年5月发表在杂志Cell Death Discovery（IF = 5.241）上的一篇文章，题目为：
A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer
基于细胞焦亡相关基因的卵巢癌预后预测
DOI: 10.1038/s41420-021-00451-x
摘要
卵巢癌（Ovarian cancer, OC）是女性最常见的恶性肿瘤之一，预后较差。细胞焦亡（Pyroptosis）是一种近期发现的，依赖于炎性小体并伴有炎症反应的细胞程序性死亡类型。然而，细胞焦亡相关基因在OC中的表达以及其与预后的关系目前尚不清楚。在本研究中，我们共筛选出31个在OC和正常卵巢组织中差异表达的细胞焦亡调节因子。根据这些差异基因，我们将全部的OC患者分为2个亚型。利用TCGA队列评估每个细胞焦亡相关基因的生存预后。使用LASSO Cox回归分析筛选得到7个基因标志物，将TCGA队列中所有OC患者分为低风险和高风险组。研究结果表明，低风险组OC患者的生存率明显高于高风险组（P<0.001）。根据TCGA队列风险评分的中值将GEO队列中的OC患者分为低风险组和高风险组，其中低风险组总生存期呈上升趋势（P=0.014）。结合临床信息，我们发现风险评分是预测OC患者总生存期的独立因素。焦亡相关基因在肿瘤免疫中发挥重要作用，可用于OCs的预后预测。
研究思路
图1 | 研究技术路线
材料与方法
1. 数据来源 在TCGA数据库中获取397个OC患者的RNA-seq数据，GTEx数据库下载88例正常人卵巢样本的RNA-seq数据，从GEO数据库下载RNA-seq数据和临床信息作为验证集（ID: GSE140082）。
2. 差异基因筛选 我们参考之前的研究（综述），共选取33个与细胞焦亡有关的基因。由于TCGA数据库中缺少正常卵巢组织数据，我们从GTEx数据库下载88个正常卵巢样本数据，并筛选正常组织和肿瘤（OC）组织之间的差异基因，这两个数据集均经归一化处理。差异基因筛选使用R包“limma”，P<0.05作为显著性筛选的标准。使用STRING数据库构建蛋白互作（PPI）网络。
3. 细胞焦亡相关基因预后模型的构建和验证 我们使用Cox回归分析来评估TCGA队列中每个基因与生存状态之间的相关性。为防止遗漏，我们将P值设为0.2，共筛选出7个生存相关的基因进行下一步分析。使用LASSO Cox回归模型（R包“glmnet”）缩小候选基因的范围，并构建预后模型。最终保留了7个基因及回归系数。TCGA表达数据标准化（R包“scale”函数）后计算风险得分，计算公式如下：Risk Score = （X: 回归系数；Y: 基因表达水平）。根据风险得分中值，我们将OC患者分为高风险组和低风险组并使用Kaplan-Meier分析比较两组间的总生存期差异。PCA分析（基于7个基因标志物）使用R包“stats”中的“prcomp”函数。ROC曲线分析使用R包“survival”、“survminer”、“time-ROC”。
4. 风险得分的独立预后分析 我们提取TCGA队列中患者的临床信息（年龄和分级），以及GEO队列患者的年龄和分期数据。这些变量与我们的回归模型中的风险评分相结合进行分析。使用单因素和多因素Cox回归分析建模。
5. 高、低风险组间差异基因富集分析 基于风险得分中值将OC患者分为高风险组和低风险组。|log2FC|≥1和FDR<0.05作为两组（高、低风险组）间差异基因的筛选标准。基于筛选出的差异基因，使用R包“clusterProfiler”进行GO和KEGG富集分析。使用ssGSEA分析（R包“gva”）计算免疫浸润细胞得分并评估免疫相关通路的活性。
6. 统计学分析 使用单因素方差分析比较OC患者正常卵巢组织和肿瘤组织间的基因表达水平，类别变量的比较使用皮尔森卡方检验。使用Kaplan-Meier法进行高、低风险组总生存期比较。使用单因素和多因素Cox回归分析进行风险评估。使用Mann-Whitney检验比较高、低风险组免疫细胞浸润和免疫通路激活情况。统计学分析均使用R来完成，研究技术路线图见图1。
研究结果
1. 正常组织和肿瘤组织之间差异基因的筛选
我们比较了OS患者379个肿瘤组织和88个正常组织样本中细胞焦亡相关基因（N=33）的表达水平，共鉴定出31个DEGs（P<0.01）。其中，13个基因（PRKACA、GSDMB、SCAF11、PJVK、CASP9、NOD1、PLCG1、NLRP1、GSDME、ELANE、TIRAP、CASP4、GSDMD）下调表达，另外18个基因（GPX4、NLRP7、NLRP2、CASP3、CASP6、TNF、IL1B、IL18、CASP8、NLRP6、GSDMA、GSDMC、PYCARD、CASP5、AIM2、NOD2、NLRC4、NLRP3）上调表达（图2）。为进一步探究33个细胞焦亡相关基因之间的相互作用，使用STRING数据库进行蛋白互作网络（PPI）分析（图3）。STRING数据库对于每一个蛋白质相互作用信息都有一个打分值，本研究选取的是打分值高于0.9的高置信度数据，以确保数据的可靠性。根据功能节点的网络属性，确定Hub基因为：CASP1、PYCARD、NLRC4、NLRP1、CASP5、NLRP3、CASP8、AIM2。其中，除CASP1外，其它基因均是肿瘤组织和正常组织之间的DEGs。细胞焦亡相关基因表达量相关性网络见图4。
图2 | 33个细胞焦亡相关基因在不同组织中（肿瘤和正常）的表达水平聚类热图。N代表正常组织、T代表肿瘤组织。红色代表高水平表达、绿色代表低水平表达。* P< 0.05、** P< 0.01、*** P< 0.001。
图3 | 细胞焦亡相关基因之间的PPI网络图。
图4 | 细胞焦亡相关基因表达相关性网络图。红色的线代表呈正相关、蓝色的线代表呈负相关，颜色越深代表相关性越强。
2. 基于DEGs的肿瘤亚型分类
为进一步探究31个细胞焦亡相关DEGs的表达水平与OC亚型之间的关系，我们对TCGA队列中全部379名OC患者样本数据进行一致性聚类分析。通过逐渐增加聚类数目k值（k=2~10），我们发现，当k=2时组内相关性最高且组间相关性较低，表明基于31个DEGs可以将379名患者人群分成2个簇（图5）。基因表达谱和临床特征（包括：肿瘤分化程度G1-G3、年龄≤60岁或>60岁、生存状态存活或死亡）以热图的形式呈现，但是我们发现这两组人群的临床特征几乎没有差别（图6）。Kaplan-Meier分析结果表明，两类患者人群总生存期（OS）不存在明显差异（图7）。
图5 | 一致性聚类分析结果将379名OC患者划分成2个簇（k=2）。
图6 | 2类（或簇）患者人群的基因表达谱热图和临床特征。G1（高分化）、G2（中分化）、G3（低分化）代表肿瘤分化程度。
图7 | 2类（或簇）患者人群Kaplan–Meier OS曲线图。
3. 基于TCGA队列的预后模型构建
使用具有完整生存信息的OS患者样本数据，共374个。使用单因素Cox回归分析初步筛选与生存相关的基因。其中，7个基因（AIM2、PLCG1、ELANE、PJVK、CASP3、CASP6、GSDMA）符合筛选标准P<0.2，用于下一步分析。其中，3个基因（PLCG1、ELANE、GSDMA）危险比HRs>1，其余4个基因（AIM2、PJVK、CASP3、CASP6）HRs>1。使用LASSO Cox回归构建预后模型，风险得分计算公式为：risk score = (−0.187*AIM2 exp.) + (0.068*PLCG1 exp.) + (0.097*ELANE exp.) + (−0.143*PJVK exp.) +(−0.086*CASP3 exp.) + (−0.033*CASP6 exp.) + (0.130*GSDMA exp.)。基于风险得分中值，我们将374个OS患者分成高、低风险组。PCA结果将不同风险得分的患者很好地分为两类。与低风险组相比，高风险组患者的死亡率高且生存时间较短。
4. 风险得分的外部验证
从GEO数据库中下载GSE140082（含380例OC患者样本数据）数据集作为验证集。首先使用“Scale”函数对基因表达数据进行标准化处理。基于TCGA队列的风险得分中值，GEO队列中203名患者被划分为低风险组，177名患者划分为高风险组。PCA结果将不同风险得分的患者很好地分为两类。研究发现，低风险组患者的生存率和死亡率均高于高风险组。此外，Kaplan-Meier分析也显示低风险组与高风险组患者的生存率存在显著差异。基于GEO队列的ROC曲线分析结果表明，我们的模型具很好的预测效果。
5. 风险模型的独立预后价值
我们使用单因素和多因素Cox回归分析来评估每个基因标签的风险得分模型是否可以作为一个独立的预后因素。单因素Cox回归分析表明，在TCGA和GEO队列中，风险得分是预测不良生存的独立因素（HR=3.285, 95% CI: 1.973-5.467; HR: 2.613, 95% CI: 1.319-5.175）（图8A, C）。多因素Cox回归分析表明，在调整其他混杂因素后，风险得分是O患者的预后因素（图8B, D）。此外，我们基于TCGA队列的临床病理特性生成一个热图，并发现低、高风险组患者人群的年龄和生存状况存在差异（P<0.05）（图9）。
图8 | 风险得分的单因素和多因素Cox回归分析。
图9 | 临床病理特征和高、低风险组关系热图。绿色代表低表达、红色代表高表达。* P<0.05。
6. 基于风险模型的功能分析
为了进一步探究基于风险模型分类的高、低风险组之间基因功能和通路的差异，我们使用R包“limma”筛选DEGs，筛选标准设置为：FDR<0.05和|log2FC|≥1。在TCGA队列的高、低风险组间共筛选出115个DEGs。其中，高风险组66个基因上调表达，其余49个基因下调表达。接下来，我们基于115个DEGs进行GO和KEGG富集分析。研究结果表明，结果表明，DEGs主要与免疫应答、趋化因子介导的信号通路和炎性细胞趋化有关（图10A, B）。
图10 | TCGA队列中高、低风险组DEGs功能富集分析。（A）GO富集气泡图（气泡大小代表Gene数；气泡的颜色代表P值高低，越红代表差异越显著）；（B）KEGG通路图（条形越长代表富集到的基因越多）。
7. 高、低风险组间免疫活性比较
我们使用单样本基因集富集分析（ssGSEA）进一步比较分析TCGA和GEO队列中低、高风险人群之间16种免疫细胞的富集分数和13种免疫相关通路的活性。在TCGA队列中（图11A），高风险人群免疫细胞的浸润程度普遍较低，尤其是CD8 +T、NK、Th（含Tfh、Th1、Th2）、TIL、Treg细胞。在TCGA队列中，除了type-2 IFN通路外，其余12种免疫通路在高风险组中的活性低于低风险组（图11B）。GEO队列中的免疫状态与ATGC类似。此外，我们发现，与高风险组相比，低风险组的DCs、iDCs和巨噬细胞显著富集，而type-2 IFN通路活性下调。
图11 | 免疫细胞和免疫相关通路ssGSEA得分比较。（A, B）TCGA队列中低（绿色箱子）、高（红色箱子）风险人群之间16种免疫细胞的富集分数和13种免疫相关通路的活性比较；GEO队列中低（蓝色箱子）、高（红色箱子）风险人群之间16种免疫细胞的富集分数和13种免疫相关通路的活性比较；ns代表无显著差异，* P<0.05，** P<0.01，*** P<0.001。
好啦~今天的内容到这里就结束啦！感谢大家的支持，我们继续精进。
生信分析咨询
欢迎扫描下方二维码联系小编
联系电话：13121367263