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
Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.
Somatic Genomics
“We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...-
Channel: Articles
05-24-2024, 01:16 PM -
ad_right_rmr
Collapse
News
Collapse
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, 06-03-2024, 06:55 AM
|
0 responses
12 views
0 likes
|
Last Post
by seqadmin
06-03-2024, 06:55 AM
|
||
Started by seqadmin, 05-30-2024, 03:16 PM
|
0 responses
26 views
0 likes
|
Last Post
by seqadmin
05-30-2024, 03:16 PM
|
||
Comprehensive Sequencing of Great Ape Sex Chromosomes Yields Insights into Evolution and Genetic Variability
by seqadmin
Started by seqadmin, 05-29-2024, 01:32 PM
|
0 responses
29 views
0 likes
|
Last Post
by seqadmin
05-29-2024, 01:32 PM
|
||
Started by seqadmin, 05-24-2024, 07:15 AM
|
0 responses
216 views
0 likes
|
Last Post
by seqadmin
05-24-2024, 07:15 AM
|
Comment