一种感知细胞中作用转录因子的方法

摘要

本发明提供的是一种感知细胞中作用转录因子的方法。设计模型计算转录因子在基因启动子上的结合值，再分别计算转录因子在基因启动子上结合值与基因表达的相关程度值；当对特定转录因子，假定甲基化抑制转录因子结合得到的相关程度值大于不考虑甲基化作用得到的相关程度值，且考虑甲基化促进转录因子结合得到的相关程度值小于不考虑甲基化作用得到的相关程度值，则为该转录因子在该细胞中存在并对基因表达进行调控。本发明可用于判定基因启动子上转录因子的结合情况，能有效解决感知细胞中作用转录因子仅能通过实验手段导致的成本高，效率低的问题。

说明书

技术领域

本发明涉及的是一种生物信息分析与处理方法，具体涉及一种对细胞中基因表达起调控作用的转录因子的感知方法。

背景技术

目前，对细胞中基因表达起调控作用转录因子的搜寻仅能通过实验手段进行，不仅成本高，而且效率低。本专利针对这一问题发明一种新的分析方法，不通过实验即可感知细胞中对基因表达起调控作用的转录因子。

众多研究成果表明，基因启动子甲基化可抑制转录因子在基因启动子区域的结合，从而导致基因表达水平的下降甚至停止表达，以达到调控基因表达的目的。本发明利用这一自然现象，通过模型来描述基因启动子甲基化对转录因子结合的作用，从而间接对细胞中作用转录因子进行感知。

发明内容

本发明的目的在于提供一种能有效解决感知细胞中作用转录因子仅能通过实验手段导致的成本高，效率低的问题的感知细胞中作用转录因子的方法。

本发明的目的是这样实现的：

包括转录因子在基因启动子上匹配值计算，基因启动子甲基化对转录因子结合的作用描述，转录因子在基因启动子上结合值计算，转录因子在基因启动子上结合值与基因表达相关程度值计算，作用转录因子分析步骤；分别在假定基因启动子甲基化不影响转录因子结合、甲基化抑制转录因子结合和甲基化促进转录因子结合的基础上，设计模型计算转录因子在基因启动子上的结合值，再分别计算转录因子在基因启动子上结合值与基因表达的相关程度值；当对特定转录因子，假定甲基化抑制转录因子结合得到的相关程度值大于不考虑甲基化作用得到的相关程度值，且考虑甲基化促进转录因子结合得到的相关程度值小于不考虑甲基化作用得到的相关程度值，则为该转录因子在该细胞中存在并对基因表达进行调控。

所述模型采用E＝e^-((M-C)/S)/(1+e^-((M-C)/S))函数形式描述基因启动子甲基化对转录因子结合抑制作用的模型；其中，M是转录因子结合位点的甲基化水平，E是该结合位点上甲基化对转录因子的抑制作用、其取值范围为[0，1]；当结合位点的甲基化水平M非常小时，E将趋近1，表示甲基化对转录因子的结合无影响；反之，当M足够大时，E将趋近0，表示在甲基化影响下，转录因子无法在该位点结合；另外，模型有两个参数，分别是函数中心值C和函数陡度S。

针对不同的转录因子，该模型参数是可以变化的。

本发明属于生物信息分析与处理领域。本发明的方法可用于判定基因启动子上转录因子的结合情况，以进一步研究转录因子的调控基因以及基因表达的调控机制。

本发明的主要贡献和特点在于：设计了描述基因启动子甲基化对转录因子结合抑制作用的模型，并在此基础上，提出了一个新方法实现了对细胞中基因表达起调控作用的转录因子的感知。从而，有效解决了目前感知细胞中作用转录因子仅能通过实验手段导致的成本高，效率低的问题。

附图说明

图1是本发明的流程图；

图2是甲基化抑制转录因子结合的通用模型曲线(C＝0，S＝0.5)；

