-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcpo_galaxy_tree.py
353 lines (317 loc) · 15.8 KB
/
cpo_galaxy_tree.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
#!/home/jjjjia/.conda/envs/py36/bin/python
#$ -S /home/jjjjia/.conda/envs/py36/bin/python
#$ -V # Pass environment variables to the job
#$ -N CPO_pipeline # Replace with a more specific job name
#$ -wd /home/jjjjia/testCases # Use the current working dir
#$ -pe smp 1 # Parallel Environment (how many cores)
#$ -l h_vmem=11G # Memory (RAM) allocation *per core*
#$ -e ./logs/$JOB_ID.err
#$ -o ./logs/$JOB_ID.log
#$ -m ea
#$ -M [email protected]
# >python cpo_galaxy_tree.py -t /path/to/tree.ph -d /path/to/distance/matrix -m /path/to/metadata
# <requirements>
# <requirement type="package" version="0.23.4">pandas</requirement>
# <requirement type="package" version="3.6">python</requirement>
# <requirement type="package" version="3.1.1">ete3</requirement>
# <requirement type="package" version="5.6.0">pyqt</requirement>
# <requirement type="package" version="5.6.2">qt</requirement>
# </requirements>
import subprocess
import pandas #conda pandas
import optparse
import os
os.environ['QT_QPA_PLATFORM']='offscreen'
import datetime
import sys
import time
import urllib.request
import gzip
import collections
import json
import numpy #conda numpy
import ete3 as e #conda ete3 3.1.1**** >requires pyqt5
#parses some parameters
parser = optparse.OptionParser("Usage: %prog [options] arg1 arg2 ...")
parser.add_option("-t", "--tree", dest="treePath", type="string", default="./pipelineTest/tree.txt", help="absolute file path to phylip tree")
parser.add_option("-d", "--distance", dest="distancePath", type="string", default="./pipelineTest/dist.tabular", help="absolute file path to distance matrix")
parser.add_option("-m", "--metadata", dest="metadataPath", type="string", default="./pipelineTest/metadata.tabular",help="absolute file path to metadata file")
parser.add_option("-o", "--output_file", dest="outputFile", type="string", default="tree.png", help="Output graphics file. Use ending 'png', 'pdf' or 'svg' to specify file format.")
# sensitive data adder
parser.add_option("-p", "--sensitive_data", dest="sensitivePath", type="string", default="", help="Spreadsheet (CSV) with sensitive metadata")
parser.add_option("-c", "--sensitive_cols", dest="sensitiveCols", type="string", default="", help="CSV list of column names from sensitive metadata spreadsheet to use as labels on dendrogram")
parser.add_option("-b", "--bcid_column", dest="bcidCol", type="string", default="BCID", help="Column name of BCID in sensitive metadata file")
parser.add_option("-n", "--missing_value", dest="naValue", type="string", default="NA", help="Value to write for missing data.")
(options,args) = parser.parse_args()
treePath = str(options.treePath).lstrip().rstrip()
distancePath = str(options.distancePath).lstrip().rstrip()
metadataPath = str(options.metadataPath).lstrip().rstrip()
sensitivePath = str(options.sensitivePath).lstrip().rstrip()
sensitiveCols = str(options.sensitiveCols).lstrip().rstrip()
outputFile = str(options.outputFile).lstrip().rstrip()
bcidCol = str( str(options.bcidCol).lstrip().rstrip() )
naValue = str( str(options.naValue).lstrip().rstrip() )
#region result objects
#define some objects to store values from results
#//TODO this is not the proper way of get/set private object variables. every value has manually assigned defaults intead of specified in init(). Also, use property(def getVar, def setVar).
class SensitiveMetadata(object):
def __init__(self):
x = pandas.read_csv( sensitivePath )
col_names = [ s for s in sensitiveCols.split(',')] # convert to 0 offset
if not bcidCol in col_names:
col_names.append( bcidCol )
all_cols = [ str(col) for col in x.columns ]
col_idxs = [ all_cols.index(col) for col in col_names ]
self.sensitive_data = x.iloc[:, col_idxs]
def get_columns(self):
cols = [ str(x) for x in self.sensitive_data.columns ]
return cols
def get_value( self, bcid, column_name ): # might be nice to get them all in single call via an input list of bcids ... for later
bcids= list( self.sensitive_data.loc[:, bcidCol ] ) # get the list of all BCIDs in sensitive metadata
if not bcid in bcids:
return naValue
else:
row_idx = bcids.index( bcid ) # lookup the row for this BCID
return self.sensitive_data.loc[ row_idx, column_name ] # return the one value based on the column (col_idx) and this row
class workflowResult(object):
def __init__(self):
self.new = False
self.ID = "?"
self.ExpectedSpecies = "?"
self.MLSTSpecies = "?"
self.SequenceType = "?"
self.MLSTScheme = "?"
self.CarbapenemResistanceGenes ="?"
self.plasmidBestMatch ="?"
self.plasmididentity =-1
self.plasmidsharedhashes ="?"
self.OtherAMRGenes="?"
self.TotalPlasmids = -1
self.plasmids = []
self.DefinitelyPlasmidContigs ="?"
self.LikelyPlasmidContigs="?"
self.row = ""
class plasmidObj(object):
def __init__(self):
self.PlasmidsID = 0
self.Num_Contigs = 0
self.PlasmidLength = 0
self.PlasmidRepType = ""
self.PlasmidMobility = ""
self.NearestReference = ""
#endregion
#region useful functions
def read(path): #read in a text file to a list
return [line.rstrip('\n') for line in open(path)]
def execute(command): #subprocess.popen call bash command
process = subprocess.Popen(command, shell=False, cwd=curDir, universal_newlines=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
# Poll process for new output until finished
while True:
nextline = process.stdout.readline()
if nextline == '' and process.poll() is not None:
break
sys.stdout.write(nextline)
sys.stdout.flush()
output = process.communicate()[0]
exitCode = process.returncode
if (exitCode == 0):
return output
else:
raise subprocess.CalledProcessError(exitCode, command)
def httpGetFile(url, filepath=""): #download a file from the web
if (filepath == ""):
return urllib.request.urlretrieve(url)
else:
urllib.request.urlretrieve(url, filepath)
return True
def gunzip(inputpath="", outputpath=""): #gunzip
if (outputpath == ""):
with gzip.open(inputpath, 'rb') as f:
gzContent = f.read()
return gzContent
else:
with gzip.open(inputpath, 'rb') as f:
gzContent = f.read()
with open(outputpath, 'wb') as out:
out.write(gzContent)
return True
def addFace(name): #function to add a facet to a tree
#if its the reference branch, populate the faces with column headers
face = e.faces.TextFace(name,fsize=10,tight_text=True)
face.border.margin = 5
face.margin_right = 5
face.margin_left = 5
return face
#endregion
#region functions to parse result files
def ParseWorkflowResults(pathToResult):
_worflowResult = {}
r = pandas.read_csv(pathToResult, delimiter='\t', header=0)
r = r.replace(numpy.nan, '', regex=True)
_naResult = workflowResult()
_worflowResult["na"] = _naResult
for i in range(len(r.index)):
_results = workflowResult()
if(str(r.loc[r.index[i], 'new']).lower() == "new"):
_results.new = True
else:
_results.new = False
_results.ID = str(r.loc[r.index[i], 'ID']).replace(".fa","")
_results.ExpectedSpecies = str(r.loc[r.index[i], 'Expected Species'])
_results.MLSTSpecies = str(r.loc[r.index[i], 'MLST Species'])
_results.SequenceType = str(r.loc[r.index[i], 'Sequence Type'])
_results.MLSTScheme = (str(r.loc[r.index[i], 'MLST Scheme']))
_results.CarbapenemResistanceGenes = (str(r.loc[r.index[i], 'Carbapenem Resistance Genes']))
_results.OtherAMRGenes = (str(r.loc[r.index[i], 'Other AMR Genes']))
_results.TotalPlasmids = int(r.loc[r.index[i], 'Total Plasmids'])
_results.plasmidBestMatch = str(r.loc[r.index[i], 'Plasmid Best Match'])
_results.plasmididentity = str(r.loc[r.index[i], 'Plasmid Identity'])
_results.plasmidsharedhashes = str(r.loc[r.index[i], 'Plasmid Shared Hash'])
for j in range(0,_results.TotalPlasmids):
_plasmid = plasmidObj()
_plasmid.PlasmidsID =(((str(r.loc[r.index[i], 'Plasmids ID'])).split(";"))[j])
_plasmid.Num_Contigs = (((str(r.loc[r.index[i], 'Num_Contigs'])).split(";"))[j])
_plasmid.PlasmidLength = (((str(r.loc[r.index[i], 'Plasmid Length'])).split(";"))[j])
_plasmid.PlasmidRepType = (((str(r.loc[r.index[i], 'Plasmid RepType'])).split(";"))[j])
_plasmid.PlasmidMobility = ((str(r.loc[r.index[i], 'Plasmid Mobility'])).split(";"))[j]
_plasmid.NearestReference = ((str(r.loc[r.index[i], 'Nearest Reference'])).split(";"))[j]
_results.plasmids.append(_plasmid)
_results.DefinitelyPlasmidContigs = (str(r.loc[r.index[i], 'Definitely Plasmid Contigs']))
_results.LikelyPlasmidContigs = (str(r.loc[r.index[i], 'Likely Plasmid Contigs']))
_results.row = "\t".join(str(x) for x in r.ix[i].tolist())
_worflowResult[_results.ID] = _results
return _worflowResult
#endregion
def Main():
if len(sensitivePath)>0:
sensitive_meta_data = SensitiveMetadata()
metadata = ParseWorkflowResults(metadataPath)
distance = read(distancePath)
treeFile = "".join(read(treePath))
distanceDict = {} #store the distance matrix as rowname:list<string>
for i in range(len(distance)):
temp = distance[i].split("\t")
distanceDict[temp[0]] = temp[1:]
#region create box tree
#region step5: tree construction
treeFile = "".join(read(treePath))
t = e.Tree(treeFile)
t.set_outgroup(t&"Reference")
#set the tree style
ts = e.TreeStyle()
ts.show_leaf_name = True
ts.show_branch_length = True
ts.scale = 2000 #pixel per branch length unit
ts.branch_vertical_margin = 15 #pixel between branches
style2 = e.NodeStyle()
style2["fgcolor"] = "#000000"
style2["shape"] = "circle"
style2["vt_line_color"] = "#0000aa"
style2["hz_line_color"] = "#0000aa"
style2["vt_line_width"] = 2
style2["hz_line_width"] = 2
style2["vt_line_type"] = 0 # 0 solid, 1 dashed, 2 dotted
style2["hz_line_type"] = 0
for n in t.traverse():
n.set_style(style2)
#find the plasmid origins
plasmidIncs = {}
for key in metadata:
for plasmid in metadata[key].plasmids:
for inc in plasmid.PlasmidRepType.split(","):
if (inc.lower().find("inc") > -1):
if not (inc in plasmidIncs):
plasmidIncs[inc] = [metadata[key].ID]
else:
if metadata[key].ID not in plasmidIncs[inc]:
plasmidIncs[inc].append(metadata[key].ID)
#plasmidIncs = sorted(plasmidIncs)
for n in t.traverse(): #loop through the nodes of a tree
if (n.is_leaf() and n.name == "Reference"):
#if its the reference branch, populate the faces with column headers
index = 0
if len(sensitivePath)>0: #sensitive metadat @ chris
for sensitive_data_column in sensitive_meta_data.get_columns():
(t&"Reference").add_face(addFace(sensitive_data_column), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("SampleID"), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("New?"), index, "aligned")
index = index + 1
for i in range(len(plasmidIncs)): #this loop adds the columns (aka the incs) to the reference node
(t&"Reference").add_face(addFace(list(plasmidIncs.keys())[i]), i + index, "aligned")
index = index + len(plasmidIncs)
(t&"Reference").add_face(addFace("MLSTScheme"), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("Sequence Type"), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("Carbapenamases"), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("Plasmid Best Match"), index, "aligned")
index = index + 1
(t&"Reference").add_face(addFace("Best Match Identity"), index, "aligned")
index = index + 1
for i in range(len(distanceDict[list(distanceDict.keys())[0]])): #this loop adds the distance matrix
(t&"Reference").add_face(addFace(distanceDict[list(distanceDict.keys())[0]][i]), index + i, "aligned")
index = index + len(distanceDict[list(distanceDict.keys())[0]])
elif (n.is_leaf() and not n.name == "Reference"):
#not reference branches, populate with metadata
index = 0
if len(sensitivePath)>0: #sensitive metadata @ chris
# pushing in sensitive data
for sensitive_data_column in sensitive_meta_data.get_columns():
# tree uses bcids like BC18A021A_S12
# while sens meta-data uses BC18A021A
# trim the "_S.*" if present
bcid = str(mData.ID)
if bcid.find( "_S" ) != -1:
bcid = bcid[ 0:bcid.find( "_S" ) ]
sens_col_val = sensitive_meta_data.get_value(bcid=bcid, column_name=sensitive_data_column )
n.add_face(addFace(sens_col_val), index, "aligned")
index = index + 1
if (n.name.replace(".fa","") in metadata.keys()):
mData = metadata[n.name.replace(".fa","")]
else:
mData = metadata["na"]
n.add_face(addFace(mData.ID), index, "aligned")
index = index + 1
if (mData.new == True): #new column
face = e.RectFace(30,30,"green","green") # TextFace("Y",fsize=10,tight_text=True)
face.border.margin = 5
face.margin_right = 5
face.margin_left = 5
face.vt_align = 1
face.ht_align = 1
n.add_face(face, index, "aligned")
index = index + 1
for incs in plasmidIncs: #this loop adds presence/absence to the sample nodes
if (n.name.replace(".fa","") in plasmidIncs[incs]):
face = e.RectFace(30,30,"black","black") # TextFace("Y",fsize=10,tight_text=True)
face.border.margin = 5
face.margin_right = 5
face.margin_left = 5
face.vt_align = 1
face.ht_align = 1
n.add_face(face, list(plasmidIncs.keys()).index(incs) + index, "aligned")
index = index + len(plasmidIncs)
n.add_face(addFace(mData.MLSTSpecies), index, "aligned")
index = index + 1
n.add_face(addFace(mData.SequenceType), index, "aligned")
index = index + 1
n.add_face(addFace(mData.CarbapenemResistanceGenes), index, "aligned")
index = index + 1
n.add_face(addFace(mData.plasmidBestMatch), index, "aligned")
index = index + 1
n.add_face(addFace(mData.plasmididentity), index, "aligned")
index = index + 1
for i in range(len(distanceDict[list(distanceDict.keys())[0]])): #this loop adds distance matrix
if (n.name in distanceDict): #make sure the column is in the distance matrice
n.add_face(addFace(list(distanceDict[n.name])[i]), index + i, "aligned")
t.render(outputFile, w=5000,units="mm", tree_style=ts) #save it as a png, pdf, svg or an phyloxml
#endregion
#endregion
start = time.time()#time the analysis
#analysis time
Main()
end = time.time()
print("Finished!\nThe analysis used: " + str(end-start) + " seconds")