如何检索标题中包含字符串列表的所有fasta记录

Question

我想要一个脚本，可以从一个文件中提取所有fasta记录，该文件包含标题中的任何字符串列表。

所以，如果我在一个文件中有这样一个列表：

2.A.1.13.5
3.A.1.208.23

像这样的fasta文件：

>gnl|TC-DB|O60645|1.F.2.1.2 Exocyst complex component 3 OS=Homo sapiens GN=EXOC3 PE=1 SV=2
MQCEDSTSFFTMKETDREAVATAVQRVAGMLQRPDQLDKVEQYRRREARKKASVEARLKA
>gnl|TC-DB|O60669|2.A.1.13.5 Monocarboxylate transporter 2 - Homo sapiens (Human).
MPPMPSAPPVHPPPDGGWGWIVVGAAFISIGFSYAFPKAVTVFFKEIQQIFHTTYSEIAW
>gnl|TC-DB|O60683|3.A.20.1.1 Peroxisome biogenesis factor 10 OS=Homo sapiens GN=PEX10 PE=1 SV=1
MAPAAASPPEVIRAAQKDEYYRGGLRSAAGGALHSLAGARKWLEWRKEVELLSDVAYFGL
>gnl|TC-DB|O60684|1.I.1.1.3 Importin subunit alpha-7 OS=Homo sapiens GN=KPNA6 PE=1 SV=1
METMASPGKDNYRMKSYKNNALNPEEMRRRREEEGIQLRKQKREQQLFKRRNVELINEEA
>gnl|TC-DB|O60706|3.A.1.208.23 ATP-binding cassette sub-family C member 9 OS=Homo sapiens GN=ABCC9 PE=1 SV=2
MSLSFCGNNISSYNINDGVLQNSCFVDALNLVPHVFLLFITFPILFIGWGSQSSKVQIHH
>gnl|TC-DB|O60721|3.A.1.208.23 Sodium/potassium/calcium exchanger 1 OS=Homo sapiens GN=SLC24A1 PE=1 SV=1
MGKLIRMGPQERWLLRTKRLHWSRLLFLLGMLIIGSTYQHLRRPRGLSSLWAAVSSHQPI
>gnl|TC-DB|O60779|2.A.1.13.5 Thiamine transporter 1 (THTR-1) (ThTr1) (Thiamine carrier 1) (TC1) - Homo sapiens (Human).
MDVPGPVSRRAAAAAATVLLRTARVRRECWFLPTALLCAYGFFASLRPSEPFLTPYLLGP

然后脚本将打印第2、第5、第6和第7张唱片，如下所示：

>gnl|TC-DB|O60669|2.A.1.13.5 Monocarboxylate transporter 2 - Homo sapiens (Human).
MPPMPSAPPVHPPPDGGWGWIVVGAAFISIGFSYAFPKAVTVFFKEIQQIFHTTYSEIAW
>gnl|TC-DB|O60706|3.A.1.208.23 ATP-binding cassette sub-family C member 9 OS=Homo sapiens GN=ABCC9 PE=1 SV=2
MSLSFCGNNISSYNINDGVLQNSCFVDALNLVPHVFLLFITFPILFIGWGSQSSKVQIHH
>gnl|TC-DB|O60721|3.A.1.208.23 Sodium/potassium/calcium exchanger 1 OS=Homo sapiens GN=SLC24A1 PE=1 SV=1
MGKLIRMGPQERWLLRTKRLHWSRLLFLLGMLIIGSTYQHLRRPRGLSSLWAAVSSHQPI
>gnl|TC-DB|O60779|2.A.1.13.5 Thiamine transporter 1 (THTR-1) (ThTr1) (Thiamine carrier 1) (TC1) - Homo sapiens (Human).
MDVPGPVSRRAAAAAATVLLRTARVRRECWFLPTALLCAYGFFASLRPSEPFLTPYLLGP

这就是我所采取的方法，但我可能在很远的地方，因为我不太清楚我在做什么。有人告诉我，BioPython很适合处理fasta格式的文件，但仍然试图掌握它。

from Bio import SeqIO
import sys

headers = []

search_list = open(sys.argv[1])
for line in search_list:
    headers.append(line.rstrip())
search_list.close()

print search_list

for seq_record in SeqIO.parse(sys.argv[2], "fasta"):
    #print seq_record
    for a in headers:
        head = str(a)
        if head in seq_record:
            print seq_record

提前为任何帮助干杯！

Answer 1

你的代码几乎什么都有。您只需检查头部是否在description中，然后打印SeqRecord元素即可。

from Bio import SeqIO
import sys

headers = list()

with open(sys.argv[1], 'r') as search_list:
    for line in search_list:
        headers.append(line.rstrip())

for seq_record in SeqIO.parse(sys.argv[2], "fasta"):
    for head in headers:
        if head in seq_record.description:
            print(seq_record.format('fasta'))
            break

注意:我在第一个found元素之后添加了一个break语句，如果找到两个标头，它会打印两次。

Answer 2

import re
list_items = '|'.join(list(open('list.txt').read().split('\n')))

with open('fasta.txt') as f:
    for line in f.read().split('\n'):
        if re.search(list_items, line ):
            print line

产出：

gnl_三磷酸腺苷结合盒亚家族C成员9 OS=Homo sapiens GN=ABCC9 PE=1 SV=2 gnl TC-DB 60721 3 A.1.208.23钠/钾/钙交换器1 OS=Homo皂甙GN=SLC24A1 PE=1 SV=1 gnl TC-DB O60779+2.A.1.13.5硫胺素转运体1 (THTR-1) (ThTr1) (硫胺素载体1) (TC1) -人。