图3是甲基化促进转录因子结合的通用模型曲线(C＝0，S＝0.5)；

图4(a)-(b)是SH-SY5Y神经瘤细胞全部基因的Z_m值分布，其中图4(a)不考虑甲基化作用，图4(b)考虑甲基化作用；

图5是SH-SY5Y神经瘤中不考虑甲基化作用下的Z值分析结果；

图6是SH-SY5Y神经瘤中考虑甲基化抑制作用下的Z_m值分析结果(S＝0.05)；

图7是SH-SY5Y神经瘤中考虑甲基化促进作用下的Z_m值分析结果。

具体实施方式

下面对本发明做更详细地描述：

1.细胞中作用转录因子的感知方法

本发明方法包括如下步骤。

1.1转录因子在基因启动子上的匹配值计算

基因启动子是一个对于基因表达调控至关重要的区域。各种转录因子在基因启动子上结合并对基因表达进行调控。为了判定转录因子是否在基因启动子上结合，最直接的方法是统计分析转录因子易于结合的碱基组合，再通过在基因启动子上进行碱基组合匹配分析来得到转录因子的结合位点。

利用转录因子与碱基结合的位置权值矩阵，可以对转录因子在基因启动子上的碱基匹配情况进行计算。该矩阵是通过统计大量实验测得的转录因子实际结合位点的碱基组合数据生成的。

本方法中，基因启动子的长度设为基因转录起始位点的前1200个碱基至后300个碱基，共1500个碱基。这样，设第i个转录因子结合位点长度为L，则在第j个基因启动子上共有1500-L+1个假定转录因子结合位点。在第k个假定结合位点上，可以根据该位点的碱基组合与转录因子位置权值矩阵的匹配程度计算出转录因子在该位点的匹配权值和A_ijk。

$A_{\mathrm{ijk}}=\underset{l=1}{\overset{L}{\Sigma }}{a}_{\mathrm{jkl}}{w}_{\mathrm{il}}^T---\left(1\right)$

式(1)中，a_jk1为描述第j个基因启动子第k个假定结合位点第1个位置上碱基种类的行向量(1000表示碱基A，0100表示碱基C，0010表示碱基G，0001表示碱基T)，w_i1为第i个转录因子位置权值矩阵中第1行与碱基相匹配的权值行向量。这样，对第j个基因，我们可以得到第i个转录因子在其启动子上所有1500-L+1个假定结合位点上的匹配权值和。那么，我们将最大的匹配权值和定义为第i个转录因子在第j个基因启动子上的匹配值S_ij。

S_ij＝max A_ijk    (2)

1.2基因启动子甲基化对转录因子结合的作用模型

在分析转录因子是否在基因启动子上结合并进行调控时，除了考虑启动子上碱基组合与转录因子位置权值矩阵的匹配程度，还应考虑启动子上甲基化对转录因子结合的作用。

这里设计了2个通用模型来描述基因启动子甲基化对转录因子结合的作用。

模型1采用如下函数形式：

$E_{\mathrm{ijk}}=\frac{e^{-\left((M_{\mathrm{jk}}-C_i)S_i\right)}}{1+e^{-((M_{\mathrm{jk}}-C_i)/S_i)}}---\left(3\right)$

其中，M_jk是第j个基因启动子上第k个假定结合位点的甲基化水平，E_ijk是其对第i个转录因子结合的抑制作用，其取值范围为[0，1]。当第k个假定结合位点的甲基化水平M_jk非常小时，E_ijk将趋近1，表示甲基化对转录因子的结合无影响；反之，当M_jk足够大时，E_ijk将趋近0，表示在甲基化影响下，转录因子无法在该位点结合。另外，模型有两个参数，分别是函数中心值C和函数陡度S。该模型曲线如图2所示。由于针对不同的转录因子，模型参数是可以变化的，因此，式(3)所描述的模型中，该两个参数分别是针对第i个转录因子的C_i和S_i。

