|
1
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import os,sys,matplotlib,zipfile,argparse,string
|
|
|
4 matplotlib.use('Agg')
|
|
|
5 from pylab import *
|
|
|
6 from lefse import *
|
|
|
7 import random as rand
|
|
|
8
|
|
|
9 colors = ['r','g','b','m','c']
|
|
|
10
|
|
|
11 def read_params(args):
|
|
|
12 parser = argparse.ArgumentParser(description='Cladoplot')
|
|
|
13 parser.add_argument('input_file_1', metavar='INPUT_FILE', type=str, help="dataset files")
|
|
|
14 parser.add_argument('input_file_2', metavar='INPUT_FILE', type=str, help="LEfSe output file")
|
|
|
15 parser.add_argument('output_file', metavar='OUTPUT_FILE', type=str, help="the file for the output (the zip file if an archive is required, the output directory otherwise)")
|
|
|
16 parser.add_argument('--width',dest="width", type=float, default=10.0 )
|
|
|
17 parser.add_argument('--height',dest="height", type=float, default=4.0)
|
|
|
18 parser.add_argument('--top',dest="top", type=float, default=-1.0, help="set maximum y limit (-1.0 means automatic limit)")
|
|
|
19 parser.add_argument('--bot',dest="bot", type=float, default=0.0, help="set minimum y limit (default 0.0, -1.0 means automatic limit)")
|
|
|
20 parser.add_argument('--title_font_size',dest="title_font_size", type=str, default="14")
|
|
|
21 parser.add_argument('--class_font_size',dest="class_font_size", type=str, default="14")
|
|
|
22 parser.add_argument('--class_label_pos',dest="class_label_pos", type=str, choices=["up","down"], default="up")
|
|
|
23 parser.add_argument('--subcl_mean',dest="subcl_mean", type=str, choices=["y","n"], default="y")
|
|
|
24 parser.add_argument('--subcl_median',dest="subcl_median", type=str, choices=["y","n"], default="y")
|
|
|
25 parser.add_argument('--font_size',dest="font_size", type=str, default="10")
|
|
|
26 parser.add_argument('-n',dest="unused", metavar="flt", type=float, default=-1.0,help="unused")
|
|
|
27 parser.add_argument('--format', dest="format", default="png", choices=["png","pdf","svg"], type=str, help="the format for the output file")
|
|
|
28 parser.add_argument('-f', dest="f", default="diff", choices=["all","diff","one"], type=str, help="wheter to plot all features (all), only those differentially abundant according to LEfSe or only one (the one given with --feature_name) ")
|
|
|
29 parser.add_argument('--feature_name', dest="feature_name", default="", type=str, help="The name of the feature to plot (levels separated by .) ")
|
|
|
30 parser.add_argument('--feature_num', dest="feature_num", default="-1", type=int, help="The number of the feature to plot ")
|
|
|
31 parser.add_argument('--archive', dest="archive", default="none", choices=["zip","none"], type=str, help="")
|
|
|
32 parser.add_argument('--background_color',dest="back_color", type=str, choices=["k","w"], default="w", help="set the color of the background")
|
|
|
33 parser.add_argument('--dpi',dest="dpi", type=int, default=72)
|
|
|
34
|
|
|
35 args = parser.parse_args()
|
|
|
36
|
|
|
37 return vars(args)
|
|
|
38
|
|
|
39 def read_data(file_data,file_feats,params):
|
|
|
40 with open(file_feats, 'r') as features:
|
|
|
41 feats_to_plot = [(f.split()[:-1],len(f.split()) == 5) for f in features.readlines()]
|
|
|
42 if not feats_to_plot:
|
|
|
43 print "No features to plot\n"
|
|
|
44 sys.exit(0)
|
|
|
45 feats,cls,class_sl,subclass_sl,class_hierarchy,params['norm_v'] = load_data(file_data, True)
|
|
|
46 if params['feature_num'] > 0:
|
|
|
47 params['feature_name'] = [line.split()[0] for line in open(params['input_file_2'])][params['feature_num']-1]
|
|
|
48 features = {}
|
|
|
49 for f in feats_to_plot:
|
|
|
50 if params['f'] == "diff" and not f[1]: continue
|
|
|
51 if params['f'] == "one" and f[0][0] != params['feature_name']: continue
|
|
|
52 features[f[0][0]] = {'dim':float(f[0][1]), 'abundances':feats[f[0][0]], 'sig':f[1], 'cls':cls, 'class_sl':class_sl, 'subclass_sl':subclass_sl, 'class_hierarchy':class_hierarchy}
|
|
|
53 if not features:
|
|
|
54 print "No features to plot\n"
|
|
|
55 sys.exit(0)
|
|
|
56 return features
|
|
|
57
|
|
|
58 def plot(name,k_n,feat,params):
|
|
|
59 fig = plt.figure(figsize=(params['width'], params['height']),edgecolor=params['fore_color'],facecolor=params['back_color'])
|
|
|
60 ax = fig.add_subplot(111,axis_bgcolor=params['back_color'])
|
|
|
61 subplots_adjust(bottom=0.15)
|
|
|
62
|
|
|
63 max_m = 0.0
|
|
|
64 norm = 1.0 if float(params['norm_v']) < 0.0 else float(params['norm_v'])
|
|
|
65 for v in feat['subclass_sl'].values():
|
|
|
66 fr,to = v[0], v[1]
|
|
|
67 median = numpy.mean(feat['abundances'][fr:to])
|
|
|
68 if median > max_m: max_m = median
|
|
|
69 max_m /= norm
|
|
|
70 max_v = max_m*3 if max_m*3 < max(feat['abundances'])*1.1/norm else max(feat['abundances'])/norm
|
|
|
71 min_v = max(0.0,min(feat['abundances'])*0.9/norm)
|
|
|
72
|
|
|
73 if params['top'] > 0.0: max_v = params['top']
|
|
|
74 if params['bot'] >= 0.0: min_v = params['bot']
|
|
|
75
|
|
|
76 if max_v == 0.0: max_v = 0.0001
|
|
|
77 if max_v == min_v: max_v = min_v*1.1
|
|
|
78
|
|
|
79 cl_sep = max(int(sum([vv[1]/norm - vv[0]/norm for vv in feat['class_sl'].values()])/150.0),1)
|
|
|
80 seps = []
|
|
|
81 xtics = []
|
|
|
82 x2tics = []
|
|
|
83 last_fr = 0.0
|
|
|
84 for i,cl in enumerate(sorted(feat['class_hierarchy'].keys())):
|
|
|
85 for j,subcl in enumerate(feat['class_hierarchy'][cl]):
|
|
|
86 fr = feat['subclass_sl'][subcl][0]
|
|
|
87 to = feat['subclass_sl'][subcl][1]
|
|
|
88 val = feat['abundances'][fr:to]
|
|
|
89 fr += cl_sep*i
|
|
|
90 to += cl_sep*i
|
|
|
91 pos = arange(fr,to)
|
|
|
92 max_x = to
|
|
|
93 col = colors[j%len(colors)]
|
|
|
94 vv = [v1/norm for v1 in val]
|
|
|
95 median = numpy.median(vv)
|
|
|
96 mean = numpy.mean(vv)
|
|
|
97 valv = [max(min(v/norm,max_v),min_v) for v in val]
|
|
|
98 ax.bar(pos,valv, align='center', color=col, edgecolor=col, linewidth=0.1 )
|
|
|
99 if params['subcl_median'] == 'y': ax.plot([fr,to-1],[median,median],"k--",linewidth=1,color=params['fore_color'])
|
|
|
100 if params['subcl_mean'] == 'y': ax.plot([fr,to-1],[mean,mean],"-",linewidth=1,color=params['fore_color'])
|
|
|
101 nna = subcl if subcl.count("_") == 0 or not subcl.startswith(cl) else "_".join(subcl.split(cl)[1:])
|
|
|
102 if nna == "subcl" or nna == "_subcl": nna = " "
|
|
|
103 xtics.append(((fr+(to-fr)/2),nna))
|
|
|
104 seps.append(float(to))
|
|
|
105 x2tics.append(((last_fr+(to-last_fr)/2),cl))
|
|
|
106 last_fr = to + float(cl_sep)
|
|
|
107 for s in seps[:-1]:
|
|
|
108 ax.plot([s,s],[min_v,max_v],"-",linewidth=5,color=params['fore_color'])
|
|
|
109 ax.set_title(k_n, size=params['title_font_size'])
|
|
|
110 xticks([x[0] for x in xtics],[x[1] for x in xtics],rotation=-30, ha = 'left', fontsize=params['font_size'], color=params['fore_color'])
|
|
|
111 yticks(fontsize=params['font_size'])
|
|
|
112
|
|
|
113 ylabel('Relative abundance')
|
|
|
114 ax.set_ylim((min_v,max_v))
|
|
|
115 a,b = ax.get_xlim()
|
|
|
116 ax.set_xlim((0-float(last_fr)/float(b-a),max_x))
|
|
|
117 ax.yaxis.grid(True)
|
|
|
118
|
|
|
119 def get_col_attr(x):
|
|
|
120 return hasattr(x, 'set_color') and not hasattr(x, 'set_facecolor')
|
|
|
121 def get_edgecol_attr(x):
|
|
|
122 return hasattr(x, 'set_edgecolor')
|
|
|
123
|
|
|
124
|
|
|
125 for o in fig.findobj(get_col_attr):
|
|
|
126 o.set_color(params['fore_color'])
|
|
|
127 for o in fig.findobj(get_edgecol_attr):
|
|
|
128 if o.get_edgecolor() == params['back_color']:
|
|
|
129 o.set_edgecolor(params['fore_color'])
|
|
|
130
|
|
|
131 for t in x2tics:
|
|
|
132 m = ax.get_ylim()[1]*0.97 if params['class_label_pos']=='up' else 0.07
|
|
|
133 plt.text(t[0],m, "class: "+t[1], ha ="center", size=params['class_font_size'], va="top", bbox = dict(boxstyle="round", ec='k', fc='y'))
|
|
|
134
|
|
|
135
|
|
|
136 plt.savefig(name,format=params['format'],facecolor=params['back_color'],edgecolor=params['fore_color'],dpi=params['dpi'])
|
|
|
137 plt.close()
|
|
|
138 return name
|
|
|
139
|
|
|
140 if __name__ == '__main__':
|
|
|
141 params = read_params(sys.argv)
|
|
|
142 params['fore_color'] = 'w' if params['back_color'] == 'k' else 'k'
|
|
|
143 features = read_data(params['input_file_1'],params['input_file_2'],params)
|
|
|
144 if params['archive'] == "zip": file = zipfile.ZipFile(params['output_file'], "w")
|
|
|
145 for k,f in features.items():
|
|
|
146 print "Exporting ", k
|
|
|
147 if params['archive'] == "zip":
|
|
|
148 of = plot("/tmp/"+str(int(f['sig']))+"_"+"-".join(k.split("."))+"."+params['format'],k,f,params)
|
|
|
149 file.write(of, os.path.basename(of), zipfile.ZIP_DEFLATED)
|
|
|
150 else:
|
|
|
151 if params['f'] == 'one': plot(params['output_file'],k,f,params)
|
|
|
152 else: plot(params['output_file']+str(int(f['sig']))+"_"+"-".join(k.split("."))+"."+params['format'],k,f,params)
|
|
|
153 if params['archive'] == "zip": file.close()
|
