代码之家 › 专栏 › 技术社区 › psoares

比较文件中字母顺序的最佳方法?

python

psoares · 技术社区 · 14 年前

我有一个文件,里面有很多字母序列。
其中一些序列可能是相等的,所以我想将它们全部与所有序列进行比较。
我在做这样的事,但这并不是我想要的:

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,等等。

4 回复 | 直到 14 年前

Kevin Jacobs 14 年前

如果目标是简单地将类似序列分组在一起,那么简单地对数据进行排序就可以做到这一点。这里有一个解决方案 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

Mark Thomas 14 年前

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

然而,你的特殊问题可以通过口述来解决,马诺吉举了一个例子。

David Webb 14 年前

比较长的字母序列会非常低效。比较序列的散列会更快。python提供了两种使用哈希的内置数据类型: set 和 dict . 最好用一下 双关语 在这里,我们可以存储所有匹配项的行号。

我假设文件在替换行上有标识符和标签,所以如果我们将文件文本拆分为新行,我们可以将一行作为ID,下一行作为要匹配的序列。

然后我们使用 双关语 以序列为键。相应的值是具有此序列的ID列表。通过使用 defaultdict from collections 我们可以很容易地处理序列不在 双关语 ;如果以前没有使用过钥匙 拖欠债务 将自动为我们创建一个值,在这种情况下为空 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

Community CDub 8 年前

下面的脚本将返回序列计数。它返回一个字典,其中单独的、不同的序列作为键,数字(每行的第一部分)出现在这些序列中。

#!/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 .