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

【问题标题】：How to count amino acids in fasta formated file?如何计算fasta格式文件中的氨基酸？
【发布时间】：2013-01-18 10:59:34
【问题描述】：

我找到了解析 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...

【问题讨论】：

您的代码只是将> 行和序列行放在两个列表中。要计算序列行中的字符，您需要collections.Counter(line)。
“我不想使用 BioPython”——为什么？如果您只是想学习 Python，那很好，否则有充分的理由使用现有的库。
你也可以问biostars：biostar.stackexchange.com

标签： python bioinformatics biopython fasta

【解决方案1】：

作为 katrielalex points out，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})

【讨论】：

是否可以计算这个序列中的序列，我的意思是：以我的例子：MRMRGRRLLPIIL 我想要这样的输出：M - 2, R - 2, RR - 2, G - 1, LL - 2, P - 1, I - 2, L - 1 ?
@mazix 这更棘手。您可以使用itertools.groupby 来拆分字符串：for aa, s in groupby('MRMRGRRLLPIILXKKK'): print ''.join(s) 会给您一个想法。但是对于重复的LL，您必须将L 的计数器增加两个，等等。

【解决方案2】：

这可以使用非常简单的bash命令来获得，我的答案如下

cat input.fasta #my input file
>gi|7290019|gb|AAF45486.1| (AE003417) EG:BACR37P7.1 gene product [Drosophila melanogaster]
    MRMRGRRLLPIIL

cat input.fasta | grep -v ">" | fold -w1 | sort | uniq -c

输出：

fold -w1 在每个字符处拆分，您对它们进行排序并计算唯一的字符

【讨论】：

【解决方案3】：

cmets 中提到的 katrielalex 的替代方法是使用另一个字典，代码如下

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

  order = []
  sequences = {}
  counts = {}

  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('*')
      for aa in sequences[name]:
        if aa in counts:
            counts[aa] = counts[aa] + 1
        else:
            counts[aa] = 1  


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

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

这个输出：

1 sequences found
{'G': 1, 'I': 2, 'M': 2, 'L': 3, 'P': 1, 'R': 4}

【讨论】：

是否可以计算这个序列中的序列，我的意思是：以我的例子：MRMRGRRLLPIILXKKK 我想要这样的输出：M - 2, R - 2, RR - 2, G - 1, LL - 2, P - 1, I - 2, L - 1, X - 1, KKK - 1 ?

【解决方案4】：

# your previous code here

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

import collections

for order in x:
  print order, ':',
  print '\n'.join('%s - %d' % (k, v) for k, v in collections.Counter(y[order]).iteritems())
  print

【讨论】：