Merge pull request #104 from LabKey/fb_merge_discvr-21.3

Merge discvr-21.3 to develop

Author: bbimber

SequenceAnalysis/build.gradle

@@ -67,3 +67,20 @@ if (project.findProject(BuildUtils.getTestProjectPath(project.gradle)) != null &
       dependsOn(createPipelineConfigTask)
    }
 }
+
+project.task("copyJars",
+        type: Copy,
+        group: "Build",
+        description: "Copy commons-math3 JAR to module's lib directory",
+        { CopySpec copy ->
+           copy.setDuplicatesStrategy(DuplicatesStrategy.EXCLUDE)
+           copy.from(project.configurations.external)
+           copy.into new File("${project.labkey.explodedModuleLibDir}")
+           copy.include {
+              "**commons-math3-**.jar"
+           }
+        }
+)
+
+project.tasks.module.dependsOn(project.tasks.copyJars)
+project.tasks.copyJars.mustRunAfter(project.tasks.populateExplodedLib)

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

@@ -131,6 +131,7 @@
 import org.labkey.sequenceanalysis.run.reference.SavedReferenceLibraryStep;
 import org.labkey.sequenceanalysis.run.reference.VirusReferenceLibraryStep;
 import org.labkey.sequenceanalysis.run.util.CombineGVCFsHandler;
+import org.labkey.sequenceanalysis.run.util.FastqcRunner;
 import org.labkey.sequenceanalysis.run.util.GenomicsDBAppendHandler;
 import org.labkey.sequenceanalysis.run.util.GenomicsDBImportHandler;
 import org.labkey.sequenceanalysis.run.variant.CombineVariantsHandler;
@@ -476,7 +477,8 @@ public Set<Class> getUnitTests()
                 ProcessVariantsHandler.TestCase.class,
                 VariantProcessingJob.TestCase.class,
                 ScatterGatherUtils.TestCase.class,
-                ChainFileValidator.TestCase.class
+                ChainFileValidator.TestCase.class,
+                FastqcRunner.TestCase.class
         );
     }
 

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

@@ -613,6 +613,11 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
             File pindelOutput = PindelAnalysis.runPindel(output, getPipelineCtx(), rs, outputDir, inputBam, referenceGenome.getWorkingFastaFile(), minFraction, minDepth, true, coverageOut, minInsertSize);
             try (CSVReader reader = new CSVReader(Readers.getReader(pindelOutput), '\t'))
             {
+                final int MAX_DEL_EVENT_COVERAGE = 20;
+                final double MIN_AF = 0.25;
+                final int MIN_LENGTH_TO_CONSIDER = 10;
+                final int MAX_DELETION_LENGTH = 5000;
+
                 String[] line;
                 while ((line = reader.readNext()) != null)
                 {
@@ -621,10 +626,50 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                         continue;
                     }
 
-                    if (Double.parseDouble(line[6]) >= 0.35)
+                    int start = Integer.parseInt(line[2]);  //1-based, coordinate prior, like VCF
+                    int end = Integer.parseInt(line[3]);  //1-based, actual coordinate, like VCF
+                    String refAllele = line[11];
+                    String altAllele = line[12];
+                    int refLength = end - start;
+                    int altLength = altAllele.length();
+
+                    // Assume LoFreq calls these well enough:
+                    if (refLength < MIN_LENGTH_TO_CONSIDER && altLength < MIN_LENGTH_TO_CONSIDER)
+                    {
+                        continue;
+                    }
+
+                    if ("D".equals(line[0]) && refLength > MAX_DELETION_LENGTH)
+                    {
+                        continue;
+                    }
+
+                    if (Double.parseDouble(line[6]) < MIN_AF)
                     {
-                        indelMap.put(line[0], indelMap.getOrDefault(line[0], 0) + 1);
+                        continue;
+                    }
+
+                    double eventCoverage = 0.0;
+                    if (StringUtils.trimToNull(line[11]) != null)
+                    {
+                        eventCoverage = Double.parseDouble(line[11]);
+                    }
+
+                    if ("D".equals(line[0]) && eventCoverage > MAX_DEL_EVENT_COVERAGE)
+                    {
+                        continue;
                     }
+
+                    indelMap.put(line[0], indelMap.getOrDefault(line[0], 0) + 1);
+
+                    VariantContextBuilder vcb = new VariantContextBuilder();
+                    vcb.start(start);
+                    vcb.stop(end);
+                    vcb.chr(line[1]);
+                    vcb.alleles(Arrays.asList(Allele.create(refAllele, true), Allele.create(altAllele)));
+                    vcb.attribute("AF", Double.parseDouble(line[6]));
+                    int dp = "I".equals(line[0]) ? Integer.parseInt(line[4]) : (int)Double.parseDouble(line[10]);
+                    vcb.attribute("DP", dp);
                 }
             }
             catch (IOException e)
