Fix handling of idxdepth/paraGRAPH coverage

Author: bbimber

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/alignment/ParagraphStep.java

@@ -1,5 +1,9 @@
 package org.labkey.sequenceanalysis.run.alignment;
 
+import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SamReader;
+import htsjdk.samtools.SamReaderFactory;
+import org.apache.commons.io.FileUtils;
 import org.json.JSONObject;
 import org.labkey.api.module.ModuleLoader;
 import org.labkey.api.pipeline.PipelineJob;
@@ -15,11 +19,14 @@
 import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
 import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper;
 import org.labkey.api.util.FileUtil;
+import org.labkey.api.writer.PrintWriters;
 import org.labkey.sequenceanalysis.SequenceAnalysisModule;
 import org.labkey.sequenceanalysis.util.SequenceUtil;
 
 import java.io.File;
 import java.io.IOException;
+import java.io.PrintWriter;
+import java.nio.charset.Charset;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.List;
@@ -95,9 +102,9 @@ else if (!svVcf.exists())
                 depthArgs.add("-b");
                 depthArgs.add(so.getFile().getPath());
 
-                File coverageFile = new File(ctx.getWorkingDirectory(), "coverage.txt");
+                File coverageJson = new File(ctx.getWorkingDirectory(), "coverage.json");
                 depthArgs.add("-o");
-                depthArgs.add(coverageFile.getPath());
+                depthArgs.add(coverageJson.getPath());
 
                 depthArgs.add("-r");
                 depthArgs.add(ctx.getSequenceSupport().getCachedGenome(so.getLibrary_id()).getWorkingFastaFile().getPath());
@@ -110,14 +117,41 @@ else if (!svVcf.exists())
 
                 new SimpleScriptWrapper(ctx.getLogger()).execute(depthArgs);
 
-                if (!coverageFile.exists())
+                if (!coverageJson.exists())
                 {
-                    throw new PipelineJobException("Missing file: " + coverageFile.getPath());
+                    throw new PipelineJobException("Missing file: " + coverageJson.getPath());
                 }
+                ctx.getFileManager().addIntermediateFile(coverageJson);
 
                 // Should produce a simple text file:
                 //    id  path    depth   read length
                 //    TNPRC-IB18  ../IB18.cram 29.77   150
+                File coverageFile = new File(ctx.getWorkingDirectory(), "coverage.txt");
+                try (PrintWriter writer = PrintWriters.getPrintWriter(coverageFile); SamReader reader = SamReaderFactory.makeDefault().open(so.getFile()))
+                {
+                    SAMFileHeader header = reader.getFileHeader();
+                    if (header.getReadGroups().size() == 0)
+                    {
+                        throw new PipelineJobException("No read groups found in input BAM");
+                    }
+                    else if (header.getReadGroups().size() > 1)
+                    {
+                        throw new PipelineJobException("More than one read group found in BAM");
+                    }
+
+                    String rgId = header.getReadGroups().get(0).getSample();
+
+                    JSONObject json = new JSONObject(FileUtils.readFileToString(coverageFile, Charset.defaultCharset()));
+                    writer.println("id\tpath\tdepth\tread length");
+                    double depth = json.getJSONObject("autosome").getDouble("depth");
+                    double readLength = json.getInt("read_length");
+                    writer.println(rgId + "\t" + so.getFile().getPath() + "\t" + depth + "\t" + readLength);
+                }
+                catch (IOException e)
+                {
+                    throw new PipelineJobException(e);
+                }
+                ctx.getFileManager().addIntermediateFile(coverageFile);
 
                 List<String> paragraphArgs = new ArrayList<>();
                 paragraphArgs.add(AbstractCommandWrapper.resolveFileInPath("multigrmpy.py", null, true).getPath());