基于数据库分析牙髓干细胞转录组差异性和功能倾向性研究*

谭思逸 许翀 周珊羽 刘军

先前的研究已经揭示了在成年牙齿中仍能再生的干细胞群，即出生后人类牙髓中的血管周围牙髓干细胞（dental pulp stem cells,DPSCs），但从老年人牙齿中分离的DPSCs的多能性明显降低［1］。由此可见，DPSCs 对形成牙本质、维持牙齿营养与美观，提高盖髓术成功率尤为重要［2］。

由上可知，人类的DPSCs 参与成牙本质细胞产生的牙本质结构的再生［3-8］。对DPSCs再生能力的基本分子机制理解并掌握，能让我们有可能利用这些功能用于再生医学和组织工程等未来细胞治疗方法。从牙髓中分离出的干细胞已经成功在体外实验分化为成脂肪、软骨、成骨牙源性谱系［9-11］。随着高通量测序技术的发展，利用高通量转录组学技术对DPSCs 的整体基因表达水平分析得以实现。机器学习方法如随机森林等算法的应用，有助于对高维数据的分析和挖掘［12］。本研究对基因表达数据库GEO的GSE123973数据集里的DPSCs和骨髓间充质干细胞（bone marrow mesenchymal stem cells,BMMSCs）的RNA-Seq 数据进行分析，找到DPSCs 的差异高表达基因，通过蛋白互作网络在差异性高表达基因集里找到起关键作用的核心基因，并通过GO和Reactome对核心基因进行功能注释，以更好地诠释DPSCs 的特性及其涉及的主要分子作用通路。

1.材料与方法

1.1 RNA-Seq 数据检索 本研究在研究初始阶段对Gene Expression Omnibus（GEO）数据库里DPSC相关的高通量基因测序数据集进行检索。

本研究数据集的纳入标准是：研究内容是关于DPSC与其他来源的干细胞比较；考虑到DPSC的异质性，我们排除了DPSC样本小于等于3个的小样本研究；使用的是基于测序技术获得的RNA-Seq数据集。

截至2021年1月，我们在GEO 数据库检索得到涉及DPSC 的数据集有9个。基于以上标准，仅数据集GSE123973 满足以上三个要求。所以本研究只对GSE123973数据集进行分析。

1.2 RNA-Seq 数据 获取本研究通过R 包GEOquery［13］从数据库中获得GSE123973 数据集里的原始RNA-Seq 数据，该数据集共10个DPSCs样本和3个BMSCs 样本。10个DPSCs 样本中，包含7个原代培养的DPSCs 样本和3个商品DPSCs样本。原代培养的DPSCs 样本从人第三磨牙中分离牙髓，提取干细胞并原代培养后进行测序。3个人源DPSCs 商业细胞样本购于LONZA。人源BMSCs 商业细胞样本购于LONZA。测序样本的基因比对参考为人类基因组GRCH37。

1.3 差异基因筛选 对基因数据分布进行初步分析后，使用R 包DESeq2［14］进行基因表达数据的标准化和差异性分析。使用Wald test得出基因表达的差异是否满足原假设。通过多重检验BH 法对P值进行了校正。

1.4 蛋白互作网络的构建 通过STRING 数据库构建DPSCs 差异表达基因的蛋白质相互作用（protein-protein interaction networks，PPI）网络［15］。并对DPSCs 差异表达基因相关的蛋白之间的相互作用网络进行研究。

1.5 核心基因的功能注释 通过蛋白互作网络获得核心基因（图3）后，对核心基因进行分子通路相关性功能注释。应用GO［16］和Reactome［17］DPSCs核心基因集的功能倾向性。

1.6 关键调控基因分析 对方法1.3 所获得的差异基因进行调控基因分析。通过对差异基因中的转录因子作为随机森林网络的调控基因，构建决策树。通过对决策树中的性能和所含差异基因的P值进行打分，分值越高的决策树被认为是越关键的调控基因。基于该排名，应用Cytoscape 3.7 对排名第一的关键调控基因网络进行可视化。

