【Rosalind】Finding a Shared Motif - all()函数的优雅使用

Rosalind习题Finding a Shared Motif

今天完成了Rosalind上这道题目，自己写的代码运行了55.5s，虽然成功得到了正确答案，但是发现讨论区的大佬们从之前的讨论答案怎么不对变成了讨论怎么能让代码更快。
我在学习了一下评论区大家的解法之后，对我的代码稍加修改，在1.4s内就跑出了和之前相同的答案，可见一个逻辑、一个运算或者是一行优雅代码的实现，就能让工作效率倍增。

---

Problem

A common substring of a collection of strings is a substring of every member of the collection. We say that a common substring is a longest common substring if there does not exist a longer common substring. For example, “CG” is a common substring of “ACGTACGT” and “AACCGTATA”, but it is not as long as possible; in this case, “CGTA” is a longest common substring of “ACGTACGT” and “AACCGTATA”.

Note that the longest common substring is not necessarily unique; for a simple example, “AA” and “CC” are both longest common substrings of “AACC” and “CCAA”.

Given: A collection of k (k≤100) DNA strings of length at most 1 kbp each in FASTA format.

Return: A longest common substring of the collection. (If multiple solutions exist, you may return any single solution.)

---

以上就是这道题的题目，大体意思就是给你一个fasta文件，然后要求找到一条最长的子序列，这条子序列出现在fasta中的每一条序列中。

思路

首先，找到最长的共有子序列我肯定想到的是从所有原始序列中最短的那一条入手，子序列的最长长度就是最短序列的长度。为了找到最短序列，我选择先对fasta文件按照序列的长短进行排序，即下面这段代码：

def sort_fasta(fasta_file):
    with open(fasta_file) as file:
        data = file.read().splitlines()
        fasta = {}
        for line in data:
            if line.startswith('>'):
                name = line
                fasta[name] = ''
            else:
                fasta[name] = fasta[name] + line
        fasta = sorted(fasta.items(), key=lambda i:len(i[1]))
        #fasta is a list made up of tuples
        return fasta

这样，fasta[0][1]就是那条最短的序列。
假设这条最短序列的长度是k，我现在想知道长度为k的子序列是否是所有fasta的序列string的substring，如果不是，我就看长度为k-1的子序列（有2条），然后如果还没有，就看k-2、k-3……的子序列，一旦找到了子序列，那么它就是我们要寻找的最长子序列。

代码实现

结合以上思路，我一开始写出了这样的代码。

def shared_motif(fasta_file):
    fasta = sort_fasta(fasta_file)
    #print(fasta)
    fasta_sequence = []
    for i in fasta:
        fasta_sequence.append(fasta[fasta.index(i)][1])
    #print(fasta_sequence)
    max_try_len = len(fasta[0][1])
    shortest = fasta[0][1]
    #print(shortest[0:10])
    #list all k-mers of the shortest sequence, then test if they are in all other sequences
    kmer_list = []
    for length in reversed(range(max_try_len+1)):
        if length == 0:
            break
        kmer_sum = max_try_len - length +1
        #print(kmer_sum)
        for i in range(kmer_sum):
            kmer_list.append(shortest[i:i+length])
            #print(kmer_list)
    #here, kmer_list should contain all the kmer from the max-mer to the 1-mer
    #print(kmer_list)


    #for loop
    for kmer in kmer_list:
        checkboard = 0
        for check in fasta_sequence:
            if kmer in check:
                checkboard = checkboard + 1
                if checkboard == len(fasta):
                    motif = kmer
                    return motif
                    break

11-18行，我提取了最短序列中的所有k-mer，比如最短序列长度为10，就从10-mer，9-mer一直提取到1-mer，并把它们全部放在一个列表里。
接着我开始遍历这个k-mer列表，从最大的k-mer开始，检查该k-mer是否存在于每个原始序列中。那么这个过程我设计了一个得分变量，每次检查都重设得分变量checkboard为0，然后只要比对上一个就让他+1，这样如果这个得分到了原始序列列表的总长度的话，说明该k-mer存在于每个原始序列中，就返回这个k-mer。

事实证明，这段代码运行没有问题，可以正确得到结果。

Better Solution

是的，就像开头说的，我这段代码运行速度很慢。我想了一下，应该是在最后一个for loop里，其实很多时候得分变量没有加分的话该循环就应该退出了，进入对下一个k-mer的检查，这个步骤应该浪费了很多计算资源。

那么有没有更好的实现检查该k-mer是否存在于每个原始序列中这一目的的方法呢？答案当然是有，用all()。

all()的参数是可迭代参数，只要可迭代参数中所有元素全为True，它就返回True，否则它就返回False。
有了这个方法，就很容易实现上述目的了~

    #all()
    for kmer in kmer_list:
        if all(kmer in check for check in fasta_sequence):
            motif = kmer
            return motif
            break

用这段代码替代刚才那段for循环，对相同的测试数据进行运算，运算时间仅需1.4s。

不过这道题目的评论区也有人提到用all()方法的运行速度可能没有传统循环方法快，有大神是这么解释的：

@Pathfinder: I believe calling all() first involves computing the boolean values for all the elements in the list. This is much slower than calling break if the substring is not in any one DNA string. Goes to show that elegant-looking code isn't always fast code.

所以，看来代码的优雅程度和它的计算效率也没有绝对的关系呀。

本文于2020年9月29日发布在https://emmettoncloud.cn/archives/44（本人博客上），现搬至CSDN博客。