线粒体功能相关基因作为结直肠癌诊断与预后标志物的应用研究

周凯翔，袁 晴，王珍妮，郭姗姗，苏丽萍，刘 洋，郭 旭，邢金良，黄启超

(1空军军医大学基础医学院生理与病理生理学教研室，陕西 西安 710032； 2西北工业大学医学研究院，陕西 西安 710072)

结直肠癌(colorectal cancer，CRC)是世界范围内发病率和死亡率最高的常见癌症之一。最新数据显示，CRC的发病率和死亡率在全球所有癌症患者中分别位列第三和第二位[1]。与大多数癌症一样，CRC的早期诊断有助于良好的预后，已有流行病学调查结果表明早期发现疾病可使五年生存率提高到90%[2-3]。然而由于缺乏敏感有效的诊断方法，CRC通常到晚期才被诊断出来，极大地影响了CRC患者的临床治疗和生存质量。尽管近年来在CRC诊断和治疗方面有了显著的改善，但被诊断为转移性CRC患者的5年生存率仍然很低，约为12年[4]。因此，迫切需要了解CRC发生发展的分子机制，并寻找新的生物标志物用于CRC的早期诊断和预后评估。

既往研究证实CRC是一种依赖于众多致癌基因和抑癌基因改变的遗传性疾病[5]。大量研究致力于鉴定CRC相关基因及其编码蛋白在包括细胞增殖、分化、凋亡和转移在内的大量生理和病理过程中的重要作用。然而，CRC的确切分子机制仍有待深入研究。线粒体作为真核细胞内重要的半自主细胞器，含有自身的线粒体DNA并编码呼吸链的13种必需蛋白质[6-7]。已有研究报道CRC细胞以线粒体氧化磷酸化作为其主要能量来源。此外，在CRC发生过程中，线粒体在调节细胞增殖和凋亡中起核心作用。除了由线粒体基因组编码的37个基因外，参与调节线粒体结构和功能的核基因超过1 600个[8]。因此，核基因组和线粒体基因组中许多基因的突变会导致线粒体功能障碍，从而导致包括肿瘤在内的许多严重疾病。然而，线粒体相关基因在CRC发生发展及预后预测中的具体作用尚不明确。因此，阐明线粒体相关基因的异常表达机制，并探索这些基因作为潜在生物标志物在CRC患者的有效诊断和预后评估中的具体作用，对揭示CRC新的诊断治疗靶点和预后策略是十分重要的。

本研究基于肿瘤基因图谱(The Cancer Genome Atlas，TCGA)数据库中CRC患者的转录组数据集，分析了线粒体功能相关基因的表达差异及其潜在的生物学功能，并将线粒体功能相关基因纳入诊断和预后评估的模型，进而建立了联合指标的CRC诊断和预后预测模型，为CRC患者的诊断和预后评估提供临床新思路。

1 材料与方法

1.1 材料

本文数据来源于TCGA数据库，分别下载了结肠癌和直肠癌的RNA-seq counts表达谱数据及对应的临床信息。其中结肠癌包括478例肿瘤样本和41例正常样本，直肠癌包括166例肿瘤样本和10例正常样本。经数据汇总后，本文分析的CRC数据包括644例肿瘤样本和51例正常样本。

线粒体功能相关基因信息来源于MitoMiner数据库[9]，综合线粒体蛋白指数数据库收录了1 626个线粒体功能相关基因的信息，其中包含了13个线粒体基因组编码基因。本文将非线粒体编码的1 613个线粒体相关基因及国家生物技术中心信息检索数据库线粒体参考基因组中的37个线粒体基因统称为线粒体功能相关基因。对数据进行质控的结果表明，本研究所用数据集中的基因对线粒体功能相关基因的覆盖率为100%。

1.2 方法

1.2.1 表达谱数据的差异基因筛选 对CRC的RNA-seq counts数据进行差异表达基因筛选，方法包括t检验和倍数变化(fold change，FC)，之后采用多重检验校正的方法Benjamini-Hochberg对P值进行校正。本文通过R软件的DESeq2程序包对RNA-seq counts数据进行标准化处理后，再分析差异基因，其中符合P.adj<0.05且|log2FC|≥1的基因在本文被定义为差异表达基因。

1.2.2 对差异表达线粒体功能相关基因进行注释 为更进一步分析差异表达线粒体功能相关基因的生物学通路，本文通过R软件的org.Hs.eg.db、clusterProfiler、DOSE、stringr等程序包进行基因本体论(gene ontology，GO)功能注释[10]及京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)功能富集分析[11]，其中GO注释为线粒体功能相关基因参与的生物学过程。采用多重检验校正的方法Benjamini-Hochberg对P值进行校正，设置满足P.adj<0.01的GO术语和KEGG通路具有显著性。