@@ -670,7 +715,9 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
             if (pangolinData != null)
             {
                 writer.writeNext(new String[]{"Pangolin", "PangolinLineage", pangolinData[1]});
-                writer.writeNext(new String[]{"Pangolin", "PangolinLineageConfidence", pangolinData[2]});
+                writer.writeNext(new String[]{"Pangolin", "PangolinConflicts", pangolinData[2]});
+                writer.writeNext(new String[]{"Pangolin", "PangolinVersions", pangolinData[3]});
+                writer.writeNext(new String[]{"Pangolin", "PangolinVersions", pangolinData[4]});
             }
             else
             {

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

@@ -3,6 +3,7 @@
 import au.com.bytecode.opencsv.CSVReader;
 import au.com.bytecode.opencsv.CSVWriter;
 import htsjdk.samtools.util.IOUtil;
+import org.apache.commons.io.FileUtils;
 import org.apache.commons.lang3.StringUtils;
 import org.apache.logging.log4j.Logger;
 import org.json.JSONObject;
@@ -30,6 +31,7 @@
 import org.labkey.api.sequenceanalysis.pipeline.SequenceOutputHandler;
 import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService;
 import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper;
+import org.labkey.api.util.FileUtil;
 import org.labkey.api.util.PageFlowUtil;
 import org.labkey.sequenceanalysis.SequenceAnalysisModule;
 import org.labkey.sequenceanalysis.SequenceAnalysisSchema;
@@ -102,6 +104,10 @@ public void processFilesOnWebserver(PipelineJob job, SequenceAnalysisJobSupport
                     row1.put("metricName", "PangolinLineage");
                     row1.put("qualvalue", line[1]);
                     row1.put("container", so.getContainer());
+                    if (StringUtils.trimToNull(line[3]) != null)
+                    {
+                        row1.put("comment", line[3]);
+                    }
                     toInsert.add(row1);
 
                     if (StringUtils.trimToNull(line[2]) != null)
@@ -111,11 +117,32 @@ public void processFilesOnWebserver(PipelineJob job, SequenceAnalysisJobSupport
                         row2.put("readset", so.getReadset());
                         row2.put("analysis_id", so.getAnalysis_id());
                         row2.put("category", "Pangolin");
-                        row2.put("metricName", "PangolinLineageConfidence");
+                        row2.put("metricName", "PangolinConflicts");
                         row2.put("value", Double.parseDouble(line[2]));
                         row2.put("container", so.getContainer());
+                        if (StringUtils.trimToNull(line[3]) != null)
+                        {
+                            row2.put("comment", line[3]);
+                        }
                         toInsert.add(row2);
                     }
+
+                    if (StringUtils.trimToNull(line[4]) != null)
+                    {
+                        Map<String, Object> row = new CaseInsensitiveHashMap<>();
+                        row.put("dataid", so.getDataId());
+                        row.put("readset", so.getReadset());
+                        row.put("analysis_id", so.getAnalysis_id());
+                        row.put("category", "Pangolin");
+                        row.put("metricName", "PangolinSummary");
+                        row.put("qualvalue", line[4]);
+                        row.put("container", so.getContainer());
+                        if (StringUtils.trimToNull(line[3]) != null)
+                        {
+                            row.put("comment", line[3]);
+                        }
+                        toInsert.add(row);
+                    }
                 }
             }
             catch (IOException e)
@@ -167,7 +194,29 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
                 for (SequenceOutputFile so : inputFiles)
                 {
                     String[] pangolinData = runPangolin(so.getFile(), ctx.getLogger(), ctx.getFileManager());
-                    writer.writeNext(new String[]{String.valueOf(so.getRowid()), (pangolinData == null ? "QC Fail" : pangolinData[1]), (pangolinData == null ? "" : pangolinData[2])});
+
+                    List<String> versions = new ArrayList<>();
+                    if (pangolinData != null)
+                    {
+                        if (StringUtils.trimToNull(pangolinData[3]) != null)
+                        {
+                            versions.add("Pangolin version: " + pangolinData[3]);
+                        }
+
+                        if (StringUtils.trimToNull(pangolinData[4]) != null)
+                        {
+                            versions.add("pangoLEARN version: " + pangolinData[4]);
+                        }
+
+                        if (StringUtils.trimToNull(pangolinData[5]) != null)
+                        {
+                            versions.add("pango version: " + pangolinData[5]);
+                        }
+                    }
+
+                    String comment = StringUtils.join(versions, ",");
+
+                    writer.writeNext(new String[]{String.valueOf(so.getRowid()), (pangolinData == null ? "QC Fail" : pangolinData[1]), (pangolinData == null ? "" : pangolinData[2]), comment, (pangolinData == null ? "" : pangolinData[7])});
                 }
             }
             catch (IOException e)
@@ -192,6 +241,11 @@ public static void updatePangolinRefs(Logger log) throws PipelineJobException
         wrapper.execute(Arrays.asList("/bin/bash", pangolin.getPath()));
     }
 
+    public static File getRenamedPangolinOutput(File consensusFasta)
+    {
+        return new File(consensusFasta.getParentFile(), FileUtil.getBaseName(consensusFasta) + ".pangolin.csv");
+    }
+
     public static String[] runPangolin(File consensusFasta, Logger log, PipelineOutputTracker tracker) throws PipelineJobException
     {
         SimpleScriptWrapper wrapper = new SimpleScriptWrapper(log);
@@ -211,17 +265,21 @@ public static String[] runPangolin(File consensusFasta, Logger log, PipelineOutp
             throw new PipelineJobException("Pangolin output not found: " + output.getPath());
         }
 
-        tracker.addIntermediateFile(output);
-        try (CSVReader reader = new CSVReader(Readers.getReader(output)))
+        try
         {
-            reader.readNext(); //header
-            String[] line = reader.readNext();
+            File outputMoved = getRenamedPangolinOutput(consensusFasta);
+            FileUtils.moveFile(output, outputMoved);
+            try (CSVReader reader = new CSVReader(Readers.getReader(outputMoved)))
+            {
+                reader.readNext(); //header
+                String[] line = reader.readNext();
 
-            return line;
+                return line;
+            }
         }
         catch (IOException e)
         {
-            throw new PipelineJobException();
+            throw new PipelineJobException(e);
         }
     }
 }