Repository: chrovis/cljam
Pull request #148 (2018-09-04 02:12:28)

Fix a bug in `cljam.algo.pileup/pileup`
#### Summary
This PR fixes a bug in `cljam.algo.pileup/pileup`.

#### Problem

IllegalArgumentException in `cljam.algo.pileup/pileup`
```
java.lang.IllegalArgumentException: Value out of range for byte: 186
```

#### Cause
The `qual` field in `cljam.io.pileup.PileupBase` is declared as a `byte` because a phred-33 encoded base quality score never exceeds `93`:

https://github.com/chrovis/cljam/blob/02addc40def7515d2f30db8032fe4cf9776c5bd8/src/cljam/io/pileup.clj#L15

The record is also constructed in `cljam.algo.pileup`:

https://github.com/chrovis/cljam/blob/02addc40def7515d2f30db8032fe4cf9776c5bd8/src/cljam/algo/pileup.clj#L85-L94

But quality scores can be `128` ~ `200` while piling alignments up with the option `:ignore-overlaps false`:

https://github.com/chrovis/cljam/blob/02addc40def7515d2f30db8032fe4cf9776c5bd8/src/cljam/algo/pileup.clj#L117-L120

#### Changes
- Use `short` primitive type hinting instead of `byte` in `cljam.io.pileup.PileupBase`
- Fix broken tests related to `pileup`

#### Tests

- `lein check` 🆗
- `lein test :all` 🆗


