- Improve subset seurat code

- move closer to including pindel in variant VCF

Author: bbimber

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

@@ -227,10 +227,51 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         //Create a BED file with all regions of coverage below MIN_COVERAGE:
         int minCoverage = getProvider().getParameterByName("minCoverage").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class);
         int positionsSkipped = 0;
+        int basesRecoveredFromDeletions = 0;
         int gapIntervals = 0;
         double avgDepth;
 
-        //TODO: rescue pindel?
+        boolean runPindel = getProvider().getParameterByName("runPindel").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class, false);
+
+        List<VariantContext> pindelConsensusVariants = new ArrayList<>();
+        int totalPindelConsensusVariants = 0;
+        int totalPindelVariants = 0;
+        if (runPindel)
+        {
+            Double minFraction = getProvider().getParameterByName("minFraction").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Double.class, 0.0);
+            int minDepth = getProvider().getParameterByName("minDepth").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 0);
+            int minInsertSize = getProvider().getParameterByName("minInsertSize").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 0);
+
+            PindelAnalysis.PindelSettings settings = new PindelAnalysis.PindelSettings();
+
+            File pindelOutput = PindelAnalysis.runPindel(output, getPipelineCtx(), rs, outputDir, inputBam, referenceGenome.getWorkingFastaFile(), minFraction, minDepth, true, coverageOut, minInsertSize);
+            File pindelVcf = PindelAnalysis.createVcf(pindelOutput, new File(pindelOutput.getParentFile(), FileUtil.getBaseName(pindelOutput) + ".all.vcf.gz"), referenceGenome, settings);
+
+            try (VCFFileReader reader = new VCFFileReader(pindelVcf);CloseableIterator<VariantContext> it = reader.iterator())
+            {
+                while (it.hasNext())
+                {
+                    VariantContext vc = it.next();
+                    if (vc.hasAttribute("IN_CONSENSUS"))
+                    {
+                        pindelConsensusVariants.add(vc);
+                        totalPindelConsensusVariants++;
+                    }
+
+                    totalPindelVariants++;
+                }
+            }
+
+            getPipelineCtx().getLogger().info("Total pindel variants: " + totalPindelVariants);
+            getPipelineCtx().getLogger().info("Total consensus pindel variants: " + totalPindelConsensusVariants);
+            if (totalPindelConsensusVariants == 0)
+            {
+                getPipelineCtx().getLogger().info("deleting empty pindel VCF: " + pindelVcf.getPath());
+                pindelVcf.delete();
+                new File(pindelVcf.getPath() + ".tbi").delete();
+            }
+        }
+
         File mask = new File(outputDir, "mask.bed");
         Map<String, Integer> gatkDepth = new HashMap<>();
         try (CSVReader reader = new CSVReader(Readers.getReader(coverageOut), '\t');CSVWriter writer = new CSVWriter(IOUtil.openFileForBufferedUtf8Writing(mask), '\t', CSVWriter.NO_QUOTE_CHARACTER))
