Fix error in GenotypeGVCF intervals for scatter job and ensure resulting VCFs respect intervals during merge

Author: bbimber

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

@@ -2,6 +2,7 @@
 
 import htsjdk.samtools.util.Interval;
 import org.apache.commons.lang3.StringUtils;
+import org.apache.logging.log4j.Logger;
 import org.jetbrains.annotations.NotNull;
 import org.labkey.api.pipeline.AbstractTaskFactory;
 import org.labkey.api.pipeline.AbstractTaskFactorySettings;
@@ -14,10 +15,13 @@
 import org.labkey.api.sequenceanalysis.SequenceOutputFile;
 import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome;
 import org.labkey.api.sequenceanalysis.pipeline.SequenceOutputHandler;
+import org.labkey.api.sequenceanalysis.run.AbstractDiscvrSeqWrapper;
 import org.labkey.api.util.FileType;
+import org.labkey.api.writer.PrintWriters;
 
 import java.io.File;
 import java.io.IOException;
+import java.io.PrintWriter;
 import java.util.ArrayList;
 import java.util.Collections;
 import java.util.HashSet;
@@ -170,7 +174,49 @@ private File runDefaultVariantMerge(TaskFileManagerImpl manager, RecordedAction
                 missing.add(vcf);
             }
 
-            toConcat.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");
+                List<Interval> expectedIntervals = jobToIntervalMap.get(name);
+
+                File intervalFile = new File(vcf.getParentFile(), "scatterIntervals.list");
+                File subsetVcf = new File(vcf.getParentFile(), SequenceAnalysisService.get().getUnzippedBaseName(vcf.getName()) + ".subset.vcf.gz");
+                File subsetVcfIdx = new File(subsetVcf.getPath() + ".tbi");
+                manager.addIntermediateFile(intervalFile);
+                manager.addIntermediateFile(subsetVcf);
+                manager.addIntermediateFile(subsetVcfIdx);
+
+                if (subsetVcfIdx.exists())
+                {
+                    getJob().getLogger().debug("Index exists, will not re-submit the VCF: " + subsetVcf.getName());
+                }
+                else
+                {
+                    try (PrintWriter writer = PrintWriters.getPrintWriter(intervalFile))
+                    {
+                        expectedIntervals.forEach(interval -> {
+                            writer.println(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd());
+                        });
+                    }
+                    catch (IOException e)
+                    {
+                        throw new PipelineJobException(e);
+                    }
+
+                    Wrapper wrapper = new Wrapper(getJob().getLogger());
+                    wrapper.execute(vcf, subsetVcf, intervalFile);
+                }
+
+                toConcat.add(subsetVcf);
+            }
+            else
+            {
+                toConcat.add(vcf);
+            }
+
             manager.addInput(action, "Input VCF", vcf);
             manager.addIntermediateFile(vcf);
             manager.addIntermediateFile(new File(vcf.getPath() + ".tbi"));
@@ -204,4 +250,33 @@ private File runDefaultVariantMerge(TaskFileManagerImpl manager, RecordedAction
 
         return combined;
     }
+
+    public static class Wrapper extends AbstractDiscvrSeqWrapper
+    {
+        public Wrapper(Logger log)
+        {
+            super(log);
+        }
+
+        public void execute(File inputVcf, File outputVcf, File intervalFile) throws PipelineJobException
+        {
+            List<String> args = new ArrayList<>(getBaseArgs());
+            args.add("OutputVariantsStartingInIntervals");
+
+            args.add("-V");
+            args.add(inputVcf.getPath());
+
+            args.add("-O");
+            args.add(outputVcf.getPath());
+
+            args.add("-L");
+            args.add(intervalFile.getPath());
+
+            execute(args);
+            if (!outputVcf.exists())
+            {
+                throw new PipelineJobException("Missing file: " + outputVcf.getPath());
+            }
+        }
+    }
 }