如何计算fasta格式文件中的氨基酸？

Question

我找到了解析fasta frmated文件的代码.我需要计算每个序列中有多少A,T,G等,例如:

>gi|7290019|gb|AAF45486.1| (AE003417) EG:BACR37P7.1 gene product [Drosophila melanogaster]
MRMRGRRLLPIIL 

在这个序列中:

M - 2
R - 4
G - 1
L - 3
I - 2
P - 1

代码很简单:

 def FASTA(filename):
  try:
    f = file(filename)
  except IOError:                     
    print "The file, %s, does not exist" % filename
    return

  order = []
  sequences = {}

  for line in f:
    if line.startswith('>'):
      name = line[1:].rstrip('\n')
      name = name.replace('_', ' ')
      order.append(name)
      sequences[name] = ''
    else:
      sequences[name] += line.rstrip('\n').rstrip('*')

  print "%d sequences found" % len(order)
  return order, sequences

x, y = FASTA("drosoph_b.fasta")

但我怎么算这些氨基酸呢？我不想使用BioPython,我想知道如何使用,例如count......

Answer 1

正如katrielalex 指出的那样,collections.Counter非常适合这项任务:

In [1]: from collections import Counter

In [2]: Counter('MRMRGRRLLPIIL')
Out[2]: Counter({'R': 4, 'L': 3, 'M': 2, 'I': 2, 'G': 1, 'P': 1})

您可以将它应用于sequences代码中的dict 值.

但是,我建议不要在现实生活中使用此代码.像BioPython这样的图书馆可以做得更好.例如,您显示的代码将生成相当庞大的数据结构.由于FASTA文件有时非常大,因此可能会耗尽内存.此外,它不会以最佳方式处理可能的异常.

最后,使用库只会节省您的时间.

使用BioPython的示例代码:

In [1]: from Bio import SeqIO

In [2]: from collections import Counter

In [3]: for entry in SeqIO.parse('1.fasta', 'fasta'):
   ...:     print Counter(entry.seq)
   ...:     
Counter({'R': 4, 'L': 3, 'I': 2, 'M': 2, 'G': 1, 'P': 1})