使用Python计数fasta文件中每个基因中特定字符的出现次数

Question

我有一个fasta文件，内容如下：

>SO_0001 
MTKIAILVGTTLGSSEYIADEMQAQLTPLGHEVHTFLHPTLDELKPYPLWILVSSTHGAGDLPDNLQPFC
KELLLNTPDLTQVKFALCAIGDSSYDTFCQGPEKLIEALEYSGAKAVVDKIQIDVQQDPVPEDPALAWLA
QWQDQI
>SO_0002  
MTTPVDAPKWPRQIPYIIASEACERFSFYGMRNILTPFLMTALLLSIPEELRGAVAKDVFHSFVIGVYFF
PLLGGWIADRFFGKYNTILWLSLIYCVGHAFLAIFEHSVQGFYTGLFLIALGSGGIKPLVSSFMGDQFDQ
>SO_0003 
MTTDTIVAQATAPGRGGVGIIRISGDKATNVAMAVLGHLPKPRYADYCYFKSASGQVIDQGIALFFKGPN
SFTGEDVLELQGHGGQIVLDMLIKRVLEVEGIRIAKPGEFSEQAFMNDKLDLTQAEAIADLIDATSEQAA
KSALQSLQGEFSKEVHELVDQVTHLRLYVEAAIDFPDEEVD

在这个文件中，">"后面的是基因ID，而紧接着">"这一行的字母则是对应的序列。我想要逐行读取这个文件，统计每个基因ID对应的序列中有多少个"C"字母。我希望输出的文件格式是用制表符分隔的，像这样：

SO_0001    Number of C's
SO_0002    Number of C's
SO_0003    Number of C's

还有其他类似的内容...

我正在使用Python，原本以为这很简单，可以把基因ID当作字典的键，但我之前只处理过用制表符分隔的文件，现在遇到麻烦了，因为每个序列的长度不同，而且它们都在基因ID的下面。有什么建议吗？

Answer 1

搜索 biopython fasta 会出现 这个页面，然后我们可以修改里面的例子：

>>> from Bio import SeqIO
>>> with open("bio.fasta") as fp:
...         record_dict = SeqIO.to_dict(SeqIO.parse(fp, "fasta"))
...     

这样会生成一个看起来像这样的数据字典：

>>> import pprint
>>> pprint.pprint(record_dict)
{'SO_0001': SeqRecord(seq=Seq('MTKIAILVGTTLGSSEYIADEMQAQLTPLGHEVHTFLHPTLDELKPYPLWILVS...DQI', SingleLetterAlphabet()), id='SO_0001', name='SO_0001', description='SO_0001', dbxrefs=[]),
 'SO_0002': SeqRecord(seq=Seq('MTTPVDAPKWPRQIPYIIASEACERFSFYGMRNILTPFLMTALLLSIPEELRGA...FDQ', SingleLetterAlphabet()), id='SO_0002', name='SO_0002', description='SO_0002', dbxrefs=[]),
 'SO_0003': SeqRecord(seq=Seq('MTTDTIVAQATAPGRGGVGIIRISGDKATNVAMAVLGHLPKPRYADYCYFKSAS...EVD', SingleLetterAlphabet()), id='SO_0003', name='SO_0003', description='SO_0003', dbxrefs=[])}

我们可以通过这个字典来访问记录并计算字符的数量：

>>> record_dict["SO_0002"]
SeqRecord(seq=Seq('MTTPVDAPKWPRQIPYIIASEACERFSFYGMRNILTPFLMTALLLSIPEELRGA...FDQ', SingleLetterAlphabet()), id='SO_0002', name='SO_0002', description='SO_0002', dbxrefs=[])
>>> record_dict["SO_0002"].seq
Seq('MTTPVDAPKWPRQIPYIIASEACERFSFYGMRNILTPFLMTALLLSIPEELRGA...FDQ', SingleLetterAlphabet())
>>> record_dict["SO_0002"].seq.count("C")
2

接下来：

>>> count = {name: record.seq.count("C") for name, record in record_dict.items()}
>>> count
{'SO_0002': 2, 'SO_0003': 1, 'SO_0001': 3}

然后：

>>> import csv
>>> with open("c_count.csv", "wb") as fp:
...     writer = csv.writer(fp, delimiter="\t")
...     for k in sorted(count):
...         writer.writerow([k, count[k]])

会生成一个像这样的文件：

localhost-2:coding $ cat c_count.csv 
SO_0001 3
SO_0002 2
SO_0003 1

建议：不要自己重新写一个 FASTA 解析器，直接使用现有的；还有，不要重新实现 csv 模块。

Answer 2

如果你已经有了你发的那种格式的数据，并且不想去研究那些专业的库，你可以试试下面这个方法。

with open('datafile.txt') as file:
  datalist = []
  for line in file:
    if line.startswith('>'):
      datalist.append([line.strip()[1:], ''])
    else:
      datalist[-1][1] += line.strip()
  for data in datalist:
    print(data[0], '   ', data[1].count('C'))