Merge pull request #6 from Zymo-Research/bacteriaOnly

Appears to be working as desired

README.md

@@ -57,6 +57,7 @@ WORKINGFOLDER	|	string	|	/data/working	|	Path to working folder for temporary fi
 OUTPUTFOLDER	|	string	|	/data/output	|	Folder within the container to write output data
 REFERENCEGENOME	|	string	|	[folderWithPackage]/reference/zrCommunityStandard.fa	|	File containing the reference sequence for the standard within the container
 FILENAMINGSTANDARD	|	string	|	ZYMO	|	How sequence files will be named (other option is "illumina")
+BACTERIAONLY | BOOL | FALSE | Analyze and calculate MIQ score on bacteria only.
 MODE | string | PE | Running mode.  PE for paired-end, SE for single end, LONG for Nanopore reads.  PE and SE use similar alignment but different analysis logic. SE and LONG use the same analysis logic, but different aligners.
 
 

analyzeStandardReads.py

@@ -27,8 +27,16 @@ def getApplicationParametersSE():
     parameters.addParameter("referenceGenome", str, default=default.referenceGenome, expectedFile=True)
     parameters.addParameter("fileNamingStandard", str, default="zymo", externalValidation=True)
     parameters.addParameter("referenceDataFile", str, default = default.referenceDataFile, expectedFile=True)
-    parameters.addParameter("goodMiqExample", str, default = default.goodMiqExample, expectedFile=True)
-    parameters.addParameter("badMiqExample", str, default=default.badMiqExample, expectedFile=True)
+    parameters.addParameter("bacteriaOnly", bool, required=False, default=False)
+    if parameters.bacteriaOnly.value:
+        print("ANALYZING BACTERIA ONLY")
+        defaultBadExample = default.badMiqExampleBacteriaOnly
+        defaultGoodExample = default.goodMiqExampleBacteriaOnly
+    else:
+        defaultBadExample = default.badMiqExample
+        defaultGoodExample = default.goodMiqExample
+    parameters.addParameter("goodMiqExample", str, default=defaultGoodExample, expectedFile=True)
+    parameters.addParameter("badMiqExample", str, default=defaultBadExample, expectedFile=True)
     if not validSampleName(parameters.sampleName.value):
         logger.error("Invalid sample name given: %s" %parameters.sampleName.value)
         raise ValueError("Invalid sample name given: %s" %parameters.sampleName.value)
@@ -48,8 +56,16 @@ def getApplicationParametersPE():
     parameters.addParameter("referenceGenome", str, default=default.referenceGenome, expectedFile=True)
     parameters.addParameter("fileNamingStandard", str, default="zymo", externalValidation=True)
     parameters.addParameter("referenceDataFile", str, default=default.referenceDataFile, expectedFile=True)
-    parameters.addParameter("goodMiqExample", str, default = default.goodMiqExample, expectedFile=True)
-    parameters.addParameter("badMiqExample", str, default=default.badMiqExample, expectedFile=True)
+    parameters.addParameter("bacteriaOnly", bool, required=False, default=False)
+    if parameters.bacteriaOnly.value:
+        print("ANALYZING BACTERIA ONLY")
+        defaultBadExample = default.badMiqExampleBacteriaOnly
+        defaultGoodExample = default.goodMiqExampleBacteriaOnly
+    else:
+        defaultBadExample = default.badMiqExample
+        defaultGoodExample = default.goodMiqExample
+    parameters.addParameter("goodMiqExample", str, default=defaultGoodExample, expectedFile=True)
+    parameters.addParameter("badMiqExample", str, default=defaultBadExample, expectedFile=True)
     if not validSampleName(parameters.sampleName.value):
         logger.error("Invalid sample name given: %s" %parameters.sampleName.value)
         raise ValueError("Invalid sample name given: %s" %parameters.sampleName.value)
@@ -68,8 +84,16 @@ def getApplicationParametersLong():
     parameters.addParameter("referenceGenome", str, default=default.referenceGenome, expectedFile=True)
     parameters.addParameter("fileNamingStandard", str, default="zymo", externalValidation=True)
     parameters.addParameter("referenceDataFile", str, default=default.referenceDataFileHMW, expectedFile=True)
-    parameters.addParameter("goodMiqExample", str, default = default.goodMiqExampleHMW, expectedFile=True)
-    parameters.addParameter("badMiqExample", str, default=default.badMiqExampleHMW, expectedFile=True)
+    parameters.addParameter("bacteriaOnly", bool, required=False, default=False)
+    if parameters.bacteriaOnly.value:
+        print("ANALYZING BACTERIA ONLY")
+        defaultBadExample = default.badMiqExampleBacteriaOnlyHMW
+        defaultGoodExample = default.goodMiqExampleBacteriaOnlyHMW
+    else:
+        defaultBadExample = default.badMiqExampleHMW
+        defaultGoodExample = default.goodMiqExampleHMW
+    parameters.addParameter("goodMiqExample", str, default=defaultGoodExample, expectedFile=True)
+    parameters.addParameter("badMiqExample", str, default=defaultBadExample, expectedFile=True)
     if not validSampleName(parameters.sampleName.value):
         logger.error("Invalid sample name given: %s" %parameters.sampleName.value)
         raise ValueError("Invalid sample name given: %s" %parameters.sampleName.value)
@@ -178,8 +202,12 @@ def getReadLengthsFromFastq(path:str):
 
 
 def analyzeStandardResult(resultTable:dict):
