Help with script, sort by coverage

[AdrianP]

Hello everyone,

I am not good with with writing optimal scripts, so I would appreciate someones help, either perl or python is fine.

I have contigs generated by assemblers in the following format:
NODE_2361_length_509_cov_1.43018_ID_745236

This is common to velvet and SPAdes. I was wondering how could I extract all contig names that in their name have a coverage above a user specified value. So I could perl selectAboveCoverage.pl 1.0 and get all above 1.0 including 1.0?

I know I have to first split by >, and then the index 0 will be contig name and index 1 will be the sequence. I need to further split index 0 by _, and index 5 will have the coverage value. I can then somehow sort by that and pick all of those above argv[0], and write them to a new file argv[1].

I know this from experience in TCL, but scripts on those languages take too long to execute.

Any help will be appreciated, thank you!

[rhinoceros]

How about

awk '{FS="_"; if($6>=1) print}' file

[AdrianP]

How about

awk '{FS="_"; if($6>=1) print}' file

You sir, are awesome.

[AdrianP]

How would I use a similar script to sort the contigs by lengths, and tell me what the top 10 or 20 lengths are?

[winsettz]

Could try adding to the above:

sort -r -n file | head

In short:

sort -r -n file

sorts the contents of file, -r to use reverse order (highest to lowest), and -n uses numeric sort. Otherwise you get 1, 10, 100, 2, 20, 200, 3, 30, 300.

head just displays the top 10 of the file, or stdout (lest it dump a pile of numbers into stdout). To adjust the number of lines displayed, use -n. For example,

head -n 20

displays the first 20 lines of a file.

[AdrianP]

Sorry but the above shows contig names, how will sort know what to sort by? Since it's given the full length of the contig?

[winsettz]

Oh, whoops.

Contents of file I wish to process

NODE_2361_length_509_cov_1.43018_ID_745236
NODE_2361_length_509_cov_5.43018_ID_745236
NODE_2361_length_509_cov_3.43018_ID_745236
NODE_2361_length_509_cov_2.43018_ID_745236
NODE_2361_length_509_cov_21.43018_ID_745236

The one-liner

sed 's/_/ /g' filename | awk '{FS=" "; print $6}' | sort -r -n | head -n 10

sed will replace _ with space, then spit out the output to stdout. | to pipe to awk.

awk will use spaces to separate columns, print $6 instructs awk to display everything from column six (which is the numbers after cov, but before ID).
You can then pipe that again to sort, with head to display just what you want to display.

The output

21.43018
5.43018
3.43018
2.43018
1.43018

Obviously not the only solution, and there are probably better ones.

--------------
Edit: If you're dealing with fasta you will have to separate the headers from the sequences. Standard format is >NODE_...; which can be captured as

grep "^>NODE" contigs.fasta > contignames

[rhinoceros]

I would just pipe the awk output to:

sort -r -g -k 4 -t _

When you have an idea what tool might work, but don't know exactly how, just google "man sort" or whatever the name of the tool happens to be. Also, with sort, it's sometimes better to use -g (general numerical value) instead of -n. For example, -n (string numerical value) doesn't work when you have exponents (1e-10 would be smaller than 2e-100 so e.g. sorting blast output by e-value would fail). Actually, I don't know if there's any reason why -g shouldn't be used always by default..