一种基于转录组数据和机器学习策略评估细胞衰老程度的方法

摘要

本发明属于计算基因组学技术领域，具体为一种基于转录组数据和机器学习策略评估细胞衰老程度的方法。本发明通过对已知的转录组数据进行模型的训练拟合得到细胞衰老评分模型，从而达到仅利用组织样本转录组测序的数据就预测其中细胞的衰老程度的目的。其细胞衰老预测的精确程度可以达到AUC值为0.87，要高于目前已有的四种细胞衰老预测工具。

说明书

技术领域

本发明属于计算基因组学技术领域，具体涉及一种基于转录组数据和机器学习策略评估细胞衰老程度的方法。

背景技术

衰老是指随着时间的推移，机体中各器官组织功能逐渐减弱的过程，通常也伴随着组织中细胞的不断衰老。细胞的衰老也会通过多种方式影响个体的衰老，增加衰老相关疾病（比如心血管疾病、神经退行性疾病及癌症等）的发病风险。因此，了解个体衰老过程中细胞层面的衰老以及准确评估细胞衰老的程度及比例对于疾病的预防与治疗具有重要的意义。同时清除个体中衰老的细胞也被证明可以延缓衰老并增加健康寿命和个体寿命，进一步提示了鉴定衰老细胞并对其进行靶向的重要性。

细胞衰老是由各种刺激的持续压力和损伤引起的，并导致细胞周期的永久性阻滞。科学家Hayflick和Moorhead最早发现细胞培养系统中人胚胎成纤维细胞的增殖潜力是有限的，培养条件下的细胞所达到的增殖极限因此也被称为海弗利克极限。这些衰老的细胞变得大而扁平，但仍有代谢活力。后来发现细胞的衰老可以分为两大类，分别是复制性衰老和诱导型早衰。培养的细胞随着代数的增加增殖速度变慢的现象被称为复制性衰老(RS)，而由于各类压力（如原癌基因诱导、阿霉素诱导、双氧水诱导、血管紧张素II诱导、高糖诱导等）所导致的类型被称为诱导型早衰。癌细胞在特定药物诱导下也可以进入细胞衰老状态，进而通过细胞凋亡或者被体内的免疫细胞识别并清除，因此不管对正常组织的细胞衰老以及癌组织中的细胞衰老状态和程度的评估均具有重要健康价值。

当前抗癌策略背后的基本原理是通过高剂量药物或辐射引起广泛的DNA损伤来杀死快速分裂的癌细胞。然而，抗癌治疗的延迟副作用，如复发、继发性癌症和化疗和放疗引起的正常组织损伤，给癌症幸存者带来了临床问题。研究表明，用化疗或放疗处理癌细胞会产生衰老状态，称为治疗诱导的细胞衰老(TIS)。后续的机制研究发现其主要涉及p53/p21和p16/pRb这两条信号通路。有意思的是，高浓度的阿霉素会诱导人类癌细胞凋亡，而低浓度则会诱导癌细胞衰老。用于治疗人类癌症的细胞衰老诱导剂具有临床意义。电离辐射也会诱导癌细胞的衰老。放射治疗研究的主要目标之一是开发更有效的方法来提高放射治疗的疗效而不会对正常组织造成毒性。因此，旨在选择性诱导细胞衰老的治疗方法可能代表一种有前途的癌症治疗新策略。综上所述，评估病人的肿瘤细胞衰老程度，对于不管是放化疗疗法，还是免疫治疗或者靶向治疗，均是评估病人预后好坏的重要指标，简单并有效的细胞衰老评估方法具有重要价值。

衰老的细胞改变了许多基因的表达，同时分泌各种细胞因子、趋化因子和酶（细胞衰老相关分泌表型），这为开发细胞衰老生物标志物创造了重要机会。之前的研究已经发现了几个与细胞衰老有关的重要特征，比如经典的细胞衰老标志物细胞衰老相关的 β-半乳糖苷酶活性增加，以SA-β-Gal染色变蓝的细胞增多为特点。由于细胞衰老是永久性的细胞周期阻滞，一些细胞周期的负调节因子，如p15、p16、p21和p27等，其RNA和蛋白质水平在细胞衰老中水平上升，可作为细胞衰老的另一个重要分子指标。细胞活力标记基因

