Python: Розалинд Консенсус и Профиль

Question

Я пытаюсь решить задачу «Консенсус и профиль» для Розалинды. Инструкции для вызова следующие:

Дано: набор из не более 10 строк ДНК одинаковой длины (не более 1 т.п.н.) в формате FASTA.

Возврат: согласованная строка и матрица профиля для коллекции. (Если существует несколько возможных согласованных строк, вы можете вернуть любую из них.)

Мой код выглядит следующим образом (большую часть его я получил от другого пользователя на этом сайте). Моя единственная проблема заключается в том, что некоторые нити ДНК разбиты на несколько отдельных строк, поэтому они добавляются в список «allstrings» как отдельные строки. Я пытаюсь выяснить, как записать каждую последовательную строку, которая не содержит ">", в виде одной строки.

import numpy as np

seq = []
allstrings = []
temp_seq = []
matrix = []
C = []
G = []
T = []
A = []
P = []
consensus = []
position = 1

file = open("C:/Users/knigh/Documents/rosalind_cons (3).txt", "r")
conout = open("C:/Users/knigh/Documents/consensus.txt", "w")

# Right now, this is reading and writing each as an individual line. Thus, it
#  is splitting each sequence into multiple small sequences. You need to figure
#  out how to read this in FASTA format to prevent this from occurring
desc = file.readlines()

for line in desc:
    allstrings.append(line)

for string in range(1, len(allstrings)):
    if ">" not in allstrings[string]:
        temp_seq.append(allstrings[string])
    else:
        seq.insert(position, temp_seq[0])
        temp_seq = []
        position += 1

# This last insertion into the sequence must be performed after the loop to empty
#  out the last remaining string from temp_seq
seq.insert(position, temp_seq[0])

for base in seq:
    matrix.append([pos for pos in base])

M = np.array(matrix).reshape(len(seq), len(seq[0]))

for base in range(len(seq[0])):
    A_count = 0
    C_count = 0
    G_count = 0
    T_count = 0
    for pos in M[:, base]:
        if pos == "A":
            A_count += 1
        elif pos == "C":
            C_count += 1
        elif pos == "G":
            G_count += 1
        elif pos == "T":
            T_count += 1
    A.append(A_count)
    C.append(C_count)
    G.append(G_count)
    T.append(T_count)

profile_matrix = {"A": A, "C": C, "G": G, "T": T}

P.append(A)
P.append(C)
P.append(G)
P.append(T)

profile = np.array(P).reshape(4, len(A))

for pos in range(len(A)):
    if max(profile[:, pos]) == profile[0, pos]:
        consensus.append("A")
    elif max(profile[:, pos]) == profile[1, pos]:
        consensus.append("C")
    elif max(profile[:, pos]) == profile[2, pos]:
        consensus.append("G")
    elif max(profile[:, pos]) == profile[3, pos]:
        consensus.append("T")

conout.write("".join(consensus) + "\n")

for k, v in profile_matrix.items():
    conout.write(k + ": " + " ".join(str(x) for x in v) + "\n")

conout.close()

Answer 1

Есть несколько способов, которыми вы можете перебирать файл FASTA как записи. Вы можете использовать готовую библиотеку или написать свою собственную.

Широко используемая библиотека для работы с данными последовательности - biopython. Этот фрагмент кода создаст список строк.

from Bio import SeqIO


file = "path/to/your/file.fa"
sequences = []

with open(file, "r") as file_handle:
    for record in SeqIO.parse(file_handle, "fasta"):
        sequences.append(record.seq)

Кроме того, вы можете написать свой собственный анализатор FASTA. Примерно так должно работать:

def read_fasta(fh):
    # Iterate to get first FASTA header        
    for line in fh:
        if line.startswith(">"):
            name = line[1:].strip()
            break

    # This list will hold the sequence lines
    fa_lines = []

    # Now iterate to find the get multiline fasta
    for line in fh:
        if line.startswith(">"):
            # When in this block we have reached 
            #  the next FASTA record

            # yield the previous record's name and
            #  sequence as tuple that we can unpack
            yield name, "".join(fa_lines)

            # Reset the sequence lines and save the
            #  name of the next record
            fa_lines = []
            name = line[1:].strip()

            # skip to next line
            continue

        fa_lines.append(line.strip())

    yield name, "".join(fa_lines)

Вы можете использовать эту функцию так:

file = "path/to/your/file.fa"
sequences = []

with open(file, "r") as file_handle:
    for name, seq in read_fasta(file_handle):
        sequences.append(seq)