如何通过基因 ID 从 Fasta 文件中检索序列答案

【问题标题】：How to retrieve sequences from a Fasta file by gene ID如何通过基因 ID 从 Fasta 文件中检索序列
【发布时间】：2021-02-14 15:38:05
【问题描述】：

我知道这个问题已经被问了一百次，但我整天都在做这个问题，我似乎无法完成这个工作。我有一个看起来像这样的 fasta 文件 ...

>BGI_novel_T016697 Solyc03g033550.3.1
CTGACGTATACAATTAAGCCGCGAAAAATCTACTTTTTTTTTAATAGATATGAATTTCTTTTGTTTCGTATAATGAAGTATTTGTTCCAACAATGTTTAATTATTAGGCATTTGGAATGTGATGGGGCAACTAACAAAGAAGCCAATATCAACATCAATTAACAAACATATGATATAATTCTAGTGAAGTGAAAGCCAAGATATGAAACTCTCCACCCACACTATCTTAAATGATCTTTTTTAAAACATTCTAATTAGGTGATAACTAAAAGCAATAATTCTACCAATTTTGAAACAAACAATATGGTCCC
>BGI_novel_T016313 Solyc03g025570.2.1
TTCAAGTGTTAGTTTCACATCATCACGTTTGGACCTACGTTTCTATATTAGAACATATTCTAACTGATCTCTAGCTGTTATTCATGGGATTGTAAGAAATTTGTATCCCTCTCCGGATTTTACTTTGATCGCCACAAAATGAACATATGCTTTCAATTTTCTATGATGAAAAATCAGCCTCTCTCAATATTGGGTTTAAA
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT
>BGI_novel_T016817 BGI_novel_G001220
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCATG
>BGI_novel_T016141 Solyc03g007600.3.1

我想从 .txt 文件中检索与基因 ID 匹配的序列：

Solyc00g256710.2.1
Solyc01g010890.3.1
Solyc01g056990.3.1
Solyc01g060050.2.1
Solyc01g081120.2.1
Solyc01g097740.3.1
Solyc01g098180.3.1
Solyc01g102320.1.1
Solyc01g106420.3.1
Solyc01g111580.3.1
Solyc01g111970.3.1
Solyc02g005530.2.1
Solyc02g031780.1.1
Solyc02g064595.1.1
Solyc02g081920.3.1
Solyc02g084220.3.1

现在，我已经尝试过 samtools 和 FaSomeRecords，但这两种方法都没有输出。我想这是因为标题还包含成绩单 ID（我可以忽略）你们对我有什么建议吗？如果您需要更多信息，请告诉我。干杯！

【问题讨论】：

标签： linux bioinformatics fasta

【解决方案1】：

使用 Perl 单行代码 grep 和 seqtk subseq 提取所需的 fasta 序列：

# Create test input:

cat > in.fasta <<EOF
>BGI_novel_T016697 Solyc03g033550.3.1
CTGACGTATACAATTAAGCCGCG
>BGI_novel_T016313 Solyc03g025570.2.1
TTCAAGTGTTAGTTTCACATCAT
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAG
>BGI_novel_T016817 BGI_novel_G001220
GCCCAAGTCATAGGTAGTGCCTG
>BGI_novel_T016141 Solyc03g007600.3.1
ACGTACGTACGTACGTACGTACG
EOF

cat > gene_ids.txt <<EOF
Solyc03g033550.3.1
Solyc03g080075.1.1
Solyc00g256710.2.1
Solyc01g010890.3.1
EOF

# Extract ids and gene ids into a tsv file:
perl -lne '@f = /^>(\S+)\s+(\S+)/ and print join "\t", @f;' in.fasta > ids_gene_ids.tsv

# Select ids that correspond to the desired gene ids:
grep -f gene_ids.txt ids_gene_ids.tsv | cut -f1 > ids.selected.txt

# Extract fasta sequence that correspond to desired gene ids:
seqtk subseq in.fasta ids.selected.txt > out.fasta                

cat out.fasta

输出：

>BGI_novel_T016697 Solyc03g033550.3.1
CTGACGTATACAATTAAGCCGCG
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAG

注意seqtk 可以安装，例如使用conda。

【讨论】：

【解决方案2】：

对我来说已经有一段时间了，但是defline
（fasta 记录中以> 开头的第一行）需要特别的格式才能被索引序列检索。

许多年（几十年）前，我会构建类似的定义

>gi|123|lcl|xyz ...

因此可以使用 'xyz' 检索序列，其中 'lcl' 是内置的爆炸“本地”命名空间，与默认索引的命名空间一起使用。

给你：

https://en.wikipedia.org/wiki/FASTA_format

将您的定义强制转换为受支持的可索引格式
重新索引，然后享受您非常快速的序列检索。

如果您的序列数据库很小，这将无济于事并且不值得索引工作。

在这种情况下，我只会使用我在某处网上找到的fasta_grep perl 脚本...

可能是这个。

http://web.mit.edu/meme_v4.11.4/share/doc/fasta-grep.html

【讨论】：

【解决方案3】：

还有一个小的awk 替代方案：

fasta 文件如下所示：