有监督模型广泛适用于基于数据点特征来预测标签（即细胞的衰老程度），但传统的有监督模型需要存在两个或更多类数据来训练模型。但在预测细胞是否衰老这种情况下，没有明确的否定类，只有一组我们希望检测的类和一些未知的类。通常，我们希望将特定亚型（衰老细胞）与所有/任何其他亚型进行对比，而不是针对某一种。因而依赖于传统的有监督方法并不能很好的解决该问题。

随着对基因的转录产物——RNA的测序技术（RNA-seq）的飞速发展，对组织、不同类型的细胞甚至单个细胞的测序费用迅速下降，测序的覆盖度也不断加深，使得利用组织或细胞的转录组数据及其分析得到的基因表达谱来评估组织中不同类细胞的衰老状态成为可能。

发明内容

为了克服现有技术中的标志物非特异性，以及检测多种指标操作繁琐且实验复杂的缺点，本发明提出了一种基于转录组数据结合机器学习策略对人的细胞或者组织的衰老状态进行评估的方法，这种方法基于公开发表的有确凿证据的多种细胞衰老模型的转录组数据来建立细胞衰老评分模型，进而从基因表达谱可靠地推断细胞的衰老程度。本发明通过创新地使用单分类逻辑回归（one-class logistic regression，OCLR）机器学习算法，从多种类型的衰老细胞表达数据集里提取细胞衰老共同的基因表达特征，从而相比于已有的细胞衰老相关标志物和计算方法获得更好的定量性能。

本发明提供的评估细胞的衰老程度的算法，采用基于机器学习算法对RNA-seq基因表达谱数据进行细胞衰老的打分预测技术，核心在于用单分类逻辑回归模型结合人的多种类型细胞衰老体系的基因表达特征，从而可以习得数据中隐藏的与细胞衰老密切相关的重要分子特征。同时考虑到由于细胞衰老是一个连续的过程，标记为非衰老的细胞在表达谱上也可以有衰老的特征，因而我们利用Spearman回归，结合基因权重向量和表达谱用来表征一个能代表细胞衰老指数的分值，该分值越大，细胞衰老程度越高。我们将这个分值定义为人类细胞衰老指数（human senescence score, hSI）。本发明方法的具体方案介绍如下。

一种基于转录组数据和机器学习策略评估细胞衰老程度的方法，具体步骤如下：

1）收集并分析公共数据库多种衰老细胞的转录组RNA-seq数据，得到衰老细胞表达谱组成的衰老细胞表达矩阵，所述衰老细胞涵盖若干种细胞类型和衰老类型；

2）以衰老细胞表达矩阵作为模型输入，先将衰老细胞表达矩阵中每个基因的所有样本的表达水平统一减去对应基因表达水平的平均值，再基于单分类逻辑回归算法训练拟合，构建细胞衰老评分模型，得到细胞衰老相关基因权重向量；

3）对待评估衰老状态的细胞进行RNA-seq测序，分析得到其基因表达谱；

4）计算待评估衰老状态的细胞的基因表达谱和细胞衰老评分模型的权重向量的Spearman相关系数，得到人类细胞衰老指数hSI，hSI值越大，细胞衰老程度越高。

本发明中，步骤1）中，细胞类型包括小鼠腹水型肝癌细胞、黑色素瘤细胞、角质形成细胞、人胚肺成纤维细胞、人皮肤成纤维细胞、星型胶质细胞、人包皮成纤维细胞、人胚肺成纤维细胞、真皮成纤维细胞、人脐静脉内皮细胞、人主动脉内皮细胞，细胞衰老类型包括电离辐射诱导的细胞衰老、复制性细胞衰老、氧化应激诱导的细胞衰老、原癌基因诱导的细胞衰老、阿霉素诱导的细胞衰老、个体老化中伴随的细胞衰老。

本发明中，步骤1）和步骤3）中，转录组RNA-seq数据的分析方法包括如下步骤：

首先使用Trim Galore程序过滤低质量的reads，并且去除3′末端低质量的碱基；接着利用STAR软件将经过质量控制的短读长比对到GRCh38人类参考基因组上去，取唯一比对的短读长进行后续分析；再使用StringTie计算每个样本的基因表达水平，表达量的标准化数值使用TPM，根据Gencode注释保留蛋白编码基因，并且去除在99%以上样本中TPM<3的低表达的基因。

