一种基于深度学习的转录因子结合位点定位方法

1.本发明属于深度学习技术领域，尤其涉及一种基于深度学习的转录因子结合位点定位方法。


背景技术：

2.蛋白质在人体内无处不在，是基因经转录和翻译等细胞活动而合成的一类由20种氨基酸组成的高分子化合物。其中有一类蛋白质会和染色体进行特异性结合——dna结合蛋白(dna binding protein,dbp)，dbp在基因的复制、重组、链裂解、转录等过程当中起着关键的作用，与细胞周期中染色质的一系列变化息息相关；转录因子(transcription factors,tfs)属于dbp中的一种，也称为反式作用因子，能够与调控区域的dna序列的非编码区发生特异性相互作用，并对基因的转录有激活或者抑制作用，形成指导基因组表达的复杂系统。理解tfs与dna的结合偏好是掌握生物体转录调控机制和构建转录调控网络的关键所在。研究tfs具有重要的作用，而其主要通过和dna序列进行特异性结合来展示作用，那么理解它们之间的结合特异性就有着广泛的意义，同时对其他转录相关的研究还有着基础性的作用。发生结合的区域称为转录因子结合位点(transcription factor binding sites,tfbss)，是dna小片段，长度为5～20个碱基，tfbss的研究成为了如今的热门研究领域，算法模型从传统的测序比对算法模型，到后来人工设计特征的机器学习算法，再到神经网络模型的发展，数据编码从one-hot编码到基于k-mer的序列编码到embedding编码，该领域一直保持着与时俱进的态势。
3.进入21世纪，高通量测序技术取得突破，能够在短时间内并行地对生物数据进行测序，产生了海量的核酸和蛋白质数据。蛋白质结合微阵列(protein binding microarrays,pbms)，为研究体外蛋白质和dna的结合偏好也就是预测tfbss提供了大量的体外结合数据。尽管有着大量的数据，但是tfbss的特性和测序技术实验本身会引入大量噪声，使得预测体外蛋白质和dna结合位点的效果依然不佳，另外传统的方法无法处理大量的数据，这又为如何利用那些大量的数据增加了许多困难。机器学习(特别是深度学习)算法的出现和计算机性能的提升为此提供了解决方案。
4.为了更加精确地预测tfbss，大量基于pbms的数据分析方法被提出，在这些方法中，直接从pbms预测出tfbss是一种简单且有效的方法，因此备受研究者的青睐。根据pbms得到的位置权重矩阵(position weight matrix,pwm)是后续建模中非常重要的一类数据，pwm中的元素表示其位置上出现某种碱基的概率大小，是一个4行l列的矩阵(其中，4对应着4种核苷酸，l表示序列长度)，这样pwm的每行对应着一种核苷酸，每列对应着tfbss中的一个位置，pwm[i,j]表示dna序列的第j个位置上出现第i种核苷酸的概率，能够表示蛋白质和dna的结合偏好。pwm数据简单易懂，因此大量基于pwm的tfbss识别方法被提出来，它们直接从pwm中学习蛋白质和dna的结合偏好。但是这样的模型存在两个不足：第一，模型中tfbss的长度是固定的；第二，模型假设tfbss各个位置之间是相互独立的。
[0005]
为了解决上面的问题，基于k-mer编码的tfbss识别方法被提出来，这类方法可以
编码出核苷酸之间的依赖关系；但是基于k-mer的方法中，基因序列仅由k-mer计数的向量表示，并没有考虑序列中每一个分段的位置，此外尽管位置特异性的序列核是存在的，但是它将序列映射到更高维度的空间，使得效率不高。
[0006]
综上，对于tfbss的定位问题，目前现有做法是先利用识别算法过滤出可能的序列，然后概率统计的方法来确定结合区域，需要遍历整条序列挑选出概率最大的位置，数据量不大时还可以应付，但是随着数据量的增长就呈现出性能差、响应时间过长的缺陷。


技术实现要素：

