This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How do I extract Fasta Sequences based on a list of IDs?

I have a fasta sequence file that looks like this: 

>sp|P31946|1433B_HUMAN 14-3-3 protein beta/alpha OS=Homo sapiens GN=YWHAB PE=1 SV=3
MTMDKSELVQKAKLAEQAERYDDMAAAMKAVTEQGHELSNEERNLLSVAYKNVVGARRSS
WRVISSIEQKTERNEKKQQMGKEYREKIEAELQDICNDVLELLDKYLIPNATQPESKVFY
LKMKGDYFRYLSEVASGDNKQTTVSNSQQAYQEAFEISKKEMQPTHPIRLGLALNFSVFY
YEILNSPEKACSLAKTAFDEAIAELDTLNEESYKDSTLIMQLLRDNLTLWTSENQGDEGD
AGEGEN
>sp|P62258|1433E_HUMAN 14-3-3 protein epsilon OS=Homo sapiens GN=YWHAE PE=1 SV=1
MDDREDLVYQAKLAEQAERYDEMVESMKKVAGMDVELTVEERNLLSVAYKNVIGARRASW
RIISSIEQKEENKGGEDKLKMIREYRQMVETELKLICCDILDVLDKHLIPAANTGESKVF
YYKMKGDYHRYLAEFATGNDRKEAAENSLVAYKAASDIAMTELPPTHPIRLGLALNFSVF
YYEILNSPDRACRLAKAAFDDAIAELDTLSEESYKDSTLIMQLLRDNLTLWTSDMQGDGE
EQNKEALQDVEDENQ 
 .........

I have another file that has a subset of headers for this file: 

>sp|P31946|1433B_HUMAN 14-3-3 protein beta/alpha OS=Homo sapiens GN=YWHAB PE=1 SV=3
>sp|P62258|1433E_HUMAN 14-3-3 protein epsilon OS=Homo sapiens GN=YWHAE PE=1 SV=1
.......

How can I use the header file to extract out my fasta sequences only for these headers?

Thank you in advance!

grep sequence
• 11,697 views • link
updated by Joe Joe • written by a.rex a.rex

2 answers

This is one of FAQ's on biostars. You should search for additional threads but here are some to get you started.

Extracting specific IDs + sequence from multifasta
extract sequences based on ids file
Extracting specific sequences from a big fasta file using ids of the sequences to be excluded

My personal favorite is a program from Jim Kent/UCSC:

faSomeRecords

Download the linux version linked (macOS available elsewhere on that site). Add execute permissions chmod a+x faSomeRecords.

 ./faSomeRecords 
faSomeRecords - Extract multiple fa records
usage:
   faSomeRecords in.fa listFile out.fa
options:
   -exclude - output sequences not in the list file.

NOTE: You should remove > from your header list file when using it as an input in place of listFile.

• 0 views • link
written by GenoMax GenoMax

Dear, is there a R way to do the same thing?

• 0 views • link
written by lessismore lessismore

Great handy tool for massive datasets.

• 0 views • link
written by sunnykevin97 sunnykevin97

You can use this:

while read line; do grep -A 5 "$line" fastaFile; done <listFile

-A in grep means lines-after the match.

• 0 views • link
written by Rashedul Islam Rashedul Islam

This will only work if the sequences are all 5 lines long which is not a generally safe assumption.

• 0 views • link
written by Joe Joe

In that case you can extract fasta sequences by ids using seqkit

while read line; do less fastaFile.fasta | seqkit grep -p $line; done < listFile >>out.fasta

• 0 views • link
written by Rashedul Islam Rashedul Islam

Dont just use commands at random from the web. Spend some time understanding why the solution might not be working for you. You have good suggestions in the other thread. This code is not well suited to the task you describe.

• 0 views • link
written by Joe Joe

Log in to answer this question.

More posts like this