Подсчет появления определенного символа для каждого гена в файле fasta с использованием python

У меня есть файл fasta следующим образом:

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

Где то, что следует за «>», является идентификатором гена, а буквы, следующие за строкой «>», являются соответствующими последовательностями. Я хочу проанализировать файл и подсчитать, сколько «C» есть в последовательности для каждого идентификатора гена. Я хотел бы, чтобы мой выходной файл был файлом с разделителями табуляции, например:

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

и так далее...

Я использую python и думал, что это будет прямолинейно, сделав ключи идентификаторов генов в словаре, но я сделал это только с файлами с разделителями табуляции, и у меня возникли проблемы, поскольку каждая последовательность имеет разную длину и ниже идентификаторы генов. Любые предложения были бы замечательными!


python counter biopython
person C9r1y    schedule 21.02.2013    source источник

Ответы (2)


arrow_upward
0
arrow_downward

Если у вас уже были данные в том формате, который вы разместили, и вы не хотите копаться в специализированных библиотеках, вы можете попробовать что-то вроде этого.

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'))
person Octipi    schedule 22.02.2013
comment
Спасибо, это сработало очень хорошо. Могу я спросить, почему вы добавляете '' в datalist.append([line.strip()[1:], '']) Какова его цель? Спасибо еще раз! - person C9r1y; 22.02.2013
comment
Я не был уверен, на сколько строк будет распределена последовательность. Он используется в качестве заполнителя для следующих строк, содержащих информацию о последовательности. - person Octipi; 22.02.2013

arrow_upward
4
arrow_downward

При поиске 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.

person DSM    schedule 22.02.2013
comment
Спасибо за шаг за шагом! Определенно тоже попробую. - person C9r1y; 22.02.2013