Perl基础实例---多序列的最长公共子序列LCS
在生物应用中, 经常需要来比较两个或者多个不同生物dna, 并且求其的最长公共子序列。
最长公共子序列问题(简称LCS),是给定N个序列X1,X2,。。。Xn 求其最长的公共子序列。
求2个序列的最长公共子序列, 一般情况下最佳的方式是动态规划, 今天我们不讲动态规划。
我们有3个序列X,Y,Z, 今天通过暴力搜索的方式来求其最大公共子序列, 该方式不一定是最佳的解决方案, 但一定是可行的解决方案。
我们定义3个序列:
my $x = 'ATCGGAGCAGT';
my $y = 'TTCCAGAGA';
my $z = 'TGAG';
my @all = sort {length $a <=> length $b} ($x, $y, $z);
my $shortest = shift @all;
my $length = length $shortest;
my $length_tmp = length $shortest;
我们先取出长度最短的一个序列, 以及计算其长度备用。
由于我们需要取最长公共子序列, 意味着这个结果$result必然同时是所有序列的子序列, 那么最长的情况也就是结果$result等于最短的序列$shortest。我们采用暴力破解的方法, 即一一枚举所有的可能, 这样肯定可以得到我们需要的公共子序列。由于我们只需要最大的子序列, 不需要所有的子序列, 因此我们从最长的可能(即$result = $shortest)开始枚举, 依次减短长度, 只要我们得到一个符合条件的, 那么其必须是最大公共子序列(也可能是同长度的几个情况之一)。
设置循环熔断条件, 当所有情况枚举完, 仍然没有结果, 那么就推出:
while($length > 0){
依次从长开始取$shortest序列的子序列,不成功则尝试下一个:
for my $index(0...$length_tmp - $length){
my $flag = 1;
我们设置一个标志位, 用来记录中间匹配结果
由于我们有3个序列, 因此需要匹配多次, 只有所有序列有成功才能算成功
通过substr函数来取子序列
my $sub = substr($shortest, $index, $length);
for my $a(@all){
只要有一个序列匹配不成功, 那么意味着不是公共子序列
则马上退出, 并尝试下一个子序列, 这样可以提高效率
unless($a =~ /$sub/){
$flag = 0;
last;
}
}
if($flag == 1){
#一旦成功就退出,因为只需要一个结果
print "$sub\n";
exit;
}
}
#同长度的子序列都失败, 则长度-1, 以此类推
$length--;
}