Невозможно проанализировать только последовательности из файла FASTA

Как я могу удалить идентификаторы типа '>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n' из последовательностей?

У меня есть такой код:

with open('sequence.fasta', 'r') as f :
    while True:
        line1=f.readline()
        line2=f.readline()
        line3=f.readline()
        if not line3:
            break
        fct([line1[i:i+100] for i in range(0, len(line1), 100)])
        fct([line2[i:i+100] for i in range(0, len(line2), 100)])
        fct([line3[i:i+100] for i in range(0, len(line3), 100)])

Вывод:

['>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n']
['CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG\n']
['AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG\n']
['CCGCCTCGGGAGCGTCCATGGCGGGTTTGAACCTCTAGCCCGGCGCAGTTTGGGCGCCAAGCCATATGAA\n']
['AGCATCACCGGCGAATGGCATTGTCTTCCCCAAAACCCGGAGCGGCGGCGTGCTGTCGCGTGCCCAATGA\n']
['ATTTTGATGACTCTCGCAAACGGGAATCTTGGCTCTTTGCATCGGATGGAAGGACGCAGCGAAATGCGAT\n']
['AAGTGGTGTGAATTGCAAGATCCCGTGAACCATCGAGTCTTTTGAACGCAAGTTGCGCCCGAGGCCATCA\n']
['GGCTAAGGGCACGCCTGCTTGGGCGTCGCGCTTCGTCTCTCTCCTGCCAATGCTTGCCCGGCATACAGCC\n']
['AGGCCGGCGTGGTGCGGATGTGAAAGATTGGCCCCTTGTGCCTAGGTGCGGCGGGTCCAAGAGCTGGTGT\n']
['TTTGATGGCCCGGAACCCGGCAAGAGGTGGACGGATGCTGGCAGCAGCTGCCGTGCGAATCCCCCATGTT\n']
['GTCGTGCTTGTCGGACAGGCAGGAGAACCCTTCCGAACCCCAATGGAGGGCGGTTGACCGCCATTCGGAT\n']
['GTGACCCCAGGTCAGGCGGGGGCACCCGCTGAGTTTACGC\n']
['\n']
...

Моя функция:

def fct(input_string):
    code={"a":0,"c":1,"g":2,"t":3}
    p=[code[i] for i in input_string]
    n=len(input_string)
    c=0

    for i, n in enumerate(range(n, 0, -1)):
        c +=p[i]*(4**(n-1))
        return c+1

fct() возвращает целое число из строки. Например, ACT дает 8, то есть: моя функция должна принимать в качестве входных строк последовательности, содержащие только следующие базы A, C, G, T

Но когда я использую свою функцию, она дает:

KeyError: '>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\n' 

Я пытаюсь удалить идентификаторы, удаляя строки, начинающиеся с >, и записывая остальное в текстовый файл, поэтому мой текстовый файл output.txt содержит только последовательности без идентификаторов, но когда я использую свою функцию fct, я обнаружил ту же ошибку:

KeyError: 'CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG\n'

Что я могу сделать?


person m28    schedule 29.07.2013    source источник
comment
С другой стороны, вам следует подумать о чтении PEP8, чтобы очистить ваше форматирование Python. Самое главное, используйте пробелы, а не табуляции!   -  person David Cain    schedule 30.07.2013
comment
Пожалуйста, не удаляйте свои комментарии и вопросы - такие вопросы могут быть полезны другим людям, у которых есть похожие проблемы.   -  person David Cain    schedule 16.08.2013


Ответы (1)


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

Анализ данных FASTA

Могу я предложить использовать отличный пакет Biopython? Он имеет отличную встроенную поддержку FASTA (чтение и запись) (см. Последовательности в Руководстве).

Чтобы проанализировать последовательности из файла FASTA:

for seq_record in SeqIO.parse("seqs.fasta", "fasta"):
    print record.description  # gi|2765658|emb|Z78533.1...
    print record.seq  # a Seq object, call str() to get a simple string

>>> print record.id
'gi|2765658|emb|Z78533.1|CIZ78533'

>>> print record.description
'gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA'

>>> print record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet())

>>> print str(record.seq)
'CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACC'  #(truncated)

Итерация по данным последовательности

В вашем коде у вас есть список строк, передаваемых в fct() (input_string на самом деле не строка, а список строк). Решение состоит в том, чтобы просто создать одну входную строку и перебрать ее.

Другие ошибки в fct:

  • Ключи к словарю нужно писать с заглавной буквы: регистр имеет значение.
  • У вас должен быть оператор return после цикла for. Сохранение вложенности означает немедленное возвращение c.
  • Зачем беспокоиться о построении p, если вы можете просто проиндексировать code при итерации по последовательности?
  • Вы записываете длину последовательности (n), используя ее в своем for цикле в качестве имени переменной.

Измененный код (с правильным форматированием PEP 8) и переменные, переименованные для большей ясности что они означают (до сих пор не знаю, что такое c):

from Bio import SeqIO


def dna_seq_score(dna_seq):
    nucleotide_code = {"A": 0, "C": 1, "G": 2, "T": 3}

    c = 0 
    for i, k in enumerate(range(len(dna_seq), 0, -1)):
        nucleotide = dna_seq[i]
        code_num = nucleotide_code[nucleotide]
        c += code_num * (4 ** (k - 1)) 
    return c + 1 


for record in SeqIO.parse("test.fasta", "fasta"):
    dna_seq_score(record.seq)
person David Cain    schedule 30.07.2013
comment
На самом деле, если вы попытаетесь запустить этот код, вы увидите другое сообщение об ошибке. Я обновил свой ответ, чтобы объяснить, что здесь происходит. - person David Cain; 30.07.2013
comment
Возможно, вам будет полезен этот вопрос: stackoverflow.com/questions/312443/ - person David Cain; 30.07.2013
comment
Хорошо. Вы пробовали преобразовать свой код для работы с объектом Seq или с последовательностью в виде строки? Если у вас все еще есть проблемы, задайте другой вопрос, и я уверен, что вы получите помощь в их решении. - person David Cain; 30.07.2013
comment
Вам нужно будет быть более конкретным. Либо отредактируйте свой вопрос, указав точный код, который вы используете, либо задать новый вопрос. - person David Cain; 30.07.2013