本发明中，步骤1）中，RNA-seq数据的分析方法还包括最后采用ComBat工具处理表达数据，以减少批次效应对模型训练的影响的步骤。

本发明利用目前广泛使用的RNA-seq测序技术，对基因表达进行定量后可快速的推断样本的细胞衰老情况。基于目前开源的各种细胞类型的细胞衰老的表达数据，采用单分类逻辑回归模型进行参数拟合，最终得到了较为可靠的评估样本细胞衰老分值的细胞衰老评分模型。和现有技术相比，本发明的有益效果在于：

（1）利用机器学习模型和海量的开源数据，可以不依赖于细胞分选和单细胞测序，以及细胞衰老相关半乳糖苷酶染色（SA-β-Gal）等生化实验即可对组织的每种细胞类型衰老情况进行评估。

（2）本发明方法可靠、方便、快捷，可以在给定样本表达谱的时候精确地评估细胞衰老分值。

附图说明

图1是本发明方法的流程示意图。

图2是在2个独立的基于RNA-seq的衰老和非衰老细胞数据集里测试人类细胞衰老指数（hSI）的可靠性。（a）三种化合物（Adria、H

图3是细胞衰老评分模型计算的hSI与已报道的11个细胞衰老相关的标志基因在预测细胞衰老表型中的性能对比。hSI与衰老标志基因在图上的顺序按AUC的均值从高到低排列。误差线通过计算30个数据单元的AUC的均值和标准差得到。

图4是细胞衰老评分模型计算的hSI与其他四种细胞衰老评估模型（基于DNA损伤相关细胞衰老特征（DAS）[1]，基于修饰分泌细胞衰老特征（mSS）[1]，基于DNA损伤衰老和分泌衰老特征（DAS+mSS）[1]，基于复制性衰老特征（Sig.RS）[2-3]）的性能比较。每个点表示一个数据单元中某种细胞衰老评估策略的AUC值。***表示

图5是10种细胞类型在年轻和衰老皮肤组织中的细胞衰老程度分布。细胞类型包括表皮干细胞（epidermal stem cell）、红细胞（erythrocytes）、成纤维细胞（fibroblasts）、角质细胞（keratinocytes）、淋巴内皮细胞（lymphatic endothelial）、巨噬细胞（macrophage）、黑色素细胞（melanocytes）、周细胞（pericytes）、T细胞（T cell）和血管内皮细胞（vascular endothelial）。NS与***分别表示P > 0.05和P < 0.001，双端Wilcoxon秩和检验。

图6是细胞衰老评分模型在肺纤维化单细胞数据上的应用。（a）肺纤维化病变肺组织和健康组织的整体细胞衰老程度分布。***表示P<0.001, 双端Wilcoxon秩和检验。（b）整合肺纤维化病变组织和健康肺组织的单细胞转录组。（c）高斯混合模型拟合病变和健康组织所有细胞的衰老分值并预测到衰老细胞和非衰老细胞两种状态。（d）肺纤维化病变和健康肺组织中，不同细胞类型的衰老状态细胞占比。*表示P < 0.05，NS表示无显著差异，双端t检验。

图7是利用细胞衰老评分模型鉴定黑色素瘤中肿瘤细胞的衰老亚群。（a）黑色素瘤中肿瘤细胞根据细胞衰老状态分成3个亚群，使用PHATE对三个亚群的细胞进行可视化。细胞衰老标志基因

具体实施方式

下面结合附图和实施例对本发明的技术方案进行详细阐述。

所有实例中数据来源均来自美国国家生物信息中心数据库（NCBI）。

实施例中，对RNA-seq数据进行分析的流程如下：首先使用Trim Galore程序过滤低质量的reads，并且去除3′末端低质量的碱基。接下来利用STAR软件将经过质量控制的短读长比对到GRCh38人类参考基因组上去，取唯一比对的短读长进行后续分析。再使用StringTie计算每个样本的基因表达水平，表达量的标准化数值使用TPM（transcript permillion）。根据Gencode（版本号 v31）注释保留蛋白编码基因，并且去除在99%以上样本中低表达（TPM<3）的基因。 而由于收集的转录组数据来源于不同的实验室，因此采用ComBat工具处理表达数据，以减少批次效应对模型训练的影响。

