vlambda博客
学习文章列表

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--;}