sorting a multiFASTA file by DNA length

802 views Asked by At

I'm trying to sort a multiFASTA file by length. I have the alphabetical sort figured out but I can't seem to get the numerical sort. The output should be a sorted multiFASTA file. This is an option to another program. Here is the code.

sub sort {
my $length;
my $key;
my $id;
my %seqs;
my $seq;
my $action = shift;
my $match = $opts{$action};
$match =~ /[l|id]/ || die "not the right parameters\n";
my $in = Bio::SeqIO->new(-file=>"$filename", -format=>'fasta');
    while(my $seqobj = $in->next_seq()){
        my $id = $seqobj->display_id();
        my $length=$seqobj->length();
        #$seq =~s/.{1,60}\K/\n/sg;
        $seqs{$id} = $seqobj, unless $match eq 'l';
        $seqs{$length}=$seqobj, unless $match eq 'id';
    }
    if($match eq 'id'){
        foreach my $id (sort keys %seqs) {
             printf ">%-9s \n%-s\n", $id, $seqs{$id}->seq;
        }
    }
    elsif($match eq 'l'){
        foreach my $length ( sort keys %seqs){
             printf "%-10s\n%-s\n",$length, $seqs{$length}->seq;
        }
    }
}
3

There are 3 answers

5
choroba On

To sort numerically, you must provide the comparing subroutine:

sort { $a <=> $b } keys %seqs

Are you sure no two sequences can have the same length? $seqs{$length}=$seqobj overwrites the previously stored value.

0
Pierre On

A one lineer: use awk to linearize. a second awk to add a column containing the length, sort on this column, remove the column, restore the fasta sequence.

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}'  input.fa  |\
awk -F '\t' '{printf("%d\t%s\n",length($2),$0);}' |\
sort -t $'\t' -k1,1n |\
cut -f 2- |\
tr "\t" "\n"

PS: for bioinformartics questions, you should use https://www.biostars.org/, or https://bioinformatics.stackexchange.com/, etc...

0
Anicet Ebou On

You can use pyfaidx or just take a look at jim hester repos. But as @pierre said above you should ask you question on biostars for example. The answer on biostars can be found here.