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

Question

我知道这个问题已经被问过一百次了，但我已经问了一整天了，但我似乎无法解决这个问题。我有一个看起来像这样的 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（我可以忽略）你们对我有什么建议吗？如果您需要更多信息，请告诉我。干杯！

Answer 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。

Answer 2

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

很多年（几十年）以前我会构建像

这样的定义

>gi|123|lcl|xyz ...

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

给你：

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

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

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

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

可能是这个。

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

Answer 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[[=12=]]++} FNR!=NR{if($NF in a){print [=12=]; getline; print}}' a.txt fasta 
>BGI_novel_T018109 Solyc03g080075.1.1
GCAAGGGAAAGAAGTATTACTAGAGGAGATTTTCCCAACAGTTTTCATTTACACACATGGGTTAAGTATTCATAAATAAAAGAGAAAAATCTGTTTATAAGTTGGAGAGTAGTATAAATACAGGAGATTTTCCCAACAGTTTTCATTTATACACATGGGTTAAGTATTCTTAAATAGAAAATCGGAAGTATTATAAATTCTCACTCAAAGAAACCACGTTTGCTCATTTTCGTTATTCCCTTAAAAACATGGGAAGATGAAAGAAAAAAACTAACACATAAAAAGATTGTGAGTTTACTTATTCATGGAGAATTCCCCATTTAAGTTGACAATATTTTTCTATGGTCTTGAACGGCCAGAAAAGTTAATATCCACAACTATTTTCCACTCAATAAGTGTTCCGATACCGTTGAACTTTTTAATATTTTGCACGCCCTTCGTGAAATGTTTTACTCCGTTACTGTCGCGATAATGATGTTTAAAAT
>BGI_novel_T016141 Solyc03g007600.3.1
GCCCAAGTCATAGGTAGTGCCTGTGCGGGTTGACACTCAACATGTGACCGCCACCACATTTTGGCATTTCCCTGAAACTGATAGGTTACAAACTCAATGCCAAATCATTCCACTATGCCCATTTTATGTAGTAACTCATGACAATCAACCAGAAAATCGTAGGCATCCTCAGATTCAGCACCCTTGAAGACTGGAGGTTTCAATTTCAAGAACTTACTGAAAAATTCATGCTGATCACTTGTCATTATAGGCCCTGTAGTCAAACGAGGAAACGTGCCTATTTCCAATGAGGCATCCAT

awk节的解释： FNR==NR{a[[=17=]]++}

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

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

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

Answer 4

由于 FASTA 构造是一个两行对，由一个定义（例如 >BGI_novel_T016697 Solyc03g033550.3.1）和一个包含序列的行组成，您可能可以使用 --after-context 参数到 grep 到 return 匹配的行后跟下一行。 None 您示例中的基因 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到 return 匹配后的下一行：

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

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

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

How to retrieve sequences from a Fasta file by gene ID

linux

bioinformatics

fasta