Merge discvr-23.7 to develop

Author: bbimber

SequenceAnalysis/api-src/org/labkey/api/sequenceanalysis/pipeline/VariantProcessingStep.java

@@ -68,9 +68,9 @@ default void validateScatter(ScatterGatherMethod method, PipelineJob job) throws
 
         }
 
-        default void performAdditionalMergeTasks(SequenceOutputHandler.JobContext ctx, PipelineJob job, TaskFileManager manager, ReferenceGenome genome, List<File> orderedScatterOutputs) throws PipelineJobException
+        default void performAdditionalMergeTasks(SequenceOutputHandler.JobContext ctx, PipelineJob job, TaskFileManager manager, ReferenceGenome genome, List<File> orderedScatterOutputs, List<String> orderedJobDirs) throws PipelineJobException
         {
-
+            ctx.getLogger().debug("No additional merge tasks are implemented for: " + getClass().getName());
         }
 
         default boolean doSortAfterMerge()

SequenceAnalysis/api-src/org/labkey/api/sequenceanalysis/run/AbstractCommandWrapper.java

@@ -253,7 +253,7 @@ public Logger getLogger()
     {
         if (_log == null)
         {
-              return  LogManager.getLogger("NoOpLogger");
+              return LogManager.getLogger("NoOpLogger");
         }
 
         return _log;

SequenceAnalysis/src/org/labkey/sequenceanalysis/analysis/PicardAlignmentMetricsHandler.java

@@ -25,6 +25,7 @@
 import org.labkey.sequenceanalysis.run.util.CollectInsertSizeMetricsWrapper;
 import org.labkey.sequenceanalysis.run.util.CollectWgsMetricsWithNonZeroCoverageWrapper;
 import org.labkey.sequenceanalysis.run.util.CollectWgsMetricsWrapper;
+import org.labkey.sequenceanalysis.run.util.MarkDuplicatesWrapper;
 
 import java.io.File;
 import java.util.ArrayList;
@@ -54,7 +55,13 @@ public PicardAlignmentMetricsHandler()
                 }}, false),
                 ToolParameterDescriptor.create("collectWgsNonZero", "Run WGS Metrics Over Non-Zero Coverage", "If checked, Picard CollectWgsMetricsWithNonZeroCoverage will be run", "checkbox", new JSONObject(){{
                     put("checked", false);
-                }}, false)
+                }}, false),
+                ToolParameterDescriptor.create("markDuplicates", "Run MarkDuplicates", "If checked, Picard CollectWgsMetricsWithNonZeroCoverage will be run", "checkbox", new JSONObject(){{
+                    put("checked", false);
+                }}, false),
+                ToolParameterDescriptor.create("useOutputFileContainer", "Submit to Source File Workbook", "If checked, each job will be submitted to the same workbook as the input file, as opposed to submitting all jobs to the same workbook.  This is primarily useful if submitting a large batch of files to process separately. This only applies if 'Run Separately' is selected.", "checkbox", new JSONObject(){{
+                    put("checked", true);
+                }}, true)
         ));
     }
 
@@ -156,6 +163,13 @@ public void processFilesOnWebserver(PipelineJob job, SequenceAnalysisJobSupport
                         metricsFiles.add(mf4);
                     }
 
+                    File mf5 = new MarkDuplicatesWrapper(job.getLogger()).getMetricsFile(m.getAlignmentFileObject());
+                    if (mf5.exists())
+                    {
+                        action.addOutput(mf5, "Duplication Metrics", false);
+                        metricsFiles.add(mf5);
+                    }
+
                     TableInfo ti = SequenceAnalysisManager.get().getTable(SequenceAnalysisSchema.TABLE_QUALITY_METRICS);
                     for (File f : metricsFiles)
                     {
@@ -195,6 +209,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
             boolean collectInsertSize = params.optBoolean("collectInsertSize", false);
             boolean collectWgs = params.optBoolean("collectWgs", false);
             boolean collectWgsNonZero = params.optBoolean("collectWgsNonZero", false);
+            boolean runMarkDuplicates = params.optBoolean("markDuplicates", false);
 
             int i = 1;
             for (SequenceOutputFile o : inputFiles)
@@ -244,6 +259,28 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
                     CollectInsertSizeMetricsWrapper wrapper = new CollectInsertSizeMetricsWrapper(job.getLogger());
                     wrapper.executeCommand(o.getFile(), metricsFile, metricsHistogram);
                 }
