Как я могу получить этот вывод из файла FASTA без использования Biopython?

Question

Мне нужно получить вывод, показанный ниже, из файла FASTA, но без использования BioPython.У кого-нибудь есть идея?

Этот код использует BioPython:

from Bio import SeqIO
records = SeqIO.parse("data/assembledSeqs.fa", "fasta")
for i, seq_record in enumerate(records):
    print("Sequence %d:" % i)
    print("Number of A's: %d" % seq_record.seq.count("A"))
    print("Number of C's: %d" % seq_record.seq.count("C"))
    print("Number of G's: %d" % seq_record.seq.count("G"))
    print("Number of T's: %d" % seq_record.seq.count("T"))
    print()

Файл FASTA выглядит следующим образом:

>chr12_9180206_+:chr12_118582391_+:a1;2 total_counts: 115 Seed: 4 K:    20 length: 79
TTGGTTTCGTGGTTTTGCAAAGTATTGGCCTCCACCGCTATGTCTGGCTGGTTTACGAGC
AGGACAGGCCGCTAAAGTG
>chr12_9180206_+:chr12_118582391_+:a2;2 total_counts: 135 Seed: 4 K: 20 length: 80
CTAACCCCCTACTTCCCAGACAGCTGCTCGTACAGTTTGGGCACATAGTCATCCCACTCG
GCCTGGTAACACGTGCCAGC
>chr1_8969882_-:chr1_568670_-:a1;113 total_counts: 7600 Seed: 225 K: 20 length: 86
CACTCATGAGCTGTCCCCACATTAGGCTTAAAAACAGATGCAATTCCCGGACGTCTAAAC
CAAACCACTTTCACCGCCACACGACC
>chr1_8969882_-:chr1_568670_-:a2;69 total_counts: 6987 Seed: 197 K: 20   length: 120
TGAACCTACGACTACACCGACTACGGCGGACTAATCTTCAACTCCTACATACTTCCCCCA
TTATTCCTAGAACCAGGCGACCTGCGACTCCTTGACGTTGACAATCGAGTAGTACTCCCG

Мне нужно получить следующий вывод:

Sequence 0:
Number of A's: 14
Number of C's: 17
Number of G's: 24
Number of T's: 24

Sequence 1:
Number of A's: 17
Number of C's: 30
Number of G's: 16
Number of T's: 17

Sequence 2:
Number of A's: 27
Number of C's: 31
Number of G's: 12
Number of T's: 16

Sequence 3:
Number of A's: 31
Number of C's: 41
Number of G's: 20
Number of T's: 28

Я пробовал это, но не могу получить тот же вывод.

def count_bases (fasta_file_name):
    with open(fasta_file_name) as file_content:
        for seqs in file_content:
            if seqs.startswith('>'):
                for i, seq in enumerate('>'):
                    print("Sequence %d:" % i)
            else:
                print("Number of A's: %d" % seqs.count("A"))
                print("Number of C's: %d" % seqs.count("C"))
                print("Number of G's: %d" % seqs.count("G"))
                print("Number of T's: %d" % seqs.count("T"))
                print()
    return bases

result = count_bases('data/assembledSeqs.fa')

Answer 1

Этот код будет работать:

def count_bases (fasta_file_name):
    sequece=''
    def count():
        if len(sequece):
            print("Number of A's: %d" % sequece.count("A"))
            print("Number of C's: %d" % sequece.count("C"))
            print("Number of G's: %d" % sequece.count("G"))
            print("Number of T's: %d" % sequece.count("T"))
            print()
    with open(fasta_file_name) as file_content:
        i=0
        for seqs in file_content:
            if seqs.startswith('>'):
                count()
                print("Sequence %d:" % i)
                i=i+1
                sequece=''
            else:
                sequece=sequece+seqs.strip()
        count()

result = count_bases('data/assembledSeqs.fa')

Answer 2

Это похоже на работу.Он использует модуль регулярного выражения re для распознавания строк заголовка и получения общей длины следующей нуклеотидной последовательности.Это значение используется для определения количества следующих строк данных для совместного чтения и объединения.

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

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

import re

pattern = r""">.+?length:\s(\d+)"""
regex = re.compile(pattern)

MAX_PER_LINE = 60

def count_bases (fasta_file_name):
    bases = None  # Undefined.

    with open(fasta_file_name) as inp:
        i = 0  # Sequence counter.
        line = next(inp, None)  # Read first line.
        while line:
            match = regex.search(line)
            if match:
                length = int(match.group(1))
                nlines = (length + MAX_PER_LINE-1) // MAX_PER_LINE
                # Read and concatenate data from the required number of lines.
                seqs = ''.join(next(inp).rstrip() for _ in range(nlines))

                print("Sequence %d:" % i)
                print("Number of A's: %d" % seqs.count("A"))
                print("Number of C's: %d" % seqs.count("C"))
                print("Number of G's: %d" % seqs.count("G"))
                print("Number of T's: %d" % seqs.count("T"))
                print()
                i += 1

            line = next(inp, None)  # Read next line.

    return bases

fasta_file = 'assembledSeqs.fa'
result = count_bases(fasta_file)

Вот вариант, который использует collections.Counter для подсчета количества всех оснований одновременно.Это может быть быстрее, чем использование метода последовательности count отдельно для каждого, и поэтому полезно, если у вас есть большое число для обработки.

def count_bases (fasta_file_name):
    bases = None  # Undefined.

    with open(fasta_file_name) as inp:
        i = 0  # Sequence counter.
        line = next(inp, None)  # Read first line.
        while line:
            match = regex.search(line)
            if match:
                length = int(match.group(1))
                nlines = (length + MAX_PER_LINE-1) // MAX_PER_LINE
                # Read and concatenate data from the required number of lines.
                seqs = ''.join(next(inp).rstrip() for _ in range(nlines))

                print("Sequence %d:" % i)
                counter = Counter(seqs)  # Count the number of each base.
                for base in 'ACGT':
                    if base in counter:
                        print("Number of {}'s: {}".format(base, counter[base]))
                print()
                i += 1

            line = next(inp, None)  # Read next line.

    return bases