How To Grep Largest Contig From A Multi Fasta File

Question

[HG]

Hi everybody, can any one tell me please how to extract a largest contig from a multi-fasta file ?? using awk or grep ??

Answer 1

[Istvan Albert]

You can use Heng Li's bioawk and samtools:

Then the commands would be

# sort the sequences by length
$ cat w.fasta | bioawk -c fastx '{ print length($seq), $name }' | sort -k1,1rn | head -1
989    HR5V3UP02C00KT

# extract the sequence from the file 
$ samtools faidx w.fasta HR5V3UP02C00KT
>HR5V3UP02C00KT
TCGTACTCGTACGTAGAGGTTCGATCCTAGGGTCCTACGACGGAAGTAAAAACGGCCGGT
CCGGGCCCCGGTTCGACGTCGGACCGTAACCAACGAAAATTGGCCGGTAAAGGGGGTTCC
...

[HG]

Than you so much but could you please check the error

$ cat list4a.fasta | awk -c fastx '{ print length($seq), $name }'| sort -k1, 1rn | head -1
sort: invalid number after `,': invalid count at start of `'

[Istvan Albert]

don't type in the whole thing at first, build it one step at a time, and pipe it through a pager, that way you will notice potential errors

[HG]

Dear Istvan, I am doing according to you suggestion:

cat list4a.fasta           #working fine 
awk -c fastx '{ print length($seq), $name }'

Here is the problem

Usage: awk [POSIX or GNU style options] -f progfile [--] file ...
Usage: awk [POSIX or GNU style options] [--] 'program' file ...

could you please tell me waht is the error ??

[Istvan Albert]

did you install bioawk? that's is the version of awk you will need to use.

[HG]

Hello Istvan, I have downloaded bioawk but when i am trying to install its showing

hiren@FB11-10207:~/Desktop/spades/bioawk-master $ make 
make: `awk' is up to date.

But once i am checking bioawk like

hiren@FB11-10207:~/Desktop/spades/bioawk-master $ bioawk
bioawk: command not found

Could you please let me know any solution ????

[Pierre Lindenbaum]

There is no space after -k1,

[FadyNabil]

can you explain this code please?

Answer 2

[Prakki Rama]

cat seqs_oneline.fasta | perl -e 'while (<>) {$h=$_; $s=<>; $seqs{$h}=$s;} foreach $header (reverse sort {length($seqs{$a}) <=> length($seqs{$b})} keys %seqs) {print $header.$seqs{$header}}' | head -2

Source: sort fasta