Allow user to select seurat and/or multiseq HTO methods

Author: bbimber

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

@@ -247,6 +247,11 @@ public static interface TracksVCF
     {
         public File getScatterJobOutput(JobContext ctx) throws PipelineJobException;
 
+        default File finalizeScatterJobOutput(JobContext ctx, File primaryOutput) throws PipelineJobException
+        {
+            return primaryOutput;
+        }
+
         public SequenceOutputFile createFinalSequenceOutput(PipelineJob job, File processed, List<SequenceOutputFile> inputFiles) throws PipelineJobException;
     }
 

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

@@ -126,5 +126,5 @@ static public void setInstance(SequencePipelineService instance)
 
     abstract public PreprocessingStep.Output simpleTrimFastqPair(File fq1, File fq2, List<String> params, File outDir, Logger log) throws PipelineJobException;
 
-    abstract public File runCiteSeqCount(PipelineStepOutput output, @Nullable String outputCategory, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, boolean generateHtoCalls, boolean createOutputFiles) throws PipelineJobException;
+    abstract public File runCiteSeqCount(PipelineStepOutput output, @Nullable String outputCategory, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, boolean generateHtoCalls, boolean createOutputFiles, boolean useSeurat, boolean useMultiSeq) throws PipelineJobException;
 }

SequenceAnalysis/resources/external/scRNAseq/htoClassifier.Rmd

