Add sorting to PBSV merge, since translocations can produce an out-of-order VCF

Author: bbimber

SequenceAnalysis/api-src/org/labkey/api/sequenceanalysis/SequenceAnalysisService.java

@@ -94,6 +94,8 @@ static public void setInstance(SequenceAnalysisService instance)
 
     abstract public File combineVcfs(List<File> files, File outputGz, ReferenceGenome genome, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel) throws PipelineJobException;
 
+    abstract public File combineVcfs(List<File> files, File outputGz, ReferenceGenome genome, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel, boolean sortAfterMerge) throws PipelineJobException;
+
     abstract public String getScriptPath(String moduleName, String path) throws PipelineJobException;
 
     abstract public void sortGxf(Logger log, File input, @Nullable File output) throws PipelineJobException;

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

@@ -67,6 +67,11 @@ default void performAdditionalMergeTasks(SequenceOutputHandler.JobContext ctx, P
         {
 
         }
+
+        default boolean doSortAfterMerge()
+        {
+            return false;
+        }
     }
 
     public static interface MayRequirePrepareTask

SequenceAnalysis/src/org/labkey/sequenceanalysis/SequenceAnalysisServiceImpl.java

@@ -1,15 +1,13 @@
 package org.labkey.sequenceanalysis;
 
 import htsjdk.samtools.util.FileExtensions;
-import htsjdk.tribble.Tribble;
 import htsjdk.tribble.index.Index;
 import htsjdk.tribble.index.IndexFactory;
 import htsjdk.variant.vcf.VCFCodec;
 import org.apache.commons.lang3.StringUtils;
 import org.apache.commons.lang3.SystemUtils;
 import org.apache.logging.log4j.LogManager;
 import org.apache.logging.log4j.Logger;
-import org.jetbrains.annotations.NotNull;
 import org.jetbrains.annotations.Nullable;
 import org.labkey.api.data.CompareType;
 import org.labkey.api.data.Container;
@@ -46,7 +44,6 @@
 import org.labkey.sequenceanalysis.pipeline.ReferenceGenomeImpl;
 import org.labkey.sequenceanalysis.pipeline.ReferenceLibraryPipelineJob;
 import org.labkey.sequenceanalysis.pipeline.SequenceTaskHelper;
-import org.labkey.sequenceanalysis.run.util.BgzipRunner;
 import org.labkey.sequenceanalysis.run.util.FastaIndexer;
 import org.labkey.sequenceanalysis.run.util.GxfSorter;
 import org.labkey.sequenceanalysis.run.util.IndexFeatureFileWrapper;
@@ -60,7 +57,6 @@
 import java.sql.ResultSet;
 import java.sql.SQLException;
 import java.util.ArrayList;
-import java.util.Arrays;
 import java.util.Collection;
 import java.util.Collections;
 import java.util.HashMap;
@@ -69,7 +65,6 @@
 import java.util.Map;
 import java.util.Set;
 import java.util.function.Function;
-import java.util.function.Predicate;
 
 /**
  * User: bimber
@@ -476,7 +471,13 @@ public String createReferenceLibrary(List<Integer> sequenceIds, Container c, Use
     @Override
     public File combineVcfs(List<File> files, File outputGz, ReferenceGenome genome, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel) throws PipelineJobException
     {
-        return SequenceUtil.combineVcfs(files, genome, outputGz, log, multiThreaded, compressionLevel);
+        return combineVcfs(files, outputGz, genome, log, multiThreaded, compressionLevel, false);
+    }
+
+    @Override
+    public File combineVcfs(List<File> files, File outputGz, ReferenceGenome genome, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel, boolean sortAfterMerge) throws PipelineJobException
+    {
+        return SequenceUtil.combineVcfs(files, genome, outputGz, log, multiThreaded, compressionLevel, sortAfterMerge);
     }
 
     @Override

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

@@ -230,7 +230,8 @@ private File runDefaultVariantMerge(JobContextImpl ctx, TaskFileManagerImpl mana
                 throw new PipelineJobException("Missing one of more VCFs: " + missing.stream().map(File::getPath).collect(Collectors.joining(",")));
             }
 
-            combined = SequenceAnalysisService.get().combineVcfs(toConcat, combined, genome, getJob().getLogger(), true, null);
+            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);
 

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/analysis/PbsvJointCallingHandler.java

@@ -179,6 +179,9 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
                         ctx.getFileManager().addIntermediateFile(SequenceAnalysisService.get().ensureVcfIndex(f, ctx.getLogger(), false));
                     }
                     vcfOutGz = SequenceUtil.combineVcfs(outputs, genome, new File(ctx.getOutputDir(), outputBaseName + ".vcf.gz"), ctx.getLogger(), true, null, false);
+
+                    // NOTE: the resulting file can be out of order due to translocations
+                    SequenceUtil.sortROD(vcfOutGz, ctx.getLogger(), 2);
                 }
 
                 SequenceAnalysisService.get().ensureVcfIndex(vcfOutGz, ctx.getLogger(), true);
@@ -437,4 +440,10 @@ public void verifyAndAddMissingSamples(JobContext ctx, File input, List<File> in
             throw new PipelineJobException(e);
         }
     }
+
+    @Override
+    public boolean doSortAfterMerge()
+    {
+        return true;
+    }
 }

SequenceAnalysis/src/org/labkey/sequenceanalysis/util/SequenceUtil.java

@@ -436,12 +436,12 @@ public static void sortROD(File input, Logger log, Integer startColumnIdx) throw
         sorted.delete();
     }
 
-    public static File combineVcfs(List<File> files, ReferenceGenome genome, File outputGzip, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel) throws PipelineJobException
+    public static File combineVcfs(List<File> files, ReferenceGenome genome, File outputGzip, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel, boolean sortAfterMerge) throws PipelineJobException
     {
-        return combineVcfs(files, genome, outputGzip, log, multiThreaded, compressionLevel, true);
+        return combineVcfs(files, genome, outputGzip, log, multiThreaded, compressionLevel, true, sortAfterMerge);
     }
 
-    public static File combineVcfs(List<File> files, ReferenceGenome genome, File outputGzip, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel, boolean showTotals) throws PipelineJobException
+    public static File combineVcfs(List<File> files, ReferenceGenome genome, File outputGzip, Logger log, boolean multiThreaded, @Nullable Integer compressionLevel, boolean showTotals, boolean sortAfterMerge) throws PipelineJobException
     {
         log.info("combining VCFs: ");
 
@@ -506,6 +506,12 @@ else if (!samples.equals(header.getGenotypeSamples()))
             SimpleScriptWrapper wrapper = new SimpleScriptWrapper(log);
             wrapper.execute(Arrays.asList("/bin/bash", bashTmp.getPath()));
 
+            if (sortAfterMerge)
+            {
+                log.debug("sorting VCF");
+                sortROD(outputGzip, log, 2);
+            }
+
             SequenceAnalysisService.get().ensureVcfIndex(outputGzip, log);
 
             bashTmp.delete();