This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract sequence number "X" from multifasta file

Hi there,
I know there are different ways to extract some "wanted" sequences from multifasta file. In general, having one or more IDs, you can retrieve the sequence from multifasta file using bioperl script or another kind of software as ones pointed out here (Extracting Multiple Fasta Sequences At A Time From A File Containing Many Sequences). In my case I know that I want the sequence number 1, 5 and 10 from a multifasta file. Any idea how can I achieve this using perl/bioperl preferably?

Thanks in advance.
P.D: I'm a beginner in bioinformatics word, sorry if it is a very basic question.

fasta
• 4,779 views • link
updated by Ram Ram • written by user230613 user230613

4 answers

Here is a BioPerl example since you asked:

I put this in a script so it is easier to read, but in practice, I would just do this at the command line with a one-liner (that's how I wrote it):

perl -MBio::SeqIO -M'List::MoreUtils qw(any)' -E '@ids = (1, 5, 10); $num = 0; $seqio = Bio::SeqIO->new(-fh => \*STDIN); while ($seq = $seqio->next_seq) { $num++; say join "\n", ">".$seq->id, $seq->seq if any { $_ == $num } @ids; }' < seqs.fas

That will save you from opening an editor. It's not as readable but it works the same. :)

• 0 views • link
written by SES SES

Thank you a lot SES. One last question, could you explain me a little bit the meaning of this last two lines?:

say join "\n", ">".$seq->id, $seq->seq
if any { $_ == $num } @ids;

And another thing, using "say" is possible to print to a file like with "print?

Thanks again.

• 0 views • link
updated by Ram Ram • written by user230613 user230613

say was added in perl version 5.10, it just adds a newline to your print, so it saves you some typing. If you have something as old as 5.10, you can enable it with use feature 'say'; or you can enable all new features at the command line with -E.

join should be self-explanatory, it just adds a newline after the sequence ID and the sequence itself. In Perl, statements are terminated by a semi-colon, so that last statement is one line. That style is called post-fix, with the conditional at the end, and this makes the code much more readable (esp. when put on separate lines) by removing unnecessary braces. That statement should be easy to read, it says, "print the record if any IDs match those in the list." And that sentence is almost working code! That is a good reason to use the functions from List::Util or List::MoreUtils, they are very expressive and they make the code easy to read (and they are fast, as they are all written in C).

Note that I imported the any function from the package List::MoreUtils. All of the functions in that package are now part of List::Utils, which is shipped with Perl (the lastest version includes these functions). I used the List::MoreUtils package so the code will work with any version of Perl, not just the latest. I think it will simplify things going forward now that all the functions are in one package and in core Perl.

• 0 views • link
written by SES SES

There is a program in BBTools which will do that:

getreads.sh in=data.fasta out=filtered.fasta id=1,5,10
• 0 views • link
updated by Ram Ram • written by Brian Bushnell Brian Bushnell

I want the sequence number 1, 5 and 10 from a multifasta file

awk '/^>/ {N++;} {if(N==1 || N==5 ||N==10) print;}' in.fa
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

I might be mistaken, but I think this only prints the lines containing the sequence identifiers, correct?

• 0 views • link
written by Matt Shirley Matt Shirley

No, there is no 'next' statement in

{N++; }

so awk continues to scan the matching patterns. It would only print the headers with

{N++; next;}
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

Thank you Pierre. I'm still impressed by how concise most awk scripts can be.

• 0 views • link
written by Matt Shirley Matt Shirley

I know you are looking for a Perl solution, but you could use my pyfaidx module for this. In the latest version I've added numeric indexing, so you can do:

• 0 views • link
updated by Ram Ram • written by Matt Shirley Matt Shirley

Log in to answer this question.

More posts like this