@@ -41,7 +41,7 @@ if (nrow(barcodeData) > 0 && ncol(barcodeData) > 0){
 ```{r GenerateCalls, fig.width=12}
 
 if (nrow(barcodeData) > 0 && ncol(barcodeData) > 0){
-    dt <- GenerateCellHashingCalls(barcodeData = barcodeData, outFile = finalCallFile, allCallsOutFile = allCallsOutFile)
+    dt <- GenerateCellHashingCalls(barcodeData = barcodeData, outFile = finalCallFile, allCallsOutFile = allCallsOutFile, useSeurat = useSeurat, useMultiSeq = useMultiSeq)
 
     if (exists('whitelistFile') && !is.null(whitelistFile)){
       GenerateSummaryForExpectedBarcodes(dt, whitelistFile=whitelistFile, outputFile=metricsFile, barcodeData=barcodeData)

SequenceAnalysis/resources/external/scRNAseq/htoClassifier.sh

@@ -14,9 +14,12 @@ RAW_CALLS=$4
 DO_HTO_FILTER=$5
 MIN_READS_PER_CELL=$6
 METRICS_FILE=$7
+USE_SEURAT=$8
+USE_MULTISEQ=$9
+
 WHITELIST="whitelistFile<-NULL;"
-if [ $# -ge 8 ];then
-    WHITELIST="whitelistFile<-'"${8}"';"
+if [ $# -ge 10 ];then
+    WHITELIST="whitelistFile<-'"${10}"';"
 fi
 
 ENV_OPTS=""
@@ -33,4 +36,4 @@ fi
 
 sudo $DOCKER pull bimberlab/oosap
 
-sudo $DOCKER run --rm=true $RAM_OPTS $ENV_OPTS -v "${WD}:/work" -v "${HOME}:/homeDir" -u $UID -e USERID=$UID -w /work -e HOME=/homeDir bimberlab/oosap Rscript -e "barcodeDir <- '"${CITESEQ_COUNT_DIR}"';finalCallFile <- '"${FINAL_CALLS}"';doHtoFilter <- "${DO_HTO_FILTER}";maxValueForColSumFilter <- "${MIN_READS_PER_CELL}";allCallsOutFile <- '"${RAW_CALLS}"';metricsFile <- '"${METRICS_FILE}"';"${WHITELIST}"rmarkdown::render('htoClassifier.Rmd', output_file = '"${HTML_FILE}"')"
+sudo $DOCKER run --rm=true $RAM_OPTS $ENV_OPTS -v "${WD}:/work" -v "${HOME}:/homeDir" -u $UID -e USERID=$UID -w /work -e HOME=/homeDir bimberlab/oosap Rscript -e "barcodeDir <- '"${CITESEQ_COUNT_DIR}"';finalCallFile <- '"${FINAL_CALLS}"';doHtoFilter <- "${DO_HTO_FILTER}";maxValueForColSumFilter <- "${MIN_READS_PER_CELL}";allCallsOutFile <- '"${RAW_CALLS}"';metricsFile <- '"${METRICS_FILE}"';useSeurat <- ${USE_SEURAT};useMultiSeq <- ${USE_MULTISEQ};"${WHITELIST}"rmarkdown::render('htoClassifier.Rmd', output_file = '"${HTML_FILE}"')"

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

@@ -527,10 +527,10 @@ else if (!SequenceUtil.hasLineCount(f))
     }
 
     @Override
-    public File runCiteSeqCount(PipelineStepOutput output, @Nullable String outputCategory, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, boolean generateHtoCalls, boolean createOutputFiles) throws PipelineJobException
+    public File runCiteSeqCount(PipelineStepOutput output, @Nullable String outputCategory, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, boolean generateHtoCalls, boolean createOutputFiles, boolean useSeurat, boolean useMultiSeq) throws PipelineJobException
     {
         CellHashingHandler handler = new CellHashingHandler();
 
-        return handler.runCiteSeqCount(output, outputCategory, htoReadset, htoList, cellBarcodeList, outputDir, basename, log, extraArgs, doHtoFiltering, minCountPerCell, localPipelineDir, editDistance, scanEditDistances, parentReadset, genomeId, generateHtoCalls ? CellHashingHandler.BARCODE_TYPE.hashing : CellHashingHandler.BARCODE_TYPE.citeseq, createOutputFiles);
+        return handler.runCiteSeqCount(output, outputCategory, htoReadset, htoList, cellBarcodeList, outputDir, basename, log, extraArgs, doHtoFiltering, minCountPerCell, localPipelineDir, editDistance, scanEditDistances, parentReadset, genomeId, generateHtoCalls ? CellHashingHandler.BARCODE_TYPE.hashing : CellHashingHandler.BARCODE_TYPE.citeseq, createOutputFiles, useSeurat, useMultiSeq);
     }
 }

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

@@ -443,7 +443,7 @@ else if (line.startsWith("Percentage unmapped"))
         if (type.doGenerateCalls())
         {
             ctx.getJob().setStatus(PipelineJob.TaskStatus.running, "Generating HTO calls for edit distance: " + editDistance);
-            File htoCalls = generateFinalCalls(outputMatrix.getParentFile(), ctx.getOutputDir(), outputBasename, ctx.getLogger(), null, true, minCountPerCell, ctx.getSourceDirectory());
+            File htoCalls = generateFinalCalls(outputMatrix.getParentFile(), ctx.getOutputDir(), outputBasename, ctx.getLogger(), null, true, minCountPerCell, ctx.getSourceDirectory(), true, true);
             File html = new File(htoCalls.getParentFile(), outputBasename + ".html");
 
             if (!html.exists())
@@ -664,7 +664,7 @@ private File ensureLocalCopy(File input, File outputDir, Logger log, Set<File> t
         return input;
     }
 
-    public File generateFinalCalls(File citeSeqCountOutDir, File outputDir, String basename, Logger log, @Nullable File cellBarcodeWhitelist, boolean doHtoFiltering, Integer minCountPerCell, File localPipelineDir) throws PipelineJobException
+    private File generateFinalCalls(File citeSeqCountOutDir, File outputDir, String basename, Logger log, @Nullable File cellBarcodeWhitelist, boolean doHtoFiltering, Integer minCountPerCell, File localPipelineDir, boolean useSeurat, boolean useMultiSeq) throws PipelineJobException
     {
         log.debug("generating final calls from folder: " + citeSeqCountOutDir.getPath());
 
@@ -694,6 +694,9 @@ public File generateFinalCalls(File citeSeqCountOutDir, File outputDir, String b
         File rawCallsFile = new File(outputDir, basename + ".raw.txt");
         File metricsFile = getMetricsFile(callsFile);
         List<String> args = new ArrayList<>(Arrays.asList("/bin/bash", scriptWrapper, citeSeqCountOutDir.getName(), htmlFile.getName(), callsFile.getName(), rawCallsFile.getName(), (doHtoFiltering ? "T" : "F"), (minCountPerCell == null ? "0" : minCountPerCell.toString()), metricsFile.getName()));
+        args.add(useSeurat ? "TRUE" : "FALSE");
+        args.add(useMultiSeq ? "TRUE" : "FALSE");
+
         if (cellBarcodeWhitelist != null)
         {
             args.add(cellBarcodeWhitelist.getName());
@@ -876,7 +879,7 @@ else if (htoReadset.getReadData().size() != 1)
         return ret;
     }
 
-    public File runCiteSeqCount(PipelineStepOutput output, String category, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, BARCODE_TYPE type, boolean createOutputFiles) throws PipelineJobException
+    public File runCiteSeqCount(PipelineStepOutput output, String category, Readset htoReadset, File htoList, File cellBarcodeList, File outputDir, String basename, Logger log, List<String> extraArgs, boolean doHtoFiltering, @Nullable Integer minCountPerCell, File localPipelineDir, @Nullable Integer editDistance, boolean scanEditDistances, Readset parentReadset, @Nullable Integer genomeId, BARCODE_TYPE type, boolean createOutputFiles, boolean useSeurat, boolean useMultiSeq) throws PipelineJobException
     {
         HtoMergeResult htoFastqs = possiblyMergeHtoFastqs(htoReadset, outputDir, log);
         if (!htoFastqs.intermediateFiles.isEmpty())
@@ -963,7 +966,7 @@ public File runCiteSeqCount(PipelineStepOutput output, String category, Readset
 
             File citeSeqCountOutDir = new File(outputDir, basename + ".citeSeqCounts." + ed + "." + type.name());
             String outputBasename = basename + "." + ed + "." + type.name();
-            Map<String, Object> callMap = executeCiteSeqCountWithJobCtx(outputDir, outputBasename, citeSeqCountOutDir, htoFastqs.files.getLeft(), htoFastqs.files.getRight(), toolArgs, ed, log, cellBarcodeList, doHtoFiltering, minCountPerCell, localPipelineDir, unknownBarcodeFile, type);
+            Map<String, Object> callMap = executeCiteSeqCountWithJobCtx(outputDir, outputBasename, citeSeqCountOutDir, htoFastqs.files.getLeft(), htoFastqs.files.getRight(), toolArgs, ed, log, cellBarcodeList, doHtoFiltering, minCountPerCell, localPipelineDir, unknownBarcodeFile, type, useSeurat, useMultiSeq);
             results.put(ed, callMap);
 
             if (type.doGenerateCalls())
@@ -1064,7 +1067,7 @@ public File runCiteSeqCount(PipelineStepOutput output, String category, Readset
         }
     }
 
-    private Map<String, Object> executeCiteSeqCountWithJobCtx(File outputDir, String basename, File citeSeqCountOutDir, File fastq1, File fastq2, List<String> baseArgs, Integer ed, Logger log, File cellBarcodeList, boolean doHtoFiltering, Integer minCountPerCell, File localPipelineDir, File unknownBarcodeFile, BARCODE_TYPE type) throws PipelineJobException
+    private Map<String, Object> executeCiteSeqCountWithJobCtx(File outputDir, String basename, File citeSeqCountOutDir, File fastq1, File fastq2, List<String> baseArgs, Integer ed, Logger log, File cellBarcodeList, boolean doHtoFiltering, Integer minCountPerCell, File localPipelineDir, File unknownBarcodeFile, BARCODE_TYPE type, boolean useSeurat, boolean useMultiSeq) throws PipelineJobException
     {
         CellHashingHandler.CiteSeqCountWrapper wrapper = new CellHashingHandler.CiteSeqCountWrapper(log);
         File doneFile = new File(citeSeqCountOutDir, "citeSeqCount." + type.name() + "." + ed + ".done");
@@ -1131,7 +1134,7 @@ else if (line.startsWith("Percentage unmapped"))
 
         if (type.doGenerateCalls())
         {
-            File htoCalls = generateFinalCalls(outputMatrix.getParentFile(), outputDir, basename, log, cellBarcodeList, doHtoFiltering, minCountPerCell, localPipelineDir);
+            File htoCalls = generateFinalCalls(outputMatrix.getParentFile(), outputDir, basename, log, cellBarcodeList, doHtoFiltering, minCountPerCell, localPipelineDir, useSeurat, useMultiSeq);
             if (!htoCalls.exists())
             {
                 throw new PipelineJobException("missing expected file: " + htoCalls.getPath());

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

@@ -110,10 +110,20 @@ private VariantProcessingJob getPipelineJob()
                 File output = ((SequenceOutputHandler.TracksVCF)handler).getScatterJobOutput(ctx);
                 try
                 {
-                    //NOTE: the VCF was copied back to the source dir, so translate paths
-                    String path = _wd.getRelativePath(output);
-                    output = new File(ctx.getSourceDirectory(), path);
-                    getPipelineJob().getScatterJobOutputs().put(getPipelineJob().getIntervalSetName(), output);
+                    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()))
+                    {
+                        //NOTE: the VCF will be copied back to the source dir, so translate paths
+                        String path = _wd.getRelativePath(output);
+                        output = new File(ctx.getSourceDirectory(), path);
+                        getPipelineJob().getScatterJobOutputs().put(getPipelineJob().getIntervalSetName(), output);
+                    }
+                    else
+                    {
+                        ctx.getLogger().debug("Output has already been moved from workdir: " + output.getPath());
+                    }
                 }
                 catch (IOException e)
                 {