Patch:
diff --git a/src/cljam/io/pileup.clj b/src/cljam/io/pileup.clj
index 4c17db5..42531de 100644
--- a/src/cljam/io/pileup.clj
+++ b/src/cljam/io/pileup.clj
@@ -12,7 +12,7 @@
 (defrecord PileupBase [^boolean start?
                        mapq ;; byte?
                        ^char base
-                       ^byte qual
+                       ^short qual
                        ^boolean reverse?
                        ^boolean end?
                        insertion ;; String?


Test patch:
diff --git a/test/cljam/algo/pileup_test.clj b/test/cljam/algo/pileup_test.clj
index f975d15..269b31f 100644
--- a/test/cljam/algo/pileup_test.clj
+++ b/test/cljam/algo/pileup_test.clj
@@ -252,9 +252,9 @@
 (def ^:private reads-for-pileup
   (mapv
    p/map->SAMAlignment
-   [{:qname "R001" :flag (int 99) :rname "seq1" :pos (int 3) :end (int 10) :seq "AATTGGCCAA" :qual "AABBCCDDEE" :cigar "2S8M"
+   [{:qname "R001" :flag (int 99) :rname "seq1" :pos (int 3) :end (int 10) :seq "AATTGGCCAA" :qual "AABBCCDDz~" :cigar "2S8M"
      :rnext "=" :pnext (int 3) :tlen (int 8) :mapq (int 60)}
-    {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 3) :end (int 10) :seq "TTGGCCAATT" :qual "AABBCCDDEE" :cigar "8M2S"
+    {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 3) :end (int 10) :seq "TTGGCCAATT" :qual "AABBCCDDz~" :cigar "8M2S"
      :rnext "=" :pnext (int 3) :tlen (int -8) :mapq (int 20)}]))
 
 (defn- pileup* [region options xs]
@@ -267,69 +267,56 @@
    region
    options))
 
+(deftest pileup-default
+  (let [plps (->> reads-for-pileup
+                  (pileup* {:chr "seq1" :start 1 :end 10} {}))]
+    (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
+           (map (comp count :pile) plps)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\A] [\A]])
+           (map #(map :base (:pile %)) plps)))
+    (is (= (filter seq [[] [] [65] [65] [67] [67] [69] [69] [124] [128]])
+           (map #(map :qual (:pile %)) plps)))))
+
 (deftest overlap-correction
-  (testing "defalut"
-    (let [plps (->> reads-for-pileup
-                    (pileup* {:chr "seq1" :start 1 :end 10} {}))]
-      (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
-             (map (comp count :pile) plps)))
-      (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\A] [\A]])
-             (map #(map :base (:pile %)) plps)))
-      (is (= (filter seq [[] [] [65] [65] [67] [67] [69] [69] [71] [71]])
-             (map #(map :qual (:pile %)) plps)))))
-  (testing "ignoring overlaps"
-    (let [plps (->> reads-for-pileup
-                    (pileup* {:chr "seq1" :start 1 :end 10} {:ignore-overlaps? true}))]
-      (is (= (filter pos? [0 0 2 2 2 2 2 2 2 2])
-             (map (comp count :pile) plps)))
-      (is (= (filter seq [[] [] [\T \T] [\T \T] [\G \G] [\G \G] [\C \C] [\C \C] [\A \A] [\A \A]])
-             (map #(map :base (:pile %)) plps)))
-      (is (= (filter seq [[] [] [33 32] [33 32] [34 33] [34 33] [35 34] [35 34] [36 35] [36 35]])
-             (map #(map :qual (:pile %)) plps))))))
+  (let [plps (->> reads-for-pileup
+                  (pileup* {:chr "seq1" :start 1 :end 10} {:ignore-overlaps? true}))]
+    (is (= (filter pos? [0 0 2 2 2 2 2 2 2 2])
+           (map (comp count :pile) plps)))
+    (is (= (filter seq [[] [] [\T \T] [\T \T] [\G \G] [\G \G] [\C \C] [\C \C] [\A \A] [\A \A]])
+           (map #(map :base (:pile %)) plps)))
+    (is (= (filter seq [[] [] [33 32] [33 32] [34 33] [34 33] [35 34] [35 34] [89 35] [93 35]])
+           (map #(map :qual (:pile %)) plps)))))
 
 (deftest overlap-qual-correction
   (let [[_ aplp] (#'plp/correct-overlaps
-                  [nil [(->pbase {:base \A :qual (byte 40) :alignment (->aln {:qname "R001" :flag 99})})
-                        (->pbase {:base \T :qual (byte 32) :alignment (->aln {:qname "R001" :flag 147})})]])]
+                  [nil [(->pbase {:base \A :qual (short 40) :alignment (->aln {:qname "R001" :flag 99})})
+                        (->pbase {:base \T :qual (short 32) :alignment (->aln {:qname "R001" :flag 147})})]])]
     (is (= (map :base aplp) [\A \T]))
     (is (= (map :qual aplp) [32 0])))
   (let [[_ aplp] (#'plp/correct-overlaps
-                  [nil [(->pbase {:base \T :qual (byte 32) :alignment (->aln {:qname "R001" :flag 147})})
-                        (->pbase {:base \A :qual (byte 40) :alignment (->aln {:qname "R001" :flag 99})})]])]
+                  [nil [(->pbase {:base \T :qual (short 32) :alignment (->aln {:qname "R001" :flag 147})})
+                        (->pbase {:base \A :qual (short 40) :alignment (->aln {:qname "R001" :flag 99})})]])]
     (is (= (map :base aplp) [\T \A]))
     (is (= (map :qual aplp) [0 32]))))
 
 (deftest filter-by-base-quality
-  (testing "default"
-    (let [plps (->> reads-for-pileup
-                    (pileup* {:chr "seq1" :start 1 :end 10} {})
-                    (map (plp/filter-by-base-quality 1)))]
-      (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
-             (map (comp count :pile) plps)))
-      (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\A] [\A]])
-             (map #(map :base (:pile %)) plps)))
-      (is (= (filter seq [[] [] [65] [65] [67] [67] [69] [69] [71] [71]])
-             (map #(map :qual (:pile %)) plps)))))
-  (testing "without filtering"
-    (let [plps (->> reads-for-pileup
-                    (pileup* {:chr "seq1" :start 1 :end 10} {:min-base-quality 0})
-                    (map (plp/filter-by-base-quality 1)))]
-      (is (= (filter pos? [0 0 2 2 2 2 2 2 2 2])
-             (map (comp count :pile) plps)))
-      (is (= (filter seq [[] [] [\T \T] [\T \T] [\G \G] [\G \G] [\C \C] [\C \C] [\A \A] [\A \A]])
-             (map #(map :base (:pile %)) plps)))
-      (is (= (filter seq [[] [] [65 0] [65 0] [67 0] [67 0] [69 0] [69 0] [71 0] [71 0]])
-             (map #(map :qual (:pile %)) plps))))))
+  (let [plps (->> reads-for-pileup
+                  (pileup* {:chr "seq1" :start 1 :end 10} {:min-base-quality 0}))]
+    (is (= (filter pos? [0 0 2 2 2 2 2 2 2 2])
+           (map (comp count :pile) plps)))
+    (is (= (filter seq [[] [] [\T \T] [\T \T] [\G \G] [\G \G] [\C \C] [\C \C] [\A \A] [\A \A]])
+           (map #(map :base (:pile %)) plps)))
+    (is (= (filter seq [[] [] [65 0] [65 0] [67 0] [67 0] [69 0] [69 0] [124 0] [128 0]])
+           (map #(map :qual (:pile %)) plps)))))
 
 (deftest filter-by-map-quality
   (let [plps (->> reads-for-pileup
-                  (filter (fn [a] (<= 30 (:mapq a))))
-                  (pileup* {:chr "seq1" :start 1 :end 10} {}))]
+                  (pileup* {:chr "seq1" :start 1 :end 10} {:min-map-quality 30}))]
     (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
            (map (comp count :pile) plps)))
     (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\A] [\A]])
            (map #(map :base (:pile %)) plps)))
-    (is (= (filter seq [[] [] [33] [33] [34] [34] [35] [35] [36] [36]])
+    (is (= (filter seq [[] [] [33] [33] [34] [34] [35] [35] [89] [93]])
            (map #(map :qual (:pile %)) plps)))))
 
 (deftest about-create-mpileup