1.2.3 随机森林构建诊断模型 将差异的线粒体功能相关基因用作构建诊断模型的不同特征，随机森林的实现采用R语言caret程序包筛选与诊断相关的最优线粒体相关特征基因。十折交叉验证的方法，将数据按照9∶1的比例随机分成训练集和验证集，最后得到最优的特征基因。

1.2.4 单因素Cox回归分析 首先对临床数据中的年龄、病理分期变量进行单因素Cox回归分析，再将单因素分析中P<0.05的变量进行保留。其次，将差异表达的11个线粒体基因作为单因素分析中的变量，再将其中P<0.05的变量进行保留。选取显著的独立变量纳入模型构建并使用受试者工作特征(receiver operating characteristic，ROC)曲线下面积(area under curve，AUC)进行效能评价。

1.2.5 生存分析 将显著的独立变量纳入模型，采用R语言的predict函数对变量进行风险打分，高于中位数为高风险，低于中位数为低风险。最后用R语言的Survival和Survminer程序包绘制生存曲线。

2 结果

2.1 CRC差异表达基因筛选

在包含CRC组644例(结肠患者478例，直肠患者166例)和正常组51例的基因表达谱数据中，采用T检验和FC的差异基因分析方法设置P.adj<0.05且|log2FC|≥1的阈值，共筛选出15 115个CRC相关的差异表达基因，上调基因10 762个，下调基因4 353个。其中，差异表达的线粒体功能相关基因有257个，上调136个，下调121个(图1)。图1展示了1 650个线粒体功能相关基因在CRC组和正常组中基因表达值的情况。差异表达的11个线粒体基因的详细信息见表1，其中MT-ND6基因为下调基因，剩余10个基因为上调基因。

图1 线粒体功能相关基因火山图

表1 11个差异表达线粒体基因

2.2 CRC差异表达线粒体功能相关基因的功能分析

对257个差异表达线粒体功能相关基因进行GO功能注释和KEGG功能富集分析，结果分别富集到2 910条GO生物学通路和205条KEGG通路。符合阈值P.adj<0.01的生物学术语154条(图2展示了最显著的15条)，它们主要与小分子分解代谢过程、含嘌呤化合物代谢过程和有机酸分解代谢过程等有关；符合阈值P.adj<0.01的KEGG通路13条(图3)，它们主要与脂肪酸代谢、精氨酸和脯氨酸代谢及脂肪酸降解等有关。

纵坐标表示不同的GO术语，横坐标表示注释到某条GO术语的基因比例。图中不同P.adj值表示为不同色阶，红色表示P.adj值较小，蓝色表示P.adj值较大。图2 差异表达线粒体功能相关基因的GO功能注释

纵坐标表示不同的KEGG通路，横坐标表示注释到某条KEGG通路的基因比例。图中不同P.adj值表示为不同色阶，红色表示P.adj值较小，蓝色表示P.adj值较大。图3 差异表达线粒体功能相关基因的KEGG功能富集

2.3 诊断特征筛选

十折交叉验证将数据按照9∶1的比例随机分成训练集和验证集，得到8个最优特征基因(图4A)。随后对其进行评估，绘制出ROC曲线，其AUC值为0.919(图4B)。

A：随机森林筛选后的特征基因；B：特征基因的ROC曲线。图4 诊断特征筛选结果

2.4 生存率的独立预测因素

644例肿瘤样本中纳入生存分析的CRC患者600例，选取P<0.05的因素作为模型构建，通过单因素Cox回归分析，发现年龄 ≥ 67岁、Stage Ⅲ & Ⅳ和MT-TL1基因是生存曲线的独立危险因素(表2)。

表2 CRC患者预后单因素分析

2.5 预后评估模型的构建及评估

通过对单因素分析显著的变量，年龄、Stage分期、MT-TL1基因表达值模型的建立，采用R语言的predict函数对多因素进行风险打分，高于中位数为高风险组，低于中位数为低风险组。通过风险得分来绘制3年和5年的ROC曲线(图5)，结果表明3年生存的AUC值为0.749，5年生存的AUC值为0.721。

图5 3年生存和5年生存的ROC曲线

2.6 生存曲线绘制

经过模型构建的评估，并结合临床信息将600例CRC划分高风险组(221例)和低风险组(379例)，绘制总生存曲线(P<0.000 1，图6)。此结果表明高风险组生存较差，低风险组生存较好。

图6 总生存曲线

3 讨论

