DNA限制性图谱的绘制(数学建模论文)

DNA限制性图谱的绘制(数学建模论文) DNA限制性图谱的绘制 摘要 本文是关于用SPDP法(即简化的部分消化方法)来确定DNA链的基本的片段排列顺序的问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 。 针对问题一，我们运用穷举法，并根据实际DNA链的排列实际将穷举法进行限定，即用简化穷举法。并以此为基础建立二叉树数据结构模型，通过MATLAB软件计算，对问题一中的两个实例进行求解。由于实例一，由于所给数据量非常小，我们可以通过笔算进行试值的办法或者通过MATLAB程序得出其排列顺序为2-6-1-4-3。对于实例二，由于存在可能解一共有八组，我们不再在这里表示出来。 SPDP算法是基于二叉树数据结构的穷举法产生的，首先经过了标准化的过程，然后利用 b-b表与b-c表两个表逻辑、拓扑的关系相互约束，以此对穷举法进行限定，使穷举的次数大大降低，达到实用的效果。但是当所给实验数据较多是任然相当费时费力，而且SPDP所得出的结果也不唯一，所给实验数据越多则实验所得可能的X也越多，还需要根据实际遗传学DNA序列的规律进行排除，选择最终唯一正确的值。所以该方法在实验数据大时实用性较低。 针对问题二，随着实验数据的增多不能分辨的元素也越来越多，使得可能的X也增多。 X的增多意味着任然需要很大的后续工作量。在误差存在的条件下，一组数据的相互分辨能力决定于两个因素，一是误差限的大小，二是数据值分布的疏密程度，即数据的方差(或标准差)。所以，当数据的平均绝对误差限的数量级达到与数据标准差相当的或更大程度时，数据之间的分辨就会变得非常困难，问题的求解也会变得基本无法进行了。 在模型建立后，我们对其进行了推广及应用，并取得了积极有效的成果。 一、问题重述 绘制DNA限制性图谱(restriction mapping)是遗传生物学中的重要问题。由于DNA分子很长，目前的实验技术无法对其进行直接测量，所以生物学家们需要把DNA分子切开，一段一段的来测量。在切开的过程中，DNA片段在原先DNA分子上的排列顺序丢失了，如何找回这些片段的排列顺序是一个关键问题。 为了构造一张限制性图谱，生物学家用不同的生化技术获得关于图谱的间接的信息，然后采用组合方法用这些数据重构图谱。一种方法是用限制性酶(restriction enzyme)来消化DNA分子。这些酶在限制性位点(restriction sites)把DNA链切开，每种酶对应的限制性位点不一样。对于每一种酶，每个DNA分子可能有多个限制性位点，此时可以按照需要来选择切开某几个位点(不一定连续)。DNA分子被切开后，得到的每个片段的长度就是重构这些片段的原始顺序的基本信息。在多种获取这种信息的实验方法中，有一种广泛采用的方法:部分消化(the partial digest, PDP)方法。 在PDP中，采用一种酶，通过实验得到任意两个限制性位点之间片段的长度。假设与使用的酶对应的限制性位点有n个， 通过大量实验，可得到n+2个点(n个位点加上两个端点)中任意 22CCn，2n，2两点之间的距离，共个值。然后用这个距离来重构n个限制性位点的位置(解不一定 ,,,唯一，两个端点对应于最长的距离)。若是线段上的点集中所有点之间距离的集合，PDP ,,,就是给定求。下图给出了一个例子。 图1. A,B是DNA分子的两个端点。 a，b，c和d是限制性位点。 通过实验可以得到 ,,,,,={2,3,4,5,2,5,9,14,16,7,12,14,9,11,7}. 再通过来求，对应于上图的,={0,2,5,9,14,16}是一种解。 上述方法要把DNA分子在任意的两个限制性位点处切开，这对于当前的实验技术来说有相当难度，而且，还要对实验数据进行处理，也很复杂。最近研究人员提出了一种新的方法，称为简化的部分消化方法(SPDP)。这个方法与PDP的不同就在于它避免了在任意两个位点切开DNA分子的难题和处理重复数据的困难。仍假设与使用的酶对应的限制性位点有n个。首先DNA分子被复制成n+1份，前n个复制品中的每一个在一个限制性位点处被切开，最后一个复制品在所有的限制性位点处被切开。这样我们分别得到2n个片段长度(称为第一组数据)和n+1个片段长度(称为第二组数据)。在没有误差的前提下，第一组数据中2n个长度可以分成n对，每对的和都等于DNA分子的总长度;第二组数据中n+1个长度的和也等于DNA分子的总长度。 SPDP问题是如何利用这两组数据重构出这n+1个片段在DNA分子上的排列，使得这个排列在n个位点切开后得到的2n个片段长度与实验得到的2n个长度相等。下图给出了一个例子。 图2. 这个例子对应的位点有4个。(a) 就是我们希望重构的顺序。 (b)中的前4对为第一组数据，它通过切开一个位点得到，每对长度的和都是16，剩下的为第二组数据，含5个片段长度，它通过切开所有位点得到，它们的长度总和也是16， 但实验结果只告知每段的长度，不知道它们在DNA分子上的排列顺序。 现对上述SPDP问题，建立数学模型，并研究以下问题: 1.设计求解该问题的算法， 并评估该算法的效率和效果。对下述2个实例给出答案: 实例1: 第一组数据:2，14，8，8，9，7，13，3 第二组数据:2，1，4，3，6 实例2: 第一组数据:1，14，12，3，7，8，9，6，11，4，12，3，13，2，5，10 第二组数据:1，1，2，1，2，2，1，2，3 2.讨论在实验中测量片段长度时的误差，将在多大程度上影响算法的效果，当误差到多大程度时，限制性图谱的重构将无法进行。 二、模型假设及符号说明 1、假设DNA在用限制性内切酶剪切的时候，每一个限制性位点处被切开。 2、假设每一段DNA都能被准确顺利的计数。 3、在实验中测量片段长度时的误差近似服从正态分布。 4、 记录 混凝土 养护记录下载土方回填监理旁站记录免费下载集备记录下载集备记录下载集备记录下载 第一组数据时在各个限制性切点切得准确，不出现漏切现象 符号说明 符号 符号表示量 备注 L 被研究的基因链总长度 第i个被切且仅切为两段的基因 b i ck 第k个基因片段的长度 bs 每对b中的较小者 C 第二组数组成的集合 cn C集合中的数 bj 与bi 相匹配的位于基因后半段的数 第二组数中与b1相匹配的数，即基因的第一个片段长度 c , cj 基因的最后一个片段长度 p c值重复的个数 c中的某个元素 c , c 每对b相减的差 k p kc数据c中的个数 k c数据中最大的为 c m nc的所有的组合的个数 两个b数相减可以等于 kk 三、 问题分析 SPDP 问题的目标是求出 DNA 限制性图谱，由于 SPDP 提供的两组数据相对于原有的 DNA 限制性图谱已经不完整，通常情况下有多个解，当解不唯一的情况下，必须解出所有可能解，然后利用生物学方法判断这些可能接中哪一个是原 DNA 限制性图谱，如果算法求出的解不全，有可能真实的 DNA 限制性图谱没有被求出，那么对下一步的科学工作会导致错误，所以 SPDP 算法是要求出所有符合第一组和第二组数据的全部可能解。 在解决该实际问题时，我们选择使用穷举法来求解，再通过不断优化改进，得到一个高效简明的穷举算法 。核心思想是通过对第一组数据与第二组数据分析，不仅考虑组内数据的关系，同时也分析组间的数据关系，得到对于SPDP解的约束，在搜索中，利用约束，可以去掉绝大部分的不可能解，只需要在极少剩余的可能解搜索即可得到全部SPDP问题的解。先从第一组数据入手，进行一定的组合，而利用第二组数据对第一组数据的组合进行检验。这里，我们首先要对第一组数据进行一定的处理。在第一组数据中，DNA都是被切成两段，当DNA总长度L知道的时候，其实只要知道每对b中的一个值就得到了全部的信息。现在我们将每对b中较小的一 bb'ii个挑出来(若与相等任取一个就可以)。这样我们就得到了n个数据，再对这n个数据进 b，b，b？？b123n行升序排列，得到一个不严格递增的b序列，在这里不妨记作，它们满足 L,,,？？,,bbbb123n2，它包含了第一组数据的所有的信息。 为了方便计算，我们不妨先对基因图谱的起点和终点做一个 规定 关于下班后关闭电源的规定党章中关于入党时间的规定公务员考核规定下载规定办法文件下载宁波关于闷顶的规定 ，假设从起点算起，第一 c,ckk1n，1个位点到起点的距离小于或者等于从终点算起第一个位点到终点的距离，即。 bbii对于重新组织后的第一组数据，其含义为:对于某一个，它表示离这段基因起点或者终点 bi距离的位置有一个位点。则对于第二种穷举方法也就由此而得出，对于每个它有两种状态， bi一种是靠近起点(不妨用“0”表示)，另一种是靠近终点(不妨用“1”来表示)。而当全部的的状态确定后，DNA的顺序也就确定了，但此时要对这个顺序的合理性进行检验，即检验这个顺序是否符合第二组数据c的要求。具体 措施 《全国民用建筑工程设计技术措施》规划•建筑•景观全国民用建筑工程设计技术措施》规划•建筑•景观软件质量保证措施下载工地伤害及预防措施下载关于贯彻落实的具体措施 如下: b，b，b？？bbbssssii1231l将“0”状态的全部挑出来按照升序排列。再将“1”状态的全部挑 b，b，b？？bssssl，1l，2l，3n出来按照升序排列。记第二组数据c组成的集合为C,{ccc？？c}123n，1。然后检验是否满足以下条件: 当全部检验条件符合的时候，这组组合即为合理的顺序，从DNA的起点算起，在前半段的基因 b,bb,0sssii,10中，第i段长度为 (其中)，而对于后半段则有相同的处理，其示意图如图5表示。流程图如图4表示。这里值得说明一下的就是对与这个循环的结束条件，因为对于基因的排序问题，必须要求出所有的可能解，因此对于这个循环就必须要把所有的可能组合全部进 n2行一遍以后才可以停止，即需要循环次。 组合b 不合 条 件 输出答案 用c检验 符合条件 图4 bb bbbbssssss12lml，2l，1„„ „„ b 图5 s1 四、模型建立与求解方法 问题一 通过查找文献，我们发现SPDP问题有它其中的一些特有的性质与规律，我们希望能找到这样的性质与规律来减少循环的次数，并且将复杂的问题化约为标准的情形便于算法的实现。基于上述这种考虑，我们寻找了一些性质，并将比较重要的总结为一些基本的定理，以便后面进行算法设计与优化。 在给出定理之前，先对一些术语做一定的规定。首先规定“同侧”与“异侧”。我们将整段基因 b与bb与bijij分为前半段与后半段，以基因正中为分界。在第二种模型中，称为同侧，如果所 b与bb与bijij表示的位点的位置在同一段中(同属于前半段或后半段)。称为异侧，如果所表示的位点不在同一段中。而对于整段基因的每个片段我们都将认为它是一个元素。那么从起点开始的第一个片段我们称之为“首元素”。而从终点开始的第一个片段我们称之为“尾元素”。而包含基因中点位置的基因片段我们称之为“中间元素”。 下面我们给出SPDP问题的一些基本定理 定理一(首元素定理) L,,,？？,,bbbb123nb，b，b？？bb2123n1 在模型中，对于第二组数据满足。则为 cb,bc，c，c？？cb,1s123n，111首元素片段的长度，。即在中存在与相等，并且从起点算 cb,1起第一个基因片段的长度就为，或者说该片段就是。 定理一说明:此定理唯一的确定了首元素，节省了排列这个元素的时间，它更重要的意义在于把问题进行标准化。 定理二(尾元素定理) cb,1 在确定了首元素之后，我们把第二组数据中的除去，把与其对应的除去。在剩下的b cb,cbjiji与c中，我们寻找相等的b与c，比如，那么这个就是可能的尾元素，为它的对 bi应元。即从整段基因的终点算起，第一个片段(实际上是这段基因最后一个片段)的长度为， cj或者说最后一段基因就是。 定理二说明:此定理给出了尾元素(最后一个片段)可能的情况，与首元素一起可以将整段基因的首尾确定下来。当尾元素的可能性超过一个是我们每次只考虑其中一种首、尾元素的情况，然后对首、尾元素的组合进行穷举。可以看出这个循环最多是n-1次。 定理三(首尾排序定理) bb,b,b？？bbs234s,11 当首、尾元素确定后，若首元素长度为，尾元素长度为，则所表示的位点的位置必在同侧并且在前半段，并且它们所确定的位点在整段基因中彼此相邻，即b与bii,1所确定的两个同侧位点之间没有别的位点，其中i=2,3,4„„s-1。 bb,bsss,1定理三说明:首先，对尾元素要进行一个补充，实际问题中有可能出现的情况， bbb,bss,1ss,1此时我们的脚标s选取两者脚表较大的一个，取而不取。而且对于相等的情况，bbss,1与一定在异侧，表示在与中点对称的两个位置有两个位点。而且由此观点我们还可以知 bbbbb,bb,b,b,bs，1ss,1sss,1ss，1ss,1道如果，则必然有。因此严格大于，严格小于。其实这个结论也非常重要，但很直观，就不在这里做为定理。而定理三一方面可以节省b,b,b？？b234s,1在组合中的状态选择，更主要的是当首尾元素确定后，它可以将SPDP问题进行标准化，即通过首尾排序定理后，我们可以把最靠近起点的s-1段基因看作一个整体(因 bb,bs,1s,1s为这s-1段的顺序已经知道了)，当作一个基因片段，它的长度为，并且。而b 中剩下的元素都比首尾元素大。 定理四(尾元素约束定理) bb,bs,1ss,1 在定理三的说明中，我们知道b中的元素有可能两两相等。如果，我们记录为 bbb？？bkkkk123l重值代表。我们找出所有的重值的情况并记录其重值代表，其中角标满足 b,bbk,k,？？,ks,ksks12l11。若尾元素长度为，则有约束，即。 定理四说明:此定理是对尾元素定理的一个补充，通过这个约束限制了尾元素的可能情况，降低计算的复杂程度。 定理五(中间元素对偶定理) bn 由前面的定义知道中间元素是包含中点的片段。那么此片段的两端有一端的位置必在所 bj确定的位点，而另一段的位置若为所确定的位点，那么一定有在第二组数据c中有元素与 bb,bbjnjn相等。则这类就为另一端点。并且若将视为从整段基因中点往两端算起的首位点 bj(首元素)，它是距离中点最近的位点，而将视为中点另一侧的第一个出现的位点，那么中 bbjn间元素的两端可以等效为基因的起点和终点，等效为首元素，等效为尾元素，那么定理一到定理四的性质在此处都将成立，这就是中间元素的两端与起点、终点的对偶性质。 定理五说明:前四个定理实际上是从基因的两端出发并不断的向中间考察，而定理五刚好相反，是先确定最中间的情况，再想两边考察。对于定理五的应用在3.5进行讨论。 定理六(相邻重值等效方法) b与bii，2 若均为重值代表，即 b,b,b,bii，1i，2i，3。那么在第二组数 据c中一定有两个相等的 c,c,b,b,,i，2i，并且在整段基因 中，它们的位置如图6所示。因此这两段 基因片段的位置就已经确定下来，因此我图6 们可以将A点与B点接起来，C点与D点 接起来。这样，原来的基因长度就缩短了 2(b,b)b,b？？bi，2ii，4i，5n。在数据中除去上述几个元素，在剩余的b中我们要对进行一定 b'bb',b,(b,b)kkkki，2i的处理，对每个数处理后的值为，其中k=i+4, i+5„„n，为处理后的值。当如果有更多的重值代表两两相邻，我们可以用同样的方法在一开始就把这些位置已 经确定的片段去除，等效出新的一整段基因以及它的初始数据。当最后计算出全部的结果之后再将这些已知的片段插回去。 定理六说明:此定理最重要的作用就是可以将问题进行标准化，使要解决的问题不含有那些已经知道基因片段位置的这些特殊情况，从另一个角度看这也是通过对数据特征的观察降低计算的复杂度。 标准化SPDP问题的建立 有了上面一些基本规律、定理的准备，我们就可以建立标准化的SPDP问题。建立过程如下: 第一步:利用定理六(相邻重值等效方法)去除由于相邻对称的位点所确定的基因片段， 调整两组数据b与c的值。 第二步:利用定理一(首元素定理)确定首元素，除去b中的首元素元与c中的对应元。 第三步:利用定理二(尾元素定理)与定理四(尾元素约束定理)确定尾元素的可能值。 并在计算中对每个有可能的尾元素进行循环，当确定其中一个可能值作为尾元素 时，去除b中的尾元素元与c中的对应元。 第四步:利用定理三(首尾排序定理)，并将前s-1个片段看为等效首元素(尾元素为bs,1 )。这里要注意:若发现前s-1个片段中存在一个片段的长度在c中没有对应bs 相等的值时，说明所选的尾元素并非真正的尾元素，此时返回第二步。 经过上述四步处理后，SPDP的原模型已被标准化。对此时的第一组数据b重新按递增排序， L,,？？,,并重新给予脚标，将重值赋予同样的脚标。这样就可以得到新的，并bbb12l2记录有重值的b，仍用{}这样一组数来表示。由定理三说明中的分析知道b若有重值最多只bki 有两个元素等相等。而对于第二组数据c也同样重新赋予脚标并排序得到，c,c,？？,c12m并用记录每一个c值重复的个数。 p,p,p？？p123m 通过上面五步的操作，我们得到标准化的SPDP问题。为了保持定理的习惯，以及后面算法设计的方便。我们对于第一组数据b不去处首尾元素，则为首元素， 为尾元素，b数据满足bb12 L,,？？,,bbb.而在c数据中我们却已经把首尾元素对应相等的值已经去除。则对于12l2 bb，b？？bi34l标准化问题，我们要做的事情就是要确定的状态组合(判定每个属于前半段还是后半段)。而我们的数据结构与算法就需要针对这样的标准SPDP问题进行选取与设计。 SPDP数据结构的选取(利用二叉树) b，b？？b34l 当我们已经将原模型转化为SPDP的标准问题后，我们就开始需要对的状态 bi进行确定。因为在确定它们的状态的时候是逐一进行的，因此每轮到确定某个，它就有两个可能的选择(属于前半段或者属 于后半段)。根据这样的一个特 点我们选取二叉树作为SPDP标 准问题的主体数据结构。 接下来的问题就是要确定我们 bi如何的逐一的对进行选择和图7 建树。这里，我们选择的顺序是 bb33从开始从小到大的对b数据进行判定并建立二叉树。我们以为例来看下如何建树。如图bb,bb,b331327所示，要不然在前半段，要不然在后半段，并且我们考察与,这两个差值之 cbbb,312中至少有一个可以与c中的某个元素相等。其原因为:是首元素，为尾元素，并且是 bb12数据b中除了、最小的元素，它所决定的位点应该是距离前半段第一个位点或者后半段第 a或者a12一个位点距离最近的点，一定是c中若干片段中的一个，不然在这种首尾元素组合下 ccbb,bb,b,,31323没有合理解。这样，我们就可以通过判定能否在c中找到=或者=而判定 cbb,b,331状态的可能性。倘若=，那么就有可能在前半段，建立二叉树时现在所处在的结点 ccb,b,,32就有左树枝，而左树枝结点的数据中，c中的数据就要去除。当=，则现在所处的 cb,b,32,结点有右树枝，并且对右树枝所在的结点做同样的处理。如果，则右树枝为空， bi并且不再被考虑，因为那是一条不可能的路径。依次类推，我们从小到大依次考虑每一个， bbii并通过上述的判定来建立二叉树。最后需要对重值的b做一个解释。如果有重值，当考虑到 bi时，应该对两边同时进行判断，是否两个能同时在两段都成立，若有其中一边不成立则说明 bi这条支路无解。如此运行，直到所有的都被考虑完，循环停止并将最后合理的组合所确定的 开始 判断该段b能 进入下一个 否放入前半段节点进行搜是 索 DNA图谱 否 该节点子树判断该段b能 没有解，搜否放入后半段否 是 索其他节点 DNA图谱 图8 DNA顺序输出。关于这个数据结构的流程如图 对SPDP算法进行进一步的优化，建立b-b表与b-c表 我们建立二叉树数据结构的思想来源于上面建立的数学模型，而这个数学模型建立的立足点是通过对第一组数据b进行组合，而拿第二组数据c来进行检验。不论是最初的穷举法还是二叉树数据结构的算法都体现了第一组数据与第二组数据这样的一种关系。然而这样的关系并非完美，因为尽管c数据比b数据包含的信息要少，但对b进行组合的时候，数据c仍然是非常重要的一个约束，并且时刻的在约束着b的每一步组合。基于这样的想法，我们希望在二叉树的结点计算时可以加入c对b的约束。为了实现这样的想法，我们在SPDP的标准问题中建立了b-b表与b-c表。 b-b表的建立 首先我们建立b-b表。在二叉树数据结构算法中，我们每步对结点树枝的是否存在的判断都是通过两个b数的差是否在这个结点中数据c(可以将它看作一个集合)中。因此判b,bij 断是否等于至关重要。由此，我们就对数据b的所有元素两两作差，建立b-b表。我cb,bijk 们可以让横坐标表示被减数的脚标i,而纵坐标表示减数j，(j，则在(i,j)的位置记为-1，这样区cb,bcmmij 分的好处在后面的算法里可以体现。那么，二叉树结构中每步=的判断就转变成为判b,bckij 断b-b表(i,j)位置上的值。 表2.b-b表 减数角标 14 403 1002 2301 被减数角标345 b-c表的建立与b-c表约束 b,bcb-b表是等式=中i与j之间的关系，而下面我们建立的b-c表是j与k之间的kij b关系。如表3所示，表3是最基本b-c表，横坐标为c的脚标k，而纵坐标是被减数的脚标j b,bcb,bcj。如果=，那么在(k,j)的位置记i，若，则在(k,j)的位置记0。,ijkijk 这样，这张b-c表就有了它自身的含义，它表示着对于每一个，有可能组合出这个的两cckk个b数的差的所有可能情况。这里就有一个非常重要的约束，对于每一个，两个b数相减可ck以等于的所有的组合的个数与数据c中的个数满足+1。这里要注意，之cncpnpp,kkkkkkk所以要加1是因为“中间元素”并非为两个b值相减而得到的，所以有可能在这里的个数cpkk比可以小1。但一旦在出现了这种情况，在别的c数据的列中就不能出现这种情况，即对nckk 于其它的c数据要求。倘若出现了与约束相矛盾的情况时，再往后计算肯定得不到合np,,, 理解了，需要换下一个可能的尾元素。 表3. b-c初始表 b50000 b45000 b30455 b20304 b13003 c1c2c3c4 p1p2p3p4….. 得到了b-c表的这个约束后，下一步需要解决的就是如何在二叉树建立的过程中在不改变开始时候做的b-c表而对每个结点的进行计算并与该节点的相比较。首先，我们发现nppkkk在每次节点中的变化是因为每次计算需要除去c数据中的一个值而导致的变化。对于每个节pk 点，都需要记录所有的的值。然而，对于,我们并不需要每次都改变b-c表，而是直接从pnkk b-c表中计算出当前情况n的值。这样，我们就需要对b-c表进行一个补充。如表4所示，对k 于每一列c,我们又增加了两个子列，一列为状态数(J函数)，另一列为排序(R函数)。我们k 现在来考察二叉树在建立过程中其中的一个节点的情况，如表4所示，假设该 表4. 完整b-c表 J函数被减数R函数J函数被减数R函数J函数被减数R函数J函数被减数R函数 b5002002001003 b4152002001003 b3001142151153 b2001131000142 b1131000000131 c1c2c3c4 bbb节点的前半段已经考虑到，后半段已经考虑到。当前需要考虑到底属于前半段还是属j,i b于后半段。倘若被放到了前半段，生成该节点的左树枝及其节点，我们发现，对于这个新生, bb成的节点，对于任何，-在这里已经是不可能的组合，同理我们发现，纵坐标比i,,,,i 小的与i,j之间的那些可能与c相等的b的差值的组合在这里已经没有意义了，需要在所有初始时b的差的组合总数中减去这些在这个状态下没有意义的组合个数。因此我们这样来定义我们的J函数与R函数:由于在3.3.1标准化SPDP问题的过程中应用了定理六(相邻重值等效方 法)，如果对于某个所对应的一列，若(j,k)上的值不为0(即=)，那么J(j)=1，ccb,bkijk j,1 否则为0。而R(j)=。那么R(j)表示以从到中为被减数的所有有可能被别的bbbJ(s)j,1,1s,1 相减等于这个的个数。那么很容易，对于新生成的节点，=。对这个cnR(l),R(,)，J(j)kk公式举一个具体的实例加以说明，如表4，当我们已经计算到前半段考虑到b4后半段考虑到b3 时，，这表示当前情况下只有一种情况可以得到c1,即n,R(5),R(4)，J(3),2,2，1,11 b5-b4=c1。 经过上述的改进，我们就可以通过一个静态的b-c表对c中所有的元素在建立树的过程中 方便的检查它的约束。如果与约束矛盾的情形前面已经讨论过，而b-c表还有非常好的一个特 性就是它在某些情况会有“临界”发生。比如对于所对应的一列，当在某个计算过程中，发ck现或者(中间元素的情况)说明可以组合出的b的差值的个数已经和ccp,np,n,1kkkkkk 的个数相等，那么这时已经发生了“临界”的情况，要求对于当前考虑建立树的节点以及它b, 以后的中凡是满足J(s)=1 (s),与(s,k)所对应的一定要在同侧，而bbb,,,ss 一定要在与的另一侧。这样大大的提高了计算的速度。我们称这个约bbb,b？？b,s，1s，2,,1s 束为b-c表约束，它的本质是数据c对数据b组合的约束。 对b-b表、b-c表的进一步讨论 在上文中我们分别建立了b-b表，b-c 表这两张非常有用的表。对于b-c表的约束 与作用在3.4.1中已经做了详细的说明。我 们发现，其实b-b表，b-c表本质上是等价 的，它们都来自于b,b是否等于c这个kij 逻辑关系。而从b-b表中，我们可以从图的 角度发现对这个问题更深刻的理解.我们将 一个充满了数据的问题由此全部转化成为 一个图的问题，通过对这个图的逻辑关系的 分析与讨论，可以得到更加好的优化方法以 及解决带有误差问题的思路。在2.2的问题 分析中所举出的那么一组合理解具有指数 个数的特例也是通过对b-b表的图的关系分 析所构造出来的。下面就先着重对b-b表的 图的逻辑关系进行一个分析。 图9 D G E A 图10 在b-b表中。每一个点(i,j)都表示一个的操作，而(i,j)上所记录的值是b,bb,bijij这个操作与的逻辑关系。而这于前面我们所讲的二叉树结构的建立是完全一样的。那么我们ck 就尝试将二叉树建立的过程就放到b-b表中，使b-b表中的每一个路径都对应二叉树中的一个分支。下面我们就以一个简单的例子来看看b-b表中的路径是如何的与二叉树中的分支相对应。如图9与图10,首先为了看两个图的对应关系，我们先不考虑二叉树建立过程中每步对c的判断，也不考虑b-b表中每个点(i,j)上记录的值。那么在图9中，我们可以看到二叉树的一个分支A-D-E-G,而b-b表中的A-D-E-G刚好与之完全对应。我们来看这个过程:在初始状态b、b12分别是首元素与尾元素。有两个选择，一个是放在一边，而另一个选择是放在一边，在bbb312b-b表中，它的选择就是A点(A点在图中表示(3，1)点)与B点。而对于我们这个分支它选择了A点。下一个考虑的是，同样道理它也有两个选择D点与F点，容易的看出，每次选择b4 时，b-b表中的图的这一列的最上面的点(如点B、D、F、G)都是可以选择的一个路径，而对于另一个路径在取决于前面的形成过程。比如考虑时，由于B点是上一次的一个选择确没有b4 选择，则将B点平移到了C点，成为了这次考虑的一个路径的选择。同样的，在这个分支选b4择了D点，当考虑b时，它的选择就有E与F; 5 b当考虑时，它的选择有G与H，而最后选择了G，最后形成了A-D-E-G这个路径。而在实际6 b,b建立二叉树时的每一步都需要检验是否等于ij c，则在b-b表中就需要经过的点(i,j)上记录k 的值不能是0，必须为某一个k，并且还要要求这个 c在已经到达这个点时仍然在数据c之中(c中的k 元素是不断减少的)。这样一来，建立二叉树的过程 就全部可以在b-b表的图的关系中得以表现。下面 我们再来观察b-b表具有的更多的特性。首先，如 b,b果b-b表中每一点的值均记为，则我们发现ij 在同一列中，下面的点的值比上面的大;同一行中， 右边的点的值比左边的大，则我们用箭头表示值增 大的方向，我们可以得出图11，这是点的值的大小图11 分布图。下面我们就用这样的关系来构造一组合理 N 32解的个数在个以上的数据(即问题分析中所举 bb出的实例)。首先我们希望构造一组数据，使得为首元素一定在前半段，为尾元素一定在12 bbbb后半段，而一定可以既能在前半段也能在后半段;在前半段，在后半段，一定既能3564 在前半段又能在后半段，如此一直类推下去。如果能构造出这样的解，它解的个数一定不少于N 3，而我们现在将刚才所说的规律所要求经过的点画在b-b图上，得到图12，这样我们就发2 现了规律，它很有规律的六个一组的进行重复。现在的任务就是需要对这些点赋上c的值，使上述对b-b表的点的值的分布可以看出，如果我们规定了，则对于这个问c,c,c,？？123 题最好的构造如图12所示。 图12 在这张图的指导下我们得到了我们构造的数据的递归式: 比如我们可以利用这个递归公式得到这样的初始条件: 表5. b数据1234591314243435577980126172173c数据1122345910112122234546479394 问题分析中的数据测试就是以这组特殊的实例来说明了由于SPDP问题要求要把所有的合理解都要计算出来，而这组数的合理解就已经有指数个，因此不论如何优化算法都不可能得到多项式的计算方法。 图12 图14 b-b表中对在二叉树的建立过程中也可以起到与3.4.2类似的效果。根据前面所述的b-b表中路径选取的规律有图13-16这几个图中任何一个情况发生，该图没有一个合理解。其中，图13是图14的特殊情况。有了以上的结论，我们又发现，当图13-16这些图中有一个点并非0(或者-1)时，不论采取何种路径，在这个图中这个点必然要被经过，比如图17-18所示的情况。这个情况实际上是与b-c表中临界后的点的情况是一样的，统称这些为临界点(注:但它们临界的原因不同)。而在b-b表中，倘若A点为临界点，则A点所在的行与列上其它所有的点都不可能被经过，此时将这些本来不为0的点改为0，称这样的操作为b-b表的约束。这样，我们就可以成功的将b-c表约束与b-b表约束结合在一起，大大优化算法的复杂度，具体过程如下: 1. 建立b-b表，进行b-b表约束，如果出现了无解状态就停止。 2. 通过b-b表建立b-c表，再进行b-c表约束。如果出现了无解状态就停止。 3. b-c表约束后产生的临界点再带入b-b表进行b-b表约束。b-b表约束后修改b-c表， 如果b-c表没有产生新的临界点，将进入计算;否则继续进行3。 下面我们将讨论一下b-b表能否也如b-c表的约束一样，在二叉树的建立过程中也进行约束。理论上说是可以的，只要将b-b表中未考虑的元素的区域取出来，根据已经到达的这一步的实际的c数据修改b-b表(或重新建立b-b表)，再进行上述的算法优化。然而对于每一步都做这样的事情，实际上也在增加了时间复杂度。b-c表约束也存在这样的问题，只是它的复杂度比b-b表约束要小的多。而我们在这里就可以采取一个策略，就是并不是在每一步都进行b-c表约束与b-b表约束，而是在某些步的时候进行抽查，比如说当树支多到一定的程度的时候就采取一次优化。这两个表在实际的算法中起到了巨大的作用。 最后，在这里想讨论一下对于这个问题的终极算法设计的可能性。首先解释终极算法:终极算法就是在建立二叉树的过程中可以在最早的时间预测到某个支路继续走下去将是“死胡同”。这种想法与b-b表约束的道理类似，然而，b-b表只是用了部分的信息，它并没有将c所包含的所有的信息用完全。如图19所示，A、B、C、D点的值同为，然而，如果的个数只有两个cc11的时候，就相当于在路径的选择过程中A、B、C、D四个点选择两个点经过。然而如果将这样的约束引入到b-b表约束，实现这个约束的复杂度就已经到达了指数的难度。因此终极算法的实现几乎是不可能的。然而幸运的是对于实际的数据，利用已有的两个约束已经可以大大的降低了复杂度，使它相当接近终极算法。 最终SPDP问题算法 经过了以上的准备，我们可以得到了最终对SPDP的算法。 1.对具体的数据进行标准化处理，得到SPDP的标准问题 2.建立b-b表与b-c表进行3.4.3中b-b表、b-c表共同约束，如果出现无解情况返回1的循环中;否则进入3 3.建立二叉树数据结构，并在b-b表、b-c表约束及其共同约束的操作下进行有效穷举，输出答案。之否再返回到1的循环，直到1的循环结束为止。 最后在本节小结SPDP问题算法:SPDP算法是基于二叉树数据结构的穷举法产生的，首先经过了标准化的过程，然后利用 b-b表与b-c表两个表逻辑、拓扑的关系相互约束，使穷举的次数大大降低，达到实用的效果。 问题二 显而易见，实际的SPDP问题中因为生物化学手段的限制，实际测量中不可能得到准确的DNA片断长度，测量时的误差是不可避免的。所以我们必须考虑在有误差的情况下，如何在原算法的基础上实现SPDP的算法: 确定b-b 表和b-c表 b-b 表和b-c表b-b表在程序的匹配判断当中起到了相当重要的作用。所以我们应该首先在误差存在的条件下将b-b 表和b-c表确定下来，这样在这种情况下程序的算法也就相对确定了。为了确定这两张表, 我们需要了解误差对确定变量值之间的关系。于是我们首先需要明确下列问题: 长度测量值的概率分布 由于片断长度的测量值误差是由很多微小的因素引起的分量叠加而形成的，在统计学上我们认为这些微小的分量满足李雅普诺夫(Liapunov)定理的条件。于是由该定理可知:同一个物理量的测量值的误差近似服从正态分布(其数学期望为0)。[2]也就是说，每个物理量的测量值满足数学期望为其真值的正态分布。 22b~N(,,,)c~N(,,,)bcibbicciiiiii 记为，，其中和分别为第一组数据和第二组数据中的长度测量值。 由3σ原理可知，其绝对误差限可以认为[注]是: r,3, bbii r,3, ccii 测量量的误差限可以由仪器精度和方法精度进行估计，也可以通过样本方差进行估计，所以可以当作是已知量。 两个片段长度的不可分辨判据 22c,c~N(,,,,,，,)ijccccijij 由概率论原理可知， 故可以认为 2222(c,c),(,,,),3,，,,r，rijccccccijijij ,,,c和cccijij 故当时，即的真值相等时，其测量值之差的可能范围满足 22c,c,r，rijccij 所以一旦其测量值之差满足了这个条件，我们就不能判定这两个值的真值是否相等。于是我们称这两个测量值不可分辨。上式的条件称为不可分辨判据。 b-b表的确定 由于b-b表和b-c表是完全等价的，故我们只需确定b-b表，就可以推算出b-c表。 b,bcijk 如前面所述，b-b表的某个域T[i][j]的值决定于是否等于某一个。在无误差的 b,bccijkk计算中，需要严格等于。但是在包含误差的计算当中，很可能没有一个能够严格的 b,bij等于，而我们需要判定是否他们可以相等。 222b,b,c~N(,,,，，),,,,,,ijkbbcbbcijkijk 由于 b,bcijk 故当和的真值相等时，可以认为 222222b,b,c,3,，,，,,r，r，rijkbbcbbcijkijk b，b，cb,bcijkijk 当满足上式条件时, 和的真值就可能相等。在这种情况下，为了不漏掉可行解，我们就认为他们在误差允许条件下相等。 b,bccijkk于是就有可能出现多于一个的和 相等。这些也是不可相互分辨的，即单凭长度的信息不能将他们区分开(在DNA链上位置互易后但从长度上讲仍然满足所取得的数据要求), b,bccijkk这跟他们相等的情况是等效的。所以认为这些等效相等，所有满足和 相等的的集 AC,{c|1,k,N，1}ijk合记为,称为等效集合。记(其中N为限制性定位点的个数)，故易知 A,C或者C,{c}:ijkc1,i,N,1,j,Nk,其中为上面的算法中最后出现在中间的元素。我们需要将这个 Aij并集进行划分，但是并不一定两两交集为空。所以我们要对这个并集进行如下处理: Aij ?1与任意其他集合相交为空的集合不变。 AA,Ai’j'i'j'ij ?2如果,去除集合。 A:A,,i'j'ij ?3除?2的情况外，如果，则 maxA:A,minA:AA,A:Ai'j'iji'j'ijiji'j'ijA. 如果 满足不可分辨判据，则令，而去除 Ai’j'集合。 maxA,maxAminA,,minAi'j'iji'j'ijB. 否则，不妨设, 则必有。 AAminAijiji'j'令 为原集合的元素中，与不满足不可分辨判据的元素的集合。 AAmaxAi'j'i'j'ij令 为原集合的元素中，与不满足不可分辨判据的元素的集合。 AAA:A,x,A,,,x,AAleftbigi'j'ijbigsmallsmall其他中元素的集合可以作划分{，},使。 AAA:AA:Aiji’j'ijbigi'j'small令 = ， = 任意满足条件的划分都要考虑，因为每一种划分将对应不同的c向量，也就将对应不同的解。 AAiji’j' ?4在?3情况中将和划分开后，如果各自仍与其他集合有非空交集，则继续按照?2和?3的原则进行划分。 这样确定下来的集合就形成了原并集的一个划分，每个集合内的元素在程序判断时当作相等处理，将它们归并在一起，元素的值取划分当中元素的均值。 b,bcijk 每一种c集合划分，表中域T[i][j]的值则为离最近的新值的下标。 下面以下表数据为例对上述方法进行说明。 表10. c[1]—初始对应的Aij c[n+1]b[1]--b[n] 113100A31A86A97 118101A31A86A97 214104A53A86A97 230105A53A86A97 319109A64A86 340110A52A64 429110A52A64 444112A10A52A64 529113A10A52A64 118A20A41 (假设所有数据相对误差限为1%,并且为了形式统一和算法设计方便设b[0] = 0 ，无误差) 按照?1?2?3?4的原则进行划分，最终得到c的两组划分 {c~c},{c~c},{c}1459101、 {c~c},{c~c},{c}1569102、 所对应的b-b表均为: 表11. 00 001 001-1 002-1-1 001-1-1-1 000-1-1-1-1 002-1-1-1-1-1 0013-11-1-1-1 023-1-1-1-1-1-1-1 从这组数据可以看出，当误差限相对于数据的标准差来说较大的时候，会产生很多不可分辨元 素，对问题的求解带来了很坏的影响。 五、计算机结果及分析 问题1 根据上文中的模型和算法，运用Matble软件编程计算除了结果，结果如下: 实例1 第一组数据:2，14，8，8，9，7，13，3 第二组数据:2，1，4，3，6 利用上述算法:计算出相应的可能DNA限制性图谱: 2，6，1，4，3 实例2 第一组数据:1，14，12，3，7，8，9，6，11，4，12，3，13，2，5，10 第二组数据:1，1，2，1，2，2，1，2，3 利用上述算法:计算出相应的可能DNA限制性图谱: 第一组: 1，1，1，1，2，2，2，2，3 第二组: 1，1，1，2，2，2，2，1，3 第三组: 1，2，1，1，2，2，3，1，2 第四组: 1，2，3，2，2，1，1，1，2 第五组: 1，2，1，2，1，3，2，1，2 第六组: 1，2，2，1，2，3，1，1，2 第七组: 1，2，2，3，1，2，1，1，2 第八组: 1，2，1，3，2，1，2，1，2 问题2 从前面模型中给出数据以及带误差的算法过程可以看出，随着误差的增长，会出现两种问题: 1、不能分辨的元素越来越多，这样大大地减少的解的确定度。每一组n个不能分辨的元素， n倍。 都要使可行解的个数增加Pn 2、对C划分的可能性数量增多。C划分的可能性数量是随着A之间重叠元素的个数呈指ij 数增长的。 在误差存在的条件下，一组数据的相互分辨能力决定于两个因素，一是误差限的大小，二是数据值分布的疏密程度，即数据的方差(或标准差)。所以，当数据的平均绝对误差限的数量级达到与数据标准差相当的或更大程度时，数据之间的分辨就会变得非常困难，问题的求解也会变得基本无法进行了。以前面所给数据为例，这组C的标准差为5.652，而b, c，的平均绝对误差限为2.009。数据已经较难分辨了。如果相对误差限变为2%， 那么数据就会变得基本无法分辨了。 六、参考文献 [1] 蒲俊，数学建模方法与训练[M]，四川:西华大学数学建模 培训 焊锡培训资料ppt免费下载焊接培训教程 ppt 下载特设培训下载班长管理培训下载培训时间表下载 教材，2016 [2] 盛骤等概率论与数理统计，p149-150 北京:高等教育出版社，2001.12 [3] 刘若鹏，DNA限制性图谱绘制SPDP问题的解决算法[J]，2010 七、附录 %data1=[1 14 12 3 7 8 9 6 11 4 12 3 13 2 5 10];data2=[1 1 2 1 2 2 1 2 3]; function DNAreset(data1,data2) matrix=[]; DNA_sequence=[]; sequence=1; for i=1:length(data1) for j=i:length(data1)-1 if judge(data1(i,:),data1(j+1,:)) temp=data1(i,1); data1(i,1)=data1(i,2); data1(i,2)=temp; end end end for i=1:length(data1) matrix=data1; matrix(i,:)=exchange(matrix(i,:)); for j=1:length(data1) matrix(j,:)=exchange(matrix(j,:)); for h=1:length(data1) matrix(h,:)=exchange(matrix(h,:)); for g=1:length(data1) matrix(g,:)=exchange(matrix(g,:)); temp=search(matrix,data1,data2); if length(temp)~=0 DNA_sequence(sequence,:)=temp; sequence=sequence+1; end end end end end len=length(DNA_sequence); del_number=0; for i=1:len-del_number j=i+1; while j1&j