python - 使用来自单独文本文件的 ID 从 7 GB fasta 文件中提取多个序列

Question

在这个网站上搜索了几个小时并尝试了很多不同的方法都没有用之后，我决定发布我自己的问题。我目前有一个文本文件 (id.txt)，其中包含大约 100 行以下形式的 IDS:

5377-P3-D5-MSITS2a_R1reads1_1125821

5377-P3-D5-MSITS2a_R1reads1_1126992

我有一个 7 GB 的 fasta 文件，其中包含表单中的条目

>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT  

>5377-P3-D5-MSITS2a_R1reads1_1129826 M00532:203:000000000-BKM3D:1:1110:14480:9405 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAAACTCTCGAGGTTACAGCCTTGCTGAATTATTAACCCTTGTCGTTCGCGTACTTCTTGTTTCCTTGGTGTGTTCGCCCACCACAAGTAAAAACATAAACCTTTTGTAA

所有来自 id.text 的 ID 都可以在 seq.fasta 中找到。预期的输出将从 id.text 文件中的 fasta 文件中找到匹配的 ID 号并产生:

>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA  

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

目前，我可以在 bash 中使用 grep 一次从 fasta 文件中提取一个序列，只需从文件中复制并粘贴一个 ID。

例如:grep 5377-P3-D5-MSITS2a_R1reads1_1126992 seq.fasta -A 1

结果:

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

但是，我有多个文本文件，每个文件包含 50-300 个 ID，我想用它们从 FASTA 文件中提取序列，单独提取序列似乎不必要地耗时。我想找到一种方法来查找和输出来自位于单独文本文件中的多个 ID 的 fasta 文件的序列。我主要在 bash 中试验了 awk 和 grep 命令，主要基于该站点上的其他答案，几乎我尝试的每个命令都没有结果，也没有错误消息。

我试过的例子:

awk -F '>' 'NR==FNR{ids[$0]; next} NF>1{f=($2 in ids)}f' id.txt seq.fasta

awk 'NR==FNR{ids[$0];next} /^>/{f=($1 in ids)} f' id.txt seq.fasta

grep -Fwf id.txt seq.fasta

grep -Ff id.txt seq.fasta

我觉得我已经尝试了这两个命令的许多变体(基于其他堆栈溢出和 biostar 的建议)并且在 bash 中，没有任何反应，没有结果或没有错误消息。我在编码方面也是一个相对初学者，所以我无法准确指出出了什么问题。我也对任何可以使用的 python 或其他代码持开放态度。任何帮助或建议将不胜感激。谢谢!

Answer 1

grep 对我来说似乎是最好的主意。我认为您可能需要从搜索字符串中删除 * 字符，因为它们与文件中的内容不匹配。有了这个改变，当我尝试你的提取物时它似乎起作用了:

$ cat fasta 
*>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0   
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA* 

*>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT* 

*>5377-P3-D5-MSITS2a_R1reads1_1129826 M00532:203:000000000-BKM3D:1:1110:14480:9405 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAAACTCTCGAGGTTACAGCCTTGCTGAATTATTAACCCTTGTCGTTCGCGTACTTCTTGTTTCCTTGGTGTGTTCGCCCACCACAAGTAAAAACATAAACCTTTTGTAA*
$ cat ids.txt 
5377-P3-D5-MSITS2a_R1reads1_1125821
5377-P3-D5-MSITS2a_R1reads1_1126992
$ grep -Ff ids.txt fasta 
*>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0   
*>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
$