[0007]
为解决上述技术问题，本发明提出了一种基于深度学习的转录因子结合位点定位方法，将其转换成二分类的问题，也就是将属于结合位点区域和不属于结合位点区域看作是两种类别，使用全卷积网络的编码器-解码器结构来学习数据模式，构建模型参数，使用全局池化层提取全局语义信息，细化残差块结构来丰富跳跃结构的空间位置信息，外加保守信息作为辅助输入数据，最后得到一个能定位出tfbss的算法模型fcnarrb 。
[0008]
为实现上述目的，本发明提供了一种基于深度学习的转录因子结合位点定位方法，包括：
[0009]
将与转录因子结合的dna序列进行one-hot编码，获取数据集，基于k折交叉验证的方法将所述数据集划分训练集和测试集；
[0010]
基于全卷积网络构建fcnarrb 模型，设置损失函数与评价指标；
[0011]
基于所述训练集与所述损失函数对所述fcnarrb 模型进行训练，训练后的所述fcnarrb 模型用于定位转录因子结合位点，并通过所述测试集与所述评价指标对训练后的所述fcnarrb 模型的定位结果进行测试评价。
[0012]
可选地，将与转录因子结合的dna序列进行one-hot编码的方式为：对所述dna序列中的碱基{a,c,g,t}按照one-hot编码进行数据编码，所述dna序列中的保守信息数据选取对应位置上的数据，两者构成所述dna序列的编码信息。
[0013]
可选地，基于全卷积网络构建所述fcnarrb 模型包括：
[0014]
将dna序列中的保守信息数据作为输入数据加入fcnarrb模型中构成所述fcnarrb 模型；
[0015]
所述fcnarrb模型由fcna模型与rrb构成，所述rrb用来提取空间信息到解码器中；
[0016]
所述fcna模型为在所述全卷积网络构的编码器最后加上一个全局池化层获得。
[0017]
可选地，所述全卷积网络结构包括编码器与解码器；
[0018]
所述编码器包括包含三个卷积块，由一个卷积层、relu层、一个最大池化层与dropout层；
[0019]
所述解码器包括三个反卷积块和四个跳跃结构。
[0020]
可选地，所述fcnarrb模型包括：第一卷积块、第二卷积块、第三卷积块、全局池化层、第一反卷积块、第二反卷积块、第三反卷积块、第四反卷积块、四个跳跃连接以及sigmoid层。
[0021]
可选地，在所述fcnarrb模型中，所述第一卷积块与第二卷积块的参数设置均为：卷积层卷积核数量64，大小15，步长1，填充0；池化层卷积核大小4，步长4，填充0；
[0022]
所述第三卷积块的参数设置为：卷积层卷积核数量64，大小3，步长1，填充0；池化
层卷积核大小2，步长2，填充0；
[0023]
所述全局池化层的参数设置为：输出数据尺寸为1；
[0024]
所述第一反卷积块、第二反卷积块、第三反卷积块与第四反卷积块中的上采样层全部采用线性扩展；
[0025]
所述第一反卷积块与第二反卷积块的参数设置均为：上采样层采用线性扩展，卷积层卷积核数量64，大小3，步长1，填充1；
[0026]
所述第三反卷积块的参数设置为：卷积层卷积核数量4，大小3，步长1，填充1；
[0027]
所述第四反卷积块的参数设置为：卷积层卷积核数量1，大小3，步长1，填充1；
[0028]
所述sigmoid层的参数设置为：将encoder的值转换成概率。
[0029]
可选地，对所述fcnarrb模型进行训练还包括：
[0030]
在所述fcnarrb模型训练过程中，选择使用adam优化算法来更新参数，同时对于dropout层的丢弃率参数、l2正则化权重衰减参数和adam算法当中使用的动量参数，采用网格搜索的策略来寻找最优组合参数，所述dropout层为全卷积网络构中。
[0031]
可选地，所述损失函数的表达式如下：
[0032][0033][0034]
l＝l
p
ln[0035]
其中，p为正例样本，k为选择的反例样本中计入损失函数中的反例个数，n为反例样本集合，pi为正例样本i被预测为正例的概率，n
p
为正例样本个数，li和lj分别为遍历正例样本集合和反例样本集合的交叉熵损失函数，l
p
为正例样本的损失函数，ln为反例样本的损失函数，topk为前k个损失函数值的反例样本集合。
[0036]
可选地，所述评价指标采用iou与prauc；
[0037]
所述iou用于评价模型预测的dense label和原始label的匹配程度，所述prauc用于综合评价模型的准确率和召回率。
[0038]
与现有技术相比，本发明具有如下优点和技术效果：
[0039]
本发明将与转录因子结合的dna序列进行one-hot编码，将tfbss的定位问题转换为碱基级别上的二分类问题，除了对dna序列进行编码外，还为其生成一个dense label数组，数组上属于tfbss区域的位置值为1，其余位置为0，然后使用一个编码器-解码器结构来学习，构建模型参数。使用fcn可以解决可变长度的序列，而且还使用了保守信息作为dna序列的附加输入数据，为了学习全局语义特征在编码器的最后附加上一个全局池化层，为了获取足够的空间位置信息在编码器和解码器之间的跳跃连接加上细化残差块结构来传递信息，构成fcnarrb模型，将在fcnarrb模型中加入dna序列中的保守信息数据作为输入数据之一的，记为fcnarrb 模型，通过本发明所提出的fcnarrb 模型能够实现转录因子结合位点的精准预测定位，具有性能好、响应时间快的优点。
附图说明
[0040]
构成本技术的一部分的附图用来提供对本技术的进一步理解，本技术的示意性实
施例及其说明用于解释本技术，并不构成对本技术的不当限定。在附图中：
[0041]
图1为本发明实施例的一种基于深度学习的转录因子结合位点定位方法流程示意图；
[0042]
图2为本发明实施例的全卷积网络结构fcn示意图；
[0043]
图3为本发明实施例的fcna模型结构示意图；
[0044]
图4为本发明实施例的fcnarrb模型结构示意图；
[0045]
图5为本发明实施例的模型fcn、fcna、fcnarrb和fcnarrb 的iou(左)prauc(右)示意图；
[0046]
图6为本发明实施例的模型fcn、fcna、fcnarrb和fcnarrb 的iou(左)和prauc(右)平均值示意图。
具体实施方式
[0047]
需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本技术。
[0048]
需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。
[0049]
实施例
[0050]
如图1所示，本实施例提供了一种基于深度学习的转录因子结合位点定位方法，包括：
[0051]
将与转录因子结合的dna序列进行one-hot编码，获取数据集，基于k折交叉验证的方法将所述数据集划分训练集和测试集；
[0052]
进一步的，将与转录因子结合的dna序列进行one-hot编码的方式为：对所述dna序列中的碱基{a,c,g,t}按照one-hot编码进行数据编码，所述dna序列中的保守信息数据选取对应位置上的数据，两者构成所述dna序列的编码信息。
[0053]
在本实施例中，将给定的dna序列(假设长度为l)进行one-hot编码作为输入数据，同时还会根据tfbss的位置得到dense label标签(长度为l，tfbss所在的位置元素值为1，其余为0)，将tfbss的定位问题转变成每个位置的二分类问题，最后得到的序列上1所在的位置就表示其可能属于tfbss，根据连续1的情况来划定tfbss的范围，也就是完成了tfbss的定位问题。
[0054]
从encode(encyclopedia of dna elements)项目选择了gm12878细胞系的21个转录因子，并且从hocomoco数据库中下载了高质量的pcms(position count matrices)，转录因子对应的pcm表示tfbs的各个碱基出现的频率。
[0055]
对于每一个转录因子数据集，选择以顶峰位置为中心附近的500个bp为一条序列，对序列中的碱基{a,c,g,t}按照one-hot编码进行数据编码，保守信息数据就选取对应位置上的数值即可，这两者构成了序列的编码信息。至于dense label信息，按照如下表1确定序列对应的标签信息的算法来确定标签信息。
[0056]
表1
[0057][0058]
gm12878细胞系上选择的21个tfs的正例和反例样本数量如下表2所示。为了准确地评估提出模型的性能，采用了k折交叉验证的策略，也就是将数据集划分成k等份，其中k-1份作为训练数据集，另外的一份作为测试数据集验证模型效果，然后重复k次。
[0059]
表2
[0060][0061][0062]
基于全卷积网络构建fcnarrb 模型，设置损失函数与评价指标；
[0063]
在本实施例中，由于全卷积网络结构fcn能从碱基级别来识别tfbss，这样就可以确定tfbss在给定序列中的位置信息。但是全卷积网络结构fcn提取的特征在整体上可能无法获取全局范围的语义信息，所以在卷积层提取的特征图后再加上一个全局池化层的结构来从整体上获取语义信息，这样就得到了带全局池化层的全卷积网络fcna。加上全局池化层之后，网络模型能够很好的得到全局的语义信息，这对于分类任务是很有帮助的，但是有时会缺少空间信息，使得有时候难以定位到tfbs的位置，细化残差块rrb的出现为此提出了一个方向，可以在编码器到解码器的跳跃连接中加入rrb块来获取空间信息，由此提出了带有rrb模块的fcnarrb模型(fcna with rrb)；在fcnarrb模型加入dna序列中的保守信息数据作为输入数据之一即获得了fcnarrb 模型。下面对模型展开介绍。
[0064]
fcn模型
[0065]
全卷积网络结构fcn如图2所示，包含左边的编码器和右边的解码器两个部分。自顶向下的编码器部分包含三个卷积块，由一个卷积层、relu层、一个最大池化层和dropout层组成。卷积块的计算公式如下：
[0066]
x＝conv(i,w)
[0067]
y＝relu(0,x)
[0068]
z＝maxpool(y)
[0069]
o＝dropout(z)
[0070]
式中，i是输入数据，w是卷积核参数，conv表示卷积运算，relu是激活函数，maxpool表示最大池化操作，dropout表示丢弃操作，就是让部分神经元不向后传播数据，o是经过卷积块的输出数据。
[0071]
卷积层用于扫描序列，能够给序列的片段进行打分。最大池化层用来减少输入数据的计算复杂度，并选出局部相邻区域内最好的代表点。dropout层是经常用来防止过拟合的技巧。
[0072]
自下向上的解码器包含三个反卷积块和四个跳跃结构，每个反卷积块由一个上采样层、一个归一层、一个relu层和一个卷积层组成，跳跃结构是直接相加运算。反卷积模块的计算公式如下：
[0073]
自下向上的解码器包含三个反卷积块和四个跳跃结构，每个反卷积块由一个上采样层、一个归一层、一个relu层和一个卷积层组成，跳跃结构是直接相加运算。反卷积模块的计算公式如下：
[0074]y
＝upsample(z

)
[0075]y
＝y