实施例2-4中，对目标组织中待测细胞进行RNA-seq测序的方法如下：

对目标组织进行取样，裂解以后，总RNA跟带着oligo d(T)探针的磁珠结合，洗脱并获取结合的mRNA，用镁离子溶液高温打碎mRNA，随机引物反转第一条cDNA，之后再合成第二条cDNA，获得双链cDNA，对双链cDNA末端修复，加A加接头，选择特定大小DNA片段，并通过聚合酶链式反应扩增并纯化，获得最终RNA-seq文库。再进行高通量测序。

实施例1：结合大量公开数据和机器学习模型构建人类细胞衰老预测模型

利用表1所示的公开发表的9项不同细胞类型、不同衰老类型的转录组测序RNA-seq数据，结合单分类逻辑回归模型，进行模型拟合，在训练集内部通过逐个保留单个衰老细胞样本，采用留一交叉验证法（leave-one-out cross-validation，LOOCV）检验细胞衰老评分模型的有效性。

表1、细胞衰老评分模型训练数据来源汇总表。

表中：Accession number代表数据集所在数据库的检索号（GEO及EBI），PMID代表数据集所属文章的在PubMed数据库的检索号。RS、OIS、OSIS、IRIS分别表示复制性衰老、癌基因诱导的细胞衰老、氧化压力诱导的细胞衰老和辐射诱导的细胞衰老，healthy aging表示正常个体衰老，Dox-induced表示阿霉素诱导的细胞衰老。HCA-2：小鼠腹水型肝癌细胞；Melanocytes：黑色素瘤细胞；Keratinocyte：角质形成细胞；IMR90：人胚肺成纤维细胞；BJ：人皮肤成纤维细胞；Astrocytes：星型胶质细胞；HFF：人包皮成纤维细胞；MRC-5：人胚肺成纤维细胞；Dermal fibroblast：真皮成纤维细胞；HUVEC：人脐静脉内皮细胞；HAEC：人主动脉内皮细胞；WI-38：人胚肺成纤维细胞。

每一次的检验将剩下的衰老细胞样本训练单分类逻辑回归模型，然后对保留的单个衰老细胞样本以及其他非衰老细胞样本计算人类细胞衰老指数（hSI）。

最后使用AUC（area under ROC curve, 位于受试者工作特征曲线下方的面积）对模型性能进行预评估，AUC即正样本得分高于其他负样本得分的概率。在训练集进行LOOCV检验后得到该细胞衰老评估模型的得分为AUC=0.95，表现很好。

实施例中，进一步利用本发明方法对两种基因表达谱类型来源的数据（基于测序的RNA-seq数据和基于杂交的基因芯片数据 [4-5]，这些数据均已知细胞的衰老和非衰老状态）进行细胞衰老打分评估，结果表明与非衰老样本相比，所有的衰老细胞样本都显示更高的衰老打分（图2），证明了这一方法的可靠性。由于我们测试的数据集中涉及的细胞类型和细胞衰老的诱导方式与训练数据集存在较大差异，因此也提示了本发明细胞衰老评分模型适用于不同的转录组测量平台和不同类型的衰老细胞转录组，具有较强的普适性。

实施例中，同时本发明在衰老和非衰老细胞共培养的单细胞转录组数据 [6]里测试细胞衰老评分模型的性能，该数据集包括体外癌基因诱导的初级衰老细胞和受到旁分泌影响而引发的次级衰老细胞，并且这些细胞的衰老表型都通过经典的SA-β-Gal染色验证。基于AUC（位于受试者工作特征曲线下方的面积）作为方法的性能评估指标，结果表明本方法在基于单细胞测序技术上的数据，评分性能可以达到AUC=0.87，并且在与11种细胞衰老相关标志物（图3）和另外四种对细胞衰老的评估方法（图4）的比较中，本方法均展现了最好的性能。以上结果证实了人类细胞衰老指数hSI方法是一种可靠的、方便的计算方法，可以在给定样本表达谱的时候精确地评估细胞衰老分值。

实施例2：利用人类细胞衰老指数hSI揭示衰老皮肤组织中细胞衰老状态的异质性