图1 研究分析流程图

2.结果

2.1 RNA-Seq数据预处理 通过GEOquery共获得13个样本，DPSCs组10个样本和BMMSCs组3个样本。本研究分析使用的参照基因组序列为Ensembl GRCh37，每个样本包含63677个基因数据。对中位数计数低于3的基因进行过滤，获得20394个基因。DPSCs的差异表达将在这20394个基因中进行分析得出。

2.2 DPSCs大部分基因表达模式异于BMMSCs通过过滤低表达基因后，对基因数据进行差异表达分析，经过多重检验后，获得8200个差异表达基因（FDR＜0.01）。取其中组间差异性最大的100个基因进行热图聚类和可视化，发现DPSCs和BMMSCs各自得到明确的聚类（图2）。

图2 DPSCs与BMMSCs差异表达基因热图

2.3 差异基因对应的蛋白互作网络投射出核心基因 由于差异表达基因较多，在进行基因的功能注释前，需要找到关键的调控基因，以减少过多基因造成的背景噪音影响。把差异基因输入STRING数据库进行蛋白互作网络的构建，得到差异基因相关的蛋白互作网络。选取高互作得分的节点（confidence＞0.7）和应用MCL聚类（inflation parameter为5）结合，过滤无任何连接的节点，得到核心基因网络（图3）。

图3 DPSCs核心基因的蛋白互作网络

2.4 功能注释显示DPSCs 干细胞功能差异性基于前期蛋白互作网络鉴定到的核心基因集，通过Fish精确试验从GO数据库中得到其显著相关的注释信息，从而对核心基因集进行功能分类注释。把核心基因改变的倍数作为参数，得到主要GO条目及其相关基因的关系图，体现了其干细胞功能差异性（图4）。本研究也应用Reactome数据库进行通路研究，包括经典的代谢通路、信号转导、基因转录调控、细胞凋亡与疾病。把核心基因改变的倍数作为参数，得到主要Reactome条目及其相关基因的关系图（图5）。

图4 DPSC核心基因集参与的GO条目及网络

图5 DPSC核心基因集参与的Reactome条目及网络

2.5 DPSCs 干细胞核心基因调控网络 进一步对以上重要差异基因进行调控网络的分析，获得排名第一的调控基因DMBX1，由图6 左所示，与BMMSCs 样本相比，这个关键调控基因在DPSCs 样本中有明显上调的差异表达。进一步应用网络可视化软件cytoscape 对这个关键调控基因网络的连接关系进行可视化，得到其调控网络（图6 右）：

图6 DPSC核心调控基因DMBX1的差异表达及其调控网络

3.讨论

不少前期研究关注人类脱落乳牙干细胞（SHEDs）与DPSCs 的功能比较。在血管生成功能方面，Zhu SY等的研究发现SHEDs和DPSCs显示出与周细胞相似分化模式，为DPSCs替代周细胞的临床应用奠定了理论基础[18]。Wang HH等的研究发现，在早期和晚期传代时，SHEDs比DPSCs表现出更高的增殖率和成骨分化能力，其结果可能为选择安全有效的用于再生医学和治疗的SHEDs和DPSCs提供指导[19]。尽管这些研究让口腔科学研究受益匪浅，但系统地比较DPSCs和其他MSCs的研究还相对较少。前期研究表明，与DPSCs相似，BMMSCs的再生能力也随着年龄的增长会下降。老年人的骨髓间充质干细胞的增殖和分化能力与年轻人相比显著下降[20]。另外，BMMSCs可作为不同环境背景的干细胞对照，因此，BMMSCs是研究DPSCs差异性特别是口腔组织分化功能的良好对照。

