· 9 years ago · Nov 30, 2016, 01:14 PM
1#!/usr/bin/env python
2
3__author__ = "William Walters"
4__copyright__ = "NA"
5__credits__ = ["William Walters"]
6__license__ = "GPL"
7__version__ = "1.0"
8__maintainer__ = "William Walters"
9__email__ = "william.a.walters@gmail.com"
10
11""" The purpose of this script is to automate the collapsing of rare taxa in an OTU table
12into a category of 'Other'. Situations may arise where one desires a limited number of
13taxa present in a graph, such as 15 total taxa, to allow for easier recognition of
14colors matching taxa. This script will take an OTU table, find the X most abundant taxa
15(specified by the taxa_count parameter), and collapse all but the X-1 most abundant taxa
16into a group (by default named "Other"). The OTU table metadata 'taxonomy' is checked for,
17and if it exists, the script assumes that that the table is an OTU table and converts the
18metadata 'taxonomy' value to "Other" for the grouped rare taxa. Otherwise, the OTU ID is
19used to group the data into the abundant and "Other" taxa.
20
21This script was written using these dependencies:
22biom version 2.1.5
23QIIME 1.9.1
24"""
25
26from operator import itemgetter
27from os.path import split, splitext, join, basename
28
29from numpy import array, arange, asarray
30from biom import load_table
31from biom.table import Table
32
33from qiime.util import (parse_command_line_parameters, get_options_lookup, create_dir,
34 write_biom_table, make_option)
35
36options_lookup = get_options_lookup()
37
38script_info = {}
39script_info['brief_description'] = """Collapse rare taxa into 'Other' category"""
40script_info['script_description'] = """Collapse rare taxa into 'Other' category"""
41script_info['script_usage'] = []
42script_info['script_usage'].append(("""Standard Example:""",
43 """Use the default settings""",
44 """%prog -i otu_table.biom"""))
45
46script_info['output_description'] = """OTU table with OTUs/taxonomy collapsed, text files
47 with IDs that were collapsed/not collapsed."""
48
49script_info['required_options'] = [
50 make_option('-i', '--input_otu_table', dest='input_otu_table',
51 type='existing_filepath', help='Input OTU table. Can be OTU level, or' +
52 ' be the output of summarize_taxa.py from QIIME.')
53]
54
55script_info['optional_options'] = [
56
57 make_option('-t', '--taxa_count', dest='taxa_count',
58 type='int', default=10,
59 help='Count of taxa to include in output. [default: %default]'),
60
61 make_option('-c', '--collapse_string', default='Other', type='string',
62 help=
63 'Taxa string to use for the collapsed category. [default: %default]'),
64
65 make_option('-o', '--dir_prefix', default='.', type='new_dirpath',
66 help='directory prefix for output files [default: %default]')
67]
68
69script_info['version'] = __version__
70
71
72def main():
73 option_parser, opts, args = parse_command_line_parameters(**script_info)
74
75 taxonomy_metadata_label = 'taxonomy'
76
77 output_dir = opts.dir_prefix
78
79 create_dir(output_dir)
80 otu_table = load_table(opts.input_otu_table)
81
82 # Check for presence of 'taxonomy' metadata, simple test for None data
83 if not otu_table.metadata(axis='observation'):
84 taxonomy_present = False
85 else:
86 for (otu_val, otu_id, otu_metadata) in otu_table.iter(axis='observation'):
87 if taxonomy_metadata_label not in otu_metadata:
88 taxonomy_present = False
89 break
90 else:
91 taxonomy_present = True
92
93 taxa_preserved = []
94 taxa_collapsed = []
95
96 otu_abundance = []
97 sample_metadata = {}
98
99 ids = otu_table.ids()
100
101 # Get the N-1 top abundant taxa, based on taxa_count variable
102 taxa_count = opts.taxa_count - 1
103
104
105 if taxonomy_present:
106 for (otu_val, otu_id, otu_metadata) in otu_table.iter(axis='observation'):
107 otu_abundance.append((array(otu_val).sum(), otu_id,
108 otu_metadata[taxonomy_metadata_label], array(otu_val)))
109 else:
110 for (otu_val, otu_id, metadata) in otu_table.iter(axis='observation'):
111 otu_abundance.append((array(otu_val).sum(), otu_id, array(otu_val)))
112
113 otu_abundance.sort(key=itemgetter(0), reverse=True)
114
115 excluded_ids = []
116
117
118 top_taxa = otu_abundance[0:taxa_count]
119 bottom_taxa_abundance_vals = [0] * len(otu_abundance[taxa_count:][0][-1])
120
121 if taxonomy_present:
122 for (sum, otu_id, taxa, abundance_vals) in otu_abundance[taxa_count:]:
123 abundance_vals = list(abundance_vals)
124 bottom_taxa_abundance_vals =\
125 [bottom_taxa_abundance_vals[i] +\
126 abundance_vals[i] for i in xrange(len(bottom_taxa_abundance_vals))]
127 excluded_ids.append("%s\t%s\n" % (otu_id, "".join(taxa)))
128 else:
129 for (sum, otu_id, abundance_vals) in otu_abundance[taxa_count:]:
130 abundance_vals = list(abundance_vals)
131 bottom_taxa_abundance_vals =\
132 [bottom_taxa_abundance_vals[i] +\
133 abundance_vals[i] for i in xrange(len(bottom_taxa_abundance_vals))]
134 excluded_ids.append("%s\n" % otu_id)
135
136 final_data = []
137 final_otus = []
138 final_taxa = []
139
140 if taxonomy_present:
141 for (sum, otu_id, taxa, abundance_vals) in otu_abundance[0:taxa_count]:
142 final_data.append(list(abundance_vals))
143 final_otus.append(otu_id)
144 final_taxa.append({'taxonomy':taxa})
145 else:
146 for (sum, otu_id, abundance_vals) in otu_abundance[0:taxa_count]:
147 final_data.append(list(abundance_vals))
148 final_otus.append(otu_id)
149
150 # Add in collapsed "other" category
151
152 final_data.append(bottom_taxa_abundance_vals)
153 final_otus.append(unicode(opts.collapse_string))
154 final_taxa.append({'taxonomy':[unicode(opts.collapse_string)]})
155
156 # Create table for writing
157
158 if taxonomy_present:
159 table = Table(asarray(final_data), final_otus, ids, final_taxa)
160 else:
161 table = Table(asarray(final_data), final_otus, ids)
162
163
164 output_name = join(output_dir,
165 basename(opts.input_otu_table).replace(".biom", "_rare_collapsed.biom"))
166
167 write_biom_table(table, output_name)
168
169 output_excluded_name = join(output_dir, "rare_merged_ids.txt")
170 f = open(output_excluded_name, "w")
171
172 for id in excluded_ids:
173 f.write(id)
174
175
176
177
178if __name__ == "__main__":
179 main()