@@ -259,40 +300,63 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
 
                 if (depth < minCoverage)
                 {
-                    positionsSkipped++;
+                    //Check if within pindel SNV calls:
+                    Interval currentPosition = new Interval(tokens[0], Integer.parseInt(tokens[1]), Integer.parseInt(tokens[1]));
+                    boolean withinDeletion = false;
+                    for (VariantContext vc : pindelConsensusVariants)
+                    {
+                        if (vc.isIndel())
+                        {
+                            // If we overlap the SNV (which should be just deletions (note: the first position of the variant is the non-indel REF base), allow low coverage:
+                            if (currentPosition.overlaps(vc) && currentPosition.getStart() != vc.getStart())
+                            {
+                                basesRecoveredFromDeletions++;
+
+                                //TODO: enable this once sure:
+                                //withinDeletion = true;
+                                break;
+                            }
+                        }
+                    }
 
-                    if (intervalOfCurrentGap != null)
+                    if (!withinDeletion)
                     {
-                        if (intervalOfCurrentGap.getContig().equals(tokens[0]))
+                        positionsSkipped++;
+
+                        if (intervalOfCurrentGap != null)
                         {
-                            //extend
-                            intervalOfCurrentGap = new Interval(intervalOfCurrentGap.getContig(), intervalOfCurrentGap.getStart(), Integer.parseInt(tokens[1]));
+                            if (intervalOfCurrentGap.getContig().equals(tokens[0]))
+                            {
+                                //extend
+                                intervalOfCurrentGap = new Interval(intervalOfCurrentGap.getContig(), intervalOfCurrentGap.getStart(), Integer.parseInt(tokens[1]));
+                            }
+                            else
+                            {
+                                //switched contigs, write and make new:
+                                writer.writeNext(new String[]{intervalOfCurrentGap.getContig(), String.valueOf(intervalOfCurrentGap.getStart() - 1), String.valueOf(intervalOfCurrentGap.getEnd())});
+                                gapIntervals++;
+                                intervalOfCurrentGap = currentPosition;
+                            }
                         }
                         else
                         {
-                            //switched contigs, write and make new:
-                            writer.writeNext(new String[]{intervalOfCurrentGap.getContig(), String.valueOf(intervalOfCurrentGap.getStart()-1), String.valueOf(intervalOfCurrentGap.getEnd())});
-                            gapIntervals++;
-                            intervalOfCurrentGap = new Interval(tokens[0], Integer.parseInt(tokens[1]), Integer.parseInt(tokens[1]));
+                            //Not existing gap, just start one:
+                            intervalOfCurrentGap = currentPosition;
                         }
-                    }
-                    else
-                    {
-                        //Not existing gap, just start one:
-                        intervalOfCurrentGap = new Interval(tokens[0], Integer.parseInt(tokens[1]), Integer.parseInt(tokens[1]));
+
+                        continue;
                     }
                 }
-                else
-                {
-                    //We just existed a gap, so write:
-                    if (intervalOfCurrentGap != null)
-                    {
-                        writer.writeNext(new String[]{intervalOfCurrentGap.getContig(), String.valueOf(intervalOfCurrentGap.getStart()-1), String.valueOf(intervalOfCurrentGap.getEnd())});
-                        gapIntervals++;
-                    }
 
-                    intervalOfCurrentGap = null;
+                //Otherwise this is a valid position to report:
+                //We just exited a gap, so write:
+                if (intervalOfCurrentGap != null)
+                {
+                    writer.writeNext(new String[]{intervalOfCurrentGap.getContig(), String.valueOf(intervalOfCurrentGap.getStart()-1), String.valueOf(intervalOfCurrentGap.getEnd())});
+                    gapIntervals++;
                 }
+
+                intervalOfCurrentGap = null;
             }
 
             //Ensure we count final gap
@@ -422,43 +486,6 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         Set<String> variantsBcftools = runBcftools(inputBam, referenceGenome, mask, minCoverage);
         int variantsBcftoolsTotal = variantsBcftools.size();
 
-        boolean runPindel = getProvider().getParameterByName("runPindel").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Boolean.class, false);
-
-        List<VariantContext> pindelConsensusVariants = new ArrayList<>();
-        int totalPindelConsensusVariants = 0;
-        if (runPindel)
-        {
-            Double minFraction = getProvider().getParameterByName("minFraction").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Double.class, 0.0);
-            int minDepth = getProvider().getParameterByName("minDepth").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 0);
-            int minInsertSize = getProvider().getParameterByName("minInsertSize").extractValue(getPipelineCtx().getJob(), getProvider(), getStepIdx(), Integer.class, 0);
-
-            PindelAnalysis.PindelSettings settings = new PindelAnalysis.PindelSettings();
-
-            File pindelOutput = PindelAnalysis.runPindel(output, getPipelineCtx(), rs, outputDir, inputBam, referenceGenome.getWorkingFastaFile(), minFraction, minDepth, true, coverageOut, minInsertSize);
-            File pindelVcf = PindelAnalysis.createVcf(pindelOutput, new File(pindelOutput.getParentFile(), FileUtil.getBaseName(pindelOutput) + ".vcf.gz"), referenceGenome, settings);
-
-            try (VCFFileReader reader = new VCFFileReader(pindelVcf);CloseableIterator<VariantContext> it = reader.iterator())
-            {
-                while (it.hasNext())
-                {
-                    VariantContext vc = it.next();
-                    if (vc.hasAttribute("IN_CONSENSUS"))
-                    {
-                        pindelConsensusVariants.add(vc);
-                        totalPindelConsensusVariants++;
-                    }
-                }
-            }
-
-            getPipelineCtx().getLogger().info("Total consensus pindel variants: " + totalPindelConsensusVariants);
-            if (totalPindelConsensusVariants == 0)
-            {
-                getPipelineCtx().getLogger().info("deleting empty pindel VCF: " + pindelVcf.getPath());
-                pindelVcf.delete();
-                new File(pindelVcf.getPath() + ".tbi").delete();
-            }
-        }
-
         //Create final VCF:
         int totalVariants = 0;
         int totalGT1 = 0;