$ cat fasta 
>BGI_novel_T016697 Solyc03g033550.3.1
CTGACGTATACAATTAAGCCGCGAAAAATCTACTTTTTTTTTAATAGATATGAATTTCTTTTGTTTCGTATAATGAAGTATTTGTTCCAACAATGTTTAATTATTAGGCATTTGGAATGTGATGGGGCAACTAACAAAGAAGCCAATATCAACATCAATTAACAAACATATGATATAATTCTAGTGAAGTGAAAGCCAAGATATGAAACTCTCCACCCACACTATCTTAAATGATCTTTTTTAAAACATTCTAATTAGGTGATAACTAAAAGCAATAATTCTACCAATTTTGAAACAAACAATATGGTCCC
>BGI_novel_T016313 Solyc03g025570.2.1
TTCAAGTGTTAGTTTCACATCATCACGTTTGGACCTACGTTTCTATATTAGAACATATTCTAACTGATCTCTAGCTGTTATTCATGGGATTGTAAGAAATTTGTATCCCTCTCCGGATTTTACTTTGATCGCCACAAAATGAACATATGCTTTCAATTTTCTATGATGAAAAATCAGCCTCTCTCAATATTGGGTTTAAA
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT
>BGI_novel_T016817 BGI_novel_G001220
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCATG
>BGI_novel_T016141 Solyc03g007600.3.1
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCATG

ID文件a.txt这样：

$ cat a.txt 
Solyc00g256710.2.1
Solyc01g010890.3.1
Solyc01g056990.3.1
Solyc01g060050.2.1
Solyc01g081120.2.1
Solyc01g097740.3.1
Solyc01g098180.3.1
Solyc01g102320.1.1
Solyc01g106420.3.1
Solyc01g111580.3.1
Solyc01g111970.3.1
Solyc02g005530.2.1
Solyc02g031780.1.1
Solyc02g064595.1.1
Solyc02g081920.3.1
Solyc02g084220.3.1
Solyc03g080075.1.1
Solyc03g007600.3.1

我们可以这样做：

$ awk 'FNR==NR{a[$0]++} FNR!=NR{if($NF in a){print $0; getline; print}}' a.txt fasta 
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT
>BGI_novel_T016141 Solyc03g007600.3.1
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCAT

awk-stanza 的解释： FNR==NR{a[$0]++}

当我们对第一个文件 ID 进行操作时，我们将每个 ID 添加到一个数组中；

FNR!=NR{if($NF in a){print $0; getline; print}}'

当我们到达快速文件时，我们获取每行的最后一个字段$NF 并检查它是否存在于数组中；如果是，我们打印该行和它后面的行。

【讨论】：

【解决方案4】：

由于 FASTA 构造是由定义线（例如 >BGI_novel_T016697 Solyc03g033550.3.1）后跟包含序列的行组成的两行对，您可能可以使用 --after-context 参数到 grep 来返回匹配的行由下一行。您的示例中的所有基因 ID 都不在您提供的示例 FASTA 文件中，但如果我将匹配的基因 ID 添加到该文件进行测试，您可以看到它是如何工作的。

基因查找文件：

$ cat lookup.txt
cat lookup.txt
Solyc00g256710.2.1
Solyc01g010890.3.1
Solyc01g056990.3.1
Solyc01g060050.2.1
Solyc01g081120.2.1
Solyc01g097740.3.1
Solyc01g098180.3.1
Solyc01g102320.1.1
Solyc01g106420.3.1
Solyc01g111580.3.1
Solyc01g111970.3.1
Solyc02g005530.2.1
Solyc02g031780.1.1
Solyc02g064595.1.1
Solyc02g081920.3.1
Solyc02g084220.3.1
Solyc03g080075.1.1

测序 FASTA 文件：

$ cat seqs.fa
>BGI_novel_T016697 Solyc03g033550.3.1
CTGACGTATACAATTAAGCCGCGAAAAATCTACTTTTTTTTTAATAGATATGAATTTCTTTTGTTTCGTATAATGAAGTATTTGTTCCAACAATGTTTAATTATTAGGCATTTGGAATGTGATGGGGCAACTAACAAAGAAGCCAATATCAACATCAATTAACAAACATATGATATAATTCTAGTGAAGTGAAAGCCAAGATATGAAACTCTCCACCCACACTATCTTAAATGATCTTTTTTAAAACATTCTAATTAGGTGATAACTAAAAGCAATAATTCTACCAATTTTGAAACAAACAATATGGTCCC
>BGI_novel_T016313 Solyc03g025570.2.1
TTCAAGTGTTAGTTTCACATCATCACGTTTGGACCTACGTTTCTATATTAGAACATATTCTAACTGATCTCTAGCTGTTATTCATGGGATTGTAAGAAATTTGTATCCCTCTCCGGATTTTACTTTGATCGCCACAAAATGAACATATGCTTTCAATTTTCTATGATGAAAAATCAGCCTCTCTCAATATTGGGTTTAAA
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT
>BGI_novel_T016817 BGI_novel_G001220
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCATG

带有-w 的Grep 命令用于整个单词匹配（以防您可能得到部分结果），-f 使用术语/模式文件进行搜索，-A1 在匹配后返回下一行：

$ grep -w -f lookup.txt -A1 seqs.fa
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT

除非你有一些复杂的模式匹配和/或一些需要完成的后处理，我认为一个简单的 grep 应该让你更简单一些。

【讨论】：