Repository: chrovis/cljam
Pull request #146 (2018-08-31 06:54:25)

Fix options of cljam.algo.pileup/pileup
#### Summary
- Add an `ignore-overlaps?` option to `cljam.algo.pileup/pileup`
- Skip filtering bases when `min-base-quality` is ~zero~ not positive

#### Tests

- `lein check` 🆗
- `lein test :all` 🆗


Patch:
diff --git a/src/cljam/algo/pileup.clj b/src/cljam/algo/pileup.clj
index 2e9ecac..a808cbb 100644
--- a/src/cljam/algo/pileup.clj
+++ b/src/cljam/algo/pileup.clj
@@ -145,13 +145,18 @@
 
 (defn pileup
   "Piles up alignments in given region and returns a lazy sequence of
-  `cljam.io.pileup.LocusPile`s."
+  `cljam.io.pileup.LocusPile`s.
+
+  The following options are available:
+  - `min-base-quality` Minimum quality of called bases [13]
+  - `min-map-quality` Minimum quality of alignments [0]
+  - `ignore-overlaps?` Disable detecting overlapped bases of PE reads [false]"
   ([sam-reader region]
    (pileup sam-reader region {}))
   ([sam-reader
     {:keys [chr start end] :or {start 1 end Integer/MAX_VALUE}}
-    {:keys [min-base-quality min-map-quality]
-     :or {min-base-quality 13 min-map-quality 0}}]
+    {:keys [min-base-quality min-map-quality ignore-overlaps?]
+     :or {min-base-quality 13 min-map-quality 0 ignore-overlaps? false}}]
    (when-let [len (:len (refs/ref-by-name (sam/read-refs sam-reader) chr))]
      (let [s (max 1 start)
            e (min len end)]
@@ -164,8 +169,12 @@
             (pileup-seq s e)
             (sequence
              (comp (map resolve-bases)
-                   (map correct-overlaps)
-                   (map (filter-by-base-quality min-base-quality))
+                   (if ignore-overlaps?
+                     identity
+                     (map correct-overlaps))
+                   (if (pos? min-base-quality)
+                     (map (filter-by-base-quality min-base-quality))
+                     identity)
                    (keep (partial ->locus-pile chr)))))))))
 
 (defn align-pileup-seqs
@@ -186,17 +195,29 @@
       (step xs))))
 
 (defn mpileup
-  "Pile up alignments from multiple sources."
-  [region & sam-readers]
-  (apply align-pileup-seqs (map #(pileup % region) sam-readers)))
+  "Pile up alignments from multiple sources.
+
+  The following `options` are available:
+  - `min-base-quality` Minimum quality of called bases [13]
+  - `min-map-quality` Minimum quality of alignments [0]
+  - `ignore-overlaps?` Disable detecting overlapped bases of PE reads [false]"
+  [region options & sam-readers]
+  (apply align-pileup-seqs (map #(pileup % region options) sam-readers)))
 
 (defn create-mpileup
-  "Creates a mpileup file from the BAM file."
+  "Creates a mpileup file from the BAM file.
+
+  The following `options` are available:
+  - `min-base-quality` Minimum quality of called bases [13]
+  - `min-map-quality` Minimum quality of alignments [0]
+  - `ignore-overlaps?` Disable detecting overlapped bases of PE reads [false]"
   ([in-sam out-mplp]
    (create-mpileup in-sam nil out-mplp))
   ([in-sam in-ref out-mplp]
    (create-mpileup in-sam in-ref out-mplp nil))
   ([in-sam in-ref out-mplp region]
+   (create-mpileup in-sam in-ref out-mplp region nil))
+  ([in-sam in-ref out-mplp region options]
    (with-open [s (sam/reader in-sam)
                w (plpio/writer out-mplp in-ref)]
      (let [regs (if region
@@ -206,4 +227,4 @@
                      {:chr name :start 1 :end len})
                    (sam/read-refs s)))]
        (doseq [reg regs]
-         (plpio/write-piles w (pileup s reg)))))))
+         (plpio/write-piles w (pileup s reg options)))))))


Test patch:
diff --git a/test/cljam/algo/pileup_test.clj b/test/cljam/algo/pileup_test.clj
index bbc5aa6..f975d15 100644
--- a/test/cljam/algo/pileup_test.clj
+++ b/test/cljam/algo/pileup_test.clj
@@ -220,7 +220,7 @@
 (deftest mpileup
   (with-open [r1 (sam/reader test-sorted-bam-file)
               r2 (sam/reader r1)]
-    (let [mplp (plp/mpileup {:chr "ref"} r1 r2)]
+    (let [mplp (plp/mpileup {:chr "ref"} {} r1 r2)]
       (is (= (keep-indexed #(when (pos? %2) (inc %1)) test-bam-pileup-ref)
              (map first mplp))))))
 
@@ -257,24 +257,35 @@
     {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 3) :end (int 10) :seq "TTGGCCAATT" :qual "AABBCCDDEE" :cigar "8M2S"
      :rnext "=" :pnext (int 3) :tlen (int -8) :mapq (int 20)}]))
 
-(defn- pileup* [region xs]
+(defn- pileup* [region options xs]
   (plp/pileup
    (reify p/IAlignmentReader
      (p/read-refs [_]
-       [{:name "seq1",  :len 249250621}])
+       [{:name "seq1", :len 249250621}])
      (p/read-alignments [_ _]
        xs))
-   region))
+   region
+   options))
 
 (deftest overlap-correction
-  (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 "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))))))
 
 (deftest overlap-qual-correction
   (let [[_ aplp] (#'plp/correct-overlaps
@@ -289,20 +300,31 @@
     (is (= (map :qual aplp) [0 32]))))
 
 (deftest filter-by-base-quality
-  (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 "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))))))
 
 (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} {}))]
     (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]])