+    if not parameters.bacteriaOnly.value:
+        analysisMethod = "Genomic"
+    else:
+        analysisMethod = "GenomicBacteriaOnly"
     referenceData = miqScoreNGSReadCountPublic.referenceHandler.StandardReference(parameters.referenceDataFile.value)
-    calculator = miqScoreNGSReadCountPublic.MiqScoreCalculator(referenceData, analysisMethod="Genomic", floor=0)
+    calculator = miqScoreNGSReadCountPublic.MiqScoreCalculator(referenceData, analysisMethod=analysisMethod, floor=0)
     miqScoreResult = calculator.calculateMiq(resultTable, parameters.sampleName.value)
     miqScoreResult.makeReadFateChart(readFatePrintNames=readFatePrintNames)
     miqScoreResult.makeRadarPlots()
@@ -219,7 +247,11 @@ def generateReport(result:miqScoreNGSReadCountPublic.MiqScoreData):
     templateFile = open(templateFilePath, 'r')
     template = templateFile.read()
     templateFile.close()
-    goodMiq, badMiq = miqScoreNGSReadCountPublic.loadExampleData(parameters.goodMiqExample.value, parameters.badMiqExample.value, referenceData, "Genomic")
+    if not parameters.bacteriaOnly.value:
+        analysisMethod = "Genomic"
+    else:
+        analysisMethod = "GenomicBacteriaOnly"
+    goodMiq, badMiq = miqScoreNGSReadCountPublic.loadExampleData(parameters.goodMiqExample.value, parameters.badMiqExample.value, referenceData, analysisMethod)
     replacementTable = generateReportReplacementTable(result, goodMiq, badMiq, readFatePrintNames=readFatePrintNames)
     report = miqScoreNGSReadCountPublic.reportGeneration.generateReport(template, replacementTable)
     reportFilePath = os.path.join(parameters.outputFolder.value, "%s.html" % parameters.sampleName.value)
@@ -243,6 +275,8 @@ def generateReport(result:miqScoreNGSReadCountPublic.MiqScoreData):
         parameters = getApplicationParametersSE()
     elif applicationMode == "LONG":
         parameters = getApplicationParametersLong()
+    else:
+        raise RuntimeError("This should never be reached as the application mode checker should have caught an invalid mode. Please investigate how this happened. Application mode: %s" %applicationMode)
     logger.debug("Starting analysis")
     bamFilePath = os.path.join(parameters.outputFolder.value, "%s.bam" %parameters.sampleName.value)
     if applicationMode == "PE":
@@ -254,6 +288,8 @@ def generateReport(result:miqScoreNGSReadCountPublic.MiqScoreData):
     elif applicationMode == "LONG":
         miqScoreShotgunPublicSupport.alignmentAnalysis.minimap2.minimapAlign(parameters.reads.value, parameters.workingFolder.value, bamFilePath, parameters.referenceGenome.value)
         readTable = miqScoreShotgunPublicSupport.alignmentAnalysis.alignmentAnalysisSE.bamFileProcessor(bamFilePath)
+    else:
+        raise RuntimeError("This should never be reached as the application mode checker should have caught an invalid mode. Please investigate how this happened. Application mode: %s" %applicationMode)
     standardAnalysisResults = analyzeStandardResult(readTable)
     saveResult(standardAnalysisResults)
     generateReport(standardAnalysisResults)

defaults/environment.py

@@ -19,4 +19,8 @@
 badMiqExample = os.path.join(referenceFolder, "badMiq.json")
 goodMiqExampleHMW = os.path.join(referenceFolder, "goodMiqHMW.json")
 badMiqExampleHMW = os.path.join(referenceFolder, "badMiqHMW.json")
+goodMiqExampleBacteriaOnly = os.path.join(referenceFolder, "goodMiqBacteriaOnly.json")
+badMiqExampleBacteriaOnly = os.path.join(referenceFolder, "badMiqBacteriaOnly.json")
+goodMiqExampleBacteriaOnlyHMW = os.path.join(referenceFolder, "goodMiqBacteriaOnlyHMW.json")
+badMiqExampleBacteriaOnlyHMW = os.path.join(referenceFolder, "badMiqBacteriaOnlyHMW.json")
 logFile = os.path.join(outputFolder, "dada2.%s.log" %timestamp)

reference/zrCommunityStandard.json

@@ -142,6 +142,16 @@
             "scerevisiae": 2,
             "cneoformans": 2
         },
+        "GenomicBacteriaOnly": {
+            "paeruginosa": 12.5,
+            "ecoli": 12.5,
+            "senterica": 12.5,
+            "lfermentum": 12.5,
+            "efaecalis": 12.5,
+            "saureus": 12.5,
+            "lmonocytogenes": 12.5,
+            "bsubtilis": 12.5
+        },
         "16s": {
             "paeruginosa": 4.2,
             "ecoli": 10.1,

reference/zrCommunityStandardHMW.json

@@ -131,6 +131,15 @@
             "lmonocytogenes": 14,
             "bsubtilis": 14,
             "scerevisiae": 2
+        },
+         "GenomicBacteriaOnly": {
+            "paeruginosa": 14.29,
+            "ecoli": 14.29,
+            "senterica": 14.29,
+            "efaecalis": 14.29,
+            "saureus": 14.29,
+            "lmonocytogenes": 14.29,
+            "bsubtilis": 14.29
         }
     }
 }