A dirty workaround for some rare blast "stochastic" alignments

mikessh · Sep 16, 2015 · d5170d1 · d5170d1
1 parent b5e2ca9
commit d5170d1
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 11 deletions.
diff --git a/src/main/groovy/com/antigenomics/higblast/clonotype/Clonotype.groovy b/src/main/groovy/com/antigenomics/higblast/clonotype/Clonotype.groovy
@@ -43,7 +43,6 @@ class Clonotype implements Comparable<Clonotype> {
     final List<Mutation> mutations
     final long count
     final double freq
-    final byte minQual
     final byte[] cdrInsertQual, mutationQual
     final boolean hasCdr3, complete, inFrame, noStop
     final Cdr3Markup cdr3Markup
@@ -60,20 +59,16 @@ class Clonotype implements Comparable<Clonotype> {
         this.count = data.count.longValue()
         this.freq = (double) count / total
         this.cdrInsertQual = new byte[data.cdrInsertQual.length()]
-        def minQual = Util.MAX_QUAL
+        def cdrQualCount = data.cdrQualCount.longValue()
         (0..<cdrInsertQual.length).each {
-            def qual = (byte) ((double) data.cdrInsertQual.get(it) / count)
+            def qual = (byte) ((double) data.cdrInsertQual.get(it) / cdrQualCount)
             cdrInsertQual[it] = qual
-            minQual = Math.min(qual, minQual)
         }
         this.mutationQual = new byte[data.mutationQual.length()]
         (0..<mutationQual.length).each {
             def qual = (byte) ((double) data.mutationQual.get(it) / count)
             mutationQual[it] = qual
-            minQual = Math.min(qual, minQual)
         }
-        this.minQual = minQual
-
         def representativeMapping = key.representativeMapping
 
         this.cdr3Markup = representativeMapping.mapping.cdr3Markup

diff --git a/src/main/groovy/com/antigenomics/higblast/clonotype/ClonotypeData.groovy b/src/main/groovy/com/antigenomics/higblast/clonotype/ClonotypeData.groovy
@@ -35,11 +35,12 @@ import java.util.concurrent.atomic.AtomicLong
 import java.util.concurrent.atomic.AtomicLongArray
 
 class ClonotypeData {
-    final AtomicLong count
+    final AtomicLong count, cdrQualCount
     final AtomicLongArray mutationQual, cdrInsertQual
 
     ClonotypeData(ReadMapping readMapping) {
         this.count = new AtomicLong(0)
+        this.cdrQualCount = new AtomicLong(0)
         this.mutationQual = new AtomicLongArray(readMapping.mutationQual.length)
         this.cdrInsertQual = new AtomicLongArray(readMapping.cdrInsertQual.length)
 
@@ -49,6 +50,11 @@ class ClonotypeData {
     void update(ReadMapping readMapping) {
         count.incrementAndGet()
         readMapping.mutationQual.eachWithIndex { it, i -> mutationQual.addAndGet(i, it) }
-        readMapping.cdrInsertQual.eachWithIndex { it, i -> cdrInsertQual.addAndGet(i, it) }
+        if (readMapping.cdrInsertQual.length == cdrInsertQual.length()) {
+            // protect against very rare cases of ambiguous D alignment
+            // not a big deal as this quality is used just for display
+            readMapping.cdrInsertQual.eachWithIndex { it, i -> cdrInsertQual.addAndGet(i, it) }
+            cdrQualCount.incrementAndGet()
+        }
     }
 }
diff --git a/src/main/groovy/com/antigenomics/higblast/clonotype/ClonotypeKey.groovy b/src/main/groovy/com/antigenomics/higblast/clonotype/ClonotypeKey.groovy
@@ -68,7 +68,6 @@ class ClonotypeKey {
 
         cdr3nt == that.cdr3nt &&
                 vSegment == that.vSegment &&
-                dSegment == that.dSegment &&
                 jSegment == that.jSegment &&
                 mutations == that.mutations
     }
@@ -78,7 +77,6 @@ class ClonotypeKey {
         int result
         result = cdr3nt.hashCode()
         result = 31 * result + vSegment.hashCode()
-        result = 31 * result + dSegment.hashCode()
         result = 31 * result + jSegment.hashCode()
         31 * result + mutations.hashCode()
     }