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