搜索
楼主: Panda姐

【Panda姐-perl练习题2】找出两条序列的最长共有序列

[复制链接]

4

主题

51

帖子

327

积分

中级会员

Rank: 3Rank: 3

积分
327
发表于 2016-9-27 21:04:13 | 显示全部楼层
提供一种常规思路,考虑利用start位置和length的变化来做嵌套循环,然后选取length最大的公共字符串输出。
rosalind有一个寻找多条序列(非相同序列)最长的公有序列(Finding a Shared Motif),相当于Panda姐这一题的扩展,链接:http://rosalind.info/problems/lcsm/,大家可以看一下。
——————
以下为搬运内容:
那个网站的参考答案(python)时间复杂度O(log(len)∗len∗len∗k),arr是一个存放各个序列的列表,具体如下:
[Python] 纯文本查看 复制代码
def substr_in_all(arr, part):
  for dna in arr:
    if dna.find(part)==-1:
      return False
  return True

def common_substr(arr, l):
  first = arr[0]
  for i in range(len(first)-l+1):
    part = first[i:i+l]
    if substr_in_all(arr, part):
      return part
  return ""

def longest_common_substr(arr):
  l = 0; r = len(arr[0])

  while l+1<r:
    mid = (l+r) / 2
    if common_substr(arr, mid)!="":
      l=mid
    else:
      r=mid

  return common_substr(arr, l)

本人本科阶段Perl应该不会学,不过看到了有人提供了Perl的solution,顺便就搬运过来了:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;