+
+                if (runMarkDuplicates)
+                {
+                    job.getLogger().info("running MarkDuplicates");
+                    job.setStatus(PipelineJob.TaskStatus.running, "RUNNING MARKDUPLICATES");
+                    MarkDuplicatesWrapper wrapper = new MarkDuplicatesWrapper(job.getLogger());
+                    File metricsFile = wrapper.getMetricsFile(o.getFile());
+                    File tempBam = new File(ctx.getOutputDir(), FileUtil.getBaseName(o.getFile()) + ".markDuplicates.bam");
+                    ctx.getFileManager().addIntermediateFile(tempBam);
+                    ctx.getFileManager().addIntermediateFile(new File(tempBam.getPath() + ".bai"));
+
+                    if (tempBam.exists())
+                    {
+                        tempBam.delete();
+                    }
+
+                    wrapper.executeCommand(o.getFile(), tempBam, null);
+                    if (!metricsFile.exists())
+                    {
+                        throw new PipelineJobException("Unable to find file: " + metricsFile);
+                    }
+                }
             }
 
             action.setEndTime(new Date());

SequenceAnalysis/src/org/labkey/sequenceanalysis/pipeline/PicardMetricsUtil.java

@@ -4,12 +4,11 @@
 import htsjdk.samtools.metrics.MetricBase;
 import htsjdk.samtools.metrics.MetricsFile;
 import org.apache.commons.beanutils.ConversionException;
-import org.apache.commons.io.Charsets;
 import org.apache.commons.lang3.StringUtils;
 import org.apache.logging.log4j.Logger;
-import org.apache.logging.log4j.LogManager;
 import org.labkey.api.data.ConvertHelper;
 import org.labkey.api.pipeline.PipelineJobException;
+import org.labkey.api.reader.Readers;
 import picard.analysis.AlignmentSummaryMetrics;
 import picard.analysis.CollectWgsMetricsWithNonZeroCoverage;
 import picard.analysis.InsertSizeMetrics;
@@ -18,9 +17,7 @@
 
 import java.io.BufferedReader;
 import java.io.File;
-import java.io.FileInputStream;
 import java.io.IOException;
-import java.io.InputStreamReader;
 import java.util.ArrayList;
 import java.util.HashMap;
 import java.util.List;
@@ -38,9 +35,9 @@ public static List<Map<String, Object>> processFile(File f, Logger log) throws P
             throw new PipelineJobException("Unable to find file: " + f.getPath());
         }
 
