Как выполнить выборку файла FASTA на основе заголовков, если заголовки содержат определенные строки?

Question

У меня есть такой файл:

>gi|373248686|emb|HE586118.1| Streptomyces albus subsp. albus salinomycin biosynthesis cluster, strain DSM 41398
GGATGCGAAGGACGCGCTGCGCAAGGCGCTGTCGATGGGTGCGGACAAGGGCATCCACGT
CGAGGACGACGATCTGCACGGCACCGACGCCGTGGGTACCTCGCTGGTGCTGGCCAAGGC
>gi|1139489917|gb|KX622588.1| Hyalangium minutum strain DSM 14724 myxochromide D subtype 1 biosynthetic gene cluster and tRNA-Thr gene, complete sequence
ATGCGCAAGCTCGTCATCACGGTGGGGATTCTGGTGGGGTTGGGGCTCGTGGTCCTTTGG
TTCTGGAGCCCGGGAGGCCCAGTCCCCTCCACGGACACGGAGGGGGAAGGGCGGAGTCAG
CGCCGGCAGGCCATGGCCCGGCCCGGCTCCGCGCAGCTGGAGAGTCCCGAGGACATGGGG
>gi|930076459|gb|KR364704.1| Streptomyces sioyaensis strain BCCO10_981 putative annimycin-type biosynthetic gene cluster, partial sequence
GCCGGCAGGTGGGCCGCGGTCAGCTTCAGGACCGTGGCCGTCGCGCCCGCCAGCACCACG
GAGGCCCCCACGGCCAGCGCCGGGCCCGTGCCCGTGCCGTACGCGAGGTCCGTGCTGAAC

, и у меня есть текстовый файл, содержащий список чисел:

373248686
930076459
296280703
......

Я хочу сделать следующее:

if the header in the fasta file contains the numbers in the text file:
        write all the matches(header+sequence) to a new output.fasta file.

Как это сделать в python?Это кажется простым, некоторые циклы for могут выполнять эту работу, но почему-то я не могу этого сделать, и если мои файлы действительно большие, цикл в другом цикле может занять много времени.Вот что я попробовал:

from Bio import SeqIO                                                               
import sys                                                                          

wanted = []
for line in open(sys.argv[2]):
    titles = line.strip() 
    wanted.append(titles)


seqiter = SeqIO.parse(open(sys.argv[1]), 'fasta')      
sys.stdout = open('output.fasta', 'w')                               

new_seq = []

for seq in seqiter:
    new_seq.append(seq if i in seq.id for i in wanted)


SeqIO.write(new_seq, sys.stdout, "fasta")
sys.stdout.close()            

Получил эту ошибку:

new_seq.append(seq if i in seq.id for i in wanted)
                                        ^
SyntaxError: invalid syntax

Есть ли лучший способ сделать это?

Спасибо!

Answer 1

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

keepers = {}
with open("ids.txt", "r") as id_handle:
    for curr_id in id_handle:
        keepers[curr_id] = True

Понимание списка генерируетсписок, поэтому вам не нужно добавлять его в другой список.

keeper_seqs = [x for x in seqiter if x.id in keepers]

При работе с большими файлами вам нужно будет перебирать seqiter и записывать записи по одному, чтобы избежать проблем с памятью.

Вы также никогда не должны присваивать sys.stdout без веской причины, если вы хотите вывести на STDOUT, просто используйте print или sys.stdout.write().

Answer 2

Используйте программу, подобную этой

from Bio import SeqIO
import sys

# read in the text file
numbersInTxtFile = set()
# hint: use with, then you don't need to
# program file closing. Furhtermore error
# handling is comming along with this too
with open(sys.argv[2], "r") as inF:
    for line in inF:
        line = line.strip()
        if line == "": continue
        numbersInTxtFile.add(int(line))

# read in the fasta file
with open(sys.argv[1], "r") as inF:
    for record in SeqIO.parse(inF, "fasta"):
        # now check if this record in the fasta file 
        # has an id we are searching for
        name = record.description
        id = int(name.split("|")[1])
        print(id, numbersInTxtFile, id in numbersInTxtFile)
        if id in numbersInTxtFile:   
            # we need to output
            print(">%s" % name)
            print(record.seq)

, который затем можно вызывать из командной строки

python3 nameOfProg.py inputFastaFile.fa idsToSearch.txt > outputFastaFile.fa