How to retrieve sequences from a Fasta file by gene ID

Question

I know this question has been asked a hundred times but I've been at it all day and I can't seem to make this work. I have a fasta file that looks like this ...

>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

and I want to retrieve the sequences that match the gene IDs from a .txt file:

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

Now, I have already tried samtools and FaSomeRecords, but both these methods produce no output. I guess it's because the title also contains transcript IDs (which I can ignore) Do you guys have any suggestions for me? please let me know if you require additional information. Cheers!

Answer 1

Use a Perl one-liner, grep and seqtk subseq to extract the desired fasta sequences:

# 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

Output:

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

Note that seqtk can be installed, for example, using conda.