qual range warnings added

broadinstitute · Jan 31, 2024 · 9db48a1 · 9db48a1
1 parent 5d10053
commit 9db48a1
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 1 deletion.
diff --git a/src/main/java/picard/analysis/CollectQualityYieldMetricsFlow.java b/src/main/java/picard/analysis/CollectQualityYieldMetricsFlow.java
@@ -30,6 +30,7 @@
 import htsjdk.samtools.reference.ReferenceSequence;
 import htsjdk.samtools.util.Histogram;
 import htsjdk.samtools.util.IOUtil;
+import htsjdk.samtools.util.Log;
 import org.broadinstitute.barclay.argparser.Argument;
 import org.broadinstitute.barclay.argparser.ArgumentCollection;
 import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
@@ -60,6 +61,11 @@
 )
 @ExperimentalFeature
 public class CollectQualityYieldMetricsFlow extends SinglePassSamProgram {
+
+    private static final byte MIN_QUAL = 0;
+    private static final byte MAX_QUAL = 100;
+    private final Log log = Log.getInstance(CollectQualityYieldMetricsFlow.class);
+
     private static final int CYCLE_SIZE = 4;
     private QualityYieldMetricsCollectorFlow collector = null;
     public Histogram<Integer> qualityHistogram = new Histogram<>("KEY", "QUAL_COUNT");
@@ -216,7 +222,19 @@ private byte[] getFlowQualities(FlowBasedRead fread) {
             double[] errorProbs = computeErrorProb(fread);
             byte[] quals = new byte[errorProbs.length];
             for ( int i = 0 ; i < errorProbs.length ; i++ ) {
-                quals[i] = (byte)Math.round(-10 * Math.log10(errorProbs[i]));
+                if ( errorProbs[i] == 0.0 ) {
+                    // this is a special case that should not happen
+                    log.warn(fread.getReadName() + ": zero errorProb on flow: " + i);
+                    quals[i] = MAX_QUAL;
+                } else {
+                    long q = Math.round(-10 * Math.log10(errorProbs[i]));
+                    if ( q < MIN_QUAL || q > MAX_QUAL ) {
+                        // this is an out-of-range condition. should not happen as well
+                        log.warn(fread.getReadName() + ": qual " + q + " is out of range on flow: " + i);
+                        q = Math.max(MIN_QUAL, Math.min(MAX_QUAL, q));
+                    }
+                    quals[i] = (byte)q;
+                }
             }
             return quals;
         }

diff --git a/src/main/java/picard/flow/FlowBasedRead.java b/src/main/java/picard/flow/FlowBasedRead.java
@@ -318,6 +318,10 @@ private void clipProbs() {
             }
         }
     }
+
+    public String getReadName() {
+        return samRecord.getReadName();
+    }
 }