Is there a good one-liner for calculating the N50 size of a list of contig sizes? Average size is easy with awk but I didn't manage to find an easy way for the N50 size so far. Thanks!
Seqanswers Leaderboard Ad
Collapse
Announcement
Collapse
No announcement yet.
X
-
If you have the mean, you can clearly grep/awk your relevant file for the correct data to compute N50. Why not just write a little generic statistics program that computes descriptive stats for a file or for data passed via stdin? Here's a Python function for median (or you could use numpy) that you could dump into a little script. This would compute the N50...
PHP Code:def median (numlist):
"""
Abstract: Returns the median value of the passed list of numbers.
Usage: median(numlist)
"""
numlist.sort()
# take the mean of the two middle elements if there are an even number
# of elements. otherwise, take the middle element
if len(numlist) % 2 == 0:
medianpos = len(numlist)/2
median = float(numlist[medianpos] + numlist[medianpos-1]) /2
else:
medianpos = len(numlist)/2
median = numlist[medianpos]
return median
Good luck,
Aaron
-
You can call Python scripts from the command line, and pipe stdin and stout too - great for use at the Unix command line.
Are you planning to supply the list of lengths as space separated integers? e.g. "120 2500 745 63 ..."? That would be easy to parse in Python from a string into a list of integers.
Comment
-
It comes up in several occasions, e.g. a file containing
contig1 length=1000 reads=582
contig2 length=1500 reads=900
contig3 length=1400 reads=766
...
I want to do something like:
$ cut -f2 file | cut -d"=" -f2 | sort -nr | awk-or-whatever-for-the-N50-size
The files can have different formats, I just have a sorted list of contig sizes at some point and it would be useful to have a command to get the N50 size.
Comment
-
OK - so the stdin is one integer per line. How about a python script like this,
see also http://seqanswers.com/forums/showthread.php?t=2332
Code:#!/usr/bin/python import sys def N50(numlist): """ Abstract: Returns the N50 value of the passed list of numbers. Usage: N50(numlist) Based on the Broad Institute definition: https://www.broad.harvard.edu/crd/wiki/index.php/N50 """ numlist.sort() newlist = [] for x in numlist : newlist += [x]*x # take the mean of the two middle elements if there are an even number # of elements. otherwise, take the middle element if len(newlist) % 2 == 0: medianpos = len(newlist)/2 return float(newlist[medianpos] + newlist[medianpos-1]) /2 else: medianpos = len(newlist)/2 return newlist[medianpos] assert N50([2, 2, 2, 3, 3, 4, 8, 8]) == 6 lengths = [] for line in sys.stdin : lengths.append(int(line)) print N50(lengths)
Code:$ grep "^>" 454AllContigs.fna | cut -d"=" -f2 | cut -d" " -f1 | ./stdin_N50.py 386
Comment
-
What about the awk experts out there? Any idea for the direct command-line version without script? (you have to store the script somewhere, and if you are using another computer you don't have it at hand... I actually really prefer code in the command line instead of applying a script - if possible..)
Comment
-
Originally posted by seqseq View PostIt comes up in several occasions, e.g. a file containing
contig1 length=1000 reads=582
contig2 length=1500 reads=900
contig3 length=1400 reads=766
...
I want to do something like:
$ cut -f2 file | cut -d"=" -f2 | sort -nr | awk-or-whatever-for-the-N50-size
The files can have different formats, I just have a sorted list of contig sizes at some point and it would be useful to have a command to get the N50 size.
Code:perl -lanF'[\s=]' -e 'push(@contigs,$F[2]);$total+=$F[2];END{foreach(sort{$b<=>$a}@contigs){$sum+=$_;$L=$_;if($sum>=$total*0.5){print "TOTAL: $total\nN50 : $L";exit;} ;}}' file
cheers,
Sven
Comment
-
Here's a stab at a one-liner :
sort -rn contig_lengths.txt | awk 'BEGIN {sum=0} {sum += $1; print $1, sum}' | tac - | awk 'NR==1 {halftot=$2/2} lastsize>halftot && $2<halftot {print} {lastsize=$2}'
For a file "contig_lengths.txt" with one contig size per line, this prints out the length of the N50 contig, along with the number of bases in contigs this size or greater.
Comment
-
Impressive
Assuming no errors, these perl and awk solutions do qualify as "Unix one line" solution, but they are too long to be typed by hand, and cut-and-pasting is asking for errors.
I personally would go a short script (in user's language of choice), probably reading the FASTA file directly.
Comment
-
Originally posted by nhansen View Postsort -rn contig_lengths.txt | awk 'BEGIN {sum=0} {sum += $1; print $1, sum}' | tac - | awk 'NR==1 {halftot=$2/2} lastsize>halftot && $2<halftot {print} {lastsize=$2}'
So in total:
sort -rn contig_lengths.txt | awk '{sum += $0; print $0, sum}' | tac | awk 'NR==1 {halftot=$2/2} lastsize>halftot && $2<halftot {print} {lastsize=$2}'
Comment
Latest Articles
Collapse
-
by seqadmin
Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...-
Channel: Articles
04-04-2024, 04:25 PM -
-
by seqadmin
Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...-
Channel: Articles
03-22-2024, 06:39 AM -
ad_right_rmr
Collapse
News
Collapse
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, 04-11-2024, 12:08 PM
|
0 responses
25 views
0 likes
|
Last Post
by seqadmin
04-11-2024, 12:08 PM
|
||
Started by seqadmin, 04-10-2024, 10:19 PM
|
0 responses
29 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 10:19 PM
|
||
Started by seqadmin, 04-10-2024, 09:21 AM
|
0 responses
24 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 09:21 AM
|
||
Started by seqadmin, 04-04-2024, 09:00 AM
|
0 responses
52 views
0 likes
|
Last Post
by seqadmin
04-04-2024, 09:00 AM
|
Comment