Allow bbmap to retain unmapped reads

Author: bbimber

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/alignment/BBMapWrapper.java

@@ -20,6 +20,7 @@
 import org.labkey.api.sequenceanalysis.pipeline.ToolParameterDescriptor;
 import org.labkey.api.sequenceanalysis.run.AbstractAlignmentPipelineStep;
 import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
+import org.labkey.api.util.FileUtil;
 
 import java.io.File;
 import java.io.IOException;
@@ -136,7 +137,8 @@ public AlignmentOutput performAlignment(Readset rs, List<File> inputFastqs1, @Nu
                 params.add("minid=" + val);
             }
 
-            File bam = getWrapper().doAlignment(inputFastq1, inputFastq2, outputDirectory, basename, params);
+            boolean retainUnmapped = getProvider().getParameterByName("retainUnmapped").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class, false);
+            File bam = getWrapper().doAlignment(inputFastq1, inputFastq2, outputDirectory, basename, params, retainUnmapped);
             if (!bam.exists())
             {
                 throw new PipelineJobException("BAM not created, expected: " + bam.getPath());
@@ -145,6 +147,17 @@ public AlignmentOutput performAlignment(Readset rs, List<File> inputFastqs1, @Nu
             output.setBAM(bam);
             output.addCommandsExecuted(getWrapper().getCommandsExecuted());
 
+            if (retainUnmapped)
+            {
+                File unmappedBam = getWrapper().getUnmappedFilename(bam);
+                if (!unmappedBam.exists())
+                {
+                    throw new PipelineJobException("Unable to find file: " + unmappedBam.getPath());
+                }
+
+                output.addSequenceOutput(unmappedBam, rs.getName() + ": BBmap unmapped reads", "Alignment", rs.getReadsetId(), null, referenceGenome.getGenomeId(), "BBMap Unmapped Reads");
+            }
+
             return output;
         }
 
@@ -184,14 +197,18 @@ public Provider()
                     }}, true),
                     ToolParameterDescriptor.create("semiperfectmode", "Semi-perfectmode", "Allow only perfect and semiperfect (perfect except for N's in the reference) mappings", "checkbox", new JSONObject()
                     {{
-                        put("checked", true);
-                    }}, true),
+                        put("checked", false);
+                    }}, false),
                     ToolParameterDescriptor.create("minid", "Minimum Identity", "Approximate minimum alignment identity to look for. Higher is faster and less sensitive", "ldk-numberfield", new JSONObject()
                     {{
                         put("minValue", 0);
                         put("maxValue", 1);
                         put("decimalPrecision", 2);
-                    }}, 0.95)
+                    }}, 0.95),
+                    ToolParameterDescriptor.create("retainUnmapped", "Retain Unmapped", "If checked, unmapped reads are written to a separate BAM file", "checkbox", new JSONObject()
+                    {{
+                        put("checked", false);
+                    }}, false)
                ), null, "https://prost.readthedocs.io/en/latest/bbmap.html", true, true);
         }
 
@@ -202,12 +219,18 @@ public BBMapAlignmentStep create(PipelineContext context)
         }
     }
 
+    protected File getUnmappedFilename(File mappedSamOrBam)
+    {
+        String fn = FileUtil.getBaseName(mappedSamOrBam.getName()) + ".unmapped." + FileUtil.getExtension(mappedSamOrBam.getName());
+        return new File(mappedSamOrBam.getParentFile(), fn);
+    }
+
     protected File getExe()
     {
         return SequencePipelineService.get().getExeForPackage("BBMAPPATH", "bbmap.sh");
     }
 
-    public File doAlignment(File inputFastq1, @Nullable File inputFastq2, File outputDirectory, String basename, List<String> options) throws PipelineJobException
+    public File doAlignment(File inputFastq1, @Nullable File inputFastq2, File outputDirectory, String basename, List<String> options, boolean retainUnmaped) throws PipelineJobException
     {
         List<String> args = new ArrayList<>();
         args.add(getExe().getPath());
@@ -243,6 +266,18 @@ public File doAlignment(File inputFastq1, @Nullable File inputFastq2, File outpu
 
         args.add("outm=" + outputSam.getPath());
 
+        File outputUnmappedSam = null;
+        if (retainUnmaped)
+        {
+            outputUnmappedSam = getUnmappedFilename(outputSam);
+            if (outputUnmappedSam.exists())
+            {
+                outputUnmappedSam.delete();
+            }
+
+            args.add("outu=" + outputUnmappedSam.getPath());
+        }
+
         Integer maxRam = SequencePipelineService.get().getMaxRam();
         if (maxRam != null)
         {
@@ -284,6 +319,32 @@ public File doAlignment(File inputFastq1, @Nullable File inputFastq2, File outpu
 
         outputSam.delete();
 
+        // repeat for unmapped:
+        if (outputUnmappedSam != null)
+        {
+            File outputUnmappedBam = getUnmappedFilename(outputBam);
+            if (outputUnmappedBam != null && outputUnmappedBam.exists())
+            {
+                outputUnmappedBam.delete();
+            }
+
+            samtoolsRunner = new SamtoolsRunner(getLogger());
+            stArgs = new ArrayList<>();
+            stArgs.add(samtoolsRunner.getSamtoolsPath().getPath());
+            stArgs.add("view");
+            stArgs.add("-o");
+            stArgs.add(outputUnmappedBam.getPath());
+            stArgs.add(outputUnmappedSam.getPath());
+            samtoolsRunner.execute(stArgs);
+
+            if (!outputUnmappedBam.exists())
+            {
+                throw new PipelineJobException("File not found: " + outputBam.getPath());
+            }
+
+            outputUnmappedSam.delete();
+        }
+
         return outputBam;
     }