@@ -602,10 +629,11 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                 }
             }
 
-            //TODO: add pindel
             if (!pindelConsensusVariants.isEmpty())
             {
-
+                //TODO: enable once sure:
+                //pindelConsensusVariants.stream().forEach(consensusVariants::add);
+                getPipelineCtx().getLogger().info("total consensus variants that would be added: " + pindelConsensusVariants.size());
             }
 
             try (CloseableIterator<VariantContext> iterator = allVariants.iterator())
@@ -634,6 +662,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         double pctNoCover = positionsSkipped / (double)dict.getReferenceLength();
         getPipelineCtx().getLogger().info("Total positions with coverage below threshold (" + minCoverage + "): " + positionsSkipped + "(" + fmt.format(pctNoCover) + ")");
         getPipelineCtx().getLogger().info("Total # gap intervals: " + gapIntervals);
+        getPipelineCtx().getLogger().info("Total # of low coverage bases recovered within consensus deletions: " + basesRecoveredFromDeletions);
 
         String description = String.format("Total Variants: %s\nTotal GT 1 PCT: %s\nTotal GT 50 PCT: %s\nTotal Indel GT 1 PCT: %s\nPositions Below Coverage: %s\nTotal In LoFreq Consensus: %s\nTotal Indel In LoFreq Consensus: %s\nTotal Consensus Variant in PBS: %s", totalVariants, totalGT1, totalGT50, totalIndelGT1, positionsSkipped, totalGTThreshold, totalIndelGTThreshold, totalConsensusInPBS);
         description += "\n" + "Strand Bias Recovered: " + filteredVariantsRecovered;

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

