比較文件內部字母序列的最佳方法？

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