sub kmp {
    my @T = (-1, 0);
    my ($pos, $cnd) = (2, 0);
    while ( $pos <= @_ ) {
        if ($_[$pos-1] eq $_[$cnd]) { $T[$pos++] = ++$cnd }
        elsif ($cnd > 0) { $cnd = $T[$cnd] // 0 }
        else { $T[$pos++] = 0 }
    }
    shift @T;
    return @T;
}

use List::Util qw(max min);
sub lcs {
    my @first = split //, my $first = shift;
    my $n = @first;
    my @string = map [ split //, $_ ], @_;
    my ($pos, $length) = (0, 0);
    SUFFIX: for my $p ( 0 .. $n - 1 ) {
        shift @first if $p;
        my @max;
        for (@string) {
            my @kmp = kmp @first, '$', @$_;
            my $max = max @kmp;
            next SUFFIX if $max < $length;
            push @max, $max;
        }
        if ((my $min = min @max) > $length) {
            ($pos, $length) = ($p, $min);
        }
    }
    return substr $first, $pos, $length;
}

open my $f, '< rosalind_lcs.txt' or die "could not open: $!";
say lcs <$f>;
回复 支持 1 反对 0

使用道具 举报

0

主题

9

帖子

97

积分

注册会员

Rank: 2

积分
97
发表于 2016-10-7 23:30:35 | 显示全部楼层
本帖最后由 musicyiyi 于 2016-10-8 05:54 编辑

[Perl] 纯文本查看 复制代码
#! usr/bin/perl
use warnings;
use strict;


my $s1="AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGCATGAGCTACCTC";
my $s2="CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCAATGAGCTACCTC";

my $size2=length($s2);

my $len=$size2;
my $getmatch=0;
my $offset=0;

LABEL:for ($len=$size2; $len>0; $len--){
              for($offset=0; $offset<=$size2-$len; $offset++){
                   if($s1 =~ substr($s2,$offset,$len)){last LABEL;}
             }
          }

print "the length of longest common substring is $len\n";
print "one of the longest common subsring is: ".substr($s2,$offset,$len)."\n";


my %common=();
$common{substr($s2,$offset,$len)}=1; #all=1

for(my $newoffset=$offset+$len; $newoffset<=$size2-$len; $newoffset++){
         if($s1 =~ substr($s2,$newoffset,$len)){
             if (!exists $common{substr($s2,$newoffset,$len)} ){                  
                 print substr($s2,$newoffset,$len)."\n";
                 $common{substr($s2,$newoffset,$len)}=1;
             }
         }
     }
     
回复 支持 反对

使用道具 举报

11

主题

53

帖子

236

积分

中级会员

Rank: 3Rank: 3

积分
236
QQ
发表于 2016-10-8 18:25:12 | 显示全部楼层
厉害,楼主师兄更厉害!赞一个!
苛求远离完美
回复 支持 反对

使用道具 举报

2

主题

9

帖子

233

积分

中级会员

Rank: 3Rank: 3

积分
233
发表于 2016-12-29 07:01:45 | 显示全部楼层
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/env python
a='AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC';
b='CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA';
match_list= [];
number_list=[];
for i in range(len(b)) :   
	for j in range(i):
		l=b[j:len(b)-i+1+j];
		if a.find(l)>0:
			match_list.append(l);
			number_list.append(len(l));
m=max(number_list);
for i in match_list:
	if len(i)==m:
		print i;
回复 支持 反对

使用道具 举报

0

主题

2

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-4-8 17:44:08 | 显示全部楼层
本帖最后由 守夜人 于 2017-4-8 17:47 编辑

[AppleScript] 纯文本查看 复制代码
  1 #!/usr/bin/env python
  2 # encoding: utf-8
  3      
  4 def quest_max_com_seq(seq1,seq2):
  5     com_seq_list = []
  6     for i in range(len(seq1)):
  7         win_size=i
  8         for n in range(len(seq1)-win_size):
  9             if seq1[n:n+win_size] in seq2:
 10                 com_seq_list.append(seq1[n:n+win_size])
 11     print(com_seq_list)
 12     max_seq = com_seq_list[-1] #因为列表是有序的,最后添加入列表的元素是最大共有序列
 13     print("the max common seq is:%s"%max_seq)                           
 14      
 15 if __name__=='__main__':
 16     seq1='AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGT    CGATCAGCTACCTCGC'
 17     seq2='CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGT    CA'
 18     quest_max_com_seq(seq1,seq2)
回复 支持 反对

使用道具 举报

0

主题

6

帖子

83

积分

注册会员

Rank: 2

积分
83
发表于 2017-9-28 16:48:49 | 显示全部楼层
这道题感觉好难啊,不会
回复 支持 反对

使用道具 举报

5

主题

10

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2018-1-30 12:18:18 | 显示全部楼层
对于数量不多的可以,如果是十万行规模的,又什么方法可以解决?
回复 支持 反对

使用道具 举报

0

主题

3

帖子

29

积分

新手上路

Rank: 1

积分
29
发表于 2018-4-9 15:56:44 | 显示全部楼层
本帖最后由 立秋之后 于 2018-4-9 16:00 编辑

写了几种方案,方案1和2好像有问题,方案3应该是比较好的一个方案,十万级的数据应该问题也不大.
[Python] 纯文本查看 复制代码
#-*- coding:utf8 -*-
from random import randint
import time

def mkSeqs(num,minLen,maxLen):
    #随机生成序列
    seqs = [mkseq(randint(minLen,maxLen)) for i in range(num)]
    return  seqs  

def mkseq(len):
    #生成单条序列
    base = ['A','T','C','G']
    lseq = [base[randint(0,3)] for i in range(len)]
    seq = ''.join(lseq)
    return seq
    
def rndgoal():
    #随机生成位点分值
    goalsDict = {}
    base = ['A','T','C','G']
    goals = [i for i in range(16)]
    for bef in base:
        for aft in base:
            rnd  = randint(0,len(goals)-1)
            goalsDict[str(bef)+str(aft)] = goals[rnd] ** 3
            del goals[rnd]
    return goalsDict
    
def sortSeq(seqs,glist):
    #计算序列的分值并排序
    seqsscore = []
    for seq in seqs:
        seqScore = 0
        for i in range(len(seq) - 1):
            dbbase = seq[i:i+2]
            seqScore += glist[dbbase]
        seqsscore.append((seqScore/len(seq),seq))
    return(sorted(seqsscore,key = lambda seqsscore : seqsscore[0]))
    
def seqSame(seq1,seq2):
    #两条序列对比,返回最大相似序列的长度和序列
    #确保seq1大于seq2
    len1 = len(seq1)
    len2 = len(seq2)
    if len1<len2:
        seq1, seq2 = seq2, seq1
        len1, len2 = len2, len1
    maxsame = 0
    seq = ''
    for i in range(len2):
        samebase = 0
        tmpseq = ''
        for t in range(i,len2):
            if seq2[i:t] in seq1:
                if maxsame < t - i:
                    maxsame = t - i
                    seq = seq2[i:t]
            
    return maxsame,seq

def borderMatch(sortseq):
    #比对相邻的序列
    maxlen = 0
    maxseq = 0
    maxbef = ''
    maxnex = ''
    for n in range(len(sortseq)):
        for i in range(1,3):
            try:
                bef = sortseq[n][1]
                nex = sortseq[n + 1][1]
                samelen,sameseq = seqSame(bef,nex)
                if maxlen < samelen:
                    maxlen = samelen
                    maxseq = sameseq
                    maxbef = bef
                    maxnex = nex
            except:
                continue
    return maxlen,maxseq,maxbef,maxnex

def countLen(seqs):
    lenth = 0
    for seq in seqs:
        lenth += len(seq)
    return lenth 

def countMinSeq(len):
    #在一定 的总长度下的最短片段长度
    countLen = 1
    minlen = 1
    maxlen = 4
    while True:
        if minlen <= len < maxlen:
            return countLen
        else:
            minlen *= 4
            maxlen *= 4
            countLen += 1

def mkdict(shortSeq):
    #根据最短序列生成字典
    seqDict = {}
    seqList = ['A','T','C','G']
    for seq in seqList:
        if len(seq) < shortSeq:
            seqList.append(seq+'A')
            seqList.append(seq+'T')
            seqList.append(seq+'C')
            seqList.append(seq+'G')
        else:
            seqDict[seq] = []
    return seqDict

def matchDict(seq,index,shortlen,seqDict):
    for i in range(len(seq) - shortlen +1):
        short = seq[i:i+ shortlen]
        seqDict[short].append((index,i))

def mkNewDict(oldSeqDict,seqs,shortlen):
    SeqDict = {}
    for s in oldSeqDict.keys():
        if  len(oldSeqDict[s])<2:
            continue
        else:
            for i in oldSeqDict[s]:
                seqId = i[0]
                index = i[1]
                newseq = seqs[seqId][index:index+shortlen+1]
                if len(newseq) < shortlen + 1:
                    continue
                try:
                    SeqDict[newseq]
                except:
                    SeqDict[newseq] = []
                SeqDict[newseq].append((seqId,index))
    return SeqDict
            
            
def path1(seqs):
    '''
    给相邻的碱基不同的分数,排序后只对相邻的两个进行判定
    '''
    sta = time.time()
    glist = rndgoal()
    sortseq = sortSeq(seqs,glist)
    maxlen,maxseq,bef,nex = borderMatch(sortseq)
    print(time.time() - sta)
    print(maxlen,maxseq,bef,nex)

def path2(seqs):
    '''
    直接遍历所有的序列取最大相似
    '''
    sta = time.time()
    num = len(seqs)
    maxlen = 0
    maxseq = ''
    bef = ''
    nex = ''
    for i in range(num):
        for t in range(i + 1,num):
            samelen,sameseq = seqSame(seqs[i],seqs[t]) 
        if maxlen < samelen:
            maxlen = samelen
            maxseq = sameseq
            bef = seqs[i]
            nex = seqs[t]
    print(time.time() - sta)
    print (maxlen,maxseq,bef,nex)
        
def path3(seqs):
    '''
    建立搜索所引来对比
    '''
    sta = time.time()
    seqslen = countLen(seqs)
    shortSeq = countMinSeq(seqslen)#可能出现最多重复的最短序列长
    seqDict = mkdict(shortSeq)
    for i in range(len(seqs)):
        matchDict(seqs[i],i,shortSeq,seqDict)
    shortSeq += 1
    while True: 
        seqDict = mkNewDict(seqDict,seqs,shortSeq)
        shortSeq += 1
        matchFinish = True
        for i in seqDict.keys():
            if len(seqDict[i]) > 1:
                record = seqDict
                matchFinish = False
                break
        if matchFinish:
            break
    print(time.time() - sta)
    for i in record.keys():
        if len(record[i])> 1:
            print(record[i],i)

if __name__ == '__main__':
    seqs = mkSeqs(200,300,400)
    path1(seqs)
    path2(seqs)
    path3(seqs)    



回复 支持 1 反对 0

使用道具 举报

7

主题

13

帖子

172

积分

注册会员

Rank: 2

积分
172
发表于 2018-5-31 10:03:21 | 显示全部楼层
[Python] 纯文本查看 复制代码
a = 'AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGCCGTAGCTACGAT'
b = 'CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA'
c=[]
d={}
for i in range(len(a)):
    for j in range(i+1,len(a)):
        if a[i:j+1] in b:
            c.append(a[i:j+1])
for i in c:
    d[i]=len(i)
e=sorted(zip(d.values(),d.keys()))
for i in e:
    if i[0]==e[-1][0]:
        print(i)
回复 支持 反对

使用道具 举报

0

主题

2

帖子

29

积分

新手上路

Rank: 1

积分
29
发表于 2018-10-21 00:09:53 | 显示全部楼层
xuehzh95 发表于 2016-9-27 21:04
提供一种常规思路,考虑利用start位置和length的变化来做嵌套循环,然后选取length最大的公共字符串输出。
...

你搬过来的程序确实不错, 效率比较高,但有一点小瑕疵:
1)对于str1 = AT, str2 = CAT这样的简单问题,显然最长共有序列是AT,但程序确给出错误答案
2)如果有多条同样长度的最长共有序列,无法全部列出。
针对这两个问题,我做了点小修改,并加了些注释,python程序如下:
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 13 10:41:27 2018
程序用于寻找多个字符串中的最长的共同子串,通过调用longest_common(strList)实现,其中strList是所有序列字符串的列表,返回最长共同子序列
@author: leaf
"""
import numpy as np
import time

def substr_in_all(strList, substr):
    """判断字符串列表strList中所有字符串是否都包含子串substr
    """
    for each in strList:
        if each.find(substr) == -1: return False
    return True
        
def have_substr_of_len_n(strList, n):
    """判断字符串列表strList中是否含有长度为n的共同子串
    """
    shortest_str = strList[0]
    substr_num = len(shortest_str) - n + 1  #长度为n的子串数目
    for i in range(substr_num):
        substr = shortest_str[i : i + n]
        if substr_in_all(strList[1:], substr): return True
    return False

def get_allstr_of_len_n(strList, n):
    """获得字符串列表中,长度为n的所有共同子串
    """
    allStr = []
    shortest_str = strList[0]
    substr_num = len(shortest_str) - n + 1
    for i in range(substr_num):
        substr = shortest_str[i : i + n]
        if substr_in_all(strList[1:], substr): allStr.append(substr)
    return allStr
   

def longest_common(strList):
    """获得字符串列表中的最长共同子串,利用二分法来确定最长子串的长度。
    search_region_left: 为搜索区间左边界,一定含有这个长度的子串
    search_region_right: 为搜索区间右边界,一定不包含这个长度的子串
    区间缩小到宽度为1,则左边界即为最长子串的长度
    """
    short_ind = np.argmin([len(each) for each in strList])
    if short_ind != 0:
        strList[0], strList[short_ind] = strList[short_ind], strList[0] #将最短序列放到strList[0]

    shortest_str_length = len(strList[0])
    if have_substr_of_len_n(strList, shortest_str_length): return strList[0]
    search_region_left, search_region_right = 0, shortest_str_length
    while search_region_left + 1 < search_region_right:
        middle = (search_region_left + search_region_right) // 2
        if have_substr_of_len_n(strList, middle): search_region_left = middle
        else: search_region_right = middle
    if search_region_left == 0: return 'no common substring'
    else: return get_allstr_of_len_n(strList, search_region_left)


str1 = 'AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC'*2000
str2 = 'CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA'*2000
t1 = time.clock()
print(longest_common([str1, str2]))
print('time_needed:', time.clock() - t1)

在我的笔记本上跑了一下,10万级别的两条序列,需要几分钟的样子。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-23 13:16 , Processed in 0.105423 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.