從大的fasta文件中提取特定的fasta序列

Question

我想使用以下腳本從大的fasta文件中提取特定的fasta序列，但輸出爲空。

transcripts.txt文件包含我想從assembly.fasta到selected_transcripts.fasta導出的列表轉錄本ID（ID和序列）。 例如：

1. transcripts.txt：

    
   Transcript_00004|5601 
   Transcript_00005|5352

2. assembly.fasta：

    
   >Transcript_00004|5601 
   GATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT 
   >Transcript_00004|5360 
   CGATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT 
   

的ID由>符號開頭：>Transcripts_00004|5601。

我要讀的assembly.fasta文件，如果在assembly.fasta成績單ID是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() 

Answer 1

根據你的例子，有幾個小問題可以解釋你爲什麼沒有得到任何東西。您transcripts.txt在一行中的多個條目，因此my_list將在my_line[0]的first_line的所有項目，在你的循環您可以通過線通過my_list迭代，所以你的第一個項目將是

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

此外，如果assembly.fasta在標題行中沒有>，您將無法取回任何帶有ID和序列的記錄。假設您將>添加到頭中，並且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() 

Reading被改變，因此所有的ID分別追加到my_list。此外，每個項目都被刪除了空白區域，以避免在與record.id相比時在字符串中出現換行符。