Извлечь последовательности из файла FASTA в несколько файлов, файл на основе header_IDs в отдельном файле

Я ищу решение python для извлечения нескольких последовательностей из файла FASTA в несколько файлов на основе соответствия списку идентификаторов заголовков. в отдельном файле.

Это немного более сложная версия проблемы, размещенной на Извлечь последовательности из файла FASTA на основе записей в отдельном файле и https://www.biostars.org/p/2822/, который выводит только один файл для всех совпадений.

Я новичок в Python и пытаюсь найти способ:

  • Возьмите файл со строками, который будет в заголовках фаста
  • Все записи, которые соответствуют строке, записаны в отдельный файл fasta

Файл header_ID_strings выглядит так:

CAP357_2030_09WPI, CAP357_2040_11WPI, CAP357_2050_13WPI и т. Д.

образец моего фаст-файла выглядит так:

> CAP357_2030_009wpi_v1v3

#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
056_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGG
> CAP357_2040_011wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
008_00006_001.1
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
> CAP357_2040_011wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
030_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
> CAP357_2040_011wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
004_00001_000.2
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
> CAP357_2050_013wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
047_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT

ожидаемый результат

file1: CAP357_2030_009wpi_v1v3.fasta
> CAP357_2030_009wpi_v1v3

#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
056_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGG

file2: CAP357_2040_011wpi_v1v3.fasta
> CAP357_2040_011wpi_v1v3

#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
008_00006_001.1
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
> CAP357_2040_011wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
030_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
> CAP357_2040_011wpi_v1v3
#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)
004_00001_000.2
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
и т.д. ...

Этот код взят из приведенной выше ссылки, но я хочу, чтобы:
* совпадения записывались в отдельные аутфайлы
* Мне не нужно указывать каждый аутфайл отдельно, если это возможно (у меня будет до 30 аутфайлов) < Br>

#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)

Пока это то, что я придумал (сценарий ниже), но не знаю, как обойтись вручную, указав все файлы и переменные (я включил здесь только три)

from Bio import SeqIO
import pandas as pd
import sys

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file2020 = sys.argv[3]
output_file2030 = sys.argv[4]
output_file2040 = sys.argv[5]

colnames = ["2020", "2030", "2040"]
headerlist = pd.read_csv(id_file, names = colnames, header = None)
infile = list(SeqIO.parse(input_file, "fasta"))
2020_seq = tuple(headerlist.2020)
2030_seq = tuple(headerlist.2030)
2040_seq = tuple(headerlist.2040)

count2020 = 0
count2030 = 0
count2040 = 0
for record in infile:
    if record.id in 2020_seq:
        SeqIO.write([record], output_file2020, "fasta")
        countSU += 1
    elif record.id in PI_seq:
        SeqIO.write([record], output_file2030, "fasta")
        countPI += 1
    elif record.id in REC_seq:
        SeqIO.write([record], output_file2040, "fasta")
        countREC += 1
    else:
        print("no matches found")

print("number of SU is", count2020)
print("number of PI is", count2030)
print("number of REC is", count2040)

python regex biopython fasta
person Colin Anthony    schedule 14.12.2014    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Пара кратких предложений:

Если все ваши заголовки следуют одному шаблону, вы можете извлечь уникальные элементы:

record.description.split("_")[1]

(дает "2040" из "CAP357_2040_011wpi_v1v3

record.description.split("_")[1]
008_00006_001.1")

Если вы используете dict, вы можете собирать коллекции записей:

collected = {}
for record in records:
    descr = record.description.split("_")[1]
    try:
        collected[descr].append(record)
    except KeyError:
        collected[descr] = [record ,]

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

file_name = "outfile%s" 
for (descr, records) in collected.items():   # iteritems in python2
    with open(os.path.join(file_path, file_name % descr), 'w') as f:
        SeqIO.write(records, f, 'fasta')
person iayork    schedule 14.12.2014
comment
Спасибо за совет! Я получаю сообщение об ошибке Нет такого файла или каталога, которое, как я понимаю, связано с открытием файла, который еще не был создан? Когда я вручную создаю пустой файл outfile4260 в каталоге out (чтобы проверить, что не так), я получаю, что файл не открыт для записи ошибки. Нужно ли мне сохранять вывод как переменную перед записью с помощью SeqIO.write? Извините, если это действительно очевидно, это мой первый скрипт на Python. - person Colin Anthony; 16.12.2014
comment
PS: я установил file_path = os.getcwd () - person Colin Anthony; 16.12.2014
comment
Ошибка в моем коде: дескриптор файла должен иметь форму открытого (путь / к / файлу, 'w') ('w' для 'записи' - будет 'r' для чтения файла). Попробуйте с буквой "w" (ответ исправил) - person iayork; 16.12.2014

arrow_upward
2
arrow_downward

Для полноты картины вот «окончательный» сценарий:

#!/usr/bin/env python
# a script to extract fasta records from a fasta file to multiple separate fasta files based on a particular ID (time point) in a particular field, for a given delimiter
# to run, navigate to file location with command prompt and enter: python split_fasta_by_collections.py infile.fasta
from Bio import SeqIO
import os
import sys

records = SeqIO.parse(sys.argv[1], "fasta")
collected = {}
for record in records:
    descr = record.description.split("_")[1] # "_" sets the delimeter, "1" sets the field where counting starts at 0 for the first field
    try:
        collected[descr].append(record)
    except KeyError:
        collected[descr] = [record ,]

file_name = "outfile%s.fasta" 
file_path = os.getcwd() #sets the output file path to your current working directory

for (descr, records) in collected.items():  
    with open(os.path.join(file_path, file_name % descr), 'w') as f:
        SeqIO.write(records, f, 'fasta')
person Colin Anthony    schedule 17.03.2015