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 に答える