比较文件内部字母序列的最佳方法？

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

如果目标是简单地组样序列一起，然后简单地排序的数据就可以了。下面是一个使用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 

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'] 

更新

更改代码，使用dafaultdict和for循环。谢谢@KennyTM。

更新2

更改代码使用append而非+。谢谢@Dave Webb。

Answer 3

一般来说，对于这种类型的工作，您可能需要调查Biopython，它具有许多解析和处理序列的功能。

但是，您可以使用字典来解决您的特定问题，这是Manoj向您提供的一个示例。

Answer 4

比较长的字母序列将是相当低效的。比较序列的散列会更快。 Python提供了两种使用散列的内置数据类型：set和dict。这里最好使用dict，因为我们可以存储所有匹配的行号。

我认为该文件对备用线标识和标签，所以如果我们分裂的新行文件的文本，我们可以把一个行的id和未来的序列相匹配。

然后我们使用一个dict，序列作为关键字。相应的值是具有该序列的ID列表。通过使用defaultdict from collections，我们可以轻松处理不在dict中的序列的情况;如果之前没有使用密钥defaultdict会自动为我们创建一个值，在这种情况下为空list。

因此，当我们完成文件的工作时，dict的值将实际上是list的list s，每个条目包含共享序列的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