@@ -468,6 +468,7 @@ public static File createVcf(File pindelOutput, File vcfOutput, ReferenceGenome
             String[] line;
             while ((line = reader.readNext()) != null)
             {
+                int inConsensus = 1;
                 if (!("D".equals(line[0]) || "I".equals(line[0]) || "S".equals(line[0])))
                 {
                     continue;
@@ -483,17 +484,17 @@ public static File createVcf(File pindelOutput, File vcfOutput, ReferenceGenome
                 // Assume LoFreq calls these well enough:
                 if (refLength < settings.MIN_LENGTH_TO_CONSIDER && altLength < settings.MIN_LENGTH_TO_CONSIDER)
                 {
-                    continue;
+                    inConsensus = 0;
                 }
 
                 if (("D".equals(line[0]) || "S".equals(line[0])) && refLength > settings.MAX_DELETION_LENGTH)
                 {
-                    continue;
+                    inConsensus = 0;
                 }
 
                 if (Double.parseDouble(line[6]) < settings.MIN_AF)
                 {
-                    continue;
+                    inConsensus = 0;
                 }
 
                 double eventCoverage = 0.0;
@@ -504,7 +505,7 @@ public static File createVcf(File pindelOutput, File vcfOutput, ReferenceGenome
 
                 if (("D".equals(line[0]) || "S".equals(line[0])) && eventCoverage > settings.MAX_DEL_EVENT_COVERAGE)
                 {
-                    continue;
+                    inConsensus = 0;
                 }
 
                 VariantContextBuilder vcb = new VariantContextBuilder();
@@ -514,7 +515,7 @@ public static File createVcf(File pindelOutput, File vcfOutput, ReferenceGenome
                 vcb.alleles(Arrays.asList(Allele.create(refAllele, true), Allele.create(altAllele)));
                 vcb.attribute(VCFConstants.ALLELE_FREQUENCY_KEY, Double.parseDouble(line[6]));
 
-                vcb.attribute("IN_CONSENSUS", 1);
+                vcb.attribute("IN_CONSENSUS", inConsensus);
                 vcb.attribute("GATK_DP", (int)Double.parseDouble(line[7]));
 
                 int dp = "I".equals(line[0]) ? Integer.parseInt(line[4]) : (int) Double.parseDouble(line[10]);

singlecell/api-src/org/labkey/api/singlecell/pipeline/AbstractSingleCellPipelineStep.java

@@ -323,7 +323,7 @@ else if (NumberUtils.isCreatable(val))
         }
         else if ("sequenceanalysis-trimmingtextarea".equals(pd.getFieldXtype()))
         {
-            val = val.replace("'", "\'");
+            val = val.replace("'", "\\\'");
             String[] vals = val.split(",");
             return "c('" + StringUtils.join(vals, "','") + "')";
         }

singlecell/resources/chunks/SubsetSeurat.R

@@ -1,17 +1,19 @@
 for (datasetId in names(seuratObjects)) {
-    seuratObj <- seuratObjects[[datasetId]]
-    seuratObjects[[datasetId]] <- NULL
+	seuratObj <- seuratObjects[[datasetId]]
+	seuratObjects[[datasetId]] <- NULL
 
-    #TODO: this is a stopgap for a former bug in RunCellHashing. Retain until all existing seurat objects lacking this field are removed
-    if (!'HTO.Classification' %in% names(seuratObj@meta.data) && 'consensuscall.global' %in% names(seuratObj@meta.data)) {
-        seuratObj$HTO.Classification <- seuratObj$consensuscall.global
-    }
+	#TODO: this is a stopgap for a former bug in RunCellHashing. Retain until all existing seurat objects lacking this field are removed
+	if (!'HTO.Classification' %in% names(seuratObj@meta.data) && 'consensuscall.global' %in% names(seuratObj@meta.data)) {
+		seuratObj$HTO.Classification <- seuratObj$consensuscall.global
+	}
 
-    <SUBSETS>
+	<SUBSET_CODE>
 
-    newSeuratObjects[[datasetId]] <- seuratObj
+	# Cleanup
+	rm(seuratObj)
+	gc()
+}
 
-    # Cleanup
-    rm(seuratObj)
-    gc()
+if (length(newSeuratObjects) == 0) {
+	stop('No cells remained in any seurat objects after subsetting')
 }

singlecell/src/org/labkey/singlecell/pipeline/singlecell/SubsetSeurat.java

@@ -55,7 +55,7 @@ public Collection<String> getRLibraries()
         return ret;
     }
 
-    final static String EXPRESSION = "<SUBSETS>";
+    final static String EXPRESSION = "<SUBSET_CODE>";
     final static String DELIM = "<>";
 
     @Override
@@ -77,9 +77,25 @@ protected List<String> loadChunkFromFile() throws PipelineJobException
             {
                 for (String subset : values)
                 {
-                    String toSub = "seuratObj <- subset(seuratObj, subset = " + subset + ")";
-                    ret.add(line.replaceAll(EXPRESSION, toSub));
-                    ret.add(line.replaceAll(EXPRESSION, "print(paste0('Cells after subset: ', ncol(seuratObj)))"));
+                    String subsetEscaped = subset.replace("'", "\\\'");
+
+                    ret.add("\tcells <- c()");
+                    ret.add("\ttryCatch({");
+                    ret.add("\t\tcells <- WhichCells(so, expression = " + subset + ")");
+                    ret.add("\t}, error = function(e){");
+                    ret.add("\t\tif (!is.null(e) && e$message == 'Cannot find cells provided') {");
+                    ret.add("\t\t\tprint(paste0('There were no cells remaining after the subset: ', '" + subsetEscaped + "'))");
+                    ret.add("\t\t}");
+                    ret.add("\t})");
+                    ret.add("");
+                    ret.add("\tif (length(cells) == 0) {");
+                    ret.add("\t\tprint(paste0('There were no cells after subsetting for dataset: ', datasetId, ', with subset: ', '" + subsetEscaped + "'))");
+                    ret.add("\t} else {");
+                    ret.add("\t\tseuratObj <- subset(seuratObj, cells = cells)");
+                    ret.add("\t\tprint(paste0('Cells after subset: ', ncol(seuratObj)))");
+                    ret.add("\t\tnewSeuratObjects[[datasetId]] <- seuratObj");
+                    ret.add("\t}");
+                    ret.add("");
                 }
             }
             else