前期研究多聚焦在DPSCs和其他来源的干细胞的差异表达基因的数量和特征上。本研究主要以基因互作网络为切入点，通过对数据库中的DPSCs和BMMSCs转录组数据进行比较分析，获得DPSCs的差异表达基因。这些差异基因在DPSCs中大部分已有表达，而在BMMSCs却处于静止表达阶段，说明相对于BMMSCs，DPSCs更像是一群处于分化阶段的细胞。应用差异表达基因构建蛋白互作网络，筛选核心基因。通过选取高得分的节点和应用MCL聚类分析结合，有助于挖掘核心的调控基因，这些基因及相关表达蛋白可能是DPSCs生长、分化过程中受影响或者发挥激活、抑制作用的关键位点。其中获得关键节点如细胞分裂周期蛋白7(cell division cyck 7，CDC7)，苯甲亚胺醇1（budding uninhibited by benzimidazoles 1,BUB1）等。对核心基因进行功能注释，获得与BMMSCs相比DPSCs的基因功能倾向性。GO功能注释显示DPSCs差异性GO分类主要集中在染色体区（GO:0098687，chromosomal region；GO:0000775,chromosome，centromeric region；GO:0000793，condensed chromosome），生长特异性（GO:0006260，DNA replication；GO:0006261，DNA-dependent，DNA replication），以及分化特异性（GO:0001503，ossification）。Reactome数据库比对结果显示，DNA复制相关的基因（R-HSA-68962，Activation of the pre-replicative complex；RHSA-69306，DNA Replication）仍是DPSCs核心基因集参与的通路。另外，细胞周期的改变也是Reactome数据库提示的DPSCs核心基因集参与的一系列通路。通过对分析挖掘，我们发现，DPSCs的高表达基因主要富集在GO和Reactome功能注释中的细胞周期（包括DNA复制）以及分化特性相关（骨细胞分化）。骨分化的影响表现了DPSCs的特性，证明了其具有干细胞多向分化能力的特征,在再生牙髓治疗中具有重要潜力[21,22]。最后，我们应用随机森林算法，构建基因调控网络，获得DPSCs的核心调控基因。本研究从差异基因整体表达情况和相关通路入手，通过对数据库层层深入的分析，引入无偏估计的随机森林算法处理高维度数据获得最有代表性的转录因子网络，模型泛化能力强。

在DPSCs和BMMSCs差异表达基因的蛋白互作网络的研究中，我们发现了网络关键节点如CDC7，BUB1等（图4，图5）。CDC7参与有丝分裂细胞周期G1/S转变和DNA复制启动。BUB1是纺锤体组装检查点最主要的组件蛋白之一，它作为监控机制的“感受器蛋白”，影响细胞的增殖、分化等过程。这些结果提示虽然DPSCs和BMMSCs都是干细胞，但其二者的细胞周期有所不同，主要表现在G1期。细胞周期功能是细胞进行细胞分裂、分化的基础，其功能的核心包括细胞周期蛋白、细胞周期蛋白依赖性激酶等[23]。越来越多的证据表明[24]，不同的细胞有不同的细胞周期运行机制，尤其体现在不同的干细胞上，不同的干细胞具备其独特细胞周期模式。并且，干细胞的自我更新和多向分化潜能等特性和功能都与其独特的细胞周期运行方式有关[25]。

最后，有别于以往研究从基因表达差异为依据的基因重要性排序方式，我们应用随机森林算法计算核心转录因子及其调控网络在DPSCs和BMMSCs之间的整体差异性。我们发现转录因子DMBX1是DPSCs最关键的转录调节因子。DMBX1是中枢神经系统发育相关的调控基因，还能调节细胞周期[26]。其在DPSCs的高表达考虑与牙髓神经的发育相关。另外，我们还发现其目标基因不仅包括细胞周期基因如周期蛋白依赖性激酶样2 (CDKL2)，还包括HLA II类组织相容性抗原HLA-DOB，说明该调控基因在DPSCs中很有可能具有免疫调控功能。

综上，本研究通过DPSCs与BMMSCs的差异基因调控网络分析，以转录因子和受调控基因作为节点，以调控关系作为边的有向网络，对转录基因与目标基因之间的关系网络进行精确的量化，从而更深入地阐述DPSCs与BMMSCs间的差异生物学整体特征。