Add option to count variants using bcftools index, which is much faster than parsing the file itself

Author: bbimber

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

@@ -0,0 +1,25 @@
+package org.labkey.api.sequenceanalysis.pipeline;
+
+import org.apache.logging.log4j.Logger;
+import org.jetbrains.annotations.Nullable;
+import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
+
+import java.io.File;
+
+/**
+ * User: bimber
+ * Date: 12/15/12
+ * Time: 9:11 PM
+ */
+public class BcftoolsRunner extends AbstractCommandWrapper
+{
+    public BcftoolsRunner(@Nullable Logger logger)
+    {
+        super(logger);
+    }
+
+    public File getSamtoolsPath()
+    {
+        return SequencePipelineService.get().getExeForPackage("SAMTOOLSPATH", "samtools");
+    }
+}

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

@@ -433,7 +433,7 @@ private void appendParents(Set<String> distinctSubjects, PedigreeRecord p, List<
     @Override
     public String getVCFLineCount(File vcf, Logger log, boolean passOnly) throws PipelineJobException
     {
-        return ProcessVariantsHandler.getVCFLineCount(vcf, log, passOnly);
+        return passOnly ? ProcessVariantsHandler.getVCFLineCount(vcf, log, passOnly, false) : ProcessVariantsHandler.getVCFLineCount(vcf, log, passOnly, true);
     }
 
     @Override

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

@@ -321,20 +321,20 @@ public void liftOverVcf(JobContext ctx, ReferenceGenome targetGenome, ReferenceG
         Long mapped = null;
         if (output.exists())
         {
-            String mappedStr = ProcessVariantsHandler.getVCFLineCount(output, job.getLogger(), false);
+            String mappedStr = ProcessVariantsHandler.getVCFLineCount(output, job.getLogger(), false, true);
             mapped = StringUtils.trimToNull(mappedStr) == null ? 0L : Long.parseLong(mappedStr);
             job.getLogger().info("total variants mapped: " + mappedStr);
-            job.getLogger().info("passing variants mapped: " + ProcessVariantsHandler.getVCFLineCount(output, job.getLogger(), true));
+            job.getLogger().info("passing variants mapped: " + ProcessVariantsHandler.getVCFLineCount(output, job.getLogger(), true, false));
             SequenceAnalysisService.get().ensureVcfIndex(output, job.getLogger());
         }
 
         Long unmapped = 0L;
         if (unmappedOutput != null && unmappedOutput.exists())
         {
-            String unmappedStr = ProcessVariantsHandler.getVCFLineCount(unmappedOutput, job.getLogger(), false);
+            String unmappedStr = ProcessVariantsHandler.getVCFLineCount(unmappedOutput, job.getLogger(), false, true);
             unmapped = StringUtils.trimToNull(unmappedStr) == null ? 0L : Long.parseLong(unmappedStr);
             job.getLogger().info("total unmapped variants: " + unmappedStr);
-            job.getLogger().info("passing unmapped variants: " + ProcessVariantsHandler.getVCFLineCount(unmappedOutput, job.getLogger(), true));
+            job.getLogger().info("passing unmapped variants: " + ProcessVariantsHandler.getVCFLineCount(unmappedOutput, job.getLogger(), true, false));
             SequenceAnalysisService.get().ensureVcfIndex(unmappedOutput, job.getLogger());
         }
 

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

@@ -27,6 +27,7 @@
 import org.labkey.api.sequenceanalysis.SequenceAnalysisService;
 import org.labkey.api.sequenceanalysis.SequenceOutputFile;
 import org.labkey.api.sequenceanalysis.pipeline.AbstractResumer;
+import org.labkey.api.sequenceanalysis.pipeline.BcftoolsRunner;
 import org.labkey.api.sequenceanalysis.pipeline.PipelineStepCtx;
 import org.labkey.api.sequenceanalysis.pipeline.PipelineStepProvider;
 import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome;
@@ -430,8 +431,8 @@ public static File processVCF(File input, Integer libraryId, JobContext ctx, Res
             {
                 currentVCF = output.getVCF();
 
-                ctx.getJob().getLogger().info("total variants: " + getVCFLineCount(currentVCF, ctx.getJob().getLogger(), false));
-                ctx.getJob().getLogger().info("passing variants: " + getVCFLineCount(currentVCF, ctx.getJob().getLogger(), true));
+                ctx.getJob().getLogger().info("total variants: " + getVCFLineCount(currentVCF, ctx.getJob().getLogger(), false, true));
+                ctx.getJob().getLogger().info("passing variants: " + getVCFLineCount(currentVCF, ctx.getJob().getLogger(), true, false));
                 ctx.getJob().getLogger().debug("index exists: " + (new File(currentVCF.getPath() + ".tbi")).exists());
 
                 try
@@ -468,8 +469,23 @@ public static File processVCF(File input, Integer libraryId, JobContext ctx, Res
         return null;
     }
 
-    public static String getVCFLineCount(File vcf, Logger log, boolean passOnly) throws PipelineJobException
+    private static String countUsingBcfTools(File vcf, Logger log) throws PipelineJobException
     {
+        return new BcftoolsRunner(log).executeWithOutput(Arrays.asList("index", "-n", vcf.getPath()));
+    }
+
+    public static String getVCFLineCount(File vcf, Logger log, boolean passOnly, boolean useBcfTools) throws PipelineJobException
+    {
+        if (useBcfTools)
+        {
+            if (passOnly)
+            {
+                throw new PipelineJobException("bcftools VCF count cannot be used with passOnly");
+            }
+
+            return countUsingBcfTools(vcf, log);
+        }
+
         String cat = vcf.getName().endsWith(".gz") ? "zcat" : "cat";
         SimpleScriptWrapper wrapper = new SimpleScriptWrapper(null);