Bugfix to sorting of indel tsv

Author: bbimber

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

@@ -29,6 +29,7 @@
 import org.labkey.api.sequenceanalysis.pipeline.ReferenceGenome;
 import org.labkey.api.sequenceanalysis.pipeline.SequenceAnalysisJobSupport;
 import org.labkey.api.sequenceanalysis.pipeline.SequenceOutputHandler;
+import org.labkey.api.sequenceanalysis.pipeline.SequencePipelineService;
 import org.labkey.api.sequenceanalysis.pipeline.ToolParameterDescriptor;
 import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper;
 import org.labkey.api.util.PageFlowUtil;
@@ -64,6 +65,7 @@
 public class MergeLoFreqVcfHandler extends AbstractParameterizedOutputHandler<SequenceOutputHandler.SequenceOutputProcessor>
 {
     private static final String MIN_AF = "minAfThreshold";
+    private static final String MIN_AF_INDEL = "minIndelAfThreshold";
     private static final String MIN_COVERAGE = "minCoverage";
 
     public MergeLoFreqVcfHandler()
@@ -77,6 +79,11 @@ public MergeLoFreqVcfHandler()
                     put("maxValue", 1);
                     put("decimalPrecision", 2);
                 }}, 0.10),
+                ToolParameterDescriptor.create(MIN_AF, "Min AF to Include (Indels)", "An indel will be included in the indel output if there is a passing variant above this AF in at least one sample.", "ldk-numberfield", new JSONObject(){{
+                    put("minValue", 0);
+                    put("maxValue", 1);
+                    put("decimalPrecision", 2);
+                }}, 0.05),
                 ToolParameterDescriptor.create(MIN_COVERAGE, "Min Coverage To Include", "A site will be reported as ND, unless coverage is above this threshold", "ldk-integerfield", new JSONObject(){{
                     put("minValue", 0);
                 }}, 25),
@@ -221,6 +228,7 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
             basename = basename + (basename.endsWith(".") ? "" : ".");
 
             final double minAfThreshold = ctx.getParams().optDouble(MIN_AF, 0.0);
+            final double minIndelAfThreshold = ctx.getParams().optDouble(MIN_AF_INDEL, 0.05);
             final int minDepth = ctx.getParams().optInt(MIN_COVERAGE, 0);
 
             Set<String> errors = new HashSet<>();
@@ -279,16 +287,6 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
                             double af = vc.getAttributeAsDouble("AF", 0.0);
                             if (af >= minAfThreshold)
                             {
-                                if (vc.isIndel())
-                                {
-                                    uniqueIndels.add(vc.getContig() + "<>" + vc.getStart() + "<>" + vc.getReference().getBaseString() + "<>" + vc.getAlternateAlleles().get(0).getBaseString());
-                                    List<Integer> depths = vc.getAttributeAsIntList("DP4", 0);
-                                    int alleleDepth = depths.get(2) + depths.get(3);
-
-                                    int length = vc.getLengthOnReference() - 1;  //NOTE: indels are left-padded including one ref base
-                                    writer.writeNext(new String[]{ctx.getSequenceSupport().getCachedReadset(so.getReadset()).getName(), String.valueOf(so.getRowid()), String.valueOf(so.getReadset()), "LoFreq", vc.getContig(), String.valueOf(vc.getStart()), String.valueOf(vc.getEnd()), String.valueOf(length), vc.getReference().getBaseString(), vc.getAlternateAlleles().get(0).getBaseString(), vc.getAttributeAsString("GATK_DP", "ND"), vc.getAttributeAsString("DP", "ND"), String.valueOf(alleleDepth), vc.getAttributeAsString("AF", "ND")});
-                                }
-
                                 for (int i = 0; i < vc.getLengthOnReference(); i++)
                                 {
                                     int effectiveStart = vc.getStart() + i;
@@ -302,6 +300,16 @@ public void processFilesRemote(List<SequenceOutputFile> inputFiles, JobContext c
                                     siteToAllele.put(key, site);
                                 }
                             }
+
+                            if (af > minIndelAfThreshold && vc.isIndel())
+                            {
+                                uniqueIndels.add(vc.getContig() + "<>" + vc.getStart() + "<>" + vc.getReference().getBaseString() + "<>" + vc.getAlternateAlleles().get(0).getBaseString());
+                                List<Integer> depths = vc.getAttributeAsIntList("DP4", 0);
+                                int alleleDepth = depths.get(2) + depths.get(3);
+
+                                int length = vc.getLengthOnReference() - 1;  //NOTE: indels are left-padded including one ref base
+                                writer.writeNext(new String[]{ctx.getSequenceSupport().getCachedReadset(so.getReadset()).getName(), String.valueOf(so.getRowid()), String.valueOf(so.getReadset()), "LoFreq", vc.getContig(), String.valueOf(vc.getStart()), String.valueOf(vc.getEnd()), String.valueOf(length), vc.getReference().getBaseString(), vc.getAlternateAlleles().get(0).getBaseString(), vc.getAttributeAsString("GATK_DP", "ND"), vc.getAttributeAsString("DP", "ND"), String.valueOf(alleleDepth), vc.getAttributeAsString("AF", "ND")});
+                            }
                         }
                     }
 
@@ -711,7 +719,7 @@ private void sortTsvFile(JobContext ctx, File input) throws PipelineJobException
                 writer.println("set -e");
                 writer.println("{");
                 writer.println("zcat " + input.getPath() +" | head -n 1;");
-                writer.println("zcat " + input.getPath() +" | tail +2 | sort -k5,5 -k6,6n -k7,7n;");
+                writer.println("zcat " + input.getPath() +" | tail -n +2 | sort -k5,5 -k6,6n -k7,7n;");
                 writer.println("} | gzip -c > " + output.getPath() + "\n");
             }
             catch (IOException e)
@@ -723,6 +731,13 @@ private void sortTsvFile(JobContext ctx, File input) throws PipelineJobException
             wrapper.setWorkingDir(script.getParentFile());
             wrapper.execute(Arrays.asList("/bin/bash", script.getPath()));
 
+            long lineCountIn = SequencePipelineService.get().getLineCount(input);
+            long lineCountOut = SequencePipelineService.get().getLineCount(output);
+            if (lineCountIn != lineCountOut)
+            {
+                throw new PipelineJobException("Error sorting file");
+            }
+
             input.delete();
             script.delete();
 

SequenceAnalysis/src/org/labkey/sequenceanalysis/util/SequenceUtil.java

@@ -113,7 +113,7 @@ public String getPrimaryExtension()
     public static long getLineCount(File f) throws PipelineJobException
     {
         FileType gz = new FileType(".gz");
-        try (InputStream is = gz.isType(f) ? new GZIPInputStream(new FileInputStream(f)) : new FileInputStream(f);BufferedReader reader = new BufferedReader(new InputStreamReader(is));)
+        try (InputStream is = gz.isType(f) ? new GZIPInputStream(new FileInputStream(f)) : new FileInputStream(f);BufferedReader reader = new BufferedReader(new InputStreamReader(is, StringUtilsLabKey.DEFAULT_CHARSET));)
         {
             long i = 0;
             while (reader.readLine() != null)