Sequences Collection Based On Length

Question

[Mat]

Hello, I made my own database of protein sequences with formatdb. I would like to know if it is possible and how to retrieve all the sequences of a certain length (e.g all the sequences smaller than 200 aa).

Thanks

Matteo

Answer 1

[Neilfws]

What Pierre said. Assuming you have dumped to fasta file myseqs.fa, here's a Bioperl approach:

#!/usr/bin/perl -w

use strict;
use Bio::SeqIO;

my $inseq  = Bio::SeqIO->new(-file => "myseqs.fa", -format => "fasta");
my $outseq = Bio::SeqIO->new(-file => ">myseqs200.fa", -format => "fasta");

# write to file myseqs200.fa
while(my $seq = $inseq->next_seq) {
  if($seq->length <= 200) {
    $outseq->write_seq($seq);
  }
}

Answer 2

[Pierre Lindenbaum]

Use fastacmd to dump your database as fasta (option -D 1 ) and filter the result using your favorite tool (perl, awk, etc... )