图2中，横轴是甲基化水平，纵轴则表示甲基化对转录因子结合的作用。针对不同的转录因子，通过取不同模型参数，通用模型曲线是不同的，图2中模型参数假选为C＝0，S＝0.1。

用该通用模型来描述基因启动子甲基化对转录因子结合的抑制作用是基于以下两点生物特性：首先是结合位点的甲基化水平越高，其对转录因子结合的抑制作用越强。其次是结合位点的甲基化水平对转录因子结合的抑制作用并不随甲基化水平的变化而成线性关系。当甲基化水平较高或较低时，其对转录因子结合的抑制作用将出现饱和现象，对甲基化水平的变化也将不敏感；而在中间值时，其抑制作用对甲基化水平的变化将较敏感。图2所示的通用模型曲线较好地反映了上述两点生物特性。

为了描述特定细胞环境下基因启动子甲基化对不同转录因子的抑制作用，该通用模型设有两个可变参数，即函数中心值C和函数陡度S。这样，特定细胞环境下，通过选取不同的模型参数，每个转录因子将具有不同的模型来描述其受基因启动子甲基化的影响。

与模型1相对应，模型2采用如下函数形式：

$E_{\mathrm{ijk}}=\frac{1}{1+e^{-((M_{\mathrm{jk}}-C_i)/S_i)}}---\left(4\right)$

该模型曲线如图3所示。与模型1所描述的甲基化对转录因子结合起抑制作用相反，模型2所描述的是甲基化对转录因子结合起促进作用。

1.4转录因子在基因启动子上的结合值计算

考虑基因启动子甲基化对转录因子结合的作用后，可以按式(5)计算第i个转录因子在第j个基因启动子上的结合值：

S_ij＝max(A_ijk×E_ijk)        (5)

其中，E_ijk表示第j个基因启动子上第k个假定结合位点的甲基化水平对第i个转录因子结合的作用。

这样，当不考虑基因启动子甲基化对转录因子结合的作用(相当于E_ijk＝1)，假定基因启动子对转录因子结合起抑制作用，假定基因启动子对转录因子结合起促进作用，可以计算得到第i个转录因子在第j个基因启动子上的3个结合值。

1.5转录因子在基因启动子上的结合值与基因表达相关程度值计算

计算得到转录因子在基因启动子上的结合值后，并不意味着该转录因子一定在其启动子上结合并对该基因表达进行调控。只有结合值足够大，才能认为转录因子在基因启动子上结合。

为了分析转录因子是否在基因启动子上结合并对该基因表达进行调控，需进行结合值的阈值分析，即当转录因子在基因启动子上的结合值大于某一阈值时，才可认为该转录因子在该基因启动子上结合并对该基因表达进行调控。

转录因子结合阈值的确定方法有很多。这里采用一种根据基因表达和转录因子结合值相关性分析的阈值确定方法。设有n个基因，g₁，...，g_n是这些基因表达值的对数且遵循正态分布，其均值为μ，均方差为σ。设定一个阈值，根据转录因子在基因启动子上结合值是否大于该阈值得到一个基因子集Gⁱ。设该基因子集Gⁱ有k个基因，其序号为i₁，...i_k。则可如式(6)计算一个Z值：

$Z(\mathrm{TF},G^i)=(\frac{1}{k}\underset{j=1}{\overset{k}{\Sigma }}{g}_{\mathrm{ij}}-\mu )/(\sigma /\sqrt{k})---\left(6\right)$

由式(6)分析可得，Z值反映着基因子集Gⁱ中基因表达平均值与全部基因表达平均值的差异程度。Z值的绝对值越大，差异程度越大，表明基因子集Gⁱ中的基因越可能被该转录因子结合并调控。而Z值的符号则反映着该转录因子的生物功能是促进基因表达还是抑制基因表达。这样，通过调整阈值，我们可得到一系列不同的Z值。其中，当Z值的绝对值最大时，其对应的阈值即为最大似然阈值，该Z值称为Z_m。可以看出，上述过程是在分析基因表达和转录因子在基因启动子上结合值的相关性。相关性越大，无疑Z_m的绝对值也越大。

