一种分析细胞命运转换关键转录因子的方法及系统

1.本发明涉及生物技术领域，尤其涉及一种分析细胞命运转换关键转录因子的方法及系统。


背景技术：

2.通过过度表达特定的转录因子来转化细胞命运已经是再生医学中的一种常用手段。在这之前，确定细胞命运转化的关键因子就需要大量的实验，也正是因为如此仍有大量具备临床应用价值的细胞类型的转化尚未实现。以上仅是说明未发现的部分，而在许多已发表的案例中，转换后的细胞虽然已显示了目标细胞的表型和标志蛋白，但这种转换过程往往是不完整的，一些重要基因的表达往往未能得到有效的抑制或激活。
3.因此，如何鉴定细胞实现完整转换所需的关键转录因子是将该领域技术进一步推向临床应用的关键。在鉴定过程中，目前，集成多组学分析方法已经具备一定掌握细胞命运转化过程中基因调控机制的能力，由此延伸的内容在于如何通过分析由该种方法获取的多种组学数据获得更准确的预测结果，是当前急需解决的技术问题。


技术实现要素：

4.为至少解决现有技术中存在的技术问题之一，本发明的目的在于提供一种分析细胞命运转换关键转录因子的方法及系统。
5.根据本发明实施例的第一方面，一种分析细胞命运转换关键转录因子的方法，包括以下步骤：
6.获取初始细胞和目标细胞；
7.对所述初始细胞进行转录组测定，得到第一基因表达谱；对所述目标细胞进行转录组测定，得到第二基因表达谱；将所述第一基因表达谱和所述第二基因表达谱进行比对，得到差异基因；
8.对所述目标细胞进行扰动，得到扰动结果，对所述扰动结果进行转录组测定，得到所述目标细胞在扰动下的基因表达谱矩阵；
9.从所述基因表达谱矩阵中提取转录因子表达谱生成第一矩阵；
10.根据所述差异基因从所述基因表达谱矩阵中提取差异基因表达谱生成第二矩阵；
11.获取所述目标细胞的转录因子结合位点信息，根据所述转录因子结合位点信息生成第三矩阵；
12.根据所述第一矩阵、所述第二矩阵和所述第三矩阵得到正则化模型，根据所述正则化模型计算得到关键转录因子预测结果。
13.进一步，所述对所述目标细胞进行扰动，得到扰动结果，对所述扰动结果进行转录组测定，得到所述目标细胞的基因表达谱矩阵这一步骤中，所述扰动包括第一扰动和第二扰动，所述第一扰动是对同一种基因的扰动过程，所述第二扰动是对多种基因的并行扰动过程。
14.进一步，所述第二扰动得到的所述扰动结果通过单细胞转录组测定得到所述目标细胞的所述基因表达谱矩阵。
15.进一步，所述第一扰动和所述第二扰动至少存在1次。
16.进一步，所述获取所述目标细胞的转录因子结合位点信息，根据所述转录因子结合位点信息生成第三矩阵这一步骤，还包括：
17.判断是否存在所述目标细胞的超级增强子区域信息，当所述超级增强子区域信息存在时，将所述目标细胞的所述转录因子结合位点信息和所述超级增强子区域信息结合生成所述第三矩阵；当所述超级增强子区域信息不存在时，根据所述目标细胞的所述转录因子结合位点信息生成所述第三矩阵。
18.进一步，当存在所述超级增强子区域信息时，通过所述超级增强子区域信息对所述转录因子结合位点信息进行过滤，从而生成所述第三矩阵。
19.进一步，所述正则化模型中的所述第三矩阵作为初始值进行计算，通过对所述第三矩阵进行优化计算，得到最优解，所述最优解即是所述关键转录因子预测结果。
20.根据本发明实施例的第二方面，一种分析细胞命运转换关键转录因子的系统，包括以下模块：
21.细胞获取模块，用于获取初始细胞和目标细胞；
22.差异基因比对模块，用于对所述初始细胞进行转录组测定，得到第一基因表达谱；对所述目标细胞进行转录组测定，得到第二基因表达谱；将所述第一基因表达谱和所述第二基因表达谱进行比对，得到差异基因；
23.基因表达谱矩阵测定模块，用于对所述目标细胞进行扰动，得到扰动结果，对所述扰动结果进行转录组测定，得到所述目标细胞在扰动下的基因表达谱矩阵；
24.第一矩阵生成模块，用于从所述基因表达谱矩阵中提取转录因子表达谱生成第一矩阵；
25.第二矩阵生成模块，用于根据所述差异基因从所述基因表达谱矩阵中提取差异基因表达谱生成第二矩阵；
26.第三矩阵生成模块，用于获取所述目标细胞的转录因子结合位点信息，根据所述转录因子结合位点信息生成第三矩阵；
27.正则化模型计算模块，用于根据所述第一矩阵、所述第二矩阵和所述第三矩阵得到正则化模型，根据所述正则化模型计算得到关键转录因子预测结果。
28.根据本发明实施例的第三方面，一种分析细胞命运转换关键转录因子的系统，包括以下装置：
29.至少一个处理器；
30.至少一个存储器，用于存储至少一个程序；
31.当所述至少一个程序被所述至少一个处理器执行，使得所述至少一个处理器实现如第一方面所述的方法。
32.本发明的有益效果是：通过与目标细胞相关的转录因子信息、差异基因信息及转录因子结合位点信息构建正则化模型，并通过正则化模型对影响细胞命运转换的关键转录因子进行预测，相较于现有技术，本发明的分析方法可以有效提高对关键转录因子的预测准确度。
附图说明
33.为了更清楚地说明本发明实施例或者现有技术中的技术方案，下面对本发明实施例或者现有技术中的相关技术方案附图作以下介绍，应当理解的是，下面介绍中的附图仅仅为了方便清晰表述本发明的技术方案中的部分实施例，对于本领域的技术人员而言，在无需付出创造性劳动的前提下，还可以根据这些附图获取到其他附图。
34.图1是本发明实施例提供的方法流程图；
35.图2是本发明实施例提供的模块连接图；
36.图3是本发明实施例提供的装置连接图。
具体实施方式
37.以下将结合实施例和附图对本发明的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本发明的目的、方案和效果。
38.为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
39.本发明的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤或单元。
40.在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本发明的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
41.首先，对本发明实施例中涉及的相关名词术语进行介绍和说明：
42.细胞命运：也称之为细胞命运决定，即细胞在发生可识别的形态变化之前，因受到约束而向特定方向分化，确定其未来发育命运的行为。
43.关键转录因子：已有大量案例表明，我们可以通过改变少量转录因子的表达，来获得某种细胞的特异性表达谱，从而改变细胞类型，这种转录因子即称为关键转录因子，例如：胚胎干细胞中的pou5f1(oct
‑
4)、sox2、nanog、klf4，肌肉细胞中的pou5f1(myod1、myog和myf5)，血细胞中的gata1、cebpa、sfpi1，t细胞中的tbx21、foxp3，肠细胞中的cdx2；关键转录因子可以在转录层次的最高层起作用并调节多个组织特异性基因。
44.参照图1，所示为根据本发明实施例提供的一种分析细胞命运转换关键转录因子的方法，包括以下步骤s100～s700：
45.s100、获取初始细胞和目标细胞；
46.s200、对所述初始细胞进行转录组测定，得到第一基因表达谱；对所述目标细胞进行转录组测定，得到第二基因表达谱；将所述第一基因表达谱和所述第二基因表达谱进行
比对，得到差异基因；
47.s300、对目标细胞进行扰动，得到扰动结果，对扰动结果进行转录组测定，得到目标细胞在扰动下的基因表达谱矩阵；需要说明的是，扰动包括第一扰动和第二扰动，第一扰动是对同一种基因的扰动过程，第二扰动是对多种基因的并行扰动过程；第二扰动得到的扰动结果需要通过单细胞转录组测定得到目标细胞的基因表达谱矩阵；第一扰动和第二扰动至少存在1次；一般地，为了满足预测需求，第一扰动或第二扰动可能进行数十次甚至上百次，因此获取的数十次或者上百次的转录组测定结果经过汇总、筛选及比对，最终形成可揭示目标细胞特异性基因调控网络的基因表达谱矩阵；
48.我们在该步骤中的扰动可以理解为，我们为了得到目标细胞特异性的基因调控信息，可对目标细胞扰动后进行转录组测序，一般来说，存在两类方式可获得所需的扰动后的转录组信息，第一类是每次实验实施单一的扰动，例如过表达、敲降或敲除某一个基因、改变培养基成分、用一种药物或者化学小分子处理细胞等，然后对扰动前后的细胞进行转录组测序，每次实验得到一组转录组数据，分别针对多个不同基因、成分、药物或小分子实施扰动，可得到多组转录组数据信息；第二类是每次实验可同时对多个目标进行扰动，比如pooled crispr screen对大量基因实施敲除或者pooled overexpression screen对大量基因过表达，然后通过单细胞转录组测序得到大量细胞扰动后的转录组；
49.s400、从所述基因表达谱矩阵中提取转录因子表达谱生成第一矩阵；
50.s500、根据所述差异基因从所述基因表达谱矩阵中提取差异基因表达谱生成第二矩阵；
51.s600、获取目标细胞的转录因子结合位点信息，根据转录因子结合位点信息生成第三矩阵；
52.可选地，步骤s600还存在以下情况，即：
53.首先判断是否存在目标细胞的超级增强子区域信息，当超级增强子区域信息存在时，将目标细胞的转录因子结合位点信息和超级增强子区域信息结合生成第三矩阵；当超级增强子区域信息不存在时，根据目标细胞的转录因子结合位点信息生成第三矩阵；
54.当存在超级增强子区域信息时，通过超级增强子区域信息对转录因子结合位点信息进行过滤，从而生成第三矩阵；
55.需要说明的是，在结合转录因子结合位点信息和超级增强子区域信息时，我们可以通过如下过程实现：
56.将转录因子绑定信息转换为矩阵x0，用作初始猜测以指导近似求解，每个基因启动子(基因转录起始位点的
‑
10kbp至+10kbp区间内)中的转录因子结合位点均来自chip
‑
seq/chip数据。对于转录因子i和基因j，如果转录因子i的结合位点存在于基因j的启动子上，则被指定为转录因子i(a，i)和基因j的表达谱的皮尔逊相关系数(pcc)(b，j)在所有样本中；否则，将分配为0。合并超级增强子信息时，将使用超级增强子区域来过滤转录因子结合位点。当转录因子结合位点在超增强器区域之外时，该转录因子结合位点定义的重置为0。
57.s700、根据第一矩阵、第二矩阵和第三矩阵得到正则化模型，根据正则化模型计算得到关键转录因子预测结果；需要说明的是，正则化模型中的第三矩阵作为初始值进行计
算，在对第三矩阵进行优化计算的过程中，直到得到满足要求的解，即是关键转录因子的预测结果。
58.对于正则化模型的构建过程表述如下：
59.一般来说，转录因子与目标细胞之间的调节关系可假设为大致由线性系统构成，我们可以构成以下计算公式，即公式(1)，
60.ax＝b+ε
ꢀꢀꢀꢀ
(1)
61.其中a代表第一矩阵，a∈r
m
×
r
，即m个样品或细胞中r个转录因子的表达矩阵；b代表第二矩阵，b∈r
m
×
n
，即m个样本或细胞中初始细胞和目标细胞之间n个差异基因的表达矩阵；ε代表噪声矩阵，ε∈r
m
×
n
；x代表根据第三矩阵x0形成的初始矩阵进行运算求解的未知矩阵，x∈r
r
×
n
，即r个转录因子和n个差异基因之间调节关系的调节矩阵；如前文，a和b是通过扰动实验得出的。
62.在识别起始细胞和目的细胞之间的差异表达基因之后，关键转录因子预测是找出调节目标涵盖大多数差异表达基因的少量转录因子，也就是说可以将其描述为一种优化问题，即找到一个x使得仅使用少量选定的关键转录因子将ax和b之间的差异最小化。
63.需要注意的是，关键转录因子的数量由x的l
2,0
范数衡量，该范数由x的非零行数(即x的组稀疏度)定义，因此，可以通过组稀疏约束优化问题来表达关键转录因子的推断过程，即公式(2)和公式(3)。
64.min‖ax
‑
b‖
f
ꢀꢀꢀꢀ
(2)
65.s.t.‖x‖
2,0
≤k
ꢀꢀꢀꢀ
(3)
66.参照如上公式(2)，f代表frobenius范式由公式(4)定义，‖x‖
2,0
定义的是组稀疏促进惩罚，k代表给定的组稀疏水平，且组稀疏度k可以根据实际需求进行人工设置，它代表我们希望作为关键转录因子的个数，并限制了‖x‖
2,0
的值，较小的‖x‖
2,0
值意味着x的组稀疏度较高，这也表明预计可以调节差异表达基因的关键转录因子数量较少。
[0067][0068]
通常，解决约束优化问题，即公式(1)是比较困难的，为避免这种障碍，一种用于求解线性逆问题的组稀疏近似解的方法是公式(1)正则化为无约束的组稀疏优化问题，其表达过程参见公式(5)。
[0069][0070]
其中λ代表正则化参数，且λ＞0，可在准确性和稀疏度之间取得权衡。l
2,0
范式可以将x的每一行(即每行代表一个转录因子与所有差异基因之间的调节关系)作为一组处理。因此，通过诱导组稀疏性，本发明涉及的方法可以选择量化关键转录因子对所有目标差异基因的调节强度的组。
[0071]
参照表1，所示为不同方法对小鼠成纤维细胞(即前文所提及的初始细胞)转换成胚胎干细胞(即前文所提及的目标细胞)所需关键转录因子的预测结果比对。
[0072]
表1小鼠成纤维细胞转换胚胎干细胞关键转录因子预测对比表
[0073]
[0074][0075]
需要说明的是，表格中的0表示该转录因子并没有被预测为前20的关键转录因子。
[0076]
在将小鼠成纤维细胞转换为胚胎干细胞的过程中，关键转录因子包括pou5f1、sox2、klf4、myc、nanog、esrrb、sall4、lin28、tetl、prdm14、jarid2、zic3以及glis，表1中特别列入设置与某类关键转录因子相关的文献数量，根据相关文献数量可从侧面佐证这种关键转录因子在小鼠成纤维细胞转换为胚胎干细胞这一转换过程中的重要程度。
[0077]
涉及的预测方法包括mogrify、cellnet、chip_superenh、hubtf、toptf、ptfbs_superenh、glasso、gso以及gso_superenh；其中，mogrify、cellnet、chip_superenh、hubtf、toptf、ptfbs_superenh以及glasso均为现有技术，gso以及gso_superenh为本发明涉及的方法；mogrify代表一种基于网络的方法；cellne代表一种根据转录组相关性分析构建的特定细胞网络分析方法；chip_superenh代表转录因子的chip
‑
seq/chip结合位点在超级增强子中的富集分析；hubtf代表小鼠胚胎干细胞转录调控网络中高度连接的核心转录因子；toptf代表小鼠胚胎干细胞转录调控网络中处于调控层级最顶端的转录因子；ptfbs_superenh代表转录因子的预测结合位点在超级增强子中的富集分析；glasso代表用group lasso整合转录组和转录因子结合位点数据进行分析；gso代表用组稀疏优化l
0,2
正则化模型整合转录因子结合位点信息进行分析；gso_superenh代表用组稀疏优化l
0，2
正则化模型整合转录因子结合位点信息和超级增强子区域信息进行分析。
[0078]
需要说明的是，现有技术中ptfbs_superenh、glasso这两种方法，因其在所有转录因子的预测中呈现的预测优先度均未呈现在前20，因此不在表1中展示；虚线框中数据表示某种转录因子在某种预测方法下的预测优先度，预测优先度的数值范围为1到20，数值越小代表预测优先度越高；表1下面两行分别代表某种转录因子在某种预测方法下的评估得分和权重得分；评估得分的数值范围为1到100，数值越大代表预测效果越好；权重得分的数值范围为1到600，数值越大代表预测效果越好。
[0079]
结合上述对表1的说明内容，在应用本发明后，对小鼠成纤维细胞转换为胚胎干细胞的关键转录因子的预测效果明显优于其他预测方法，这也表明本发明相较于现有技术可以实现更好的预测效果。
[0080]
本发明通过整合细胞特异性的多组学数据来提取高质量的细胞特异性的基因调控信息；同步进行基因调控网络构建和关键转录因子预测，且将细胞命运转换过程中需要改变的所有基因视为一个完整的目标基因组，确保了所有必需基因都被作为推断的关键转录因子的靶标；同时还整合了超级增强子的基因组区域，以缩小关键转录因子候选范围，最终实现针对关键转录因子对目标基因的调控作用的量化过程；有效地提高了对关键转录因子的预测准确度。
[0081]
参照图2，本发明还提供了一种分析细胞命运转换关键转录因子的系统，包括以下模块：
[0082]
细胞获取模块201，用于获取初始细胞和目标细胞；
[0083]
差异基因比对模块202，与细胞获取模块201连接实现交互，用于对初始细胞进行转录组测定，得到第一基因表达谱；对目标细胞进行转录组测定，得到第二基因表达谱；将第一基因表达谱和第二基因表达谱进行比对，得到差异基因；
[0084]
基因表达谱矩阵测定模块203，与细胞获取模块201连接实现交互，用于对目标细胞进行扰动，得到扰动结果，对扰动结果进行转录组测定，得到目标细胞在扰动下的基因表达谱矩阵；
[0085]
第一矩阵生成模块204，分别与差异基因比对模块202和基因表达谱矩阵测定模块203连接实现交互，用于从基因表达谱矩阵中提取转录因子表达谱生成第一矩阵；
[0086]
第二矩阵生成模块205，与基因表达谱矩阵测定模块203连接实现交互，用于根据差异基因从基因表达谱矩阵中提取差异基因表达谱生成第二矩阵；
[0087]
第三矩阵生成模块206，与细胞获取模块201连接实现交互，用于获取目标细胞的转录因子结合位点信息，根据转录因子结合位点信息生成第三矩阵；
[0088]
正则化模型计算模块207，分别与第一矩阵生成模块204、第二矩阵生成模块205和第三矩阵生成模块206连接实现交互，用于根据第一矩阵、第二矩阵和第三矩阵得到正则化模型，根据正则化模型计算得到关键转录因子预测结果
[0089]
图2所示的实施例中的内容均适用于本系统实施例中，本系统实施例所具体实现的功能与图2所示的方法实施例相同，并且达到的有益效果与图2所示的方法实施例所达到的有益效果也相同。
[0090]
参照图3，本发明还提供了一种分析细胞命运转换关键转录因子的系统，包括以下设备：
[0091]
至少一个处理器301；
[0092]
至少一个存储器302，用于存储至少一个程序；
[0093]
当至少一个程序被至少一个处理器301执行，使得至少一个处理器301实现如图1所示的方法。
[0094]
图1所示的方法实施例中的内容均适用于本系统实施例中，本系统实施例所具体实现的功能与图1所示的方法实施例相同，并且达到的有益效果与图1所示的方法实施例所达到的有益效果也相同。
[0095]
本发明实施例还公开了一种计算机程序产品或计算机程序，该计算机程序产品或计算机程序包括计算机指令，该计算机指令存储在计算机可读存介质中。计算机设备的处理器可以从计算机可读存储介质读取该计算机指令，处理器执行该计算机指令，使得该计
算机设备执行图1所示的方法。
[0096]
图1所示的方法实施例中的内容均适用于本存储介质实施例中，本存储介质实施例所具体实现的功能与图1所示的方法实施例相同，并且达到的有益效果与图1所示的方法实施例所达到的有益效果也相同。
[0097]
可以理解的是，上文中所公开方法中的全部或某些步骤、系统可以被实施为软件、固件、硬件及其适当的组合。某些物理组件或所有物理组件可以被实施为由处理器，如中央处理器、数字信号处理器或微处理器执行的软件，或者被实施为硬件，或者被实施为集成电路，如专用集成电路。这样的软件可以分布在计算机可读介质上，计算机可读介质可以包括计算机存储介质(或非暂时性介质)和通信介质(或暂时性介质)。如本领域普通技术人员公知的，术语计算机存储介质包括在用于存储信息(诸如计算机可读指令、数据结构、程序模块或其他数据)的任何方法或技术中实施的易失性和非易失性、可移除和不可移除介质。计算机存储介质包括但不限于ram、rom、eeprom、闪存或其他存储器技术、cd
‑
rom、数字多功能盘(dvd)或其他光盘存储、磁盒、磁带、磁盘存储或其他磁存储装置、或者可以用于存储期望的信息并且可以被计算机访问的任何其他的介质。此外，本领域普通技术人员公知的是，通信介质通常包含计算机可读指令、数据结构、程序模块或者诸如载波或其他传输机制之类的调制数据信号中的其他数据，并且可包括任何信息递送介质。
[0098]
上面结合附图对本发明实施例作了详细说明，但是本发明不限于上述实施例，在技术领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。