搜索
123
返回列表 发新帖
楼主: Panda姐

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

[复制链接]

0

主题

2

帖子

29

积分

新手上路

Rank: 1

积分
29
发表于 2018-10-21 00:13:30 | 显示全部楼层
刚才的程序代码格式有点问题,重新发一个:
[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)
回复 支持 反对

使用道具 举报

1

主题

4

帖子

171

积分

注册会员

Rank: 2

积分
171
发表于 2019-4-30 11:56:21 | 显示全部楼层
本帖最后由 Aerthbun 于 2019-4-30 11:59 编辑



[AppleScript] 纯文本查看 复制代码
############################
###找出两条序列的最长公共序列###
############################
print("###只适合1w以内的数据###\n")
from random import randint

def mkseq(len):
	#生成单条序列
    base = ['A','T','C','G']
    lseq = [base[randint(0,3)] for i in range(len)]
    seq = ''.join(lseq)
    return seq
   
def maxnum(list,seq_1,seq_2):
	#二维列表最大值 
	one=list[0][0]
	for i in range(0,len(seq_1)):
		for j in range(0,len(seq_2)):
			if list[j][i] > one:
				one = list[j][i]
	print("最长公共序列长度:\n  "+str(one))
	return one

def lscore(seq_1,seq_2):
	#最长公共序列打分列表
	lscore_list=[[0 for i in range(len(seq_1))]for i in range(len(seq_2))]
	for i in range(1,len(seq_1)):
		if seq_1[i]==seq_2[0]:
			lscore_list[0][i]=lscore_list[0][i]+1
	for i in range(1,len(seq_2)):
		if seq_2[i]==seq_1[0]:
			lscore_list[i][0]=lscore_list[i][0]+1
	for i in range(1,len(seq_1)):
		for j in range(1,len(seq_2)):
			if seq_1[i]==seq_2[j]:
				lscore_list[j][i]=lscore_list[j-1][i-1]+1
	return lscore_list	

def seqprint(list,seq_1,seq_2,max):
	print("最长公共序列:")
	for i in range(0,len(seq_1)):
		for j in range(0,len(seq_2)):
			if list[j][i]==max:
				print("  "+str(seq_1[j-1-max:j-1]))
	
print("###开始生成序列###")
seq1=mkseq(5000)
seq2=mkseq(5000)
print("###开始打分###")
ls=lscore(seq1,seq2)
print("###开始计算最大长度###")
maxscore=maxnum(ls,seq1,seq2)
print("###开始生成最长序列###")
seqprint(ls,seq1,seq2,maxscore)



回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-20 21:05 , Processed in 0.027914 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.