z
[0076]y
＝bn(y

)
[0077]o
＝conv(x

,w

)
[0078]
式中，z

是输入数据，上采样操作需要记录编码器的数据尺寸信息，bn操作用来进行归一化，将不同范围的数据归一化到0～1之间，relu是激活函数，w

是卷积核参数，conv表示卷积运算，o

是经过卷积块的输出数据。
[0079]
fcna模型
[0080]
全局池化层能够捕捉全局的上下文信息，也就是语义信息，对于dna序列来说，全局的语义信息(也就是侧翼区域信息)会影响转录因子和dna序列的结合行为。这样在最后一个卷积块后面加上一个全局池化层来捕捉转录因子结合位点的全局信息。
[0081]
对比于fcn模型，fcna模型在编码器最后加上了一个全局池化层，这样就相应地就会多一个反卷积块和一个跳跃连接。fcna的具体结构如图3所示。
[0082]
fcnarrb模型
[0083]
fcn模型和fcna模型采用的跳跃结构是直接将编码器对应层数据加到解码器的对应层当中，这样的操作对于空间信息的提取就显得有些粗糙了。细化残差块rrb可以用来提取空间信息到解码器中，所以就在跳跃结构当中增加了rrb先提取编码器的特征再加到解码器对应的层中。模型结构如图4所示，rrb模块能够过滤数据中存在的噪声，其计算公式可以简单的如下所示：
[0084]
x＝x f(x)
[0085]
式中，x是输入数据，f(x)表示经过卷积运算得到的结果。
[0086]
在fcnarrb模型中，为了检验保守信息的作用，可以将保守信息和序列的one-hot编码结合在一起进行消融试验来检验保守信息对模型效果的影响。由此，将在fcnarrb模型中加入dna序列中的保守信息数据作为输入数据之一的；记为fcnarrb 模型；
[0087]
其输入数据是(b,501,5)，其中b表示一次性处理的数据数量，501表示序列长度，5表示输入数据的通道数，其中4个通道是dna序列进行one-hot编码，另外1个是保守信息，fcnarrb模型中各层的具体参数设置及输出数据尺寸见下表3。
[0088]
表3
[0089][0090]
这样给定一条长度为501的序列，按照表3的算法得到这条序列的dense label，它标识了tfbss的位置信息。对序列进行one-hot编码输入到自顶向下的编码器和自下向上的解码器进行学习，最后输出一条长度为501的数字序列，其中的每个数字都经过sigmoid函数处理成当前位置属于tfbss的概率，如何设置了阈值，就可以得到一条0和1构成的长度为501的序列。最后通过对比这条输出序列和dense label序列就可以使用指标来评估模型的效果。
[0091]
设置损失函数与评价指标；
[0092]
由于选择的序列长度为501bp，而转录因子结合位点的长度一般为5～20bp，所以denselabel中标签为1和为0的比例就在1:95～1:24，显然这是分布极不均衡的情况，这和目标检测领域hard negative mining(难分样本挖掘)问题很像，为此提出hard sampling方法，从所有样本中选择一定量的正样本和负样本，只有被选择的样本才计算损失，一般会
倾向选择一些难负例样本，选择hard sampling方法来处理本发明中定位tfbss中的0-1分布不均衡问题，损失函数计算公式如下所示：
[0093][0094][0095]
l＝l
p
ln[0096]
其中，p为正例样本，k为选择的反例样本中计入损失函数中的反例个数，n为反例样本集合，pi为正例样本i被预测为正例的概率，n
p
为正例样本个数，li和lj分别为遍历正例样本集合和反例样本集合的交叉熵损失函数，l
p
为正例样本的损失函数，ln为反例样本的损失函数，topk为前k个损失函数值的反例样本集合。
[0097]
hard sampling损失函数计算算法如下表4所示：
[0098]
表4
[0099][0100][0101]
评价指标：
[0102]
要定位tfbss，这涉及到两个方面：(1)被预测为tfbss的位置要尽可能的和原始序列中的区域对齐，交合比(intersection overunion,iou)在目标检测领域能够反映预测检测框和真实检测框的检测效果，具有尺度不变性，将交合比引入到tfbss的定位问题当中，用来评价tfbss预测的区间和真实区间的检测效果；(2)对于序列中0和1的预测要尽可能的准确，引入了prauc指标，prauc表示pr曲线、x＝1直线以及横轴之间的面积，0～1之间，值越大表示模型效果越好。
[0103]
iou应用到tfbss的定位问题中，关注1的位置，同时也关注0的位置，iou的计算公式如下：
[0104][0105]
假如真实序列是{0,0,0,1,1,1,0}，模型预测的序列是{0,0,1,1,1,0,0}，那么分别计算niou(negative iou)和piou(positive iou)。对于niou，两条序列中至少有一个0的位置有{1,2,3,5,6}，而两个都为0的位置有{1,2,6}，所以niou是0.6；对于piou，两条序列中至少有一个1的位置有{3,4,5,6}，而两个都为1的位置有{4,5}，所以piou是0.5；最后iou就取niou和piou的平均值，也就是0.55。
[0106]
prauc是pr曲线下的面积，p指精确度(precision rate)，表示所有预测为正例的
样本中真实类别为正例的比例，r指召回率(recall rate)，表示所有真实类别为正例的样本中被预测为正例的比例，计算公式如下：
[0107][0108][0109]
式中，tp、fp和fn分别表示真正例(true positive)、假正例(false positive)和假反例(false negative)，具体如下表5的混淆矩阵所示。
[0110]
表5
[0111][0112]
pr曲线的绘制介绍如下：首先，按照预测的概率从大到小排序；其次，不断减小阈值，将样本预测为正例，依次计算准确率p和召回率r；最后，将第二步得到的一系列点标识在xoy坐标系上，连接成线，即为pr曲线，曲线下面积就是prauc。auc整体上均衡的关注正例样本和反例样本，prauc则偏向于关注正例样本。
[0113]
基于训练集与损失函数对fcnarrb 模型进行训练，训练后的fcnarrb 模型用于定位转录因子结合位点，并通过测试集与评价指标对训练后的fcnarrb 模型的定位结果进行测试评价。
[0114]
在本实施例的模型训练过程当中，还选择使用adam优化算法来更新参数。同时对于dropout层的丢弃率参数、l2正则化权重衰减参数和adam算法当中使用的动量参数，采用网格搜索的策略来寻找最优组合参数。具体内容见下表6。
[0115]
表6
[0116][0117]
本实施例通过测试集与评价指标对训练后的所述fcnarrb 模型的定位结果进行测试评价。评价如下：
[0118]
为了验证全局池化层(global average pooling layer,gapl)、细化残差块(refinement residual block,rrb)和进化信息(conservation information,cons)的效果，进行了一系列的消融试验。
[0119]
为了描述方便，将仅使用全卷积神经网络的模型记为fcn，在全卷积网络的基础上加上gapl的记为fcna，在fcna的基础上附加上rrb的记为fcnarrb，最后将cons数据作为输
入数据之一的fcnarrb记为fcnarrb 。结合前面介绍的iou和prauc两个评价指标，结果如下图5和图6所示。
[0120]
从图中可以看出，gapl、rrb和cons均有一定的效果提升。加上gapl的fcna模型在iou和prauc两个指标上实现了对fcn模型的全包围，且在gm12878细胞系选择的21个tfs上平均提升9.32％(iou)和7.61％(prauc)，足以可见gapl能够在全局范围上获取语义信息，从而显著提升模型的效果。加上rrb的fcnarrb模型在iou和prauc两个指标上实现了对fcna模型的全包围，且在gm12878细胞系选择的21个tfs上平均提升7.40％(iou)和4.64％(prauc)，这也表明rrb对于跳跃连接来说能够很好地提取空间信息和过滤噪声信息，从而也能显著提升模型的效果。至于cons，对于模型的效果提升没有gapl和rrb显着，在选择的21个tfs中，iou指标有提升的概率为47.6％，平均提升0.09％，prauc指标有提升的概率是76.2％，平均提升0.27％。
[0121]
在本实施例中，还使用了fcn、fcna、fcnarrb和fcnarrb 模型来查找motif logo。很明显，结果表明，在hocomoco数据库中fcna优于fcn，fcnarrb优于fcna，fcnarrb 优于fcnarrb，可知gapl、rrb和cons对于motif定位是有效的。
[0122]
以上，仅为本技术较佳的具体实施方式，但本技术的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本技术揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本技术的保护范围之内。因此，本技术的保护范围应该以权利要求的保护范围为准。