2016-05-13 6 views
3

Я хочу извлечь конкретные последовательности fasta из большого файла fasta, используя следующий скрипт, но выход пуст.Извлечь конкретные последовательности fasta из большого файла fasta

Файл transcripts.txt содержит идентификаторы транскриптов списка, которые я хочу экспортировать (как идентификаторы, так и последовательности) от assembly.fasta до selected_transcripts.fasta. Например:

  1. transcripts.txt:
     
    Transcript_00004|5601 
    Transcript_00005|5352
  2. assembly.fasta:
     
    >Transcript_00004|5601 
    GATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT 
    >Transcript_00004|5360 
    CGATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT 
    

Идентификаторы предшествует символ >: >Transcripts_00004|5601.

Я должен прочитать файл assembly.fasta, если запись ID в assembly.fasta тот же этой записи в transcripts.txt, я должен написать эту расшифровку ID и его последовательность в selected_transcripts.fasta. Итак, в приведенном выше примере мне нужно написать только первый расшифрованный стенограмму.

Любые предложения? Спасибо.

from Bio import SeqIO 

my_list = [line.split(',') for line in open("/home/universita/transcripts.txt")] 

fin = open('/home/universita/assembly.fasta', 'r') 
fout = open('/home/universita/selected_transcripts.fasta', 'w') 

for record in SeqIO.parse(fin,'fasta'): 
    for item in my_list: 
     if item == record.id: 
      fout.write(">" + record.id + "\n") 
      fout.write(record.seq + "\n") 

fin.close() 
fout.close() 
+1

см. Https://www.biostars.org/p/68718/ – Pierre

+0

Можете ли вы отредактировать свой вопрос и включить некоторые из 'transcripts.txt', а также часть' assembly.fasta', поэтому у нас есть некоторые данные для работы? – MattDMo

+0

Вы разделили строку транскриптов после каждого двоеточия, но это пространство разделено. Это специально? –

ответ

1

На основании ваших примеров есть несколько небольших проблем, которые могут объяснить, почему вы ничего не получаете. Ваш transcripts.txt имеет несколько записей в одной строке, поэтому my_list будет иметь все элементы на first_line в my_line[0], в цикле вы перебирать my_list линиями, так что ваш первый пункт будет

['Transcript_00004|5601', 'Transcript_00005|5352']

Кроме того, если assembly.fasta не имеет > в строке заголовка, вы не получите никаких записей с идентификаторами и последовательностями. Следующий код должен позаботиться об этих проблемах, предполагая, что вы добавили > в заголовки, а функция split теперь использует пробел, а не двоеточие.

from Bio import SeqIO 

my_list = [] 
with open("transcripts.txt") as transcripts: 
    for line in transcripts: 
     my_list.extend(line.split(' ')) 

fin = open('assembly.fasta', 'r') 
fout = open('selected_transcripts.fasta', 'w') 

for record in SeqIO.parse(fin,'fasta'): 
    for item in my_list: 
     if item.strip() == record.id: 
      fout.write(">" + record.id + "\n") 
      fout.write(record.seq + "\n") 


fin.close() 
fout.close() 

Чтение стенограммы было изменено таким образом, что все идентификаторы добавляются my_list отдельно. Также каждый элемент лишен белого пространства, чтобы избежать разрыва строки в строке при сравнении с record.id.