此外，在计算转录因子在基因启动子结合值时，需要对该转录因子受甲基化作用模型中的两个参数进行确定。这可以通过变动模型参数值C和S，在转录因子结合阈值分析中得到一系列不同的Z_m值。当Z_m的绝对值最大时，此时的模型参数值C和S即为最大似然模型参数值。而由最大似然模型参数值C和S构造的模型即可用于计算最终的基因表达和基因启动子上的结合值与基因表达相关程度值。

这样，当不考虑基因启动子甲基化对转录因子结合的作用，假定基因启动子对转录因子结合起抑制作用，假定基因启动子对转录因子结合起促进作用，对一个特定转录因子，可以计算得到相应的3个基因表达和基因启动子上的结合值的相关程度值。

1.6作用转录因子分析

由于细胞中基因启动子甲基化对转录因子结合能力具有抑制作用已经是公认的一种自然现象，因此，当对特定转录因子，考虑基因启动子甲基化抑制作用得到的基因表达和基因启动子上结合值的相关程度值大于不考虑甲基化作用得到的相关程度值，且考虑甲基化促进作用得到的相关程度值小于不考虑甲基化作用得到的相关程度值，则可认为该转录因子在该细胞中存在并对基因表达进行调控。

2.方法测试

为了验证本发明方法的有效性，通过具体生物数据对本方法进行了测试。

2.1数据

本发明方法测试选用的生物实验数据来自国际公开的大型生物信息数据库GEO(Gene Expression Omnibus)，分别是针对SH-SY5Y神经瘤细胞所作的生物实验数据(GSE4600，GSE9568)。该检测数据包括如下部分：

2.1.1基因表达数据

本发明测试所需的SH-SY5Y神经瘤细胞基因表达数据从美国国家生物信息中心管理的GEO数据库获得(GSE4600)，该组数据由Affymetrix公司的HG U133 plus 2.0基因芯片生成。

2.1.2基因启动子甲基化数据

本发明测试所需的SH-SY5Y神经瘤细胞基因启动子甲基化数据也从美国国家生物信息中心管理的GEO数据库获得(GSE9568)，该组数据由NimbleGen公司的1500b_Pomoter_Array基因芯片生成。

此外，我们还从Transfac7.0数据库获得459个人类转录因子的位置权值矩阵。而人体基因启动子序列信息和人体基因描述信息则从美国UCSC网站上获取。

根据人体基因描述信息，从上述SH-SY5Y神经瘤细胞生物实验数据中共整理出10065个基因的表达数据和启动子甲基化水平数据。

2.2基因启动子甲基化对转录因子结合的作用

首先，不考虑基因启动子甲基化对转录因子结合的作用，通过转录因子在基因启动子的结合值和基因表达相关性计算Z_m值。其分布直方图如图4(a)所示。图中，横轴为Z_m值，纵轴为频次。

与此对应，考虑基因启动子甲基化对转录因子结合的抑制作用，通过转录因子在基因启动子的结合值和基因表达相关性计算Z_m值。其分布直方图如图4(b)所示。

由图4(a)，图4(b)可见，不考虑甲基化作用时，正Z_m集合和负Z_m集合的均值分别为5.33和-3.24；而考虑甲基化抑制作用时，正Z_m集合和负Z_m集合的均值分别为7.23和-4.66。可见，考虑基因启动子甲基化抑制作用后，Z_m值的绝对值明显变大，这表明转录因子结合值和基因表达相关性程度明显提高，证明基因启动子甲基化确实对转录因子的结合起抑制作用。

