动态规划实现生物碱基序列全局匹配

问题

1. Find the best globe alignment of sequences TTCGGGGATG and TATATGGGTCG using daynamic programming.
2. Find the local alignments of sequences ATGGTTCCTTGGTA and GGAGTATATTTATGTAC using dynamic programming.

思路分析：

匹配的优劣主要有匹配的得分来衡量，如匹配正确‘GG’为1分，匹配错误如‘GC’为-1分，而一方缺失如‘G_{’为-2分，‘}~’非法。总得分为最后匹配完成后所有碱基对的总和。

直接使用暴力搜索显然是时间复杂性不允许的，所以可以考虑使用动态规划法，而匹配的结果可能不止一个，这又涉及到图论中有向图两点间路径的遍历问题，可以考虑采用回溯法，数据结构方面，使用MATLAB结构体数组表示图，使用栈来实现回溯法。

而输入和结果输出方面使用简单的MATLAB图形界面编程即可，不必过分在意。

实现步骤：

总函数为Align，可以将所有函数文件放在一个文件Align.m里，Align函数只要包含输入部分程序，然后调用GlobeAlignment和LocalAlignment函数即可。

​ GlobeAlignment的功能是实现碱基链的总体匹配，调用了下面四个函数：

Vmat = getStruct(Squence_1,Squence_2);
Vmat= globeStruct(Vmat,Squence_1,Squence_2);
[GlobeResult,Vmat] = globeSolve(Vmat);
printAlignResult(Vmat,GlobeResult);

其中，getStruct函数主要构建二维结构体数组，每个数组元素的结构体元素如下：

alignValue:分值，如x(1,1).alignValue = 1;

alignChar:匹配元素,类型为二维数组，如[A,C]；

maxDrec:可联通元素位置,初始为[],如对于第(1,1)个元素，maxDrec可能是[ ];

isGetted:是否访问过的记录，但不是只要访问过就改记录，而是出栈之后记录，而只要当连续入栈出栈时才依据其作判断。

globeStruct函数主要是对矩阵的前三个量进行赋值，alignValue和alignDrec的赋值规则及最优轨迹的寻找如下图1至图6所示：

GlobeSolve函数为核心程序，是在以上步骤完成后，进行的回溯法遍历起点到终点的轨迹，从而输出所有的最佳匹配。整个算法流程如下（以上图为例）：

1．栈内元素（坐标）为（1,1）（栈顶元素），（2,2），（3,3），（4,3），（5,4）；

2．（1,1）出栈，检查（2,2）点除（1,1）点外的可达点，没有则退栈；

3．现在栈顶为（3,3），再次退栈，栈顶为（4,3）；

4.（3,2）入栈，为栈顶，两可通点任选一个比如（2,2）入栈；

5.（1,1）入栈，输出整个栈，（1,1）再次出栈；

6.（2,1）入栈,(1,1)入栈，输出栈；

7.这一步很关键，（1,1），（1,2）出栈后，按照上述规则（2,2）又会入栈，这样不只是影响起点出栈和遍历结束，加入初始最优路径是（1,1），（2,1）,(3,2),(4,3),(5,4),则程序将反复在这两条轨迹之间遍历，不能遍历所以路径，所以这一步要加控制语句，即当紧接出栈后的入栈时，要检查一下准备入栈的元素的是否入过栈的记录，比如图中（2,1）出栈后，检查（2,2）已入过栈，则不再入栈，直接在退栈。这样做的依据是程序的递归性决定了某个访问过的分支起点（图中为（2,1））不可能再次访问。

​ 其他部分程序思想较为简单，且程序中已有注释，不再解释。

​ 由于LocalAlignment的思想与GlobeAligment类似，由于时间因素，没有实现。

结果表现：

直接运行Align命令或者函数程序得到下面选择对话框：

选择GlobeAlignment，得到输入对话框，输入两个碱基链：

之后得到