-        try (BufferedReader reader = new BufferedReader(new InputStreamReader(new FileInputStream(f), Charsets.UTF_8)))
+        try (BufferedReader reader = Readers.getReader(f))
         {
-            MetricsFile metricsFile = new MetricsFile();
+            MetricsFile<MetricBase,?> metricsFile = new MetricsFile<>();
             metricsFile.read(reader);
             List<MetricBase> metrics = metricsFile.getMetrics();
             if (metrics.get(0).getClass() == DuplicationMetrics.class)

SequenceAnalysis/src/org/labkey/sequenceanalysis/pipeline/ProcessVariantsHandler.java

@@ -209,7 +209,8 @@ public static File getScatterOutputByCategory(JobContext ctx, final String categ
 
         if (scatterOutputs.isEmpty())
         {
-            throw new PipelineJobException("Unable to find final VCF");
+            ctx.getLogger().info("No outputs of category: " + category + " were found");
+            return null;
         }
         else if (scatterOutputs.size() > 1)
         {
@@ -411,6 +412,7 @@ public static File processVCF(File input, Integer libraryId, JobContext ctx, Res
             resumer.getFileManager().addIntermediateFile(outputFileIdx);
         }
 
+        File effectiveInput = currentVCF; //this will be tested at the end to determine if a new file was actually created
         for (PipelineStepCtx<VariantProcessingStep> stepCtx : providers)
         {
             ctx.getLogger().info("Starting to run: " + stepCtx.getProvider().getLabel());
@@ -484,7 +486,7 @@ public static File processVCF(File input, Integer libraryId, JobContext ctx, Res
             resumer.setStepComplete(stepIdx, input.getPath(), action, currentVCF);
         }
 
-        if (currentVCF.exists())
+        if (currentVCF != null && currentVCF.exists() && !currentVCF.equals(effectiveInput))
         {
             resumer.getFileManager().removeIntermediateFile(currentVCF);
             resumer.getFileManager().removeIntermediateFile(new File(currentVCF.getPath() + ".tbi"));
@@ -872,14 +874,15 @@ else if (AbstractGenomicsDBImportHandler.TILE_DB_FILETYPE.isType(input))
     }
 
     @Override
-    public void performAdditionalMergeTasks(JobContext ctx, PipelineJob job, TaskFileManager manager, ReferenceGenome genome, List<File> orderedScatterOutputs) throws PipelineJobException
+    public void performAdditionalMergeTasks(JobContext ctx, PipelineJob job, TaskFileManager manager, ReferenceGenome genome, List<File> orderedScatterOutputs, List<String> orderedJobDirs) throws PipelineJobException
     {
         List<PipelineStepCtx<VariantProcessingStep>> providers = SequencePipelineService.get().getSteps(job, VariantProcessingStep.class);
         for (PipelineStepCtx<VariantProcessingStep> stepCtx : providers)
         {
-            if (stepCtx.getProvider() instanceof VariantProcessingStep.SupportsScatterGather ssg)
+            VariantProcessingStep vps = stepCtx.getProvider().create(ctx);
+            if (vps instanceof VariantProcessingStep.SupportsScatterGather ssg)
             {
-                ssg.performAdditionalMergeTasks(ctx, job, manager, genome, orderedScatterOutputs);
+                ssg.performAdditionalMergeTasks(ctx, job, manager, genome, orderedScatterOutputs, orderedJobDirs);
             }
         }
     }

SequenceAnalysis/src/org/labkey/sequenceanalysis/pipeline/VariantProcessingRemoteMergeTask.java

@@ -3,6 +3,7 @@
 import htsjdk.samtools.util.Interval;
 import org.apache.commons.lang3.StringUtils;
 import org.jetbrains.annotations.NotNull;
+import org.jetbrains.annotations.Nullable;
 import org.labkey.api.pipeline.AbstractTaskFactory;
 import org.labkey.api.pipeline.AbstractTaskFactorySettings;
 import org.labkey.api.pipeline.PipelineJob;
@@ -127,13 +128,16 @@ private VariantProcessingJob getPipelineJob()
         if (handler instanceof SequenceOutputHandler.TracksVCF)
         {
             Set<SequenceOutputFile> outputs = new HashSet<>();
-            scatterOutputs.values().forEach(f -> outputs.addAll(getPipelineJob().getOutputsToCreate().stream().filter(x -> f.equals(x.getFile())).collect(Collectors.toSet())));
+            scatterOutputs.values().forEach(f -> outputs.addAll(getPipelineJob().getOutputsToCreate().stream().filter(x -> x != null && f.equals(x.getFile())).collect(Collectors.toSet())));
             getJob().getLogger().debug("Total component outputs created: " + outputs.size());
             getPipelineJob().getOutputsToCreate().removeAll(outputs);
             getJob().getLogger().debug("Total SequenceOutputFiles on job after remove: " + getPipelineJob().getOutputsToCreate().size());
 
-            SequenceOutputFile finalOutput = ((SequenceOutputHandler.TracksVCF)getPipelineJob().getHandler()).createFinalSequenceOutput(getJob(), finalOut, getPipelineJob().getFiles());
-            manager.addSequenceOutput(finalOutput);
+            if (finalOut != null)
+            {
+                SequenceOutputFile finalOutput = ((SequenceOutputHandler.TracksVCF) getPipelineJob().getHandler()).createFinalSequenceOutput(getJob(), finalOut, getPipelineJob().getFiles());
+                manager.addSequenceOutput(finalOutput);
+            }
         }
         else
         {
@@ -152,14 +156,15 @@ private VariantProcessingJob getPipelineJob()
         return new RecordedActionSet(action);
     }
 
-    private File runDefaultVariantMerge(JobContextImpl ctx, TaskFileManagerImpl manager, RecordedAction action, SequenceOutputHandler<SequenceOutputHandler.SequenceOutputProcessor> handler) throws PipelineJobException
+    private @Nullable File runDefaultVariantMerge(JobContextImpl ctx, TaskFileManagerImpl manager, RecordedAction action, SequenceOutputHandler<SequenceOutputHandler.SequenceOutputProcessor> handler) throws PipelineJobException
     {
         Map<String, List<Interval>> jobToIntervalMap = getPipelineJob().getJobToIntervalMap();
         getJob().setStatus(PipelineJob.TaskStatus.running, "Combining Per-Contig VCFs: " + jobToIntervalMap.size());
 
         Map<String, File> scatterOutputs = getPipelineJob().getScatterJobOutputs();
         List<File> toConcat = new ArrayList<>();
         Set<File> missing = new HashSet<>();
+        int totalNull = 0;
         for (String name : jobToIntervalMap.keySet())
         {
             if (!scatterOutputs.containsKey(name))
@@ -168,45 +173,29 @@ private File runDefaultVariantMerge(JobContextImpl ctx, TaskFileManagerImpl mana
             }
 
             File vcf = scatterOutputs.get(name);
-            if (!vcf.exists())
+            if (scatterOutputs.get(name) == null)
             {
-                missing.add(vcf);
-            }
-
-            // NOTE: this was added to fix a one-time issue where -L was dropped from some upstream GenotypeGVCFs.
-            // Under normal conditions this would never be necessary.
-            boolean ensureOutputsWithinIntervals = getPipelineJob().getParameterJson().optBoolean("variantCalling.GenotypeGVCFs.ensureOutputsWithinIntervalsOnMerge", false);
-            if (ensureOutputsWithinIntervals)
-            {
-                getJob().getLogger().debug("Ensuring ensure scatter outputs respect intervals");
-
-                File subsetVcf = new File(vcf.getParentFile(), SequenceAnalysisService.get().getUnzippedBaseName(vcf.getName()) + ".subset.vcf.gz");
-                File subsetVcfIdx = new File(subsetVcf.getPath() + ".tbi");
-                manager.addIntermediateFile(subsetVcf);
-                manager.addIntermediateFile(subsetVcfIdx);
-
-                if (subsetVcfIdx.exists())
-                {
-                    getJob().getLogger().debug("Index exists, will not re-subset the VCF: " + subsetVcf.getName());
-                }
-                else
-                {
-                    OutputVariantsStartingInIntervalsStep.Wrapper wrapper = new OutputVariantsStartingInIntervalsStep.Wrapper(getJob().getLogger());
-                    wrapper.execute(vcf, subsetVcf, getPipelineJob().getIntervalsForTask());
-                }
-
-                toConcat.add(subsetVcf);
+                totalNull++;
+                continue;
             }
-            else
+            else if (!vcf.exists())
             {
-                toConcat.add(vcf);
+                missing.add(vcf);
+                continue;
             }
 
+            toConcat.add(vcf);
+
             manager.addInput(action, "Input VCF", vcf);
             manager.addIntermediateFile(vcf);
             manager.addIntermediateFile(new File(vcf.getPath() + ".tbi"));
         }
 
+        if (totalNull > 0 && !toConcat.isEmpty())
+        {
+            throw new PipelineJobException("The scatter jobs returned a mixture of null and non-null outputs");
+        }
+
         Set<Integer> genomeIds = new HashSet<>();
         getPipelineJob().getFiles().forEach(x -> genomeIds.add(x.getLibrary_id()));
         if (genomeIds.size() != 1)
@@ -216,29 +205,33 @@ private File runDefaultVariantMerge(JobContextImpl ctx, TaskFileManagerImpl mana
 
         ReferenceGenome genome = getPipelineJob().getSequenceSupport().getCachedGenome(genomeIds.iterator().next());
 
-        String basename = SequenceAnalysisService.get().getUnzippedBaseName(toConcat.get(0).getName());
-        File combined = new File(getPipelineJob().getAnalysisDirectory(), basename + ".vcf.gz");
-        File combinedIdx = new File(combined.getPath() + ".tbi");
-        if (combinedIdx.exists())
+        File combined = null;
+        if (!toConcat.isEmpty())
         {
-            getJob().getLogger().info("VCF exists, will not recreate: " + combined.getPath());
-        }
-        else
-        {
-            if (!missing.isEmpty())
+            String basename = SequenceAnalysisService.get().getUnzippedBaseName(toConcat.get(0).getName());
+            combined = new File(getPipelineJob().getAnalysisDirectory(), basename + ".vcf.gz");
+            File combinedIdx = new File(combined.getPath() + ".tbi");
+            if (combinedIdx.exists())
             {
-                throw new PipelineJobException("Missing one of more VCFs: " + missing.stream().map(File::getPath).collect(Collectors.joining(",")));
+                getJob().getLogger().info("VCF exists, will not recreate: " + combined.getPath());
             }
+            else
+            {
+                if (!missing.isEmpty())
+                {
+                    throw new PipelineJobException("Missing one of more VCFs: " + missing.stream().map(File::getPath).collect(Collectors.joining(",")));
+                }
 
-            boolean sortAfterMerge = handler instanceof VariantProcessingStep.SupportsScatterGather && ((VariantProcessingStep.SupportsScatterGather)handler).doSortAfterMerge();
-            combined = SequenceAnalysisService.get().combineVcfs(toConcat, combined, genome, getJob().getLogger(), true, null, sortAfterMerge);
+                boolean sortAfterMerge = handler instanceof VariantProcessingStep.SupportsScatterGather && ((VariantProcessingStep.SupportsScatterGather) handler).doSortAfterMerge();
+                combined = SequenceAnalysisService.get().combineVcfs(toConcat, combined, genome, getJob().getLogger(), true, null, sortAfterMerge);
+            }
+            manager.addOutput(action, "Merged VCF", combined);
         }
-        manager.addOutput(action, "Merged VCF", combined);
 
         if (handler instanceof VariantProcessingStep.SupportsScatterGather)
         {
             ctx.getLogger().debug("Running additional merge tasks");
-            ((VariantProcessingStep.SupportsScatterGather) handler).performAdditionalMergeTasks(ctx, getPipelineJob(), manager, genome, toConcat);
+            ((VariantProcessingStep.SupportsScatterGather) handler).performAdditionalMergeTasks(ctx, getPipelineJob(), manager, genome, toConcat, new ArrayList<>(jobToIntervalMap.keySet()));
         }
 
         return combined;

SequenceAnalysis/src/org/labkey/sequenceanalysis/pipeline/VariantProcessingRemoteSplitTask.java

@@ -117,8 +117,13 @@ private VariantProcessingJob getPipelineJob()
                 {
                     output = ((SequenceOutputHandler.TracksVCF)handler).finalizeScatterJobOutput(ctx, output);
 
-                    // If the output is still under the work dir, translate path. Otherwise it was already copied to the the source dir
-                    if (output.getPath().startsWith(_wd.getDir().getPath()))
+                    // If the output is still under the work dir, translate path. Otherwise it was already copied to the source dir
+                    if (output == null)
+                    {
+                        ctx.getLogger().debug("No output produced, adding null to scatter outputs");
+                        getPipelineJob().getScatterJobOutputs().put(getPipelineJob().getIntervalSetName(), null);
+                    }
+                    else if (output.getPath().startsWith(_wd.getDir().getPath()))
                     {
                         //NOTE: the VCF will be copied back to the source dir, so translate paths
                         String path = _wd.getRelativePath(output);
@@ -136,7 +141,14 @@ private VariantProcessingJob getPipelineJob()
                     throw new PipelineJobException(e);
                 }
 
-                getPipelineJob().getLogger().debug("Final scatter output: " + output.getPath());
+                if (output != null)
+                {
+                    getPipelineJob().getLogger().debug("Final scatter output: " + output.getPath());
+                }
+                else
+                {
+                    getPipelineJob().getLogger().debug("No primary scatter output produced");
+                }
             }
             else
             {

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/bampostprocessing/MarkDuplicatesStep.java

@@ -26,7 +26,7 @@
  */
 public class MarkDuplicatesStep extends AbstractCommandPipelineStep<MarkDuplicatesWrapper> implements BamProcessingStep
 {
-    public MarkDuplicatesStep(PipelineStepProvider provider, PipelineContext ctx)
+    public MarkDuplicatesStep(PipelineStepProvider<?> provider, PipelineContext ctx)
     {
         super(provider, ctx, new MarkDuplicatesWrapper(ctx.getLogger()));
     }