1. 引言
肾癌是一种起源于肾脏、肾实质和尿路上皮的恶性肿瘤 [1] ，根据美国监测、流行病学和终末结果(Surveillance, Epidemiology, and End Results, SEER)数据库(https://seer.cancer.gov/statfacts/html/kidrp.html)，肾癌的发病率逐年上升，占所有癌症死亡人数的2.4%，5年相对生存率仅为77.6%。其中，肾透明细胞癌(KIRC)是最常见、最致命的肾癌亚型，占肾癌的75% [2] 。KIRC患者的主要死因是肿瘤转移，肾癌早期无明显临床特征，25%~30%的患者在初诊时已有转移 [3] 。因此，疾病的早期诊断和治疗非常重要。近年来，随着肿瘤分子生物学的快速发展，基因组研究表明，mRNA的重要性与日俱增，筛选与疾病相关的mRNA作为疾病生物标志物用于疾病治疗的研究也有报道 [4] 。因此，有必要开展肾癌生物标志物的分子筛选，以改善肾癌的预后并降低死亡率。
细胞死亡是所有生物的基本生理过程。它的作用包括胚胎发育、器官维护和衰老，以及协调免疫反应和自身免疫 [5] 。已知程序性细胞死亡主要有三种方式：细胞焦亡、细胞凋亡和坏死。尽管细胞焦亡、细胞凋亡和坏死先前已被确定和描述为不同的独立程序性细胞死亡途径，但越来越多的证据表明这三种途径之间存在广泛的相互作用 [6] 。2009年，Han等人发现在同一细胞株内，细胞程序性坏死与细胞凋亡可以相互转变，线粒体膜的通透性变化可能间接参与这一过程 [7] 。另外在同一环境中，细胞焦亡与细胞凋亡能够共同作用 [5] 。泛凋亡是一个新发现的概念，强调了细胞焦亡、细胞凋亡和坏死之间的作用和协调 [8] [9] 。而细胞死亡与肿瘤发生有关，Hanahan等人认为抵抗细胞死亡是肿瘤的基本特征之一 [10] 。
孟德尔随机化(Mendelian randomization, MR)作为一种工具变量分析方法，近年来因其自身的特殊优势和基因组学的快速发展而受到医学研究领域的青睐。在传统的观察性研究中，潜在的混杂因素和反向因果关系会影响其因果推断能力 [11] [12] 。MR是一种工具变量(Instrumental variable, IV)分析，利用遗传变异作为IV来检测和量化因果关系 [13] 。由于MR能够克服潜在混杂因素和反向因果关系的影响，近年来越来越多地用于观察性研究 [14] [15] 。早期的MR研究通常是在小样本人群中进行的，而且只使用了少量的遗传变异 [16] ，这使得MR研究的效果不佳。然而，随着生物界发现大量与特定性状密切相关的遗传变异，以及许多大样本全基因组关联研究(Genome-Wide Association Studies, GWAS)公开发布大量暴露与遗传变异关系的集合数据 [16] ，该领域发生了一场革命。这些集合数据使研究人员能够估计大样本数据中的遗传关联，从而促进了MR研究的发展。
本研究旨在通过生物信息学方法鉴定在KIRC中发挥重要作用的泛凋亡相关基因(PANoptosis-related genes, PRG)，通过孟德尔分析揭示其与KIRC的因果关系，探讨其在KIRC中的潜在调控机制，为临床开发新的治疗手段和新药提供理论支持。
2. 资料与方法
收集数据集和 PANoptosis相关基因
从基因表达总库(GEO)中下载KIRC数据集GSE53757 (共144个样本，包括72个KIRC肿瘤组织样本和72个对照组织样本)作为训练集 [17] ，KIRC数据集GSE40435 (共202个样本，包括101个KIRC组织样本和101个对照组织样本)作为验证集。下载癌症基因组图谱(TCGA)-KIRC 数据集作为验证集，其中包含535个KIRC组织样本和72个对照组织样本。从先前的研究 [18] 中提取了包括细胞热凋亡、细胞凋亡和坏死在内的347个细胞凋亡相关基因(PANoptosis related genes, PRG)。
PANoptosis相关差异表达基因的鉴定
使用R软件包“limma” [19] 对72份KIRC肿瘤样本和72份对照组织样本的转录组数据进行差异表达基因(Differentially expressed gene, DEG)分析，以鉴定合并数据集中患病样本和对照样本之间的DEGs。DEGs的阈值设定为log2 fold change (|log2FC|)的绝对值 > 1且adj.P < 0.05，表示DEGs表达上调。同样，log2FC < −1且adj.火山图和热图显示DEG筛选结果。通过交叉鉴定PANoptosis相关基因和DEGs基因，确定了差异表达的PANoptosis相关基因(Differentially expressed PANoptosis-related genes, DE-PRGs)。
功能富集分析和蛋白质-蛋白质相互作用网络
京都基因和基因组百科全书(KEGG)是从分子水平应用的资源数据库，尤其是大型分子数据集中生成基因组测序和其他高通量实验技术。基因本体(GO)系统由三部分组成：生物过程、分子功能和细胞成分。利用R软件包“clusterProfiler” [20] ，在P < 0.05和计数 > 2的统计条件下，筛选出与DE-PRGs相关的KEGG通路和GO功能。
为了了解DE-PRGs的蛋白–蛋白相互作用(protein-protein interactions, PPI)，我们将这些基因上传到STRING 数据库进行PPI网络预测(相互作用得分大于0.7)。
孟德尔随机化分析
为获得候选基因的单核苷酸多态性(single nucleotide polymorphism, SNP)数据，首先将上述DE-PRGs转化为Ensemble ID，以获得这些基因在GWAS数据库中的暴露ID，并在此过程中剔除不存在表达量性状位点(expression quantitative trait loci, eQTL)的基因 [21] 。
然后，对存在eQTL的DE-PRGs分别进行MR分析，并筛选出标记基因。我们以固定效应方差反加权(inverse variance weighting, IVW)为主要方法，以MR Egger、加权中位数和简单中位数为次要方法，进行了四种不同的MR分析。虽然每个基因变异都满足IVs假设，但IVW方法使用元分析方法来可靠地评估暴露对结果的因果关系 [22] 。固定效应IVW结果是本研究的主要依据。
最后，为了在多变量水平上探讨暴露因素(标记基因)的因果关系，我们对标记基因进行了多变量MR分析，以观察其与KIRC的关系是否有变化。
基于DE-PRGs的诊断诺模图的构建与验证
为了评估最佳基因生物标志物区分KIRC和正常样本的能力，我们构建了受试者工作特征曲线(receiver operating characteristic, ROC)来评估整体诊断能力，并测量了曲线下面积(area under the curve, AUC)和准确率。此外，我们还使用“rms”软件包 [23] 构建了预测KIRC发病率的诺模图模型，以探讨关键DE-PRGs在KIRC诊断中的具体作用。该模型采用校准曲线和决策曲线分析法(decision curve analysis, DCA)进行拟合 [24] 。
单基因基因组富集分析
我们使用“clusterProfiler”软件包 [20] [25] 对每个特征基因进行了基因组富集分析(gene set enrichment analysis, GSEA)，进一步分析了KIRC中这些关键DE-PRGs的相关通路和潜在生物学功能。以P < 0.05作为显著富集的标准。
免疫浸润分析和靶向药物预测
CIBERSORT算法是一种通过基因表达谱分析组织细胞组成特征的方法 [26] ，该算法可估算训练集中22种免疫细胞类型的比例，从而寻找细胞生物标志物和治疗靶点。随后，通过观察正常样本和KIRC样本中22种免疫细胞浸润丰度的差异，以及所有标记基因与22种免疫细胞类型的相关性分析，探讨了标记基因与免疫微环境之间的关系。
此外，为了分析化疗药物对相关基因的影响，我们通过DGIDB (https://dgidb.org/)数据库构建了药物与基因相互作用网络，该数据库存储了大量关于药物与基因相互作用关系的文献和报告。
构建标记基因的ceRNA网络
mRNA的表达受lncRNA与miRNA竞争性结合的调控。本研究利用R软件包multiMiR分析了miRNA与靶标的相互作用及其与疾病和药物的关联 [27] ，该软件包全面收集了已预测和验证的miRNA与靶标的相互作用及其与疾病和药物的关联。仅保留了两个数据库同时预测得到的miRNA和lncRNA。并根据上述关系构建了ceRNA (mRNA-miRNA-lncRNA)调控网络。
验证标记基因的表达
为了进一步验证标记基因在正常组织和病变组织中的表达，我们在GEO和TCGA数据库中对标记基因的表达进行了验证。
实时荧光定量PCR (RT-qPCR)
采集KIRC患者的肿瘤组织和配对的邻近组织。使用Trizol试剂提取肿瘤组织和邻近组织的RNA。使用RT kit gDNA Eraser (Perfect Real Time)合成cDNA，并进行RT-qPCR。用SYBR Pre-mix Ex Taq II (TliRNaseH Plus)检测mRNA表达。该实验已获得济南市中心医院人类伦理委员会和山东大学伦理办公室的批准。
统计学分析
本研究的所有数据分析均使用R软件(4.2.2版)进行，使用wilcox.test进行两组间的比较。斯皮尔曼检验用于揭示评分基因与免疫细胞之间的关系，P < 0.05为显著。所有MR分析均由“双样本MR”和“MRPRESSO”软件包完成。双侧P值小于0.05视为具有统计学意义。
3. 结果
DE-PRGs的鉴定
为了筛选出与KIRC相关的基因，我们使用Limma软件包对72个肿瘤组织样本和72个对照组织样本的转录组数据进行了差异分析，差异筛选条件为adj.P值小于0.05且|Log2FC| > 1。经统计，共有3446个基因在疾病中与对照组相比有显著差异表达，其中上调基因1557个，下调基因1889个(图1(a)、图1(b))。我们将3446个DEGs和347个PANoptosis基因进行了交叉。如Venn图所示(图1(c))，共得到79个DE-PRG。
功能富集分析和蛋白质-蛋白质相互作用网络
对79个DE-PRGs进行了KEGG和GO富集分析，并在P < 0.05、计数 > 1的条件下筛选出富集的通路。最终结果显示，共富集了823个GO条目，包括生物过程、细胞组分和分子功能。具体条目包括细胞因子产生的正向调控、蛋白丝氨酸/苏氨酸激酶活性的调控、内部多肽酶活性的调控、半胱氨酸型内肽酶活性的调控、内部凋亡信号通路、伤口愈合、激酶活性的正向调控、肽酶活性的正向调控、半胱氨酸型内肽酶参与凋亡活性的调控等。KEGG富集了67条通路，具体包括：细胞凋亡、沙门氏菌感染、癌症中的微粒体、PI3K-Akt信号通路、坏死、癌症中的转录失调、癌症中的蛋白多糖、JAK-STAT信号通路、脂质和动脉粥样硬化通路等。
将79个DE-PRGs输入STRING数据库，PPI网络的阈值为0.7，该网络由62个节点和156条边组成。
孟德尔随机化分析
对于鉴定过程中得到的79个DE-PRGs，进行Ensemble ID转换，得到暴露因子ID。在79个DE-PRGs中，66个基因存在eQTL，其中20个基因在数据库中相应的暴露因子ID如下表1所示。
单变量MR 分析
对具有eQTL的66个基因逐一进行了MR分析，结果如下。反方差加权法的p值小于0.05，表明这些基因与KIRC疾病有因果关系；以1为分界线计算OR值，大于1为危险因子，小于1为安全因子。从森林图(图2(a))和散点图(图2(b)~(d))来看，固定效应IVW结果显示，PRF1 (P = 0.030，OR 95%置信区间[CI] = 1.599 [1.047~2.441])、TNFRSF10B (P = 0.008, OR 95% CI = 1.238 [1.058~1.447])和PDGFRB (P = 0.012, OR 95% CI = 0.546 [0.340~0.876])与KIRC存在遗传因果关系。通过筛选暴露相关的SNP并去除连锁不平衡(LD)，我们得到了4个与PDGFRB相关的SNP，3个与TNFRSF10B相关的SNP，5个与PRF1相关的SNP。结果发现，PDGFRB (P = 0.012, OR < 1)是KIRC疾病的保护因素，PRF1 (P = 0.03, OR > 1)和TNFRSF10B (P = 0.008, OR > 1)是KIRC疾病的危险因素。
此外，漏斗图(图3(a)~(c))显示，PDGFRB、PRF1和TNFRSF10B与KIRC之间的MR分析呈正态分布(P > 0.05)。Cochran’s Q统计量(MR-IVW)和Rucker’s Q统计量(MR Egger)显示，标记基因与KIRC的MR分析无异质性(P > 0.05) (表2)。MR Egger分析的截距检验表明，PDGFRB、PRF1和TNFRSF10B和KIRC的MR分析没有水平多向性(P > 0.05) (表2)。“留一法”分析表明，我们的MR分析不是由单一SNP驱动的(图3(d)~(f))。
多变量MR分析
筛选IV后，通过分析与任何暴露因子PDGFRB、TNFRSF10B和PRF1相关的IV，得出以下结果。关注每个暴露因子的P值并计算OR值。森林图(图2(e))显示，PRF1 (P = 0.922，OR 95%置信区间[CI] = 1.054 [3.686e(−01)~3.012])、TNFRSF10B (P = 0.096, OR 95% CI = 0.003 [2.888e(−06)~2.835])和PDGFRB (P = 0.016, OR 95% CI = 0.230 [6.936e(−02)~0.762])，这说明TNFRSF10B的安全/危险属性在3个暴露因子中发生了变化，由危险因子变为安全，而PRF1和PDGFRB没有变化。
Figure 1. (a) Volcano plot shows expression pattern of DEGs in the cohort of GSE53757. Green is down-regulated genes. Red is up-regulated genes. (b) Heatmap plot compares the expression density of DEGs in tumor group and normal group. (c) Venn diagram showing the intersecting 79 DE-PRGs
图1. (a) 火山图显示GSE53757队列DEGs的表达模式。绿色，下调基因。红色，上调基因。(b) 热图比较肿瘤组和正常组中DEGs的表达密度。(c) 维恩图显示79个DE-PRGs
Table 1. 66 DE-PRGs existed eQTL (20 genes demonstrated).
表1. 66个DE-PRGs存在eQTL (表格展示20个基因)
Table 2. Sensitivity analysis of 3 marker genes (PDGFRB, TNFRSF10B and PRF1)
表2. 3个标记基因(PDGFRB、TNFRSF10B和PRF1)的敏感性分析
GSEA
在本研究中，我们以三个标记基因为目标基因进行了GSEA富集分析。这3个基因大多在适应性免疫反应、线粒体内膜、线粒体基质、线粒体含蛋白复合物、呼吸链复合物、碳代谢、柠檬酸循环(TCA循环)、Th1和Th2细胞分化、TNF信号通路和其他通路条目中被显著富集。前10条明显富集的通路被可视化，筛选阈值设定为P < 0.05。
构建诊断诺模图
针对三个关键基因分别绘制了ROC曲线，以评估关键基因在区分KIRC与对照样本方面的有效性。曲线用以下两个参数绘制：真阳性率和假阳性率。AUC用于表示预测的准确性和灵敏度。图4(a)，图4(b)中，训练集(GSE53757)和验证集(GSE40435)中三个基因的AUC值均大于0.9。
最后，利用“rms”软件包构建了预测KIRC诊断的诺模图模型。校准曲线结果显示，诺模图模型与理想模型的预测结果基本一致，根据该模型做出的临床决策可能对KIRC患者最有利。
免疫微环境分析与药物预测
CIBERSORT算法得到的结果如图4(a)，图4(b)所示。B细胞记忆、浆细胞、T细胞CD8、T细胞CD4、NK细胞静止、NK细胞活化、单核细胞、巨噬细胞M1等17种免疫细胞在疾病与对照组之间的评分差异，提示17种免疫细胞可能影响疾病的进展。
进一步计算17个差异免疫细胞与3个关键基因的相关系数，分析其意义。PDGFRB、TNFRSF10B和PRF1基因与NK细胞活化、浆细胞和T细胞调节(Tregs)、嗜酸性粒细胞和B细胞记忆呈负相关，与巨噬细胞M1和树突状细胞活化呈正相关。
为了分析化疗药物对相关基因的影响，我们通过DGIDB (https://dgidb.org/)数据库构建了药物–基因相互作用网络，并将3个核心基因录入DGIDB数据库。结果发现，PDGFRB基因预测71种药物，TNFRSF10B基因预测7种药物，PRF1基因预测2种药物。
LncRNA-miRNA-mRNA网络
我们用mirwalk软件预测可能与筛选出的关键基因相关的miRNA，然后用starbase预测关键基因的miRNA。我们选择了两个数据库中预测到的miRNA进行后续分析，结果发现共有6个miRNA同时存在于两个数据库中。在这6个miRNA的基础上，我们使用starbase和mirnet (默认参数)软件预测了可能与miRNA相互作用的lncRNA，同时也只保留了两个数据库中预测的lncRNA，结果共有175个lncRNA。根据上述关系，我们保留了同时存在mRNA-miRNA-lncRNA调控关系的基因，构建了一个由6个miRNA、3个mRNA和175个lncRNA组成的、有184个节点和296条边的ceRNA (mRNA-miRNA-lncRNA)调控网络。
基因表达水平的验证
为了进一步验证3个关键基因在正常组织和病变组织中的表达情况，我们对3个关键基因的表达情况进行了分析，3个关键基因在癌症组中都有高表达。同时，我们在验证集中对这3个基因进行了验证，结果表明这3个基因在验证集中可以得到验证，并且与训练集和TCGA数据库中的表达模式相同。
通过RT-qPCR，我们利用KIRC临床样本检测了所选诊断特征的表达水平。如表3所示，KIRC肿瘤组织中含有高水平的PDGFRB (P = 0.0273)和PRF1 (P = 0.0397)，而对照组中这些标记基因的表达量较低。然而，TNFRSF10B (P = 0.6426)基因的表达结果不显著。
(a)
(b)
Figure 4. (a) Immune cell content stacking plot between tumor group and normal group. (b) CIBERSORT algorithm to evaluate the abundant difference of immune cells between tumor group and normal group
图4. (a) 肿瘤组与正常组免疫细胞含量叠加图。(b) CIBERSORT算法评价肿瘤组与正常组免疫细胞丰度差异
Table 3. RT-qPCR results of marker genes (PRF1, TNFRSF10B and PDGFRB)
表3. 标记基因(PRF1, TNFRSF10B和PDGFRB)的RT-qPCR结果
4. 讨论
KIRC是一种常见的泌尿生殖系统恶性肿瘤，早期发现和诊断相当困难，其诊断需要结合临床症状、影像学和实验室检查等多学科方法 [28] [29] [30] 。目前，我们一般采用造影剂增强计算机断层扫描(CT)作为诊断肾癌的直接依据，其对≤4厘米肾小肿瘤的预测准确率仅为79.4%，对≤1.5厘米肾小肿瘤的诊断率更低 [31] 。延迟诊断容易导致各种严重并发症和不良预后，并进一步导致对患者监测的减少 [32] 。本研究利用GEO公共数据库中的KIRC患者转录组数据和IEU OpenGWAS中的KIRC SNP数据进行生物信息学分析，全面鉴定DE-PRGs。为了探索与DE-PRGs相关的生物学功能和通路，研究人员利用功能富集分析建立了PPI网络。
对66个eQTL基因进行了MR分析，结果显示3个基因具有显著的因果关系。通过单变量和多变量MR分析，我们可以得到以下信息：1) 在单变量分析中，三个暴露因子均有显著性，PRF1和TNFRSF10B是KIRC疾病的危险因子，PDGFRB是保护因子；2) 在多变量分析中，TNFRSF10B暴露因子的安全性/危害性属性发生了变化，而PRF1、TNFRSF10B的显著性下降，变得不显著；3) 从第2点可以推断，PDGFRB是KIRC疾病的保护因子。从第2点可以推断，PDGFRB、TNFRSF10B和PRF1基因单变量显著，但在多变量MR中，校正PDGFRB的影响后，PRF1和TNFRSF10B基因从显著因子变为不显著因子。即3个暴露因素同时存在时，PDGFRB基因是一个直接因素。
诊断模型包括3个重要的生物标志物(PDGFRB、TNFRSF10B和PRF1)，所有基因的ROC曲线下面积所代表的AUC值均大于0.9，表明这3个重要的生物标志物在区分KIRC样本和健康样本方面具有一定的准确性。结合与免疫微环境变化相关的结果，PDGFRB、TNFRSF10B和PRF1可能与KIRC的诊断和免疫治疗密切相关。PDGFRB是一种跨膜受体酪氨酸激酶，与相应配体结合后会磷酸化、内化并激活下游信号转导分子的多个对接位点 [33] 。它被认为是KIRC组蛋白内酯化的关键靶基因 [34] 。此外，最新研究发现，PDGFRB信号通过激活PI3K/AKT和Ras/MAPK通路促进肿瘤细胞增殖、转移和血管生成 [35] [36] [37] ，在胃肠道肿瘤、肝细胞癌、胰腺癌等多种肿瘤中的表达均有所增加 [38] [39] 。
穿孔素1 (PRF1)编码一种在结构上与补体成分C9相似的蛋白质，可形成膜孔并释放颗粒酶，导致靶细胞的细胞溶解，在免疫中发挥重要作用 [40] 。穿孔素相关蛋白的X射线晶体结构表明 [41] ，这些蛋白的孔形成结构域与细菌胆固醇依赖性细胞溶解素(CDCs)同源。cdc最初结合并在膜表面寡聚形成一个由20~50个单体组成的预孔。在孔前组装之后，每个单体中的两个小簇α-螺旋(跨膜螺旋区1 (TMH1)和TMH2)展开并作为一对两亲性β-发夹插入膜中。因此，一个完整的孔隙由80~200股的大β桶组成 [42] 。此外，穿孔素本身的晶体结构，以及28Å分辨率的孔的低温电镜结构，为穿孔素通过类似cdc的机制形成孔的观点提供了广泛的支持 [43] 。
肿瘤坏死因子受体超家族成员10B (TNFRSF10B)是一种属于肿瘤坏死因子受体(tumor necrosis factor receptor, TNFR)超家族的蛋白质 [44] 。当细胞死亡时，其胞质在受体三聚化时形成死亡诱导信号复合体(death-inducing signaling complex, DISC) [45] 。由于其对恶性细胞的靶向选择性高于正常细胞，促凋亡激动剂的开发主要集中在TNFRSF10B上 [46] [47] 。同时，我们对验证集中的3个基因进行了验证，结果表明，验证集中的3个基因与训练集中的表达模式相同。
在RT-qPCR验证过程中，我们通过数据分析得到了一个不同的结果，即肿瘤组织中TNFRSF10B基因的表达量低于正常组织。一方面，我们认为这是一个偶然的结果，因为样本数量较少，只有3对。另一方面，在总结以往文献的过程中，我们进一步发现近年来抗TNFRSF10B抗体的表达，导致肿瘤组织中TNFRSF10B基因的表达量有下降的趋势 [48] [49] ，这可能也会对我们的研究产生一定的影响。
本研究存在一些局限性。首先，本研究选择的基因来自于之前描述的免疫相关遗传特征，因此我们可能遗漏了研究方案中未包括的具有预测作用的基因。其次，虽然我们通过RT-qPCR进行了初步的实验验证，但还需要独立的前瞻性临床研究来证实综合诊断模型的能力，并进一步研究导致肾癌分子机制改变的突变基因。
5. 结论
综上，我们的研究发现，PDGFRB、TNFRSF10B、PRF1基因与KIRC疾病之间存在因果关系，当3个标记基因同时存在时，PDGFRB基因是直接因素。基于标记基因，我们建立了一个具有良好诊断价值的KIRC诊断模型。这项研究将有助于开发KIRC的早期诊断指标和基于PANoptosis基因的免疫治疗方法。
NOTES
*通讯作者。