annotate cpo_galaxy_prediction.py @ 22:e03c8f2d2629 draft

planemo upload
author jjjjia
date Tue, 28 Aug 2018 23:21:34 -0400
parents 1543496b2db4
children 2cca036ceb91
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
1 #!/home/jjjjia/.conda/envs/py36/bin/python
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
2
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
3 #$ -S /home/jjjjia/.conda/envs/py36/bin/python
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
4 #$ -V # Pass environment variables to the job
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
5 #$ -N CPO_pipeline # Replace with a more specific job name
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
6 #$ -wd /home/jjjjia/testCases # Use the current working dir
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
7 #$ -pe smp 8 # Parallel Environment (how many cores)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
8 #$ -l h_vmem=11G # Memory (RAM) allocation *per core*
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
9 #$ -e ./logs/$JOB_ID.err
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
10 #$ -o ./logs/$JOB_ID.log
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
11 #$ -m ea
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
12 #$ -M bja20@sfu.ca
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
13
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
14 #./prediction.py -i ~/testCases/cpoResults/contigs/BC11-Kpn005_S2.fa -m ~/testCases/predictionResultsQsubTest/predictions/BC11-Kpn005_S2.mlst -c ~/testCases/predictionResultsQsubTest/predictions/BC11-Kpn005_S2.recon/contig_report.txt -f ~/testCases/predictionResultsQsubTest/predictions/BC11-Kpn005_S2.recon/mobtyper_aggregate_report.txt -a ~/testCases/predictionResultsQsubTest/predictions/BC11-Kpn005_S2.cp -r ~/testCases/predictionResultsQsubTest/predictions/BC11-Kpn005_S2.rgi.txt -e "Klebsiella"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
15 import subprocess
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
16 import pandas
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
17 import optparse
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
18 import os
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
19 import datetime
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
20 import sys
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
21 import time
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
22 import urllib.request
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
23 import gzip
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
24 import collections
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
25 import json
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
26 import numpy
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
27
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
28
6
cabceaa239e4 planemo upload
jjjjia
parents: 5
diff changeset
29 debug = False #debug skips the shell scripts and also dump out a ton of debugging messages
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
30
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
31 if not debug:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
32 #parses some parameters
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
33 parser = optparse.OptionParser("Usage: %prog [options] arg1 arg2 ...")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
34 #required
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
35 #MLSTHIT, mobsuite, resfinder, rgi, mlstscheme
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
36 parser.add_option("-i", "--id", dest="id", type="string", help="identifier of the isolate")
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
37 parser.add_option("-m", "--mlst", dest="mlst", type="string", help="absolute file path to mlst result")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
38 parser.add_option("-c", "--mobfinderContig", dest="mobfinderContig", type="string", help="absolute path to mobfinder aggregate result")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
39 parser.add_option("-f", "--mobfinderAggregate", dest="mobfinderAggregate", type="string", help="absolute path to mobfinder plasmid results")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
40 parser.add_option("-a", "--abricate", dest="abricate", type="string", help="absolute path to abricate results")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
41 parser.add_option("-r", "--rgi", dest="rgi", type="string", help="absolute path to rgi results")
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
42 parser.add_option("-e", "--expected", dest="expectedSpecies", default="NA/NA/NA", type="string", help="expected species of the isolate")
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
43 parser.add_option("-s", "--mlst-scheme", dest="mlstScheme", default= "./scheme_species_map.tab", type="string", help="absolute file path to mlst scheme")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
44 parser.add_option("-p", "--plasmidfinder", dest="plasmidfinder", type="string", help="absolute file path to plasmidfinder ")
21
1543496b2db4 planemo upload
jjjjia
parents: 18
diff changeset
45 parser.add_option("-d", "--mash", dest="mash", type="string", help="absolute file path to mash plasmiddb result")
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
46
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
47 #parallelization, useless, these are hard coded to 8cores/64G RAM
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
48 #parser.add_option("-t", "--threads", dest="threads", default=8, type="int", help="number of cpu to use")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
49 #parser.add_option("-p", "--memory", dest="memory", default=64, type="int", help="memory to use in GB")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
50
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
51 (options,args) = parser.parse_args()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
52 #if len(args) != 8:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
53 #parser.error("incorrect number of arguments, all 7 is required")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
54 curDir = os.getcwd()
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
55 ID = str(options.id).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
56 mlst = str(options.mlst).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
57 mobfindercontig = str(options.mobfinderContig).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
58 mobfinderaggregate = str(options.mobfinderAggregate).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
59 abricate = str(options.abricate).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
60 rgi = str(options.rgi).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
61 expectedSpecies = str(options.expectedSpecies).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
62 mlstScheme = str(options.mlstScheme).lstrip().rstrip()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
63 plasmidfinder = str(options.plasmidfinder).lstrip().rstrip()
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
64 mash = str(options.mash).lstrip().rstrip()
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
65 outputDir = "./"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
66 print(mlst)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
67 print(mobfindercontig)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
68 print(mobfinderaggregate)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
69 print(abricate)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
70 print(rgi)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
71 print(expectedSpecies)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
72 print(mlstScheme)
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
73 print(mash)
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
74
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
75 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
76 curDir = os.getcwd()
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
77 ID = "BC11"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
78 mlst = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.mlst"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
79 mobfindercontig = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.recon\contig_report.txt"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
80 mobfinderaggregate = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.recon\mobtyper_aggregate_report.txt"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
81 abricate = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.cp"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
82 rgi = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.rgi.txt"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
83 expectedSpecies = "Escherichia coli"
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
84 mlstScheme = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\scheme_species_map.tab"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
85 plasmidfinder = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\BC11-Kpn005_S2.origins"
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
86 mash = "D:\OneDrive\ProjectCDC\ProjectCDCInPython\ProjectCDCInPython\pipelineTest\predictions\mash.tsv"
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
87 outputDir = "./"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
88
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
89 #region result objects
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
90 #define some objects to store values from results
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
91 #//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).
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
92 class starFinders(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
93 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
94 self.file = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
95 self.sequence = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
96 self.start = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
97 self.end = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
98 self.gene = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
99 self.shortGene = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
100 self.coverage = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
101 self.coverage_map = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
102 self.gaps = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
103 self.pCoverage = 100.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
104 self.pIdentity = 100.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
105 self.database = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
106 self.accession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
107 self.product = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
108 self.source = "chromosome"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
109 self.row = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
110
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
111 class PlasFlowResult(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
112 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
113 self.sequence = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
114 self.length = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
115 self.label = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
116 self.confidence = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
117 self.usefulRow = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
118 self.row = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
119
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
120 class MlstResult(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
121 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
122 self.file = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
123 self.speciesID = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
124 self.seqType = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
125 self.scheme = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
126 self.species = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
127 self.row=""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
128
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
129 class mobsuiteResult(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
130 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
131 self.file_id = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
132 self.cluster_id = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
133 self.contig_id = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
134 self.contig_num = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
135 self.contig_length = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
136 self.circularity_status = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
137 self.rep_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
138 self.rep_type_accession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
139 self.relaxase_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
140 self.relaxase_type_accession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
141 self.mash_nearest_neighbor = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
142 self.mash_neighbor_distance = 0.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
143 self.repetitive_dna_id = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
144 self.match_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
145 self.score = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
146 self.contig_match_start = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
147 self.contig_match_end = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
148 self.row = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
149
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
150 class mobsuitePlasmids(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
151 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
152 self.file_id = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
153 self.num_contigs = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
154 self.total_length = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
155 self.gc = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
156 self.rep_types = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
157 self.rep_typeAccession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
158 self.relaxase_type= ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
159 self.relaxase_type_accession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
160 self.mpf_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
161 self.mpf_type_accession= ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
162 self.orit_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
163 self.orit_accession = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
164 self.PredictedMobility = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
165 self.mash_nearest_neighbor = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
166 self.mash_neighbor_distance = 0.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
167 self.mash_neighbor_cluster= 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
168 self.row = ""
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
169
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
170 class RGIResult(object):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
171 def __init__(self):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
172 self.ORF_ID = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
173 self.Contig = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
174 self.Start = -1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
175 self.Stop = -1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
176 self.Orientation = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
177 self.Cut_Off = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
178 self.Pass_Bitscore = 100000
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
179 self.Best_Hit_Bitscore = 0.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
180 self.Best_Hit_ARO = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
181 self.Best_Identities = 0.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
182 self.ARO = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
183 self.Model_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
184 self.SNPs_in_Best_Hit_ARO = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
185 self.Other_SNPs = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
186 self.Drug_Class = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
187 self.Resistance_Mechanism = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
188 self.AMR_Gene_Family = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
189 self.Predicted_DNA = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
190 self.Predicted_Protein = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
191 self.CARD_Protein_Sequence = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
192 self.Percentage_Length_of_Reference_Sequence = 0.00
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
193 self.ID = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
194 self.Model_ID = 0
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
195 self.source = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
196 self.row = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
197
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
198 class MashResult(object):
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
199 def __init__(self):
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
200 self.size = 0.0
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
201 self.depth = 0.0
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
202 self.identity = 0.0
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
203 self.sharedHashes = ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
204 self.medianMultiplicity = 0
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
205 self.pvalue = 0.0
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
206 self.queryID= ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
207 self.queryComment = ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
208 self.species = ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
209 self.row = ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
210 self.accession = ""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
211 self.gcf=""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
212 self.assembly=""
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
213
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
214 def toDict(self): #doesnt actually work
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
215 return dict((name, getattr(self, name)) for name in dir(self) if not name.startswith('__'))
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
216
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
217
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
218 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
219
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
220 #region useful functions
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
221 def read(path):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
222 return [line.rstrip('\n') for line in open(path)]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
223 def execute(command):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
224 process = subprocess.Popen(command, shell=False, cwd=curDir, universal_newlines=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
225
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
226 # Poll process for new output until finished
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
227 while True:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
228 nextline = process.stdout.readline()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
229 if nextline == '' and process.poll() is not None:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
230 break
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
231 sys.stdout.write(nextline)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
232 sys.stdout.flush()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
233
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
234 output = process.communicate()[0]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
235 exitCode = process.returncode
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
236
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
237 if (exitCode == 0):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
238 return output
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
239 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
240 raise subprocess.CalledProcessError(exitCode, command)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
241 def httpGetFile(url, filepath=""):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
242 if (filepath == ""):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
243 return urllib.request.urlretrieve(url)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
244 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
245 urllib.request.urlretrieve(url, filepath)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
246 return True
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
247 def gunzip(inputpath="", outputpath=""):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
248 if (outputpath == ""):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
249 with gzip.open(inputpath, 'rb') as f:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
250 gzContent = f.read()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
251 return gzContent
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
252 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
253 with gzip.open(inputpath, 'rb') as f:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
254 gzContent = f.read()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
255 with open(outputpath, 'wb') as out:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
256 out.write(gzContent)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
257 return True
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
258 def ToJson(dictObject, outputPath):
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
259 #outDir = outputDir + '/summary/' + ID + ".json/"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
260 #if not (os.path.exists(outDir)):
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
261 #os.makedirs(outDir)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
262 #with open(outputPath, 'w') as f:
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
263 #json.dump([ob.__dict__ for ob in dictObject.values()], f, ensure_ascii=False)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
264 return ""
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
265 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
266
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
267 #region functions to parse result files
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
268 def ParseMLSTResult(pathToMLSTResult, scheme):
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
269 _mlstResult = {}
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
270 scheme = pandas.read_csv(scheme, delimiter='\t', header=0)
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
271 scheme = scheme.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
272
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
273 taxon = {}
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
274 #record the scheme as a dictionary
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
275 taxon["-"] = "No MLST Match"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
276 for i in range(len(scheme.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
277 key = scheme.iloc[i,0]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
278 if (str(scheme.iloc[i,2]) == "nan"):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
279 value = str(scheme.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
280 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
281 value = str(scheme.iloc[i,1]) + " " + str(scheme.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
282
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
283 if (key in taxon.keys()):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
284 taxon[key] = taxon.get(key) + ";" + value
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
285 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
286 taxon[key] = value
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
287 #read in the mlst result
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
288 mlst = pandas.read_csv(pathToMLSTResult, delimiter='\t', header=None)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
289 _mlstHit = MlstResult()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
290
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
291 _mlstHit.file = mlst.iloc[0,0]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
292 _mlstHit.speciesID = (mlst.iloc[0,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
293 _mlstHit.seqType = str(mlst.iloc[0,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
294 for i in range(3, len(mlst.columns)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
295 _mlstHit.scheme += mlst.iloc[0,i] + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
296 _mlstHit.species = taxon[_mlstHit.speciesID]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
297 _mlstHit.row = "\t".join(str(x) for x in mlst.ix[0].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
298 _mlstResult[_mlstHit.speciesID]=_mlstHit
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
299
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
300 return _mlstResult
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
301
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
302 def ParsePlasmidFinderResult(pathToPlasmidFinderResult):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
303 #pipelineTest/contigs/BC110-Kpn005.fa contig00019 45455 45758 IncFIC(FII)_1 8-308/499 ========/=..... 8/11 59.52 75.65 plasmidfinder AP001918 IncFIC(FII)_1__AP001918
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
304 #example resfinder:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
305 #pipelineTest/contigs/BC110-Kpn005.fa contig00038 256 1053 OXA-181 1-798/798 =============== 0/0 100.00 100.00 bccdc AEP16366.1 OXA-48 family carbapenem-hydrolyzing class D beta-lactamase OXA-181
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
306
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
307 _pFinder = {} #***********************
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
308 plasmidFinder = pandas.read_csv(pathToPlasmidFinderResult, delimiter='\t', header=0)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
309 plasmidFinder = plasmidFinder.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
310
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
311
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
312 for i in range(len(plasmidFinder.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
313 pf = starFinders()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
314 pf.file = str(plasmidFinder.iloc[i,0])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
315 pf.sequence = str(plasmidFinder.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
316 pf.start = int(plasmidFinder.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
317 pf.end = int(plasmidFinder.iloc[i,3])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
318 pf.gene = str(plasmidFinder.iloc[i,4])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
319 pf.shortGene = pf.gene[:pf.gene.index("_")]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
320 pf.coverage = str(plasmidFinder.iloc[i,5])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
321 pf.coverage_map = str(plasmidFinder.iloc[i,6])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
322 pf.gaps = str(plasmidFinder.iloc[i,7])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
323 pf.pCoverage = float(plasmidFinder.iloc[i,8])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
324 pf.pIdentity = float(plasmidFinder.iloc[i,9])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
325 pf.database = str(plasmidFinder.iloc[i,10])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
326 pf.accession = str(plasmidFinder.iloc[i,11])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
327 pf.product = str(plasmidFinder.iloc[i,12])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
328 pf.source = "plasmid"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
329 pf.row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
330 _pFinder[pf.gene]=pf
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
331 #row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
332 #plasmidFinderContigs.append(str(plasmidFinder.iloc[i,1]))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
333 #origins.append(str(plasmidFinder.iloc[i,4][:plasmidFinder.iloc[i,4].index("_")]))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
334 return _pFinder
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
335
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
336 def ParseMobsuiteResult(pathToMobsuiteResult):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
337 _mobsuite = {}
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
338 mResult = pandas.read_csv(pathToMobsuiteResult, delimiter='\t', header=0)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
339 mResult = mResult.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
340
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
341 for i in range(len(mResult.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
342 mr = mobsuiteResult()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
343 mr.file_id = str(mResult.iloc[i,0])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
344 mr.cluster_id = str(mResult.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
345 if (mr.cluster_id == "chromosome"):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
346 break
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
347 mr.contig_id = str(mResult.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
348 mr.contig_num = mr.contig_id[(mr.contig_id.find("contig")+6):mr.contig_id.find("_len=")]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
349 mr.contig_length = int(mResult.iloc[i,3])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
350 mr.circularity_status = str(mResult.iloc[i,4])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
351 mr.rep_type = str(mResult.iloc[i,5])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
352 mr.rep_type_accession = str(mResult.iloc[i,6])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
353 mr.relaxase_type = str(mResult.iloc[i,7])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
354 mr.relaxase_type_accession = str(mResult.iloc[i,8])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
355 mr.mash_nearest_neighbor = str(mResult.iloc[i,9])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
356 mr.mash_neighbor_distance = float(mResult.iloc[i,10])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
357 mr.repetitive_dna_id = str(mResult.iloc[i,11])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
358 mr.match_type = str(mResult.iloc[i,12])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
359 if (mr.match_type == ""):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
360 mr.score = -1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
361 mr.contig_match_start = -1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
362 mr.contig_match_end = -1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
363 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
364 mr.score = int(mResult.iloc[i,13])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
365 mr.contig_match_start = int(mResult.iloc[i,14])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
366 mr.contig_match_end = int(mResult.iloc[i,15])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
367 mr.row = "\t".join(str(x) for x in mResult.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
368 _mobsuite[mr.contig_id]=(mr)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
369 return _mobsuite
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
370
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
371 def ParseMobsuitePlasmids(pathToMobsuiteResult):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
372 _mobsuite = {}
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
373 mResults = pandas.read_csv(pathToMobsuiteResult, delimiter='\t', header=0)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
374 mResults = mResults.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
375
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
376 for i in range(len(mResults.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
377 mr = mobsuitePlasmids()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
378 mr.file_id = str(mResults.iloc[i,0])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
379 mr.num_contigs = int(mResults.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
380 mr.total_length = int(mResults.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
381 mr.gc = int(mResults.iloc[i,3])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
382 mr.rep_types = str(mResults.iloc[i,4])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
383 mr.rep_typeAccession = str(mResults.iloc[i,5])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
384 mr.relaxase_type = str(mResults.iloc[i,6])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
385 mr.relaxase_type_accession = str(mResults.iloc[i,7])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
386 mr.mpf_type = str(mResults.iloc[i,8])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
387 mr.mpf_type_accession = str(mResults.iloc[i,9])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
388 mr.orit_type = str(mResults.iloc[i,10])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
389 mr.orit_accession = str(mResults.iloc[i,11])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
390 mr.PredictedMobility = str(mResults.iloc[i,12])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
391 mr.mash_nearest_neighbor = str(mResults.iloc[i,13])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
392 mr.mash_neighbor_distance = float(mResults.iloc[i,14])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
393 mr.mash_neighbor_cluster = int(mResults.iloc[i,15])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
394 mr.row = "\t".join(str(x) for x in mResults.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
395 _mobsuite[mr.file_id] = mr
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
396 return _mobsuite
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
397
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
398 def ParseResFinderResult(pathToResFinderResults, plasmidContigs, likelyPlasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
399 _rFinder = {}
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
400 resFinder = pandas.read_csv(pathToResFinderResults, delimiter='\t', header=0)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
401 resFinder = resFinder.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
402
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
403 for i in range(len(resFinder.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
404 rf = starFinders()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
405 rf.file = str(resFinder.iloc[i,0])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
406 rf.sequence = str(resFinder.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
407 rf.start = int(resFinder.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
408 rf.end = int(resFinder.iloc[i,3])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
409 rf.gene = str(resFinder.iloc[i,4])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
410 rf.shortGene = rf.gene
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
411 rf.coverage = str(resFinder.iloc[i,5])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
412 rf.coverage_map = str(resFinder.iloc[i,6])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
413 rf.gaps = str(resFinder.iloc[i,7])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
414 rf.pCoverage = float(resFinder.iloc[i,8])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
415 rf.pIdentity = float(resFinder.iloc[i,9])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
416 rf.database = str(resFinder.iloc[i,10])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
417 rf.accession = str(resFinder.iloc[i,11])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
418 rf.product = str(resFinder.iloc[i,12])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
419 rf.row = "\t".join(str(x) for x in resFinder.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
420 if (rf.sequence[6:] in plasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
421 rf.source = "plasmid"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
422 elif (rf.sequence[6:] in likelyPlasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
423 rf.source = "likely plasmid"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
424 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
425 rf.source = "likely chromosome"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
426 _rFinder[rf.gene]=rf
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
427 return _rFinder
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
428
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
429 def ParseRGIResult(pathToRGIResults, plasmidContigs, likelyPlasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
430 _rgiR = {}
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
431 RGI = pandas.read_csv(pathToRGIResults, delimiter='\t', header=0)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
432 RGI = RGI.replace(numpy.nan, '', regex=True)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
433
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
434 for i in range(len(RGI.index)):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
435 r = RGIResult()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
436 r.ORF_ID = str(RGI.iloc[i,0])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
437 r.Contig = str(RGI.iloc[i,1])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
438 r.Contig_Num = r.Contig[6:r.Contig.find("_")]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
439 r.Start = int(RGI.iloc[i,2])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
440 r.Stop = int(RGI.iloc[i,3])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
441 r.Orientation = str(RGI.iloc[i,4])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
442 r.Cut_Off = str(RGI.iloc[i,5])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
443 r.Pass_Bitscore = int(RGI.iloc[i,6])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
444 r.Best_Hit_Bitscore = float(RGI.iloc[i,7])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
445 r.Best_Hit_ARO = str(RGI.iloc[i,8])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
446 r.Best_Identities = float(RGI.iloc[i,9])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
447 r.ARO = int(RGI.iloc[i,10])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
448 r.Model_type = str(RGI.iloc[i,11])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
449 r.SNPs_in_Best_Hit_ARO = str(RGI.iloc[i,12])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
450 r.Other_SNPs = str(RGI.iloc[i,13])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
451 r.Drug_Class = str(RGI.iloc[i,14])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
452 r.Resistance_Mechanism = str(RGI.iloc[i,15])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
453 r.AMR_Gene_Family = str(RGI.iloc[i,16])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
454 r.Predicted_DNA = str(RGI.iloc[i,17])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
455 r.Predicted_Protein = str(RGI.iloc[i,18])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
456 r.CARD_Protein_Sequence = str(RGI.iloc[i,19])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
457 r.Percentage_Length_of_Reference_Sequence = float(RGI.iloc[i,20])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
458 r.ID = str(RGI.iloc[i,21])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
459 r.Model_ID = int(RGI.iloc[i,22])
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
460 r.row = "\t".join(str(x) for x in RGI.ix[i].tolist())
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
461 if (r.Contig_Num in plasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
462 r.source = "plasmid"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
463 elif (r.Contig_Num in likelyPlasmidContigs):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
464 r.source = "likely plasmid"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
465 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
466 r.source = "likely chromosome"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
467 _rgiR[r.Model_ID]=r
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
468 return _rgiR
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
469
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
470 def ParsePlasmidFinderResult(pathToPlasmidFinderResult):
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
471 #pipelineTest/contigs/BC110-Kpn005.fa contig00019 45455 45758 IncFIC(FII)_1 8-308/499 ========/=..... 8/11 59.52 75.65 plasmidfinder AP001918 IncFIC(FII)_1__AP001918
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
472 #example resfinder:
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
473 #pipelineTest/contigs/BC110-Kpn005.fa contig00038 256 1053 OXA-181 1-798/798 =============== 0/0 100.00 100.00 bccdc AEP16366.1 OXA-48 family carbapenem-hydrolyzing class D beta-lactamase OXA-181
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
474
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
475 _pFinder = {} #***********************
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
476 plasmidFinder = pandas.read_csv(pathToPlasmidFinderResult, delimiter='\t', header=0)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
477
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
478 for i in range(len(plasmidFinder.index)):
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
479 pf = starFinders()
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
480 pf.file = str(plasmidFinder.iloc[i,0])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
481 pf.sequence = str(plasmidFinder.iloc[i,1])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
482 pf.start = int(plasmidFinder.iloc[i,2])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
483 pf.end = int(plasmidFinder.iloc[i,3])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
484 pf.gene = str(plasmidFinder.iloc[i,4])
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
485 if (pf.gene.find("_") > -1):
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
486 pf.shortGene = pf.gene[:pf.gene.index("_")]
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
487 else:
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
488 pf.shortGene = pf.gene
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
489 pf.coverage = str(plasmidFinder.iloc[i,5])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
490 pf.coverage_map = str(plasmidFinder.iloc[i,6])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
491 pf.gaps = str(plasmidFinder.iloc[i,7])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
492 pf.pCoverage = float(plasmidFinder.iloc[i,8])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
493 pf.pIdentity = float(plasmidFinder.iloc[i,9])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
494 pf.database = str(plasmidFinder.iloc[i,10])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
495 pf.accession = str(plasmidFinder.iloc[i,11])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
496 pf.product = str(plasmidFinder.iloc[i,12])
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
497 pf.source = "plasmid"
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
498 pf.row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
499 _pFinder[pf.gene]=pf
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
500 #row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
501 #plasmidFinderContigs.append(str(plasmidFinder.iloc[i,1]))
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
502 #origins.append(str(plasmidFinder.iloc[i,4][:plasmidFinder.iloc[i,4].index("_")]))
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
503 return _pFinder
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
504
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
505 def ParseMashResult(pathToMashScreen):
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
506 mashScreen = pandas.read_csv(pathToMashScreen, delimiter='\t', header=None)
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
507
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
508 _mashPlasmidHits = {} #***********************
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
509 #parse what the species are.
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
510 for i in (range(len(mashScreen.index))):
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
511 mr = MashResult()
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
512 mr.identity = float(mashScreen.ix[i, 0])
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
513 mr.sharedHashes = mashScreen.ix[i, 1]
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
514 mr.medianMultiplicity = int(mashScreen.ix[i, 2])
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
515 mr.pvalue = float(mashScreen.ix[i, 3])
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
516 mr.name = mashScreen.ix[i, 4] #accession
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
517 mr.row = "\t".join(str(x) for x in mashScreen.ix[i].tolist())
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
518 _mashPlasmidHits[mr.name] = mr
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
519 return _mashPlasmidHits
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
520 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
521
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
522 def Main():
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
523 outputDir = "./"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
524 notes = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
525 #init the output list
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
526 output = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
527 jsonOutput = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
528
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
529 print(str(datetime.datetime.now()) + "\n\nID: " + ID + "\nAssembly: " + ID)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
530 output.append(str(datetime.datetime.now()) + "\n\nID: " + ID + "\nAssembly: " + ID)
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
531
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
532 #region parse the mlst results
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
533 print("step 3: parsing mlst, plasmid, and amr results")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
534
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
535 print("identifying MLST")
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
536 mlstHit = ParseMLSTResult(mlst, str(mlstScheme))#***********************
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
537 ToJson(mlstHit, "mlst.json") #write it to a json output
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
538 mlstHit = list(mlstHit.values())[0]
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
539
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
540 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
541
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
542 #region parse mobsuite, resfinder and rgi results
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
543 print("identifying plasmid contigs and amr genes")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
544
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
545 plasmidContigs = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
546 likelyPlasmidContigs = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
547 origins = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
548
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
549 #parse mobsuite results
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
550 mSuite = ParseMobsuiteResult(mobfindercontig) #outputDir + "/predictions/" + ID + ".recon/contig_report.txt")#*************
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
551 ToJson(mSuite, "mobsuite.json") #*************
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
552 mSuitePlasmids = ParseMobsuitePlasmids(mobfinderaggregate)#outputDir + "/predictions/" + ID + ".recon/mobtyper_aggregate_report.txt")#*************
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
553 ToJson(mSuitePlasmids, "mobsuitePlasmids.json") #*************
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
554
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
555 for key in mSuite:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
556 if mSuite[key].contig_num not in plasmidContigs and mSuite[key].contig_num not in likelyPlasmidContigs:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
557 if not (mSuite[key].rep_type == ''):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
558 plasmidContigs.append(mSuite[key].contig_num)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
559 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
560 likelyPlasmidContigs.append(mSuite[key].contig_num)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
561 for key in mSuite:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
562 if mSuite[key].rep_type not in origins:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
563 origins.append(mSuite[key].rep_type)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
564
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
565 #parse resfinder AMR results
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
566 pFinder = ParsePlasmidFinderResult(plasmidfinder)
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
567 ToJson(pFinder, "origins.json")
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
568
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
569 rFinder = ParseResFinderResult(abricate, plasmidContigs, likelyPlasmidContigs)#outputDir + "/predictions/" + ID + ".cp", plasmidContigs, likelyPlasmidContigs) #**********************
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
570 ToJson(rFinder, "resfinder.json") #*************
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
571
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
572 rgiAMR = ParseRGIResult(rgi, plasmidContigs, likelyPlasmidContigs) # outputDir + "/predictions/" + ID + ".rgi.txt", plasmidContigs, likelyPlasmidContigs)#***********************
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
573 ToJson(rgiAMR, "rgi.json") #*************
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
574
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
575 plasmidFamily = ParseMashResult(mash)
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
576 ToJson(plasmidFamily, "mash.json")
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
577
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
578 carbapenamases = []
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
579 resfinderCarbas = [] #list of rfinder objects for lindaout list
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
580 amrGenes = []
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
581 for keys in rFinder:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
582 carbapenamases.append(rFinder[keys].shortGene + "(" + rFinder[keys].source + ")")
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
583 resfinderCarbas.append(rFinder[keys])
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
584 for keys in rgiAMR:
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
585 if (rgiAMR[keys].Drug_Class.find("carbapenem") > -1 and rgiAMR[keys].AMR_Gene_Family.find("beta-lactamase") > -1):
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
586 if (rgiAMR[keys].Best_Hit_ARO not in carbapenamases):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
587 carbapenamases.append(rgiAMR[keys].Best_Hit_ARO+ "(" + rgiAMR[keys].source + ")")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
588 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
589 if (rgiAMR[keys].Best_Hit_ARO not in amrGenes):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
590 amrGenes.append(rgiAMR[keys].Best_Hit_ARO+ "(" + rgiAMR[keys].source + ")")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
591 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
592
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
593 #region output parsed mlst information
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
594 print("formatting mlst outputs")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
595 output.append("\n\n\n~~~~~~~MLST summary~~~~~~~")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
596 output.append("MLST determined species: " + mlstHit.species)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
597 output.append("\nMLST Details: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
598 output.append(mlstHit.row)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
599
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
600 output.append("\nMLST information: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
601 if (mlstHit.species == expectedSpecies):
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
602 output.append("MLST determined species is the same as expected species")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
603 #notes.append("MLST determined species is the same as expected species")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
604 else:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
605 output.append("!!!MLST determined species is NOT the same as expected species, contamination? mislabeling?")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
606 notes.append("MLST: Not expected species. Possible contamination or mislabeling")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
607
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
608 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
609
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
610 #region output the parsed plasmid/amr results
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
611 output.append("\n\n\n~~~~~~~~Plasmids~~~~~~~~\n")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
612
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
613 output.append("predicted plasmid origins: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
614 output.append(";".join(origins))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
615
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
616 output.append("\ndefinitely plasmid contigs")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
617 output.append(";".join(plasmidContigs))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
618
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
619 output.append("\nlikely plasmid contigs")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
620 output.append(";".join(likelyPlasmidContigs))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
621
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
622 output.append("\nmob-suite prediction details: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
623 for key in mSuite:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
624 output.append(mSuite[key].row)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
625
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
626 output.append("\n\n\n~~~~~~~~AMR Genes~~~~~~~~\n")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
627 output.append("predicted carbapenamase Genes: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
628 output.append(",".join(carbapenamases))
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
629 output.append("other RGI AMR Genes: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
630 for key in rgiAMR:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
631 output.append(rgiAMR[key].Best_Hit_ARO + "(" + rgiAMR[key].source + ")")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
632
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
633 output.append("\nDetails about the carbapenamase Genes: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
634 for key in rFinder:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
635 output.append(rFinder[key].row)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
636 output.append("\nDetails about the RGI AMR Genes: ")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
637 for key in rgiAMR:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
638 output.append(rgiAMR[key].row)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
639
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
640 #write summary to a file
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
641 summaryDir = outputDir + "/summary/" + ID
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
642 out = open("summary.txt", 'w')
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
643 for item in output:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
644 out.write("%s\n" % item)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
645
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
646
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
647 #TSV output
6
cabceaa239e4 planemo upload
jjjjia
parents: 5
diff changeset
648 lindaOut = []
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
649 tsvOut = []
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
650 lindaOut.append("ID\tQUALITY\tExpected Species\tMLST Scheme\tSequence Type\tMLST_ALLELE_1\tMLST_ALLELE_2\tMLST_ALLELE_3\tMLST_ALLELE_4\tMLST_ALLELE_5\tMLST_ALLELE_6\tMLST_ALLELE_7\tSEROTYPE\tK_CAPSULE\tPLASMID_2_RFLP\tPLASMID_1_FAMILY\tPLASMID_1_BEST_MATCH\tPLASMID_1_COVERAGE\tPLASMID_1_SNVS_TO_BEST_MATCH\tPLASMID_1_CARBAPENEMASE\tPLASMID_1_INC_GROUP\tPLASMID_2_RFLP\tPLASMID_2_FAMILY\tPLASMID_2_BEST_MATCH\tPLASMID_2_COVERAGE\tPLASMID_2_SNVS_TO_BEST_MATCH\tPLASMID_2_CARBAPENEMASE\tPLASMID_2_INC_GROUP")
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
651 lindaTemp = ID + "\t" #id
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
652 lindaTemp += "\t" #quality
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
653 lindaTemp += expectedSpecies + "\t" #expected
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
654 lindaTemp += mlstHit.species + "\t" #mlstscheme
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
655 lindaTemp += str(mlstHit.seqType) + "\t" #seq type
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
656 lindaTemp += "\t".join(mlstHit.scheme.split(";")) + "\t"#mlst alleles x 7
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
657 lindaTemp += "\t\t" #sero and kcap
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
658
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
659 #resfinderCarbas
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
660 index = 0
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
661 for carbs in resfinderCarbas:
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
662 if (carbs.source == "plasmid"): #
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
663 lindaTemp += "\t"
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
664 plasmid = plasmidFamily[list(plasmidFamily.keys())[index]]
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
665 lindaTemp += plasmid.name + "\t"
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
666 lindaTemp += str(plasmid.identity) + "\t"
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
667 lindaTemp += plasmid.sharedHashes + "\t"
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
668 lindaTemp += carbs.shortGene + "\t" #found an carbapenase
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
669 contig = carbs.sequence[6:] #this is the contig number
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
670 for i in mSuite.keys():
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
671 if (str(mSuite[i].contig_num) == str(contig)): #found the right plasmid
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
672 clusterid = mSuite[i].cluster_id
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
673 rep_types = mSuitePlasmids["plasmid_" + str(clusterid) + ".fasta"].rep_types
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
674 lindaTemp += rep_types
13
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
675 lindaOut.append(lindaTemp)
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
676 out = open("summary.linda.tsv", 'w')
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
677 for item in lindaOut:
a14b12a71a53 planemo upload
jjjjia
parents: 6
diff changeset
678 out.write("%s\n" % item)
6
cabceaa239e4 planemo upload
jjjjia
parents: 5
diff changeset
679
18
596bf8a792de planemo upload
jjjjia
parents: 13
diff changeset
680 tsvOut.append("new\tID\tExpected Species\tMLST Species\tSequence Type\tMLST Scheme\tCarbapenem Resistance Genes\tOther AMR Genes\tPlasmid Best Match\tTotal Plasmids\tPlasmids ID\tNum_Contigs\tPlasmid Length\tPlasmid RepType\tPlasmid Mobility\tNearest Reference\tDefinitely Plasmid Contigs\tLikely Plasmid Contigs")
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
681 #start with ID
6
cabceaa239e4 planemo upload
jjjjia
parents: 5
diff changeset
682 temp = "\t"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
683 temp += (ID + "\t")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
684 temp += expectedSpecies + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
685
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
686 #move into MLST
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
687 temp += mlstHit.species + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
688 temp += str(mlstHit.seqType) + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
689 temp += mlstHit.scheme + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
690
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
691 #now onto AMR genes
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
692 temp += ";".join(carbapenamases) + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
693 temp += ";".join(amrGenes) + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
694
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
695 #lastly plasmids
22
e03c8f2d2629 planemo upload
jjjjia
parents: 21
diff changeset
696 temp += str(plasmidFamily[list(plasmidFamily.keys())[0]].name) + "\t"
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
697 temp+= str(len(mSuitePlasmids)) + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
698 plasmidID = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
699 contigs = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
700 lengths = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
701 rep_type = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
702 mobility = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
703 neighbour = ""
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
704 for keys in mSuitePlasmids:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
705 plasmidID += str(mSuitePlasmids[keys].mash_neighbor_cluster) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
706 contigs += str(mSuitePlasmids[keys].num_contigs) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
707 lengths += str(mSuitePlasmids[keys].total_length) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
708 rep_type += str(mSuitePlasmids[keys].rep_types) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
709 mobility += str(mSuitePlasmids[keys].PredictedMobility) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
710 neighbour += str(mSuitePlasmids[keys].mash_nearest_neighbor) + ";"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
711 temp += plasmidID + "\t" + contigs + "\t" + lengths + "\t" + rep_type + "\t" + mobility + "\t" + neighbour + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
712 temp += ";".join(plasmidContigs) + "\t"
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
713 temp += ";".join(likelyPlasmidContigs)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
714 tsvOut.append(temp)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
715
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
716 summaryDir = outputDir + "/summary/" + ID
3
e6027598a35c planemo upload
jjjjia
parents: 2
diff changeset
717 out = open("summary.tsv", 'w')
1
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
718 for item in tsvOut:
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
719 out.write("%s\n" % item)
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
720 #endregion
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
721
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
722
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
723 start = time.time()#time the analysis
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
724 print("Starting workflow...")
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
725 #analysis time
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
726 Main()
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
727
fea89c4d5227 Uploaded
jjjjia
parents:
diff changeset
728 end = time.time()
5
698579246d0d planemo upload
jjjjia
parents: 4
diff changeset
729 print("Finished!\nThe analysis used: " + str(end-start) + " seconds")