关键词 序列比对; 序列拼接; 遗传算法; 启发式; 种子
对于序列拼接的问题[1]，即便在忽略了各种误差的情况下，要作出十分精确的目标序列也是极其困难的。目前，比较普遍的做法是将拼接分为三个阶段[2] ，第一，发现片段交叠；第二，构造片段排列；第三，计算表决序列。
1 若干算法
目前，对于一个典型的较小规模测序问题，处理过程可粗略地分为两步：首先用“散弹枪”法将目标分子打散，然后将打散的片段重新拼接起来。
经过实验处理后获得的片段数量多在1000-2000左右，目标分子的长度也仅为30,000bp-100,000bp。
因而将每一个片段作为图的一个顶点，构造一个完全图进行处理是合理的。目前通常的做法是将所有结点即片段两两通过局部比对，除去嵌合片及重复片段，然后采用诸如最短公超串(SCS)、重构(reconstruction)、多连叠(multicontig)等模型来构造片段排列，计算表决，如由X.Huang提出的方法[3]就是这样。
其中，针对SCS模型又有若干种算法，这些算法的传统目标多为寻找Hamilton最大权路径，这可以采用递归搜索策略，但可能需要花费指数时间[5]; 也可以采用贪心搜索策略[6], 其时间复杂度为Ο(n2)，这个是目前其他各种算法的基础。对于贪心搜索策略主要有以下两个问题，第一：由于Hamilton难题是NP完全的，故贪心算法并不总是返回最短超串的[1]如图1所示，权为4的边是第一个被考察的，从而另两条权为3的边被丢弃；第二：不适应实际情况，为了解决目标序列上的覆盖缺乏问题，通常会增加随机采样，造成所有片段长度的总和大于目标序列的7-8倍以上，而Hamilton路径将所有节点都加入，这不符合实际。
图1贪心算法实例
2 启发式算法
我们采用启发式搜索策略构造出近似的目标序列, 基本思想是：
选择两个长度较长且交叠罚分很低的结点作种子，求得两种子的最大权路径，在计算路径的权值时就判断路径长度是否小于目标序列的大致长度，若大于则停止前进，反之则选择新的不属于已有序列的片段作为新的种子进行延伸，直至长度达到目标序列的长度。
2.1 准备工作
首先，求出每对片段间的比对计分，采用的方法是半全局部比对[4]，其特点是：1.失配计绝对值较大的负分。2.控制失配碱基个数。3.比对两次，第一次A片段的首空格不计分且B片段的尾空格也不计分，第二次反之如图2所示，第一次甲的尾空格与乙的首空格罚分不计分，第二次乙的尾空格与甲的首空格罚分不计分。
第二步，构建一个有向图，将每个片段作为该图的结点，每对结点有两条边。再构建一个n×n的矩阵(n为片段总数)，将第一步得到的半全局罚分填入矩阵，这样每对片段有两条边。图中的每一个结点与矩阵的一个元组对应，每条边的权值与矩阵的一个元素对应。
图2启发式算法实例
2.2 算法描述
构造图G和以root为根结点的最长路径树T，目标结点集合D，Parent和Child作为两个临时节点，建立一个排序队列Sort_Queue，使用队列Store_Queue存放已处理过的节点，定义int型变量seed[]作为存放种子的数组存放着片段的序号，best[ID]定义为到每个节点的最高罚分与两片段表决距离的比值（初始值为-MAX，ID为片段的编号），目标种子定义为全局变量TargetSeed，TargetLen是目标序列的长度, Min是表示序列无关的罚分阙值。其中图1节点的结构为:
Typedef Struct Node *TreeNode
{
int ID; //节点编号
int len //已得序列的长度
Node *father
}
Typedef Struct Link_Node *Link
{
Tree *Node
int dis //节点与目标节点的距离
Link *next
}
long D[][]; //存储罚分的n×n矩阵
Heuristic_Path(G)
Initialize Tree T with sorce nodes and clear Sort_Queue
Compute the scores of the pieces in pair and fill in D[n][n]
While(overlap(seed[0],seed[1])>Min)
{
Random(seed[0])
Random(seed[1])//产生两个种子
root = seed[0] //作为树根
}
int id=0//用于记录其他所有节点序号
While(id∈ADJ(seed[0])){
Enter_Sort_Queue(id)
}
i=1 //为种子计数
TargetSeed = seed[i]