Mercurial > repos > george-weingart > lefse

comparison plot_features.py @ 0:e7cd19afda2e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Lefse
author george-weingart
date Tue, 13 May 2014 21:57:00 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e7cd19afda2e
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()