|
0
|
1 #!/usr/bin/env python
|
|
|
2 # Copyright (C) <2015> EMBL-European Bioinformatics Institute
|
|
|
3
|
|
|
4 # This program is free software: you can redistribute it and/or
|
|
|
5 # modify it under the terms of the GNU General Public License as
|
|
|
6 # published by the Free Software Foundation, either version 3 of
|
|
|
7 # the License, or (at your option) any later version.
|
|
|
8
|
|
|
9 # This program is distributed in the hope that it will be useful,
|
|
|
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
12 # GNU General Public License for more details.
|
|
|
13
|
|
|
14 # Neither the institution name nor the name roary_plots
|
|
|
15 # can be used to endorse or promote products derived from
|
|
|
16 # this software without prior written permission.
|
|
|
17 # For written permission, please contact <marco@ebi.ac.uk>.
|
|
|
18
|
|
|
19 # Products derived from this software may not be called roary_plots
|
|
|
20 # nor may roary_plots appear in their names without prior written
|
|
|
21 # permission of the developers. You should have received a copy
|
|
|
22 # of the GNU General Public License along with this program.
|
|
|
23 # If not, see <http://www.gnu.org/licenses/>.
|
|
|
24
|
|
|
25 __author__ = "Marco Galardini"
|
|
|
26 __version__ = '0.1.0'
|
|
|
27
|
|
|
28 def get_options():
|
|
|
29 import argparse
|
|
|
30
|
|
|
31 # create the top-level parser
|
|
|
32 description = "Create plots from roary outputs"
|
|
|
33 parser = argparse.ArgumentParser(description = description,
|
|
|
34 prog = 'roary_plots.py')
|
|
|
35
|
|
|
36 parser.add_argument('tree', action='store',
|
|
|
37 help='Newick Tree file', default='accessory_binary_genes.fa.newick')
|
|
|
38 parser.add_argument('spreadsheet', action='store',
|
|
|
39 help='Roary gene presence/absence spreadsheet', default='gene_presence_absence.csv')
|
|
|
40
|
|
|
41 parser.add_argument('--labels', action='store_true',
|
|
|
42 default=False,
|
|
|
43 help='Add node labels to the tree (up to 10 chars)')
|
|
|
44 parser.add_argument('--format',
|
|
|
45 choices=('png',
|
|
|
46 'tiff',
|
|
|
47 'pdf',
|
|
|
48 'svg'),
|
|
|
49 default='png',
|
|
|
50 help='Output format [Default: png]')
|
|
|
51 parser.add_argument('-N', '--skipped-columns', action='store',
|
|
|
52 type=int,
|
|
|
53 default=14,
|
|
|
54 help='First N columns of Roary\'s output to exclude [Default: 14]')
|
|
|
55
|
|
|
56 parser.add_argument('--version', action='version',
|
|
|
57 version='%(prog)s '+__version__)
|
|
|
58
|
|
|
59 return parser.parse_args()
|
|
|
60
|
|
|
61 if __name__ == "__main__":
|
|
|
62 options = get_options()
|
|
|
63
|
|
|
64 import matplotlib
|
|
|
65 matplotlib.use('Agg')
|
|
|
66
|
|
|
67 import matplotlib.pyplot as plt
|
|
|
68 import seaborn as sns
|
|
|
69
|
|
|
70 sns.set_style('white')
|
|
|
71
|
|
|
72 import os
|
|
|
73 import pandas as pd
|
|
|
74 import numpy as np
|
|
|
75 from Bio import Phylo
|
|
|
76
|
|
|
77 t = Phylo.read(options.tree, 'newick')
|
|
|
78
|
|
|
79 # Max distance to create better plots
|
|
|
80 mdist = max([t.distance(t.root, x) for x in t.get_terminals()])
|
|
|
81
|
|
|
82 # Load roary
|
|
|
83 roary = pd.read_csv(options.spreadsheet, low_memory=False)
|
|
|
84 # Set index (group name)
|
|
|
85 roary.set_index('Gene', inplace=True)
|
|
|
86 # Drop the other info columns
|
|
|
87 roary.drop(list(roary.columns[:options.skipped_columns-1]), axis=1, inplace=True)
|
|
|
88
|
|
|
89 # Transform it in a presence/absence matrix (1/0)
|
|
|
90 roary.replace('.{2,100}', 1, regex=True, inplace=True)
|
|
|
91 roary.replace(np.nan, 0, regex=True, inplace=True)
|
|
|
92
|
|
|
93 # Sort the matrix by the sum of strains presence
|
|
|
94 idx = roary.sum(axis=1).sort_values(ascending=False).index
|
|
|
95 roary_sorted = roary.loc[idx]
|
|
|
96
|
|
|
97 # Pangenome frequency plot
|
|
|
98 plt.figure(figsize=(7, 5))
|
|
|
99
|
|
|
100 plt.hist(roary.sum(axis=1), roary.shape[1],
|
|
|
101 histtype="stepfilled", alpha=.7)
|
|
|
102
|
|
|
103 plt.xlabel('No. of genomes')
|
|
|
104 plt.ylabel('No. of genes')
|
|
|
105
|
|
|
106 sns.despine(left=True,
|
|
|
107 bottom=True)
|
|
|
108 plt.savefig('pangenome_frequency.%s'%options.format, dpi=300)
|
|
|
109 plt.clf()
|
|
|
110
|
|
|
111 # Sort the matrix according to tip labels in the tree
|
|
|
112 roary_sorted = roary_sorted[[x.name for x in t.get_terminals()]]
|
|
|
113
|
|
|
114 # Plot presence/absence matrix against the tree
|
|
|
115 with sns.axes_style('whitegrid'):
|
|
|
116 fig = plt.figure(figsize=(17, 10))
|
|
|
117
|
|
|
118 ax1=plt.subplot2grid((1,40), (0, 10), colspan=30)
|
|
|
119 a=ax1.matshow(roary_sorted.T, cmap=plt.cm.Blues,
|
|
|
120 vmin=0, vmax=1,
|
|
|
121 aspect='auto',
|
|
|
122 interpolation='none',
|
|
|
123 )
|
|
|
124 ax1.set_yticks([])
|
|
|
125 ax1.set_xticks([])
|
|
|
126 ax1.axis('off')
|
|
|
127
|
|
|
128 ax = fig.add_subplot(1,2,1)
|
|
|
129 # matplotlib v1/2 workaround
|
|
|
130 try:
|
|
|
131 ax=plt.subplot2grid((1,40), (0, 0), colspan=10, facecolor='white')
|
|
|
132 except AttributeError:
|
|
|
133 ax=plt.subplot2grid((1,40), (0, 0), colspan=10, axisbg='white')
|
|
|
134
|
|
|
135 fig.subplots_adjust(wspace=0, hspace=0)
|
|
|
136
|
|
|
137 ax1.set_title('Roary matrix\n(%d gene clusters)'%roary.shape[0])
|
|
|
138
|
|
|
139 if options.labels:
|
|
|
140 fsize = 12 - 0.1*roary.shape[1]
|
|
|
141 if fsize < 7:
|
|
|
142 fsize = 7
|
|
|
143 with plt.rc_context({'font.size': fsize}):
|
|
|
144 Phylo.draw(t, axes=ax,
|
|
|
145 show_confidence=False,
|
|
|
146 label_func=lambda x: str(x)[:10],
|
|
|
147 xticks=([],), yticks=([],),
|
|
|
148 ylabel=('',), xlabel=('',),
|
|
|
149 xlim=(-mdist*0.1,mdist+mdist*0.45-mdist*roary.shape[1]*0.001),
|
|
|
150 axis=('off',),
|
|
|
151 title=('Tree\n(%d strains)'%roary.shape[1],),
|
|
|
152 do_show=False,
|
|
|
153 )
|
|
|
154 else:
|
|
|
155 Phylo.draw(t, axes=ax,
|
|
|
156 show_confidence=False,
|
|
|
157 label_func=lambda x: None,
|
|
|
158 xticks=([],), yticks=([],),
|
|
|
159 ylabel=('',), xlabel=('',),
|
|
|
160 xlim=(-mdist*0.1,mdist+mdist*0.1),
|
|
|
161 axis=('off',),
|
|
|
162 title=('Tree\n(%d strains)'%roary.shape[1],),
|
|
|
163 do_show=False,
|
|
|
164 )
|
|
|
165 plt.savefig('pangenome_matrix.%s'%options.format, dpi=300)
|
|
|
166 plt.clf()
|
|
|
167
|
|
|
168 # Plot the pangenome pie chart
|
|
|
169 plt.figure(figsize=(10, 10))
|
|
|
170
|
|
|
171 core = roary[(roary.sum(axis=1) >= roary.shape[1]*0.99) & (roary.sum(axis=1) <= roary.shape[1] )].shape[0]
|
|
|
172 softcore = roary[(roary.sum(axis=1) >= roary.shape[1]*0.95) & (roary.sum(axis=1) < roary.shape[1]*0.99)].shape[0]
|
|
|
173 shell = roary[(roary.sum(axis=1) >= roary.shape[1]*0.15) & (roary.sum(axis=1) < roary.shape[1]*0.95)].shape[0]
|
|
|
174 cloud = roary[roary.sum(axis=1) < roary.shape[1]*0.15].shape[0]
|
|
|
175
|
|
|
176 total = roary.shape[0]
|
|
|
177
|
|
|
178 def my_autopct(pct):
|
|
|
179 val=int(round(pct*total/100.0))
|
|
|
180 return '{v:d}'.format(v=val)
|
|
|
181
|
|
|
182 a=plt.pie([core, softcore, shell, cloud],
|
|
|
183 labels=['core\n(%d <= strains <= %d)'%(roary.shape[1]*.99,roary.shape[1]),
|
|
|
184 'soft-core\n(%d <= strains < %d)'%(roary.shape[1]*.95,roary.shape[1]*.99),
|
|
|
185 'shell\n(%d <= strains < %d)'%(roary.shape[1]*.15,roary.shape[1]*.95),
|
|
|
186 'cloud\n(strains < %d)'%(roary.shape[1]*.15)],
|
|
|
187 explode=[0.1, 0.05, 0.02, 0], radius=0.9,
|
|
|
188 colors=[(0, 0, 1, float(x)/total) for x in (core, softcore, shell, cloud)],
|
|
|
189 autopct=my_autopct)
|
|
|
190 plt.savefig('pangenome_pie.%s'%options.format, dpi=300)
|
|
|
191 plt.clf()
|
