Блог

Главная - Вопросы по программированию - Пересечение двух файлов fasta с использованием python

Пересечение двух файлов fasta с использованием python

#python #bioinformatics #intersection #biopython #seq

#python #биоинформатика #пересечение #biopython #продолжение

Вопрос:

У меня есть два больших файла fasta — их структуры разные (показаны ниже), но заголовки (начинающиеся с @) для чтения одинаковы в обоих файлах:

файл1

 >MN00153:75:000H37WNG:1:12106:12990:1333
AAAACCCC

>MN00153:75:000H37WNG:1:12106:21652:2374
AAAAGGGG

>MN00153:75:000H37WNG:1:12106:21652:2366
AGGGGGTT

файл2

 >MN00153:75:000H37WNG:1:12106:12990:1333
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCAGATCTCGCCC

>MN00153:75:000H37WNG:1:12106:21652:2374
AGATCTCGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT

>MN00153:75:000H37WNG:1:12106:21652:2366                  
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT

Я создал словарь из headers (ключа) и reads (значения) file1 с помощью скрипта:

 from Bio import SeqIO

dict={}
with open ('index2.fasta', 'r') as file1:
    for record in SeqIO.parse(file1, 'fasta'):
        dict[str(record.id)] = str(record.seq)

Что я сделал, так это зациклил чтение в file2, и если 'AGATCTCG' строка была внутри чтения, я сохранил эти заголовки чтения в списке.

Теперь у меня проблема в том, что я хочу создать вложенный файл file2 на основе dictionary и list . Если элемент моего списка существует как ключ в моем словаре, и если значение равно 'AAAACCCC' , то вывод должен быть MN00153:75:000H37WNG:1:12106:12990:1333 , НО я получаю оба MN00153:75:000H37WNG:1:12106:12990:1333 и MN00153:75:000H37WNG:1:12106:21652:2374 

 ATTACTCG_ids=[]
with open ('Read1.fasta', 'r') as file2:
    for record in SeqIO.parse(file2, 'fasta'):
        if 'AGATCTCG' in record.seq:
            ATTACTCG_ids.append(record.id)
            for i in ATTACTCG_ids:
                if dict.get(i) == 'AAAACCCC':
                    final = record.format('fasta')
                    print(final)

Может кто-нибудь помочь мне с этим?

Комментарии:

1. Ваши заголовки не начинаются с > , так что это не FASTA. @ похоже на заголовок FASTQ (но без оценки качества и т. Д.) Похоже, у вас искаженный формат файла, и я удивлен, если Biopython или любой другой нестандартный модуль правильно его анализирует

2. Я только что исправил это спасибо за предупреждение

Ответ №1:

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

 with open ('Read1.fasta', 'r') as file2:
    for record in SeqIO.parse(file2, 'fasta'):
        if 'AGATCTCG' in record.seq and dict.get(record.id) == 'AAAACCCC':
            final = record.format('fasta')
            print(final)
Метки: