-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathstampa_merge.py
117 lines (101 loc) · 4.36 KB
/
stampa_merge.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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Parse stampa results and compute last common ancestor.
"""
__author__ = "Frédéric Mahé <[email protected]>"
__date__ = "2019/09/20"
__version__ = "$Revision: 2.0"
import os
import sys
# *************************************************************************** #
# #
# Functions #
# #
# *************************************************************************** #
def last_common_ancestor(taxonomies):
"""Compute last common ancestor"""
lca = list()
if len(taxonomies) > 1:
zipped = list(zip(*taxonomies))
for level in zipped:
if len(set(level)) > 1:
level = "*"
else:
level = level[0]
lca.append(level)
else: # only one top hit
try:
lca = taxonomies[0]
except IndexError:
print(taxonomies, file=sys.stderr)
raise
return lca
def main():
"""Parse stampa results and compute last common ancestor."""
# Parse command line options and change working directory
directory = sys.argv[1]
if not os.path.exists(directory):
sys.exit("ERROR: directory %s not found!" % directory)
os.chdir(directory)
# List files
files = [f for f in os.listdir(directory) if f.startswith("hits.")]
files.sort()
# detect empty files (something went wrong upstream)
empty_files = [f for f in files if not os.stat(f).st_size]
if empty_files:
print("ERROR: empty hits file detected", file=sys.stderr)
for f in empty_files:
print(f, file=sys.stderr)
sys.exit(1)
# Parse files
for input_file in files:
previous = ("", "", "")
taxonomies = list()
accessions = list()
output_file = input_file.replace("hits.", "results.")
with open(input_file, "r") as input_file:
with open(output_file, "w") as output_file:
for line in input_file:
amplicon, identity, hit = line.strip().split("\t")
amplicon, abundance = amplicon.split("_")
if hit != "*":
try:
accession, taxonomy = hit.split(" ", 1)
except ValueError: # something's wrong
print(amplicon, identity, hit, file=sys.stderr)
raise
else: # no hit
accession = taxonomy = "No_hit"
taxonomy = taxonomy.split("|")
if previous[0] == amplicon:
taxonomies.append(taxonomy)
accessions.append(accession)
elif previous[0] == "": # deal with first item
taxonomies.append(taxonomy)
accessions.append(accession)
previous = (amplicon, abundance, identity)
elif previous[0] != amplicon:
# flush
lca = last_common_ancestor(taxonomies)
print("\t".join(previous), "|".join(lca),
",".join(accessions), sep="\t", file=output_file)
# reinitialize
taxonomies = list()
accessions = list()
taxonomies.append(taxonomy)
accessions.append(accession)
previous = (amplicon, abundance, identity)
# Deal with end of file
lca = last_common_ancestor(taxonomies)
print("\t".join(previous), "|".join(lca),
",".join(accessions), sep="\t", file=output_file)
return
# *************************************************************************** #
# #
# Body #
# #
# *************************************************************************** #
if __name__ == '__main__':
main()
sys.exit(0)