在本研究中，我们综合分析了来自TCGA的全部肠癌相关转录组测序数据。CRC组织与正常组织共检测到257个线粒体相关基因表达，其中上调基因表达136个，下调基因表达121个。功能富集分析表明，线粒体相关差异基因在小分子分解代谢、含嘌呤化合物代谢和脂肪酸代谢等生物学过程中富集。之前研究报道，将羧酸酯酶1(carboxylesterase 1，CES1)鉴定为一种必需的NF-κB调节的脂肪酶，它将肥胖相关炎症与脂肪代谢和侵袭性CRC中的能量压力适应联系起来。CES1通过促进脂肪酸氧化的细胞自主机制促进CRC细胞存活并防止三酰基甘油的毒性积聚。研究人员发现CES1表达升高与超重CRC患者的较差结果相关[12-13]。另外，胆汁酸是由肝脏中的胆固醇合成的生理产物，具有产生胆汁流动以及促进肠道对脂质、营养素和维生素的吸收和运输的作用。在CRC的癌变和发展过程中，胆汁酸似乎起着至关重要的作用。研究人员揭示了石胆酸和脱氧胆酸在N-甲基-N′-硝基-N-亚硝基胍治疗后促进大鼠CRC肿廇生长。BERNSTEIN等[14]进一步证实了脱氧胆酸在小鼠实验中的作用；他们用含2 mL/L脱氧胆酸的饮食喂养18只小鼠，8～10个月后，其中10只患上CRC。近年来，越来越多的研究发现次级胆汁酸石胆酸和脱氧胆酸是参与CRC发生和发展的主要因素。我们的研究结果表明，这些线粒体相关差异基因可能参与了CRC的不同代谢过程，从而导致肿瘤的发生和进展。

近年来，许多研究表明机器学习技术在癌症诊断领域具有广泛前景。例如，PU等[15-16]利用支持向量机方法确定了一个基于多个超甲基化CpG位点的诊断模型，准确率为0.82%。本研究以695个样本的257个差异线粒体相关基因的表达谱为原始数据，通过随机森林构建诊断模型。经过十折交叉验证的方法，结果表明8个线粒体相关基因集合(ABCA8、SLC25A34、MT-ND6、P2RY1、ABCG2、MT-TV、MT-TL1和ADHFE1)对CRC肿廇生长有较好的诊断性能，其AUC值为0.919。有趣的是，研究人员最近从一小部分患者中报道了结肠腺瘤(带蒂息肉)中ABCG2相对于相邻的正常结肠黏膜的下调。这被认为有利于致癌物质在肠道中的积累，从而促进CRC的形成。

许多研究已经基于肿瘤相关差异基因构建了一些预后模型来预测肿瘤患者的预后。例如，QIU等[17-18]基于生物信息学分析发现了口腔舌鳞状细胞癌患者的16个基因预后特征，其预测5年生存率的AUC为0.872。LIU等[19-20]通过单因素和多因素Cox比例风险回归模型研究胃癌预后相关基因，他们发现一组由9个基因组成的预后基因标记，其与胃癌患者的预后有显著相关性。在本研究中，为了评估CRC患者线粒体相关差异基因与临床结局之间的关系，我们首先进行了单因素Cox比例风险回归，结果发现了年龄、肿瘤病理分期和MT-TL1基因与生存相关。最后，我们通过风险预测模型建立了基于上述三个因素的生存预测模型。经风险评分分层后，预后特征的Kaplan-Meier曲线显示，CRC患者的临床结局存在显著差异(P<0.000 1)。ROC分析获得的AUC 3年和5年生存预测分别为0.749和0.721。这表明本CRC预后模型的特异性和敏感性相对较高。

本研究具有许多重要优势：第一，本研究纳入了TCGA数据集中全部CRC相关的测序数据，使用的数据非常全面；第二，应用基于随机森林的机器学习方法和Cox回归模型来构建线粒体相关差异基因在CRC患者中的诊断及预后模型，应用ROC曲线较为稳健地评估效能。然而我们的研究仍然存在一些局限性：第一，本研究的研究样本来自白种人、亚洲人、黑人或非裔美国人等不同人群，不同民族间基因表达谱可能存在一定差异。为了弥补这一缺陷，我们未来将继续收集CRC患者临床数据及肿瘤组织样本，利用独立的实验数据进一步验证我们的预测模型。第二，线粒体相关差异基因的诊断效率和预后价值仅在TCGA数据集中进行了分析和验证。因此，未来的研究中我们将进一步在其他数据库中验证本结果。第三，我们的研究方法主要基于生物信息学分析，需要进一步的实验来验证基于肿瘤样本和临床数据的结果。此外，体内和体外实验可以增强我们对CRC中线粒体相关差异基因功能作用的认识。

综上所述，本研究通过在大样本中分析CRC的线粒体功能相关基因的表达，阐述了线粒体功能在CRC发生发展中的可能作用，结合临床指标构建并评估了CRC诊断和预后模型，为CRC患者的预后判断提供了有效的新方法。