Печать частоты последовательностей в файле fasta (python) - PullRequest
Новая
       115

Печать частоты последовательностей в файле fasta (python)

1 голос
/ 13 июля 2020

Я пытаюсь найти разнообразие совокупности последовательностей из большого файла fasta. Конечная цель - создать гистограмму распределения.

Я написал приведенный ниже код, чтобы подсчитать, сколько раз каждая последовательность встречается в файле fasta. Я сделал это, чтобы добавить счетчик в конец идентификатора. Вместо этого форматирования я хотел бы напечатать выходной файл, в котором просто говорится, что одна последовательность встречается x раз. Последовательности y встречаются z раз и т. д. без последовательности и идентификатора. 

from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("filename.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("filename_output.fasta", 'w') as output:
    for seq, ids in sorted(dedup_records.items(), key=lambda t: len(t[1]), reverse=True):
        output.write(">{}_counts{}\n".format(ids[0], len(ids)))
        output.write(seq + "\n")

На изображении показан фрагмент выходного файла

Из этого изображения I хотел бы, чтобы вывод был напечатан: 1 последовательность встречается 1885 раз 1 последовательность встречается 1099 раз 1 последовательность встречается 280 раз.

Кроме того, когда несколько последовательностей встречаются одинаковое количество раз, каждая из них распечатывается отдельно. Я не уверен, как их объединить. EX

Пожалуйста, дайте мне знать, если у вас есть какие-либо предложения. Большое вам спасибо.

1 Ответ

0 голосов
/ 14 июля 2020

collections.Counter ()

Используйте счетчик дважды. Примерно так: 

from Bio import SeqIO
from collections import Counter

# counts the number of times each sequence occurs
sequences = SeqIO.parse("filename.fasta", "fasta")
seq_counts = Counter(str(record.seq) for record in sequences)

# counts how many repeat 2, 3, 4, ... times
count_repeats = Counter(seq_counts.values())

with open("filename_output.fasta", 'w') as output:
    for repeat, num_seqs in count_repeats.most_common():
        output.write(f">{num_seqs} sequences occur {repeat} times\n")
...