annotate cpo_galaxy_prediction.py @ 1:fea89c4d5227 draft

Uploaded

1 #!/home/jjjjia/.conda/envs/py36/bin/python
2 
3 #$ -S /home/jjjjia/.conda/envs/py36/bin/python
4 #$ -V             # Pass environment variables to the job
5 #$ -N CPO_pipeline    # Replace with a more specific job name
6 #$ -wd /home/jjjjia/testCases           # Use the current working dir
7 #$ -pe smp 8      # Parallel Environment (how many cores)
8 #$ -l h_vmem=11G  # Memory (RAM) allocation *per core*
9 #$ -e ./logs/$JOB_ID.err
10 #$ -o ./logs/$JOB_ID.log
11 #$ -m ea
12 #$ -M [e-mail]
13 
14 #~/scripts/pipeline.py -i BC11-Kpn005_S2 -f /data/jjjjia/R1/BC11-Kpn005_S2_L001_R1_001.fastq.gz -r /data/jjjjia/R2/BC11-Kpn005_S2_L001_R2_001.fastq.gz -o pipelineResultsQsub -e "Klebsiella pneumoniae" 
15 
16 import subprocess
17 import pandas
18 import optparse
19 import os
20 import datetime
21 import sys
22 import time
23 import urllib.request
24 import gzip
25 import collections
26 import json
27 import numpy
28 
29 
30 debug = True #debug skips the shell scripts and also dump out a ton of debugging messages
31 
32 if not debug:
33     #parses some parameters
34     parser = optparse.OptionParser("Usage: %prog [options] arg1 arg2 ...")
35     #required
36     parser.add_option("-i", "--id", dest="id", type="string", help="identifier of the isolate")    
37     parser.add_option("-a", "--assembly", dest="assemblyPath", type="string", help="absolute file path to contigs fasta")
38     parser.add_option("-c", "--card-db", dest="cardDB", default = "/home/jjjjia/databases/card202.json", type="string", help="absolute path to card reference database")
39     parser.add_option("-o", "--output", dest="output", default='./', type="string", help="absolute path to output folder")
40     parser.add_option("-e", "--expected", dest="expectedSpecies", default="NA/NA/NA", type="string", help="expected species of the isolate")
41 
42     #optionals
43     parser.add_option("-k", "--script-path", dest="scriptDir", default="/home/jjjjia/scripts", type="string", help="absolute file path to this script folder")
44     parser.add_option("-b", "--update-abricate-path", dest="updateAbPath", default = "", type="string", help="absolute file path to fasta sequence used for abricate database")
45     parser.add_option("-m", "--update-abricate-dbname", dest="updateAbName", default = "default", type="string", help="name of abricate database to update")
46     parser.add_option("-u", "--update-mlst", dest="updateMLST", default = "False", type="string", help="True = update MLST")
47     #used for parsing 
48     parser.add_option("-s", "--mlst-scheme", dest="mlst", default= "/home/jjjjia/databases/scheme_species_map.tab", type="string", help="absolute file path to mlst scheme")
49 
50     
51     #parallelization, useless, these are hard coded to 8cores/64G RAM
52     #parser.add_option("-t", "--threads", dest="threads", default=8, type="int", help="number of cpu to use")
53     #parser.add_option("-p", "--memory", dest="memory", default=64, type="int", help="memory to use in GB")
54 
55     (options,args) = parser.parse_args()
56     #if len(args) != 8:
57         #parser.error("incorrect number of arguments, all 7 is required")
58     curDir = os.getcwd()
59     outputDir = options.output
60     expectedSpecies = options.expectedSpecies
61     mlstScheme = options.mlst
62     tempDir = outputDir + "/shovillTemp"
63     scriptDir = options.scriptDir
64     updateAbName = options.updateAbName
65     updateAbPath = options.updateAbPath
66     updateMLST = options.updateMLST
67     cardDB=options.cardDB
68     assemblyPath=options.assemblyPath
69     ID = options.id
70 else:
71     manifestDir = ""
72     curDir = os.getcwd()
73     outputDir = "pipelineTest"
74     expectedSpecies = "Escherichia coli"
75     #threads = 8
76     #memory = 64
77     mlstScheme = outputDir + "/scheme_species_map.tab"
78     tempDir= outputDir + "/shovillTemp"    
79     scriptDir = "~/scripts"
80     updateAbName = "cpo"
81     updateAbPath = "~/database/bccdcCPO.seq"
82     updateMLST = True
83     assemblyPath = "./"
84     cardDB = "./"
85     ID = "BC11-Kpn005_S2"
86 
87 #region result objects
88 #define some objects to store values from results
89 #//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).
90 class starFinders(object):
91     def __init__(self):
92         self.file = ""
93         self.sequence = ""
94         self.start = 0
95         self.end = 0
96         self.gene = ""
97         self.shortGene = ""
98         self.coverage = ""
99         self.coverage_map = ""
00         self.gaps = ""
01         self.pCoverage = 100.00
02         self.pIdentity = 100.00
03         self.database = ""
04         self.accession = ""
05         self.product = ""
06         self.source = "chromosome"
07         self.row = ""
08 
09 class PlasFlowResult(object):
10     def __init__(self):
11         self.sequence = ""
12         self.length = 0
13         self.label = ""
14         self.confidence = 0
15         self.usefulRow = ""
16         self.row = ""
17 
18 class MlstResult(object):
19     def __init__(self):
20         self.file = ""
21         self.speciesID = ""
22         self.seqType = 0
23         self.scheme = ""
24         self.species = ""
25         self.row=""
26 
27 class mobsuiteResult(object):
28     def __init__(self):
29         self.file_id = ""
30         self.cluster_id	= ""
31         self.contig_id	= ""
32         self.contig_num = 0
33         self.contig_length	= 0
34         self.circularity_status	= ""
35         self.rep_type	= ""
36         self.rep_type_accession = ""	
37         self.relaxase_type	= ""
38         self.relaxase_type_accession = ""	
39         self.mash_nearest_neighbor	 = ""
40         self.mash_neighbor_distance	= 0.00
41         self.repetitive_dna_id	= ""
42         self.match_type	= ""
43         self.score	= 0
44         self.contig_match_start	= 0
45         self.contig_match_end = 0
46         self.row = ""
47 
48 class mobsuitePlasmids(object):
49     def __init__(self):
50         self.file_id = ""
51         self.num_contigs = 0
52         self.total_length = 0
53         self.gc = ""
54         self.rep_types = ""
55         self.rep_typeAccession = ""
56         self.relaxase_type= ""
57         self.relaxase_type_accession	= ""
58         self.mpf_type	= ""
59         self.mpf_type_accession= ""	
60         self.orit_type	= ""
61         self.orit_accession	= ""
62         self.PredictedMobility	= ""
63         self.mash_nearest_neighbor	= ""
64         self.mash_neighbor_distance	= 0.00
65         self.mash_neighbor_cluster= 0
66         self.row = ""
67 class RGIResult(object):
68     def __init__(self):
69         self.ORF_ID	= ""
70         self.Contig	= ""
71         self.Start	= -1
72         self.Stop	= -1
73         self.Orientation = ""	
74         self.Cut_Off	= ""
75         self.Pass_Bitscore	= 100000
76         self.Best_Hit_Bitscore	= 0.00
77         self.Best_Hit_ARO	= ""
78         self.Best_Identities	= 0.00
79         self.ARO = 0
80         self.Model_type	= ""
81         self.SNPs_in_Best_Hit_ARO	= ""
82         self.Other_SNPs	= ""
83         self.Drug_Class	= ""
84         self.Resistance_Mechanism	= ""
85         self.AMR_Gene_Family	= ""
86         self.Predicted_DNA	= ""
87         self.Predicted_Protein	= ""
88         self.CARD_Protein_Sequence	= ""
89         self.Percentage_Length_of_Reference_Sequence	= 0.00
90         self.ID	= ""
91         self.Model_ID = 0
92         self.source = ""
93         self.row = ""
94 
95 #endregion
96 
97 #region useful functions
98 def read(path):
99     return [line.rstrip('\n') for line in open(path)]
00 def execute(command):
01     process = subprocess.Popen(command, shell=False, cwd=curDir, universal_newlines=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
02 
03     # Poll process for new output until finished
04     while True:
05         nextline = process.stdout.readline()
06         if nextline == '' and process.poll() is not None:
07             break
08         sys.stdout.write(nextline)
09         sys.stdout.flush()
10 
11     output = process.communicate()[0]
12     exitCode = process.returncode
13 
14     if (exitCode == 0):
15         return output
16     else:
17         raise subprocess.CalledProcessError(exitCode, command)
18 def httpGetFile(url, filepath=""):
19     if (filepath == ""):
20         return urllib.request.urlretrieve(url)
21     else:
22         urllib.request.urlretrieve(url, filepath)
23         return True
24 def gunzip(inputpath="", outputpath=""):
25     if (outputpath == ""):
26         with gzip.open(inputpath, 'rb') as f:
27             gzContent = f.read()
28         return gzContent
29     else:
30         with gzip.open(inputpath, 'rb') as f:
31             gzContent = f.read()
32         with open(outputpath, 'wb') as out:
33             out.write(gzContent)
34         return True
35 def ToJson(dictObject, outputPath):
36     outDir = outputDir + '/summary/' + ID + ".json/"
37     if not (os.path.exists(outDir)):
38         os.makedirs(outDir)
39     with open(outDir + outputPath, 'w') as f:
40       json.dump([ob.__dict__ for ob in dictObject.values()], f, ensure_ascii=False)
41 #endregion
42 
43 #region functions to parse result files
44 def ParseMLSTResult(pathToMLSTResult):
45     _mlstResult = {}
46     scheme = pandas.read_csv(mlstScheme, delimiter='\t', header=0)
47     scheme = scheme.replace(numpy.nan, '', regex=True)
48 
49     taxon = {}
50     #record the scheme as a dictionary
51     taxon["-"] = "No MLST Match"
52     for i in range(len(scheme.index)):
53         key = scheme.iloc[i,0]
54         if (str(scheme.iloc[i,2]) == "nan"):
55             value = str(scheme.iloc[i,1])
56         else:
57             value = str(scheme.iloc[i,1]) + " " + str(scheme.iloc[i,2])
58         
59         if (key in taxon.keys()):
60             taxon[key] = taxon.get(key) + ";" + value
61         else:
62             taxon[key] = value
63     #read in the mlst result
64     mlst = pandas.read_csv(pathToMLSTResult, delimiter='\t', header=None)
65     _mlstHit = MlstResult()
66 
67     _mlstHit.file = mlst.iloc[0,0]
68     _mlstHit.speciesID = (mlst.iloc[0,1])
69     _mlstHit.seqType = str(mlst.iloc[0,2])
70     for i in range(3, len(mlst.columns)):
71         _mlstHit.scheme += mlst.iloc[0,i] + ";"
72     _mlstHit.species = taxon[_mlstHit.speciesID]
73     _mlstHit.row = "\t".join(str(x) for x in mlst.ix[0].tolist())
74     _mlstResult[_mlstHit.speciesID]=_mlstHit
75 
76     return _mlstResult
77 
78 def ParsePlasmidFinderResult(pathToPlasmidFinderResult):
79     #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
80     #example resfinder:
81     #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 
82 
83     _pFinder = {} #***********************
84     plasmidFinder = pandas.read_csv(pathToPlasmidFinderResult, delimiter='\t', header=0)
85     plasmidFinder = plasmidFinder.replace(numpy.nan, '', regex=True)
86 
87 
88     for i in range(len(plasmidFinder.index)):
89         pf = starFinders()
90         pf.file = str(plasmidFinder.iloc[i,0])
91         pf.sequence = str(plasmidFinder.iloc[i,1])
92         pf.start = int(plasmidFinder.iloc[i,2])
93         pf.end = int(plasmidFinder.iloc[i,3])
94         pf.gene = str(plasmidFinder.iloc[i,4])
95         pf.shortGene = pf.gene[:pf.gene.index("_")]
96         pf.coverage = str(plasmidFinder.iloc[i,5])
97         pf.coverage_map = str(plasmidFinder.iloc[i,6])
98         pf.gaps = str(plasmidFinder.iloc[i,7])
99         pf.pCoverage = float(plasmidFinder.iloc[i,8])
00         pf.pIdentity = float(plasmidFinder.iloc[i,9])
01         pf.database = str(plasmidFinder.iloc[i,10])
02         pf.accession = str(plasmidFinder.iloc[i,11])
03         pf.product = str(plasmidFinder.iloc[i,12])
04         pf.source = "plasmid"
05         pf.row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
06         _pFinder[pf.gene]=pf
07         #row = "\t".join(str(x) for x in plasmidFinder.ix[i].tolist())
08         #plasmidFinderContigs.append(str(plasmidFinder.iloc[i,1]))
09         #origins.append(str(plasmidFinder.iloc[i,4][:plasmidFinder.iloc[i,4].index("_")]))
10     return _pFinder
11 
12 def ParseMobsuiteResult(pathToMobsuiteResult):
13     _mobsuite = {}
14     mResult = pandas.read_csv(pathToMobsuiteResult, delimiter='\t', header=0)
15     mResult = mResult.replace(numpy.nan, '', regex=True)
16 
17     for i in range(len(mResult.index)):
18         mr = mobsuiteResult()
19         mr.file_id = str(mResult.iloc[i,0])
20         mr.cluster_id = str(mResult.iloc[i,1])
21         if (mr.cluster_id == "chromosome"):
22             break
23         mr.contig_id = str(mResult.iloc[i,2])
24         mr.contig_num = mr.contig_id[(mr.contig_id.find("contig")+6):mr.contig_id.find("_len=")]
25         mr.contig_length = int(mResult.iloc[i,3])
26         mr.circularity_status = str(mResult.iloc[i,4])
27         mr.rep_type = str(mResult.iloc[i,5])
28         mr.rep_type_accession = str(mResult.iloc[i,6])
29         mr.relaxase_type = str(mResult.iloc[i,7])
30         mr.relaxase_type_accession = str(mResult.iloc[i,8])
31         mr.mash_nearest_neighbor = str(mResult.iloc[i,9])
32         mr.mash_neighbor_distance = float(mResult.iloc[i,10])
33         mr.repetitive_dna_id = str(mResult.iloc[i,11])
34         mr.match_type = str(mResult.iloc[i,12])
35         if (mr.match_type == ""):
36             mr.score = -1
37             mr.contig_match_start = -1
38             mr.contig_match_end = -1
39         else:
40             mr.score = int(mResult.iloc[i,13])
41             mr.contig_match_start = int(mResult.iloc[i,14])
42             mr.contig_match_end = int(mResult.iloc[i,15])
43         mr.row = "\t".join(str(x) for x in mResult.ix[i].tolist())
44         _mobsuite[mr.contig_id]=(mr)
45     return _mobsuite
46 
47 def ParseMobsuitePlasmids(pathToMobsuiteResult):
48     _mobsuite = {}
49     mResults = pandas.read_csv(pathToMobsuiteResult, delimiter='\t', header=0)
50     mResults = mResults.replace(numpy.nan, '', regex=True)
51 
52     for i in range(len(mResults.index)):
53         mr = mobsuitePlasmids()
54         mr.file_id = str(mResults.iloc[i,0])
55         mr.num_contigs = int(mResults.iloc[i,1])
56         mr.total_length = int(mResults.iloc[i,2])
57         mr.gc = int(mResults.iloc[i,3])
58         mr.rep_types = str(mResults.iloc[i,4])
59         mr.rep_typeAccession = str(mResults.iloc[i,5])
60         mr.relaxase_type = str(mResults.iloc[i,6])
61         mr.relaxase_type_accession = str(mResults.iloc[i,7])
62         mr.mpf_type = str(mResults.iloc[i,8])
63         mr.mpf_type_accession = str(mResults.iloc[i,9])
64         mr.orit_type = str(mResults.iloc[i,10])
65         mr.orit_accession = str(mResults.iloc[i,11])
66         mr.PredictedMobility = str(mResults.iloc[i,12])
67         mr.mash_nearest_neighbor = str(mResults.iloc[i,13])
68         mr.mash_neighbor_distance = float(mResults.iloc[i,14])
69         mr.mash_neighbor_cluster = int(mResults.iloc[i,15])
70         mr.row = "\t".join(str(x) for x in mResults.ix[i].tolist())
71         _mobsuite[mr.file_id] = mr
72     return _mobsuite
73 
74 def ParseResFinderResult(pathToResFinderResults, plasmidContigs, likelyPlasmidContigs):
75     _rFinder = {}
76     resFinder = pandas.read_csv(pathToResFinderResults, delimiter='\t', header=0)
77     resFinder = resFinder.replace(numpy.nan, '', regex=True)
78 
79     for i in range(len(resFinder.index)):
80         rf = starFinders()
81         rf.file = str(resFinder.iloc[i,0])
82         rf.sequence = str(resFinder.iloc[i,1])
83         rf.start = int(resFinder.iloc[i,2])
84         rf.end = int(resFinder.iloc[i,3])
85         rf.gene = str(resFinder.iloc[i,4])
86         rf.shortGene = rf.gene
87         rf.coverage = str(resFinder.iloc[i,5])
88         rf.coverage_map = str(resFinder.iloc[i,6])
89         rf.gaps = str(resFinder.iloc[i,7])
90         rf.pCoverage = float(resFinder.iloc[i,8])
91         rf.pIdentity = float(resFinder.iloc[i,9])
92         rf.database = str(resFinder.iloc[i,10])
93         rf.accession = str(resFinder.iloc[i,11])
94         rf.product = str(resFinder.iloc[i,12])
95         rf.row = "\t".join(str(x) for x in resFinder.ix[i].tolist())
96         if (rf.sequence[6:] in plasmidContigs):
97             rf.source = "plasmid"
98         elif (rf.sequence[6:] in likelyPlasmidContigs):
99             rf.source = "likely plasmid"
00         else:
01             rf.source = "likely chromosome"
02         _rFinder[rf.gene]=rf
03     return _rFinder
04 
05 def ParseRGIResult(pathToRGIResults, plasmidContigs, likelyPlasmidContigs):
06     _rgiR = {}
07     RGI = pandas.read_csv(pathToRGIResults, delimiter='\t', header=0)
08     RGI = RGI.replace(numpy.nan, '', regex=True)
09 
10     for i in range(len(RGI.index)):
11         r = RGIResult()
12         r.ORF_ID = str(RGI.iloc[i,0])
13         r.Contig = str(RGI.iloc[i,1])
14         r.Contig_Num = r.Contig[6:r.Contig.find("_")]
15         r.Start = int(RGI.iloc[i,2])
16         r.Stop = int(RGI.iloc[i,3])
17         r.Orientation = str(RGI.iloc[i,4])
18         r.Cut_Off = str(RGI.iloc[i,5])
19         r.Pass_Bitscore = int(RGI.iloc[i,6])
20         r.Best_Hit_Bitscore = float(RGI.iloc[i,7])
21         r.Best_Hit_ARO = str(RGI.iloc[i,8])
22         r.Best_Identities = float(RGI.iloc[i,9])
23         r.ARO = int(RGI.iloc[i,10])
24         r.Model_type = str(RGI.iloc[i,11])
25         r.SNPs_in_Best_Hit_ARO = str(RGI.iloc[i,12])
26         r.Other_SNPs = str(RGI.iloc[i,13])
27         r.Drug_Class = str(RGI.iloc[i,14])
28         r.Resistance_Mechanism = str(RGI.iloc[i,15])
29         r.AMR_Gene_Family = str(RGI.iloc[i,16])
30         r.Predicted_DNA = str(RGI.iloc[i,17])
31         r.Predicted_Protein = str(RGI.iloc[i,18])
32         r.CARD_Protein_Sequence = str(RGI.iloc[i,19])
33         r.Percentage_Length_of_Reference_Sequence = float(RGI.iloc[i,20])
34         r.ID = str(RGI.iloc[i,21])
35         r.Model_ID = int(RGI.iloc[i,22])
36         r.row = "\t".join(str(x) for x in RGI.ix[i].tolist())
37         if (r.Contig_Num in plasmidContigs):
38             r.source = "plasmid"
39         elif (r.Contig_Num in likelyPlasmidContigs):
40             r.source = "likely plasmid"
41         else:
42             r.source = "likely chromosome"
43         _rgiR[r.Model_ID]=r
44     return _rgiR
45 #endregion
46 
47 def Main():
48     notes = []
49     #init the output list
50     output = []
51     jsonOutput = []
52 
53     print(str(datetime.datetime.now()) + "\n\nID: " + ID + "\nAssembly: " + assemblyPath)
54     output.append(str(datetime.datetime.now()) + "\n\nID: " + ID + "\nAssembly: " + assemblyPath)
55 
56     #region update databases if update=true
57     if not debug:
58         #update databases if necessary
59         if not (updateAbPath == "" and updateAbName == "default"):
60             print("updating abricate database: " + updateAbName + " @fasta path: " + updateAbPath)
61             cmd = [scriptDir + "/pipeline_updateAbricateDB.sh", updateAbPath, updateAbName]
62             update = execute(cmd)
63         if (updateMLST.lower() == "true"):
64             print("updating mlst database... ")
65             cmd = [scriptDir + "/pipeline_updateMLST.sh"]
66             update = execute(cmd)
67     #endregion
68 
69     print("step 3: parsing the mlst results")
70 
71     print("performing MLST")
72     #region call the mlst script
73     if not debug:
74         print("running pipeline_prediction.sh")
75         #input parameters: 1=ID 2 = assemblyPath, 3= outputdir, 4=card.json
76         cmd = [scriptDir + "/pipeline_prediction.sh", ID, assemblyPath, outputDir, cardDB]
77         result = execute(cmd)
78     #endregion
79 
80     #region parse the mlst results
81     print("step 3: parsing mlst, plasmid, and amr results")
82     
83     print("identifying MLST")    
84     pathToMLSTScheme = outputDir + "/predictions/" + ID + ".mlst"
85     mlstHit = ParseMLSTResult(pathToMLSTScheme)#***********************
86     ToJson(mlstHit, "mlst.json") #write it to a json output
87     mlstHit = list(mlstHit.values())[0]
88 
89     #endregion
90 
91     #region parse mobsuite, resfinder and rgi results
92     print("identifying plasmid contigs and amr genes")
93 
94     plasmidContigs = []
95     likelyPlasmidContigs = []
96     origins = []
97 
98     #parse mobsuite results
99     mSuite = ParseMobsuiteResult(outputDir + "/predictions/" + ID + ".recon/contig_report.txt")#*************
00     ToJson(mSuite, "mobsuite.json") #*************
01     mSuitePlasmids = ParseMobsuitePlasmids(outputDir + "/predictions/" + ID + ".recon/mobtyper_aggregate_report.txt")#*************
02     ToJson(mSuitePlasmids, "mobsuitePlasmids.json") #*************
03 
04     for key in mSuite:
05         if mSuite[key].contig_num not in plasmidContigs and mSuite[key].contig_num not in likelyPlasmidContigs:
06             if not (mSuite[key].rep_type == ''):
07                 plasmidContigs.append(mSuite[key].contig_num)
08             else:
09                 likelyPlasmidContigs.append(mSuite[key].contig_num)
10     for key in mSuite:
11         if mSuite[key].rep_type not in origins:
12             origins.append(mSuite[key].rep_type)
13 
14     #parse resfinder AMR results
15     rFinder = ParseResFinderResult(outputDir + "/predictions/" + ID + ".cp", plasmidContigs, likelyPlasmidContigs) #**********************
16     ToJson(rFinder, "resfinder.json") #*************
17 
18     rgiAMR = ParseRGIResult(outputDir + "/predictions/" + ID + ".rgi.txt", plasmidContigs, likelyPlasmidContigs)#***********************
19     ToJson(rgiAMR, "rgi.json") #*************
20 
21     carbapenamases = []
22     amrGenes = []
23     for keys in rFinder:
24         carbapenamases.append(rFinder[keys].shortGene + "(" + rFinder[keys].source + ")")
25     for keys in rgiAMR:
26         if (rgiAMR[keys].Drug_Class.find("carbapenem") > -1):
27             if (rgiAMR[keys].Best_Hit_ARO not in carbapenamases):
28                 carbapenamases.append(rgiAMR[keys].Best_Hit_ARO+ "(" + rgiAMR[keys].source + ")")
29         else:
30             if (rgiAMR[keys].Best_Hit_ARO not in amrGenes):
31                 amrGenes.append(rgiAMR[keys].Best_Hit_ARO+ "(" + rgiAMR[keys].source + ")")
32     #endregion
33 
34     #region output parsed mlst information
35     print("formatting mlst outputs")
36     output.append("\n\n\n~~~~~~~MLST summary~~~~~~~")
37     output.append("MLST determined species: " + mlstHit.species)
38     output.append("\nMLST Details: ")
39     output.append(mlstHit.row)
40 
41     output.append("\nMLST information: ")
42     if (mlstHit.species == expectedSpecies):
43         output.append("MLST determined species is the same as expected species")
44         #notes.append("MLST determined species is the same as expected species")
45     else:
46         output.append("!!!MLST determined species is NOT the same as expected species, contamination? mislabeling?")
47         notes.append("MLST: Not expected species. Possible contamination or mislabeling")
48 
49     #endregion
50 
51     #region output the parsed plasmid/amr results
52     output.append("\n\n\n~~~~~~~~Plasmids~~~~~~~~\n")
53     
54     output.append("predicted plasmid origins: ")
55     output.append(";".join(origins))
56 
57     output.append("\ndefinitely plasmid contigs")
58     output.append(";".join(plasmidContigs))
59     
60     output.append("\nlikely plasmid contigs")
61     output.append(";".join(likelyPlasmidContigs))
62 
63     output.append("\nmob-suite prediction details: ")
64     for key in mSuite:
65         output.append(mSuite[key].row)
66 
67     output.append("\n\n\n~~~~~~~~AMR Genes~~~~~~~~\n")
68     output.append("predicted carbapenamase Genes: ")
69     output.append(",".join(carbapenamases))
70     output.append("other RGI AMR Genes: ")
71     for key in rgiAMR:
72         output.append(rgiAMR[key].Best_Hit_ARO + "(" + rgiAMR[key].source + ")")
73 
74     output.append("\nDetails about the carbapenamase Genes: ")
75     for key in rFinder:
76         output.append(rFinder[key].row)
77     output.append("\nDetails about the RGI AMR Genes: ")
78     for key in rgiAMR:
79         output.append(rgiAMR[key].row)
80 
81     #write summary to a file
82     summaryDir = outputDir + "/summary/" + ID
83     out = open(summaryDir + ".txt", 'w')
84     for item in output:
85         out.write("%s\n" % item)
86 
87 
88     #TSV output
89     tsvOut = []
90     tsvOut.append("ID\tExpected Species\tMLST Species\tSequence Type\tMLST Scheme\tCarbapenem Resistance Genes\tOther AMR Genes\tTotal Plasmids\tPlasmids ID\tNum_Contigs\tPlasmid Length\tPlasmid RepType\tPlasmid Mobility\tNearest Reference\tDefinitely Plasmid Contigs\tLikely Plasmid Contigs")
91     #start with ID
92     temp = ""
93     temp += (ID + "\t")
94     temp += expectedSpecies + "\t"
95 
96     #move into MLST
97     temp += mlstHit.species + "\t"
98     temp += str(mlstHit.seqType) + "\t"
99     temp += mlstHit.scheme + "\t"
00     
01     #now onto AMR genes
02     temp += ";".join(carbapenamases) + "\t"
03     temp += ";".join(amrGenes) + "\t"
04 
05     #lastly plasmids
06     temp+= str(len(mSuitePlasmids)) + "\t"
07     plasmidID = ""
08     contigs = ""
09     lengths = ""
10     rep_type = ""
11     mobility = ""
12     neighbour = ""
13     for keys in mSuitePlasmids:
14         plasmidID += str(mSuitePlasmids[keys].mash_neighbor_cluster) + ";"
15         contigs += str(mSuitePlasmids[keys].num_contigs) + ";"
16         lengths += str(mSuitePlasmids[keys].total_length) + ";"
17         rep_type += str(mSuitePlasmids[keys].rep_types) + ";"
18         mobility += str(mSuitePlasmids[keys].PredictedMobility) + ";"
19         neighbour += str(mSuitePlasmids[keys].mash_nearest_neighbor) + ";"
20     temp += plasmidID + "\t" + contigs + "\t" + lengths + "\t" + rep_type + "\t" + mobility + "\t" + neighbour + "\t"
21     temp += ";".join(plasmidContigs) + "\t"
22     temp += ";".join(likelyPlasmidContigs)
23     tsvOut.append(temp)
24 
25     summaryDir = outputDir + "/summary/" + ID
26     out = open(summaryDir + ".tsv", 'w')
27     for item in tsvOut:
28         out.write("%s\n" % item)
29     #endregion
30 
31 
32 start = time.time()#time the analysis
33 print("Starting workflow...")
34 #analysis time
35 Main()
36 
37 end = time.time()
38 print("Finished!\nThe analysis used: " + str(end-start) + " seconds")