Извлечение последовательностей из файла fasta по номеру ID в заголовке

У меня есть файл fasta с несколькими последовательностями с заголовками, которые выглядят так:

>1016BSA34080.1
MTHSVRIITVTVNFLQHRFFIDYMSEIGLLDGEIEQMVSALQEQVHIVARARTLPEMKNLERDTHVIVKT
LKKQLTAFHSEVKKIADSTQRSRYEGKHQTYEAKVKDLEKELRTQIDPPPKSVSEKHMEDLMGEGGPDGS
GFKTTDQVLRAGIRIQNDA

>1038BSA81955.1
MQQQQARRRMEEPTAAAATASSTTSFAAQPLLSRSVAPQAASSPQASARLAESAGFRSAAVFGSAQAAVG
GRGRGGFGAPPGRGGFGAPPAAGFGAAPAFGAPPTLQAFSAAPAPGGFGAPPAPQGFGAPRAAGFGAPPA
PQAFSAVAPASSTAIPLDVTTYLGDTFGSAPTRGPP

Четырехзначное число в начале заголовка является уникальным идентификатором последовательности.

Не могли бы вы помочь мне написать сценарий Python для извлечения последовательностей по 4-значному идентификатору (в текстовом файле с одним идентификатором на строку)?

Я попытался изменить этот сценарий (на этом веб-сайте я нашел: Извлечь последовательности из файла FASTA на основе записей в отдельном файле) для моей цели (напрасно):

f2 = open('accessionids.txt','r')
f1 = open('fasta.txt','r')
f3 = open('fasta_parsed.txt','w')

AI_DICT = {}
for line in f2:
    AI_DICT[line[:-1]] = 1

skip = 0
for line in f1:
    if line[0] == '>':
        _splitline = line.split('|')
        accessorIDWithArrow = _splitline[0]
        accessorID = accessorIDWithArrow[1:-1]
        # print accessorID
        if accessorID in AI_DICT:
            f3.write(line)
            skip = 0
        else:
            skip = 1
    else:
        if not skip:
            f3.write(line)

f1.close()
f2.close()
f3.close()

Я новичок в Python, буду благодарен за любую помощь! Спасибо -Дивья


python header uniqueidentifier extract sequences
person user2662391    schedule 07.08.2013    source источник
comment
Просто чтобы убедиться, что я правильно понимаю: вы хотите извлечь полную последовательность с 4-значным идентификатором?   -  person Barranka    schedule 08.08.2013

Ответы (2)


arrow_upward
1
arrow_downward

accessionids.txt содержит только четырехзначные коды?

Если да, измените accessorID на:

accessorID = accessorIDWithArrow[1:5]

Вот несколько способов сделать это более питоническим:

Используйте набор вместо словаря для AI_DICT, используйте strip() вместо нарезки для удаления новой строки и используйте выражение генератора для создания набора

AI_SET = set((line.strip() for line in f2))

Используйте True и False вместо 0 и 1 для skip.

Я бы переделал основной цикл таким образом:

in_accession_ids = False
for line in f1:
    if line[0] == '>':
        _splitline = line.split('|')
        accessorIDWithArrow = _splitline[0]
        accessorID = accessorIDWithArrow[1:5]
        # print accessorID
        in_accession_ids = accessorID in AI_SET
    if in_accession_ids:
        f3.write(line)

Думаю, здесь логика более очевидна. Кроме того, начало с skip = 0 в оригинале или с in_accession_ids=True в моем означало бы, что вы распечатываете все до того, как найдете первый заголовок последовательности. Это может быть то, что вы хотите, а может и нет - я предположил, что не в моем переписывании.

В конце концов, вы можете захотеть заглянуть в коллекцию Biopython - это избыточно для этой конкретной задачи, но в целом довольно приятно. Среди прочего, множество инструментов для чтения файлов FASTA и связанных форматов.

http://biopython.org/wiki/Biopython

person Peter DeGlopper    schedule 07.08.2013

arrow_upward
1
arrow_downward

Используя Biopython, вы можете сделать это так (требуется установленный biopyhton):

from Bio import SeqIO

f1 = "fasta.fa"
f2 = "accessionids.txt"
f3 = "selected_seqs.fa"
selected_seqs = list()

with open(f2, "r") as seq_ids:
    accessionids = [line.rstrip("\n") for line in seq_ids]

for seq_record in SeqIO.parse(f1, "fasta")
    header = seq_record.name # (or .id or so)
    for accession_id in accessionids:
        if accession_id == header[0:4]:
            selected_seqs.append(seq_record)


SeqIO.write(selected_seqs, f3, "fasta")

Это будет проходить через ваши записи последовательности (файл fasta) и для каждой записи проверять, есть ли совпадение с идентификатором из файла accessionids.

Примечание:

  • seq_record может иметь разные теги, проверьте, в каком из них находится ваш идентификатор.
  • Если у вас есть идентификатор, который не указан в начале (и уникален для одного заголовка последовательности), вы просто используете if accession_id in header:
  • Дополнительные сведения о SeqIO см. В руководстве по биопайтону, глава 5.
person TheRibosome    schedule 13.06.2020