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>
7
+ − 21 # </requirements>
1
+ − 22
+ − 23 import subprocess
7
+ − 24 import pandas #conda pandas
1
+ − 25 import optparse
+ − 26 import os
8
+ − 27 os.environ['QT_QPA_PLATFORM']='offscreen'
1
+ − 28 import datetime
+ − 29 import sys
+ − 30 import time
+ − 31 import urllib.request
+ − 32 import gzip
+ − 33 import collections
+ − 34 import json
7
+ − 35 import numpy #conda numpy
+ − 36 import ete3 as e #conda ete3 3.1.1**** >requires pyqt5
6
+ − 37
1
+ − 38
+ − 39 #parses some parameters
+ − 40 parser = optparse.OptionParser("Usage: %prog [options] arg1 arg2 ...")
+ − 41 parser.add_option("-t", "--tree", dest="treePath", type="string", default="./pipelineTest/tree.txt", help="identifier of the isolate")
+ − 42 parser.add_option("-d", "--distance", dest="distancePath", type="string", default="./pipelineTest/distance.tab", help="absolute file path forward read (R1)")
+ − 43 parser.add_option("-m", "--metadata", dest="metadataPath", type="string", default="./pipelineTest/metadata.tsv",help="absolute file path to reverse read (R2)")
+ − 44 (options,args) = parser.parse_args()
+ − 45 treePath = str(options.treePath).lstrip().rstrip()
+ − 46 distancePath = str(options.distancePath).lstrip().rstrip()
+ − 47 metadataPath = str(options.metadataPath).lstrip().rstrip()
+ − 48
+ − 49
+ − 50 #region result objects
+ − 51 #define some objects to store values from results
+ − 52 #//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).
+ − 53 class workflowResult(object):
+ − 54 def __init__(self):
+ − 55 self.new = False
+ − 56 self.ID = ""
+ − 57 self.ExpectedSpecies = ""
+ − 58 self.MLSTSpecies = ""
+ − 59 self.SequenceType = ""
+ − 60 self.MLSTScheme = ""
+ − 61 self.CarbapenemResistanceGenes =""
+ − 62 self.OtherAMRGenes=""
+ − 63 self.TotalPlasmids = 0
+ − 64 self.plasmids = []
+ − 65 self.DefinitelyPlasmidContigs =""
+ − 66 self.LikelyPlasmidContigs=""
+ − 67 self.row = ""
+ − 68 class plasmidObj(object):
+ − 69 def __init__(self):
+ − 70 self.PlasmidsID = 0
+ − 71 self.Num_Contigs = 0
+ − 72 self.PlasmidLength = 0
+ − 73 self.PlasmidRepType = ""
+ − 74 self.PlasmidMobility = ""
+ − 75 self.NearestReference = ""
+ − 76
+ − 77 #endregion
+ − 78
+ − 79 #region useful functions
7
+ − 80 def read(path): #read in a text file to a list
1
+ − 81 return [line.rstrip('\n') for line in open(path)]
7
+ − 82 def execute(command): #subprocess.popen call bash command
1
+ − 83 process = subprocess.Popen(command, shell=False, cwd=curDir, universal_newlines=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
+ − 84
+ − 85 # Poll process for new output until finished
+ − 86 while True:
+ − 87 nextline = process.stdout.readline()
+ − 88 if nextline == '' and process.poll() is not None:
+ − 89 break
+ − 90 sys.stdout.write(nextline)
+ − 91 sys.stdout.flush()
+ − 92
+ − 93 output = process.communicate()[0]
+ − 94 exitCode = process.returncode
+ − 95
+ − 96 if (exitCode == 0):
+ − 97 return output
+ − 98 else:
+ − 99 raise subprocess.CalledProcessError(exitCode, command)
7
+ − 100 def httpGetFile(url, filepath=""): #download a file from the web
1
+ − 101 if (filepath == ""):
+ − 102 return urllib.request.urlretrieve(url)
+ − 103 else:
+ − 104 urllib.request.urlretrieve(url, filepath)
+ − 105 return True
7
+ − 106 def gunzip(inputpath="", outputpath=""): #gunzip
1
+ − 107 if (outputpath == ""):
+ − 108 with gzip.open(inputpath, 'rb') as f:
+ − 109 gzContent = f.read()
+ − 110 return gzContent
+ − 111 else:
+ − 112 with gzip.open(inputpath, 'rb') as f:
+ − 113 gzContent = f.read()
+ − 114 with open(outputpath, 'wb') as out:
+ − 115 out.write(gzContent)
+ − 116 return True
7
+ − 117 def addFace(name): #function to add a facet to a tree
6
+ − 118 #if its the reference branch, populate the faces with column headers
+ − 119 face = e.faces.TextFace(name,fsize=10,tight_text=True)
+ − 120 face.border.margin = 5
+ − 121 face.margin_right = 5
+ − 122 face.margin_left = 5
+ − 123 return face
1
+ − 124 #endregion
+ − 125
+ − 126 #region functions to parse result files
+ − 127 def ParseWorkflowResults(pathToResult):
+ − 128 _worflowResult = {}
6
+ − 129 r = pandas.read_csv(pathToResult, delimiter='\t', header=0)
1
+ − 130 r = r.replace(numpy.nan, '', regex=True)
+ − 131 for i in range(len(r.index)):
+ − 132 _results = workflowResult()
6
+ − 133 if(str(r.loc[r.index[i], 'new']).lower() == "new"):
1
+ − 134 _results.new = True
+ − 135 else:
+ − 136 _results.new = False
6
+ − 137 _results.ID = str(r.loc[r.index[i], 'ID'])
+ − 138 _results.ExpectedSpecies = str(r.loc[r.index[i], 'Expected Species'])
+ − 139 _results.MLSTSpecies = str(r.loc[r.index[i], 'MLST Species'])
+ − 140 _results.SequenceType = str(r.loc[r.index[i], 'Sequence Type'])
+ − 141 _results.MLSTScheme = (str(r.loc[r.index[i], 'MLST Scheme']))
+ − 142 _results.CarbapenemResistanceGenes = (str(r.loc[r.index[i], 'Carbapenem Resistance Genes']))
+ − 143 _results.OtherAMRGenes = (str(r.loc[r.index[i], 'Other AMR Genes']))
+ − 144 _results.TotalPlasmids = int(r.loc[r.index[i], 'Total Plasmids'])
1
+ − 145 for j in range(0,_results.TotalPlasmids):
+ − 146 _plasmid = plasmidObj()
6
+ − 147 _plasmid.PlasmidsID =(((str(r.loc[r.index[i], 'Plasmids ID'])).split(";"))[j])
+ − 148 _plasmid.Num_Contigs = (((str(r.loc[r.index[i], 'Num_Contigs'])).split(";"))[j])
+ − 149 _plasmid.PlasmidLength = (((str(r.loc[r.index[i], 'Plasmid Length'])).split(";"))[j])
+ − 150 _plasmid.PlasmidRepType = (((str(r.loc[r.index[i], 'Plasmid RepType'])).split(";"))[j])
+ − 151 _plasmid.PlasmidMobility = ((str(r.loc[r.index[i], 'Plasmid Mobility'])).split(";"))[j]
+ − 152 _plasmid.NearestReference = ((str(r.loc[r.index[i], 'Nearest Reference'])).split(";"))[j]
1
+ − 153 _results.plasmids.append(_plasmid)
6
+ − 154 _results.DefinitelyPlasmidContigs = (str(r.loc[r.index[i], 'Definitely Plasmid Contigs']))
+ − 155 _results.LikelyPlasmidContigs = (str(r.loc[r.index[i], 'Likely Plasmid Contigs']))
1
+ − 156 _results.row = "\t".join(str(x) for x in r.ix[i].tolist())
+ − 157 _worflowResult[_results.ID] = _results
+ − 158 return _worflowResult
+ − 159
+ − 160 #endregion
+ − 161
+ − 162 def Main():
+ − 163 metadata = ParseWorkflowResults(metadataPath)
+ − 164 distance = read(distancePath)
+ − 165 treeFile = "".join(read(treePath))
+ − 166
6
+ − 167 distanceDict = {} #store the distance matrix as rowname:list<string>
1
+ − 168 for i in range(len(distance)):
+ − 169 temp = distance[i].split("\t")
+ − 170 distanceDict[temp[0]] = temp[1:]
6
+ − 171
+ − 172 #region create box tree
+ − 173 #region step5: tree construction
+ − 174 treeFile = "".join(read(treePath))
+ − 175 t = e.Tree(treeFile)
1
+ − 176 t.set_outgroup(t&"Reference")
+ − 177
6
+ − 178 #set the tree style
+ − 179 ts = e.TreeStyle()
+ − 180 ts.show_leaf_name = False
1
+ − 181 ts.show_branch_length = True
+ − 182 ts.scale = 2000 #pixel per branch length unit
+ − 183 ts.branch_vertical_margin = 15 #pixel between branches
6
+ − 184 style2 = e.NodeStyle()
1
+ − 185 style2["fgcolor"] = "#000000"
+ − 186 style2["shape"] = "circle"
+ − 187 style2["vt_line_color"] = "#0000aa"
+ − 188 style2["hz_line_color"] = "#0000aa"
+ − 189 style2["vt_line_width"] = 2
+ − 190 style2["hz_line_width"] = 2
+ − 191 style2["vt_line_type"] = 0 # 0 solid, 1 dashed, 2 dotted
+ − 192 style2["hz_line_type"] = 0
+ − 193 for n in t.traverse():
+ − 194 n.set_style(style2)
+ − 195
6
+ − 196 #find the plasmid origins
1
+ − 197 plasmidIncs = {}
+ − 198 for key in metadata:
+ − 199 for plasmid in metadata[key].plasmids:
+ − 200 for inc in plasmid.PlasmidRepType.split(","):
+ − 201 if (inc.lower().find("inc") > -1):
+ − 202 if not (inc in plasmidIncs):
+ − 203 plasmidIncs[inc] = [metadata[key].ID]
+ − 204 else:
+ − 205 if metadata[key].ID not in plasmidIncs[inc]:
+ − 206 plasmidIncs[inc].append(metadata[key].ID)
+ − 207 #plasmidIncs = sorted(plasmidIncs)
6
+ − 208 for n in t.traverse(): #loop through the nodes of a tree
1
+ − 209 if (n.is_leaf() and n.name == "Reference"):
6
+ − 210 #if its the reference branch, populate the faces with column headers
+ − 211 index = 0
+ − 212 (t&"Reference").add_face(addFace("SampleID"), index, "aligned")
+ − 213 index = index + 1
+ − 214 (t&"Reference").add_face(addFace("New?"), index, "aligned")
+ − 215 index = index + 1
1
+ − 216 for i in range(len(plasmidIncs)): #this loop adds the columns (aka the incs) to the reference node
6
+ − 217 (t&"Reference").add_face(addFace(list(plasmidIncs.keys())[i]), i + index, "aligned")
+ − 218 index = index + len(plasmidIncs)
+ − 219 (t&"Reference").add_face(addFace("MLSTScheme"), index, "aligned")
+ − 220 index = index + 1
+ − 221 (t&"Reference").add_face(addFace("Sequence Type"), index, "aligned")
+ − 222 index = index + 1
+ − 223 (t&"Reference").add_face(addFace("Carbapenamases"), index, "aligned")
+ − 224 index = index + 1
+ − 225 for i in range(len(distanceDict[list(distanceDict.keys())[0]])): #this loop adds the distance matrix
+ − 226 (t&"Reference").add_face(addFace(distanceDict[list(distanceDict.keys())[0]][i]), index + i, "aligned")
+ − 227 index = index + len(distanceDict[list(distanceDict.keys())[0]])
+ − 228 elif (n.is_leaf() and not n.name == "Reference"):
+ − 229 #not reference branches, populate with metadata
+ − 230 index = 0
+ − 231 mData = metadata[n.name.replace(".fa","")]
+ − 232 n.add_face(addFace(mData.ID), index, "aligned")
+ − 233 index = index + 1
+ − 234 if (metadata[n.name.replace(".fa","")].new == True): #new column
+ − 235 face = e.RectFace(30,30,"green","green") # TextFace("Y",fsize=10,tight_text=True)
1
+ − 236 face.border.margin = 5
+ − 237 face.margin_right = 5
+ − 238 face.margin_left = 5
+ − 239 face.vt_align = 1
+ − 240 face.ht_align = 1
6
+ − 241 n.add_face(face, index, "aligned")
+ − 242 index = index + 1
1
+ − 243 for incs in plasmidIncs: #this loop adds presence/absence to the sample nodes
+ − 244 if (n.name.replace(".fa","") in plasmidIncs[incs]):
6
+ − 245 face = e.RectFace(30,30,"black","black") # TextFace("Y",fsize=10,tight_text=True)
1
+ − 246 face.border.margin = 5
+ − 247 face.margin_right = 5
+ − 248 face.margin_left = 5
+ − 249 face.vt_align = 1
+ − 250 face.ht_align = 1
6
+ − 251 n.add_face(face, list(plasmidIncs.keys()).index(incs) + index, "aligned")
+ − 252 index = index + len(plasmidIncs)
+ − 253 n.add_face(addFace(mData.MLSTSpecies), index, "aligned")
+ − 254 index = index + 1
+ − 255 n.add_face(addFace(mData.SequenceType), index, "aligned")
+ − 256 index = index + 1
+ − 257 n.add_face(addFace(mData.CarbapenemResistanceGenes), index, "aligned")
+ − 258 index = index + 1
1
+ − 259 for i in range(len(distanceDict[list(distanceDict.keys())[0]])): #this loop adds distance matrix
6
+ − 260 n.add_face(addFace(list(distanceDict[n.name])[i]), index + i, "aligned")
+ − 261
8
+ − 262 t.render("./tree.pdf", w=5000,units="mm", tree_style=ts) #save it as a png. or an phyloxml
1
+ − 263
+ − 264 #endregion
+ − 265 #endregion
+ − 266
+ − 267
+ − 268 start = time.time()#time the analysis
+ − 269
+ − 270 #analysis time
+ − 271 Main()
+ − 272
+ − 273 end = time.time()
+ − 274 print("Finished!\nThe analysis used: " + str(end-start) + " seconds")