/usr/lib/python2.7/dist-packages/pbsuite/utils/quickN50.py is in python-pbsuite-utils 15.8.24+dfsg-2.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 | #!/usr/bin/python
import sys, math
def getStats(seqLengths):
data = {}
seqLengths.sort(reverse=True)
data["numItems"] = len(seqLengths)
data["itemSum"] = sum(seqLengths)
tl = data["itemSum"]
n50_mark = data["itemSum"] * .5
n90_mark = data["itemSum"] * .90
n95_mark = data["itemSum"] * .95
data["n50"] = None
data["n50_gt_count"] = None
data["n90"] = None
data["n90_gt_count"] = None
data["n95"] = None
data["n95_gt_count"] = None
basesSeen = 0
for pos,n in enumerate(seqLengths):
basesSeen += n
if data["n50"] is None and basesSeen > n50_mark:
data["n50"] = n
data["n50_gt_count"] = pos
if data["n90"] is None and basesSeen > n90_mark:
data["n90"] = n
data["n90_gt_count"] = pos
if data["n95"] is None and basesSeen > n95_mark:
data["n95"] = n
data["n95_gt_count"] = pos
break
#may not have gaps
if data["numItems"] == 0:
return data
data["min"] = seqLengths[-1]
data["FstQu"] = seqLengths[ int(math.floor(data["numItems"]*.75)) ]
median = data["numItems"]*.50
data["median"] = int( (seqLengths[ int(math.floor(median)) ] + \
seqLengths[ int(math.floor(median)) ]) / 2)
data["mean"] = data["itemSum"]/data["numItems"]
data["TrdQu"] = seqLengths[ int(math.floor(data["numItems"]*.25)) ]
data["max"] = seqLengths[0]
return data
def run(data):
"""
list of numbers - can be a string if you want
"""
data = map(float, data)
ret = getStats(data)
outputOrder = ["itemSum",
"numItems",
"min",
"FstQu",
"mean",
"median",
"n50",
"n50_gt_count",
"TrdQu",
"n90",
"n90_gt_count",
"n95",
"n95_gt_count",
"max"]
for key in outputOrder:
print "{0}\t{1:.2f}".format(key, ret[key])
if __name__ == '__main__':
run(sys.stdin.read().strip().split('\n'))
|