Mercurial > repos > bebatut > combine_metaphlan2_humann2
comparison combine_metaphlan_humann.py @ 3:01ac9954c27f draft
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/combine_metaphlan2_humann2 commit b84cbcbffa25c55acd1a31df76e4a4f78772cbd7
author | bgruening |
---|---|
date | Thu, 20 Jul 2023 10:07:12 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
2:fdfb35745104 | 3:01ac9954c27f |
---|---|
1 #!/usr/bin/env python | |
2 # -*- coding: utf-8 -*- | |
3 | |
4 import argparse | |
5 | |
6 | |
7 def extract_clade_abundance(metaphlan_fp): | |
8 clade_abund = {} | |
9 with open(metaphlan_fp, "r") as metaphlan_f: | |
10 is_metaphlan_v4 = False | |
11 for line in metaphlan_f.readlines(): | |
12 if 'SGB' in line: | |
13 # New versions of metaphlan against a recent DB contain a header line with DB name, which contains SGB | |
14 is_metaphlan_v4 = True | |
15 if line.find("g__") == -1: | |
16 continue | |
17 | |
18 split_line = line[:-1].split("\t") | |
19 taxo = split_line[0] | |
20 if is_metaphlan_v4: | |
21 # Column order in new metaphlan versions: | |
22 # clade_name NCBI_tax_id relative_abundance additional_species | |
23 abundance = split_line[2] | |
24 else: | |
25 # Column order in the old metaphlan versions: | |
26 # clade_name relative_abundance coverage average_genome_length_in_the_clade estimated_number_of_reads_from_the_clade | |
27 abundance = split_line[1] | |
28 | |
29 genus = taxo[(taxo.find("g__") + 3):] | |
30 if genus.find("|") != -1: | |
31 genus = genus[: (genus.find("|"))] | |
32 clade_abund.setdefault(genus, {"abundance": 0, "species": {}}) | |
33 if taxo.find("t__") != -1: | |
34 continue | |
35 elif taxo.find("s__") != -1: | |
36 species = taxo[(taxo.find("s__") + 3):] | |
37 clade_abund[genus]["species"].setdefault(species, abundance) | |
38 else: | |
39 clade_abund[genus]["abundance"] = abundance | |
40 return clade_abund | |
41 | |
42 | |
43 def compute_overall_abundance(humann_fp): | |
44 overall_abundance = 0 | |
45 with open(humann_fp, "r") as humann_f: | |
46 for line in humann_f.readlines(): | |
47 if line.find("|") != -1 or line.startswith("#"): | |
48 continue | |
49 split_line = line[:-1].split("\t") | |
50 overall_abundance += float(split_line[1]) | |
51 return overall_abundance | |
52 | |
53 | |
54 def format_characteristic_name(name): | |
55 formatted_n = name | |
56 formatted_n = formatted_n.replace("/", " ") | |
57 formatted_n = formatted_n.replace("-", " ") | |
58 formatted_n = formatted_n.replace("'", "") | |
59 if formatted_n.find("(") != -1 and formatted_n.find(")") != -1: | |
60 open_bracket = formatted_n.find("(") | |
61 close_bracket = formatted_n.find(")") + 1 | |
62 formatted_n = formatted_n[:open_bracket] + formatted_n[close_bracket:] | |
63 return formatted_n | |
64 | |
65 | |
66 def combine_metaphlan_humann(args): | |
67 clade_abund = extract_clade_abundance(args.metaphlan_fp) | |
68 overall_abund = compute_overall_abundance(args.humann_fp) | |
69 | |
70 with open(args.output_fp, "w") as output_f: | |
71 s = "genus\tgenus_abundance\tspecies\tspecies_abundance\t" | |
72 s = "%s\t%s_id\t%s_name\t%s_abundance\n" % (s, args.type, args.type, args.type) | |
73 output_f.write(s) | |
74 with open(args.humann_fp, "r") as humann_f: | |
75 for line in humann_f.readlines(): | |
76 if line.find("|") == -1: | |
77 continue | |
78 | |
79 split_line = line[:-1].split("\t") | |
80 abundance = 100 * float(split_line[1]) / overall_abund | |
81 annotation = split_line[0].split("|") | |
82 charact = annotation[0].split(":") | |
83 charact_id = charact[0] | |
84 char_name = "" | |
85 if len(charact) > 1: | |
86 char_name = format_characteristic_name(charact[-1]) | |
87 taxo = annotation[1].split(".") | |
88 | |
89 if taxo[0] == "unclassified": | |
90 continue | |
91 genus = taxo[0][3:] | |
92 species = taxo[1][3:] | |
93 | |
94 if genus not in clade_abund: | |
95 print("no %s found in %s" % (genus, args.metaphlan_fp)) | |
96 continue | |
97 if species not in clade_abund[genus]["species"]: | |
98 print( | |
99 "No %s found in %s for % s" | |
100 % (species, args.metaphlan_fp, genus) | |
101 ) | |
102 continue | |
103 | |
104 s = "%s\t%s\t" % (genus, clade_abund[genus]["abundance"]) | |
105 s += "%s\t%s\t" % (species, clade_abund[genus]["species"][species]) | |
106 s += "%s\t%s\t%s\n" % (charact_id, char_name, abundance) | |
107 output_f.write(s) | |
108 | |
109 | |
110 if __name__ == "__main__": | |
111 parser = argparse.ArgumentParser() | |
112 parser.add_argument("--humann_fp", required=True) | |
113 parser.add_argument("--metaphlan_fp", required=True) | |
114 parser.add_argument("--output_fp", required=True) | |
115 parser.add_argument("--type", required=True, choices=["gene_families", "pathways"]) | |
116 args = parser.parse_args() | |
117 | |
118 combine_metaphlan_humann(args) |