0

我有一个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之下。任何建议都会很棒！

4

2 回答 2

4

搜索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模块。

于 2013-02-22T00:16:29.137 回答
0

如果您已经拥有您发布的格式的数据并且不想深入研究专门的库，您可以尝试这样的事情。

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'))
于 2013-02-22T00:20:53.280 回答