由于转录因子E2F1在SH-SY5Y神经瘤细胞中丰度较高，因此针对E2F1对本发明方法进行了详细分析。

首先，不考虑基因启动子甲基化对转录因子E2F1结合的作用，对E2F1在基因启动子的结合值进行阈值分析。不同阈值下得到的Z值如图5所示。由图5可见，Z值为正值，这表明E2F1对基因表达起促进作用，这是符合转录因子E2F1的生物功能的。当匹配阈值取9.54时，Z值取最大，即Z_m＝11.17。

图5中，横轴为转录因子结合阈值，纵轴为对应的Z值。

其次，考虑基因启动子甲基化对转录因子E2F1的结合起抑制作用，通过转录因子在基因启动子的结合值和基因表达的相关性寻找最大的Z_m值。即变动通用模型参数C和S寻找最大的Z_m值。通过搜寻，当模型参数C＝-0.15，S＝0.05时得到最大的Z_m值(Z_m＝12.531)。图6显示了当S＝0.05时，随着C值的变化，Z_m值的变化过程。图6中的横线(Z_m＝11.17)标明不考虑基因启动子甲基化影响时计算得到的Z_m值。

图6中，横轴为模型参数C值，纵轴为对应的Z_m值。

由图6可见，当模型参数C值比较小时，Z_m值非常小，表明模型参数取值不合理。随着C值的增大，Z_m值迅速增大。当C值大于-0.35时，Z_m值开始大于11.17，表明其分析结果已优于不考虑基因启动子甲基化对转录因子E2F1作用时得到的结果。Z_m值在-0.15取最大值。随后，Z_m值随C值增大而下降，并趋近于11.17。通过上述分析，证明考虑甲基化抑制作用来分析转录因子E2F1在基因启动子上的结合是合理的。

最后，考虑基因启动子甲基化对转录因子E2F1的结合起促进作用，通过转录因子在基因启动子的结合值和基因表达的相关性寻找最大的Z_m值，其结果如图7所示。

由图7可见，不论函数参数C取何值，Z_m值均远小于11.17。这表明基因启动子甲基化对转录因子E2F1的结合起促进作用这一假设是不合理的。

可见，对于在SH-SY5Y神经瘤细胞中存在的转录因子E2F1，考虑基因启动子甲基化抑制作用得到的基因表达和基因启动子上结合值的相关程度值大于不考虑甲基化作用得到的相关程度值，且考虑甲基化促进作用得到的相关程度值小于不考虑甲基化作用得到的相关程度值。

采用上述分析方法，在神经瘤细胞459个人类转录因子中，共发现91个转录因子存在该现象。其中，最明显的前10个转录因子的分析结果见表1。表1中，Z_no表示不考虑甲基化作用得到的Z_m值，Z_me表示考虑甲基化作用得到的绝对值最大的Z_m值，C，S则为针对每一个转录因子所建立的启动子甲基化作用模型参数。

表1转录因子甲基化作用分析结果

我们发现在Katarzyna采用SH-SY5Y细胞研究所作实验中，通过实验测出SH-SY5Y细胞中存在17个转录因子。该研究成果发表于PLoS ONE期刊2009年4月名为“Regulation of Clock-Controlled Genes in Mammals”的研究论文中。与之相对照，本方法发现的SH-SY5Y细胞中明显受甲基化影响的前10个转录因子中，除Whn和KROX外，有8个转录因子出现在Katarzyna的研究结果中。而Whn(亦称Foxn1)和KROX也分别在吴云舟(应用新城疫病毒治疗肿瘤的研究进展，生物工程学报，2010，26(8)：1031-1036)和Merrill(All-trans retinoic acid-responsive genes identified in the human SH-SY5Y neuroblastoma cell line and their regulated expression in the nervous system of early embryos.Biological Chemistry，2004，385(7)：605 614)采用SH-SY5Y细胞所作的研究成果中提及。这直接证明了本发明方法的有效性。