Как я могу написать скрипт для суммирования информации из файла мульти-fast без использования биопиона?

Question

У меня есть задание для класса биоинформатики, который запрашивает скрипт на Python для выполнения следующих действий для файла FASTA с несколькими последовательностями белков:

-Открытие файла .fasta, указанного пользователем

-Печать строку заголовка (то есть строку, начинающуюся с ">")

-Печать первые 10 аминокислот и в той же строке, сообщить количество аминокислот в последовательности

После нескольких часов попыток мне удалось получить только строку заголовка и первые 10 аминокислот, напечатанные только для первой последовательности. Цикл for, который я написал, похоже, не сработает (извините, если это мусор, я новичок!)

input_file = open(input("\nPlease enter the name of a FASTA file (inlcuding the .fasta extension): "))
# Opens the file specified by the user
for line in input_file:
    line = line.rstrip()
    if line[0] == ">":
        header = line
        print("\nHeader:", header)  # prints the header line
        no_lines_searched = 0  # To stop loop after first line of sequence when printing first 10 AAs
        for i in input_file:
            if ">" not in i and no_lines_searched < 1:
                print("First ten amino acid residues: ", i[0:10], "# Amino acids: ")  # Prints the first 10 AAs
                no_lines_searched = no_lines_searched+1  # To stop loop after first line of sequence

Я пытался быть умным и спроектировал второй циклтаким образом, что он вернет первые 10 аминокислот последовательности, затем остановится, пока не встретится другая последовательность (обозначаемая «>»).

Затем я планировал использовать заполнитель %s, чтобы каким-то образом рассчитать общее количество. длина последовательности для каждой последовательности в файле, но, кажется, не может пройти эту точку!

Вывод, который я получаю, выглядит следующим образом:

Header: >sp|P03378|ENV_HV1A2 Envelope glycoprotein gp160 OS=Human immunodeficiency virus type 1 group M subtype B (isolate ARV2/SF2) GN=env PE=3 SV=1
First ten amino acid residues:  MKVKGTRRNY # Amino acids:

Любая помощь будет наиболее ценной!

Answer 1

Потратив большую часть сегодняшнего дня на эту проблему, я ответил на свой вопрос. Я делал это совершенно неправильно. Вот код, который я использовал для достижения именно того, к чему стремился:

#!/usr/bin/env python3
seq = ""  # seq is a variable used to store each line of sequence data as the loop below iterates.

ten_AAs = ''  # ten_AAs is a variable used to store the first ten amino acids of each sequence (it actually stores
# the first 10 AAs from each line of the sequence, but only the first 10 residues of each sequence are printed).

seq_count = 0  # seq_count is a variable used to track the total number of sequences in the file.

infile = input("\nPlease enter the name of a FASTA file (including the .fasta extension): ")
# Stores the file specified by the user in the variable, infile.

with open(infile, "r") as fasta:  # Opens the above file such that we do not need to call filename.close() at the end.
    print("\n")
    for line in fasta:  # Initiates a for loop to summarise the file: headers, first ten residues and sequence lengths.
        line = line.rstrip()  # Removes trailing newline characters.
        if line[0] == ">":  # Indicates the start of a new FASTA sequence and is a condition for printing summary data
            if seq != "":  # Second condition for printing summary data (see lines 28-30).
                print("Header: {0}\nFirst ten amino acids: {1} \tSequence length: {2} residues\n"
                      .format(header, ten_AAs[0:10], len(seq)))  # Placeholders for summary info.
            seq_count += 1  # Increases seq_count by one for each header encountered.
            header = line  # Stores the header line for printing.
            seq = ""  # Resets seq for each iteration.
            ten_AAs = ''  # Resets ten_AAs for each iteration.
        else:
            seq += line[:]  # Appends seq with the characters of the line if it is not a header (i.e. sequence data).
            ten_AAs += line[0:10]  # Appends ten_AAs with the first ten characters of the line if it is not a header.

# After reading all lines of one sequence, the loop encounters the next FASTA sequence, starting with a ">" character,
# but 'seq' is not blank. Therefore, both 'if' statements are now TRUE and all summary data for the first sequence
# (stored in 'header', 'ten_AAs' and 'seq') are printed as per lines 18-19.

print("Header: {0}\nFirst ten amino acids: {1} \tSequence length: {2} residues"
      .format(header, ten_AAs[0:10], len(seq)))
# The final sequence is not followed by another accession number; therefore, the summary data for the final sequence
# must be printed manually, outside of the loop.

print("\nThis multi-FASTA file contains", str(seq_count), "sequences.")
# For completeness' sake, reports the number of sequences in the multi-FASTA file.

А вот вывод:

Please enter the name of a FASTA file (including the .fasta extension): multiprotein.fasta


Header: >1433G_HUMAN (P61981) 14-3-3 protein gamma (Protein kinase C inhibitor protein 1) (KCIP-1) [Homo sapiens]
First ten amino acids: VDREQLVQKA       Sequence length: 246 residues

Header: >ATP8_RAT (P11608) ATP synthase protein 8 (EC 3.6.3.14) (ATPase subunit 8) (A6L) (Chargerin II) [Rattus norvegicus]
First ten amino acids: MPQLDTSTWF       Sequence length: 67 residues

Header: >ALR_LISIN (Q92DC9) Alanine racemase (EC 5.1.1.1)
First ten amino acids: MVTGWHRPTW       Sequence length: 368 residues

Header: >CDCA4_HUMAN (Q9BXL8) Cell division cycle-associated protein 4 (Hematopoietic progenitor protein) [Homo sapiens]
First ten amino acids: MFARGLKRKC       Sequence length: 241 residues

This multi-FASTA file contains 4 sequences.