皮肤作为人体表面的保护屏障，其衰老受到内源性（如时间、遗传因素和激素）和外源性因素（如紫外线照射和污染）的共同影响。发明人选取了年轻（25和27岁）和年老（53、69和70岁）的人腹股沟皮肤的单细胞转录组测序数据，用hSI去预测在皮肤组织中不同细胞亚型在衰老和年轻个体中的变化。由图5所示，一部分细胞类型的细胞衰老程度更高，而另一些却没有观察到显著差异，提示了衰老皮肤组织中众多细胞发生衰老的进程是不同步的，这也与以往研究中发现的不同类型细胞衰老程度不同步类似。

实施例3：利用人类细胞衰老指数hSI揭示肺纤维化病变组织中参与到疾病进展过程的衰老细胞类型

肺纤维化是一种慢性和高致死性的衰老相关疾病，其特征是异常的纤维化瘢痕导致的肺功能受损。目前已知细胞衰老分泌表型（SASP）对邻近细胞会产生影响，部分地介导了肺纤维化的疾病进展，并且利用抗细胞衰老药物定向清除肺纤维化小鼠模型的衰老细胞能改善受损的肺功能。因而发明人将本方法应用到和肺纤维化相关组织的单细胞转录组数据上 [7]，以推测和肺纤维化相关的衰老病变细胞类型是哪些。从图6中可以发现肺纤维化的病变肺组织细胞（肺泡2型细胞AT2、club细胞）的衰老程度比正常组织要高。这一结果证明hSI可以应用到一些病变肺组织上。

实施例4：利用人类细胞衰老指数hSI帮助鉴定黑色素瘤中肿瘤细胞的衰老亚群

目前已有多项研究证明衰老肿瘤细胞与免疫识别清除作用之间的关系，因而黑色素瘤中衰老的肿瘤细胞可以作为免疫治疗的潜在靶标。发明人应用本发明在黑色素瘤单细胞数据中 [8]，计算了每个肿瘤细胞的衰老程度，同时利用高斯混合模型聚类算法和PHATE降维算法识别出衰老黑色素瘤细胞的连续的衰老轨迹，从图7上可以看到衰老相关的标志基因

以上利用本方法做了实施例2-4三个应用实例，分别是揭示衰老皮肤组织中细胞衰老状态的异质性（图5）、揭示肺纤维化病变组织中参与到疾病进展过程的衰老细胞类型（图6）和鉴定黑色素瘤中肿瘤细胞的衰老亚群（图7），进一步证明了本方法在评估细胞衰老中的有效性和实用性。

参考文献

[1] Lafferty-Whyte K, Bilsland A, Cairney C J, et al. Scoring ofsenescence signalling in multiple human tumour gene expression datasets,identification of a correlation between senescence score and drug toxicity inthe NCI60 panel and a pro-inflammatory signature correlating with survivaladvantage in peritoneal mesothelioma[J]. BMC genomics, 2010, 11(1): 1-16.

[2] Reyfman P A, Walter J M, Joshi N, et al. Single-celltranscriptomic analysis of human lung provides insights into the pathobiologyof pulmonary fibrosis[J]. American journal of respiratory and critical caremedicine, 2019, 199(12): 1517-1536.

[3] Barbie D A, Tamayo P, Boehm J S, et al. Systematic RNAinterference reveals that oncogenic KRAS-driven cancers require TBK1[J].Nature, 2009, 462(7269): 108-112.

[4] Purcell M, Kruger A, Tainsky M A. Gene expression profiling ofreplicative and induced senescence[J]. Cell Cycle, 2014, 13(24): 3927-3937.

[5] Sati S, Bonev B, Szabo Q, et al. 4D genome rewiring duringoncogene-induced and replicative senescence[J]. Molecular cell, 2020, 78(3):522-538. e9.

[6] Tang H, Geng A, Zhang T, et al. Single senescent cell sequencingreveals heterogeneity in senescent cells induced by telomere erosion[J].Protein & cell, 2019, 10(5): 370-375.

[7] Habermann A C, Gutierrez A J, Bui L T, et al. Single-cell RNAsequencing reveals profibrotic roles of distinct epithelial and mesenchymallineages in pulmonary fibrosis[J]. Science advances, 2020, 6(28): eaba1972.

[8] Tirosh I, Izar B, Prakadan S M, et al. Dissecting themulticellular ecosystem of metastatic melanoma by single-cell RNA-seq[J].Science, 2016, 352(6282): 189-196。