The perl stuff is faster, at least if you have a few 100,000 of contigs
Sven
Sven
You are currently viewing the SEQanswers forums as a guest, which limits your access. Click here to register now, and join the discussion
def N50(numlist):
"""
Abstract: Returns the N50 value of the passed list of numbers.
Usage: N50(numlist)
Based on the definition from this SEQanswers post
http://seqanswers.com/forums/showpost.php?p=7496&postcount=4
(modified Broad Institute's definition
https://www.broad.harvard.edu/crd/wiki/index.php/N50)
See SEQanswers threads for details:
http://seqanswers.com/forums/showthread.php?t=2857
http://seqanswers.com/forums/showthread.php?t=2332
"""
numlist.sort(reverse = True)
s = sum(numlist)
limit = s * 0.5
for l in numlist:
s -= l
if s <= limit:
return l
#!/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)
$ grep "^>" 454AllContigs.fna | cut -d"=" -f2 | cut -d" " -f1 | ./stdin_N50.py 386
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, Yesterday, 06:37 PM
|
0 responses
8 views
0 likes
|
Last Post
by seqadmin
Yesterday, 06:37 PM
|
||
Started by seqadmin, Yesterday, 06:07 PM
|
0 responses
8 views
0 likes
|
Last Post
by seqadmin
Yesterday, 06:07 PM
|
||
Started by seqadmin, 03-22-2024, 10:03 AM
|
0 responses
49 views
0 likes
|
Last Post
by seqadmin
03-22-2024, 10:03 AM
|
||
Started by seqadmin, 03-21-2024, 07:32 AM
|
0 responses
66 views
0 likes
|
Last Post
by seqadmin
03-21-2024, 07:32 AM
|
Comment