搜索
查看: 5399|回复: 20

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

[复制链接]

58

主题

103

帖子

754

积分

版主

Rank: 7Rank: 7Rank: 7

积分
754
QQ
发表于 2016-8-28 15:20:09 | 显示全部楼层 |阅读模式
本帖最后由 Panda姐 于 2016-8-30 19:01 编辑

以下2条序列:
序列一:AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC
序列二:CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA
编写脚本,找出这两条序列中最长的共有序列。

代码记录:
[Perl] 纯文本查看 复制代码
#! usr/bin/perl
use strict;
my($s1,$s2,$max,$length,$star,$common,$compare,@compare,$common,@common,$maxcom,%len,$len);
open OUTPUT,">result.txt";
$s1="AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC";
$s2="CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA";
$max=length($s2);
for($star=0;$star<=$max;$star++)
{
   for($length=2;$length<=$max-$star;$length++)
   {
     push @compare,substr($s2,$star,$length);
       
    }
}
sort @compare;
foreach $compare(@compare){
     chomp $compare;
     if($s1=~/$compare/){
     push @common,$compare;
}
}
foreach $common(@common){
$len{$common}=length($common);
}
@common=sort {$len{$b} <=>$len{$a}} keys %len;
$maxcom=shift@common;
print "$maxcom\t$len{$maxcom}\n";


师兄给的答案:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use warnings;
use strict;
my $str1 = "AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC";
my $str2 = "CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA";
my $str = $str1 . "\n" . $str2;
my @substr=sort{length($b) <=> length($a)} map{substr($str,$_)=~/^(.+)(?=.*\n.*\1)/} 0..length($str1)-1;#此>    处用到perl扩展模式匹配(?=肯定的预见匹配语法)
my @result = grep { length == length $substr[0] } @substr;
print "@result\n";








上一篇:【Panda姐-perl练习题1】统计文本中的每个单词出现次数并排序
下一篇:【Panda姐-perl练习题3】从gbk文件中提取蛋白序列
回复

使用道具 举报

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

使用道具 举报

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

使用道具 举报

4

主题

29

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2016-9-8 11:04:19 | 显示全部楼层
上边这段代码是怎么执行的啊
回复 支持 反对

使用道具 举报

7

主题

32

帖子

228

积分

版主

Rank: 7Rank: 7Rank: 7

积分
228
发表于 2016-9-16 17:23:18 | 显示全部楼层
姐姐,能不能注释一下,刚开始学,看不懂
“Nothing in Biology Makes Sense Except in the Light of Evolution”  --Theodosius Dobzhansky
回复 支持 反对

使用道具 举报

0

主题

14

帖子

225

积分

中级会员

Rank: 3Rank: 3

积分
225
发表于 2016-9-24 22:40:15 | 显示全部楼层
本帖最后由 0号菜鸟 于 2016-9-26 09:18 编辑

谢谢楼下的提醒,已修改~
分享一个python版本的,个人觉着里面有bug,但是暂时还没有发现,欢迎相互交流~
[Python] 纯文本查看 复制代码
a = 'AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGCCGTAGCTACGAT'
b = 'CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA'
match_max = '0'
match_list = []
for x in range(len(a)):
    for y in range(len(a)):
        if a[x:y+1] in b and y > x: #添加条件‘y>x’为了防止字符逆序
            match = a[x:y+1]
            if len(match) >= len(match_max):
                match_list.append(match)
                match_max = match
i = 0
match_list = match_list[::-1]
for seq in match_list:
    i+=1
    if len(seq) > len(match_list[i]):
        print match_list[:i]
        break
回复 支持 反对

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-9-25 10:33:59 | 显示全部楼层
0号菜鸟 发表于 2016-9-24 22:40
分享一个python版本的,个人觉着里面有bug,但是暂时还没有发现,欢迎相互交流~
[mw_shl_code=python,true] ...

您的这个程序写的简介易懂
如果存在两个或者多个长度一样的最长共有序列,您的这个程序只能输出其中一个最长共有序列(ps: 在这个题目没有出现这样的情况,只存在一个最长共有序列)。
我也是刚开始学python,不知道说的对不对,希望大家一起交流共同学习
回复 支持 反对

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-9-25 10:48:08 | 显示全部楼层
本帖最后由 learnyoung 于 2016-9-25 16:24 编辑

看了@0号菜鸟的脚本,自己试着写了一个,刚开始学python,程序写的一般
[Python] 纯文本查看 复制代码
# coding=utf-8
a = 'AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC'
b = 'CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA'
match =[]
for x in range(len(b)):
    for y in range(len(b)):
        if x<y and b[x:y+1] in a:#b[x:y+1]为b中所有存在的元素(以元素长度作为判断,元素长度最大为b的长度,最小为2)
            match.append(b[x:y+1])
setmatch=set(match)
#找出所有的可能存在最长共有序列(视a,b情况而定,可能存在多个共有最长序列)
t=[]
for i in setmatch:
    t.append(len(i))
for j in setmatch:
    if len(j)==max(t):
        print j
回复 支持 反对

使用道具 举报

58

主题

103

帖子

754

积分

版主

Rank: 7Rank: 7Rank: 7

积分
754
QQ
 楼主| 发表于 2016-9-25 12:41:54 | 显示全部楼层
learnyoung 发表于 2016-9-25 10:48
看了@0号菜鸟的脚本,自己试着写了一个,刚开始学python,程序写的一般
# coding=utf- ...

你把代码高亮一下,使用<>。
回复 支持 反对

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-9-25 16:16:03 | 显示全部楼层
本帖最后由 learnyoung 于 2016-9-25 16:23 编辑

[Python] 纯文本查看 复制代码
# coding=utf-8
a = 'AACTGCACGTGCATCGGATGCATCGATCGTGCGAGTAGTCGATCGATCGTAGCTAGCTCAGTCGATCAGCTACCTCGC'
b= 'CGTAGCTACGATCGCATCAGCTACCTCCGTTCGAGTCTGCGCAACGCTACGACTATCGACGTCA'
match =[]
for x in range(len(b)):
    for y in range(len(b)):
        if x<y and b[x:y+1] in a:#b[x:y+1]为b中所有存在的元素(以元素长度作为判断,元素长度最大为b的长度,最小为2)
            match.append(b[x:y+1])
setmatch=set(match)
#找出所有的可能存在最长共有序列(视a,b情况而定,可能存在多个共有最长序列)
t=[]
for i in setmatch:
    t.append(len(i))
for j in setmatch:
    if len(j)==max(t):
        print j
回复 支持 反对

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-9-25 16:25:50 | 显示全部楼层
Panda姐 发表于 2016-9-25 12:41
你把代码高亮一下,使用。

好的,版主,要达到什么级数才能取消验证码?
回复 支持 反对

使用道具 举报

0

主题

14

帖子

225

积分

中级会员

Rank: 3Rank: 3

积分
225
发表于 2016-9-25 21:12:44 | 显示全部楼层
learnyoung 发表于 2016-9-25 10:33
您的这个程序写的简介易懂
如果存在两个或者多个长度一样的最长共有序列,您的这个程序只能输出其中一个 ...

对的,我没有考虑到这个问题~
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2018-11-16 22:49 , Processed in 0.044704 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.