比较文件中字母序列的最佳方法是什么？

Question

我有一个文件，里面有很多字母的组合。
这些组合中有些可能是相同的，所以我想把它们逐一进行比较。
我正在做一些类似的事情，但这并不是我想要的结果：

for line in fl:
line = line.split()
for elem in line:
    if '>' in elem:
        pass
    else:
        for el in line:
            if elem == el:
                print elem, el

文件的例子：

>1
GTCGTCGAAGCATGCCGGGCCCGCTTCGTGTTCGCTGATA  
>2
GTCGTCGAAAGAGGTCT-GACCGCTTCGCGCCCGCTGGTA    
>3
GTCGTCGAAAGAGGCTT-GCCCGCCACGCGCCCGCTGATA  
>4
GTCGTCGAAAGAGGCTT-GCCCGCTACGCGCCCCCTGATA  
>5
GTCGTCGAAAGAGGTCT-GACCGCTTCGCGCCCGCTGGTA  
>6
GTCGTCGAAAGAGTCTGACCGCTTCTCGCCCGCTGATACG  
>7
GTCGTCGAAAGAGGTCT-GACCGCTTCTCGCCCGCTGATA

所以我想知道是否有任何组合完全等于1，或者等于2，依此类推。

Answer 1

比较很长的字母序列会非常低效。直接比较这些序列的哈希值会更快。Python 提供了两种内置的数据类型使用哈希：set 和 dict。在这里最好使用 dict，因为我们可以存储所有匹配的行号。

我假设文件的奇数行是标识符，偶数行是要匹配的序列，所以如果我们按行分割文件文本，就可以将一行作为 ID，下一行作为要匹配的序列。

接着我们使用一个 dict，将序列作为键。对应的值是一个包含所有拥有这个序列的 ID 的列表。通过使用 defaultdict from collections，我们可以轻松处理序列不在 dict 中的情况；如果这个键之前没有被使用过，defaultdict 会自动为我们创建一个值，在这种情况下是一个空的 list。

所以，当我们完成文件的处理后，dict 的值实际上会是一个包含多个列表的 list，每个条目包含共享同一序列的 ID。然后我们可以使用列表推导式来提取有趣的值，也就是那些由多个 ID 共享的序列。

from collections import defaultdict
lines = filetext.split("\n")
sequences = defaultdict(list)

while (lines):
    id = lines.pop(0)
    data = lines.pop(0)
    sequences[data].append(id)

results = [match for match in sequences.values() if len(match) > 1]
print results

Answer 2

下面这个脚本会计算序列的数量。它会返回一个字典，字典里的每个键是独特的序列，而对应的值是这些序列出现的行号（每行的第一部分）。

#!/usr/bin/python
import sys
from collections import defaultdict

def count_sequences(filename):
    result = defaultdict(list)
    with open(filename) as f:
        for index, line in enumerate(f):        
            sequence = line.replace('\n', '')
            line_number = index + 1
            result[sequence].append(line_number)
    return result

if __name__ == '__main__':
    filename = sys.argv[1]
    for sequence, occurrences in count_sequences(filename).iteritems():
        print "%s: %s, found in %s" % (sequence, len(occurrences), occurrences)

示例输出：

etc@etc:~$ python ./fasta.py /path/to/my/file
GTCGTCGAAAGAGGCTT-GCCCGCTACGCGCCCCCTGATA: 1, found in ['4']
GTCGTCGAAAGAGGCTT-GCCCGCCACGCGCCCGCTGATA: 1, found in ['3']
GTCGTCGAAAGAGGTCT-GACCGCTTCGCGCCCGCTGGTA: 2, found in ['2', '5']
GTCGTCGAAAGAGGTCT-GACCGCTTCTCGCCCGCTGATA: 1, found in ['7']
GTCGTCGAAGCATGCCGGGCCCGCTTCGTGTTCGCTGATA: 1, found in ['1']
GTCGTCGAAAGAGTCTGACCGCTTCTCGCCCGCTGATACG: 1, found in ['6']

更新

代码改成使用了 defaultdict 和 for 循环。感谢 @KennyTM。

更新 2

代码改成使用 append 而不是 +。感谢 @Dave Webb。

Answer 3

如果你的目标只是把相似的序列放在一起，那么简单地对数据进行排序就可以了。这里有一个解决方案，它使用了BioPython来解析输入的FASTA文件，先对序列进行排序，然后利用Python的标准库中的itertools.groupby函数来合并相同序列的ID，最后输出一个新的FASTA文件：

from itertools import groupby
from Bio       import SeqIO

records = list(SeqIO.parse(file('spoo.fa'),'fasta'))

def seq_getter(s): return str(s.seq)
records.sort(key=seq_getter)

for seq,equal in groupby(records, seq_getter):
  ids = ','.join(s.id for s in equal)
  print '>%s' % ids
  print seq

输出结果：

>3
GTCGTCGAAAGAGGCTT-GCCCGCCACGCGCCCGCTGATA
>4
GTCGTCGAAAGAGGCTT-GCCCGCTACGCGCCCCCTGATA
>2,5
GTCGTCGAAAGAGGTCT-GACCGCTTCGCGCCCGCTGGTA
>7
GTCGTCGAAAGAGGTCT-GACCGCTTCTCGCCCGCTGATA
>6
GTCGTCGAAAGAGTCTGACCGCTTCTCGCCCGCTGATACG
>1
GTCGTCGAAGCATGCCGGGCCCGCTTCGTGTTCGCTGATA