annotate cpo_galaxy_tree_sensitive.py @ 21:1543496b2db4 draft

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