-
Notifications
You must be signed in to change notification settings - Fork 81
/
get_failed.py
55 lines (41 loc) · 1.49 KB
/
get_failed.py
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
#!/usr/bin/env python
from sys import argv
from americangut.util import parse_mapping_file
from biom import load_table
__author__ = "Daniel McDonald"
__copyright__ = "Copyright 2015, The American Gut Project"
__credits__ = ["Daniel McDonald"]
__license__ = "BSD"
__version__ = "unversioned"
__maintainer__ = "Daniel McDonald"
__email__ = "mcdonadt@colorado.edu"
if __name__ == '__main__':
if len(argv) != 4:
print "usage: python %s mappingfile table threshold" % argv[0]
exit(1)
header, mapping_file = parse_mapping_file(open(argv[1]))
table = load_table(argv[2])
threshold = float(argv[3])
map_ids = {i[0].split('.')[0] for i in mapping_file}
table_ids = {i.split('.')[0] for i in table.ids()}
missing = {i for i in (map_ids - table_ids) if 'blank' not in i.lower()}
counts = table.sum(axis='sample')
below_threshold = set()
above_threshold = set()
for id_, count in zip([i.split('.')[0] for i in table.ids()], counts):
if 'blank' in id_.lower():
continue
if count < threshold:
below_threshold.add(id_)
else:
above_threshold.add(id_)
if missing:
print "The following IDs do not have any reads:"
for id_ in sorted(missing):
print id_
print
# remove any resequenced samples
if below_threshold - above_threshold:
print "The following have fewer than %d reads" % threshold
for id_ in sorted(below_threshold):
print id_