Improve pindel calls

Author: bbimber

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

@@ -621,7 +621,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                 String[] line;
                 while ((line = reader.readNext()) != null)
                 {
-                    if (!("D".equals(line[0]) || "I".equals(line[0])))
+                    if (!("D".equals(line[0]) || "I".equals(line[0]) || "S".equals(line[0])))
                     {
                         continue;
                     }
@@ -639,7 +639,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                         continue;
                     }
 
-                    if ("D".equals(line[0]) && refLength > MAX_DELETION_LENGTH)
+                    if (("D".equals(line[0]) || "S".equals(line[0])) && refLength > MAX_DELETION_LENGTH)
                     {
                         continue;
                     }
@@ -655,7 +655,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                         eventCoverage = Double.parseDouble(line[11]);
                     }
 
-                    if ("D".equals(line[0]) && eventCoverage > MAX_DEL_EVENT_COVERAGE)
+                    if (("D".equals(line[0]) || "S".equals(line[0])) && eventCoverage > MAX_DEL_EVENT_COVERAGE)
                     {
                         continue;
                     }

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

@@ -8,6 +8,7 @@
 import htsjdk.samtools.reference.IndexedFastaSequenceFile;
 import htsjdk.samtools.reference.ReferenceSequence;
 import htsjdk.variant.utils.SAMSequenceDictionaryExtractor;
+import org.apache.commons.lang.StringUtils;
 import org.json.JSONObject;
 import org.labkey.api.pipeline.PipelineJobException;
 import org.labkey.api.reader.Readers;
@@ -308,7 +309,7 @@ private static void parsePindelOutput(PipelineContext ctx, CSVWriter writer, Fil
 
                     //NOTE: this is the indel region itself, no flanking. so for a deletion with REF/ALT of ATTC / A--C, it reports TT. for an insertion of ATT / AGTT, it reports G
                     String pindelAllele = tokens[2].split(" ")[2];
-                    pindelAllele = pindelAllele.replaceAll("\"", "");
+                    pindelAllele = StringUtils.trimToNull(pindelAllele.replaceAll("\"", ""));
 
                     File dict = new File(fasta.getPath().replace("fasta", "dict"));
                     if (!dict.exists())
@@ -353,6 +354,10 @@ private static void parsePindelOutput(PipelineContext ctx, CSVWriter writer, Fil
 
                         eventCoverage = eventCoverage / (baseAfterEnd - basePriorToStart - 1);
                     }
+                    else if ("I".equals(type))
+                    {
+                        eventCoverage = (double)depth;
+                    }
 
                     double meanCoverage = (leadingCoverage + trailingCoverage) / 2.0;
                     double pct = (double)support / meanCoverage;
@@ -368,24 +373,31 @@ private static void parsePindelOutput(PipelineContext ctx, CSVWriter writer, Fil
                         String ref = "";
                         if ("I".equals(type))
                         {
+                            if (pindelAllele == null)
+                            {
+                                throw new IllegalArgumentException("Unexpected empty allele for insertion: " + basePriorToStart);
+                            }
+
                             ref = sequence.getBaseString().substring(basePriorToStart-1, basePriorToStart);
                             alt = ref + pindelAllele;
                         }
                         else if ("D".equals(type))
                         {
                             ref = sequence.getBaseString().substring(basePriorToStart-1, trueEnd);
                             alt = sequence.getBaseString().substring(basePriorToStart-1, basePriorToStart);
-
-                            // Pindel reports the region over the deletion. so add the leading base to match the reference (VCF-style)
-                            // Our sequence based on coordinates should match what pindel reported
-                            String predictedPindelAllele = alt + pindelAllele;
-                            if (!predictedPindelAllele.equals(ref))
+                            if (pindelAllele != null)
                             {
-                                throw new IllegalArgumentException("Unexpected pindel allele: " + ref + " / " + predictedPindelAllele + " / " + pindelAllele);
+                                type = "S";
+                                alt = alt + pindelAllele;
+
+                                if (alt.length() != ref.length())
+                                {
+                                    throw new IllegalArgumentException("Unexpected pindel allele at " + basePriorToStart + ": " + ref + " / " + alt + " / " + pindelAllele);
+                                }
                             }
                         }
 
-                        writer.writeNext(new String[]{type, contig, String.valueOf(basePriorToStart), String.valueOf(trueEnd), String.valueOf(depth), String.valueOf(support), String.valueOf(pct), String.valueOf(meanCoverage), String.valueOf(leadingCoverage), String.valueOf(trailingCoverage), (eventCoverage == null ? "" : String.valueOf(eventCoverage)), ref, alt, pindelAllele});
+                        writer.writeNext(new String[]{type, contig, String.valueOf(basePriorToStart), String.valueOf(trueEnd), String.valueOf(depth), String.valueOf(support), String.valueOf(pct), String.valueOf(meanCoverage), String.valueOf(leadingCoverage), String.valueOf(trailingCoverage), (eventCoverage == null ? "" : String.valueOf(eventCoverage)), ref, alt, (pindelAllele == null ? "" : pindelAllele)});
                         totalPassing++;
                     }
                     else

jbrowse/build.gradle

@@ -93,6 +93,11 @@ project.task("copyJBrowseCss",
 	}
 })
 
+project.tasks.register("cleanJBrowse", Delete) {
+	delete files("./resources/web/jbrowseApp")
+}
+
+project.tasks.cleanModule.dependsOn(project.tasks.named('cleanJBrowse'))
 
 project.tasks.copyJBrowse.dependsOn(project.tasks.npmInstall)
 

jbrowse/package.json

@@ -7,8 +7,8 @@
     "build": "npm run build-dev",
     "start": "cross-env NODE_ENV=development LK_MODULE_CONTAINER=DiscvrLabKeyModules LK_MODULE=jbrowse webpack-dev-server --config config/watch.config.js",
     "build-dev": "npm run clean && cross-env NODE_ENV=development LK_MODULE_CONTAINER=DiscvrLabKeyModules LK_MODULE=jbrowse webpack --config config/dev.config.js --progress --profile",
-    "build-prod": "npm run clean && cross-env NODE_ENV=production PROD_SOURCE_MAP=source-map LK_MODULE_CONTAINER=DiscvrLabKeyModules LK_MODULE=jbrowse webpack --config config/prod.config.js --progress --profile --colors",
-    "clean": "rimraf resources/web/jbrowseApp && rimraf resources/web/jbrowse/gen && rimraf resources/views/gen && rimraf resources/views/browser*",
+    "build-prod": "npm run clean && cross-env NODE_ENV=production PROD_SOURCE_MAP=source-map LK_MODULE_CONTAINER=DiscvrLabKeyModules LK_MODULE=jbrowse webpack --config config/prod.config.js --progress --profile",
+    "clean": "rimraf resources/web/jbrowse/gen && rimraf resources/views/gen && rimraf resources/views/browser*",
     "test": "cross-env NODE_ENV=test jest"
   },
   "dependencies": {