Repository: chrovis/cljam
Pull request #160 (2019-04-10 09:32:23)

Improve pileup by chunking alignments
#### Summary
Pileup alignments per chunk to improve the performance

#### Benchmark
##### Pileup many alignments in a short region
```clojure
(time (with-open [r (sam/reader "/path/to/bam/SRP006438.chr1.bam")]
  (dorun (pileup/pileup r {:chr "chr1", :start 152210000, :end 152214999}))))
```

- master: 8523 ms
- PR: 6013 ms (-29.4%)

##### Pileup-bench
- master: `time: 56.355322 sec, sd: 4.858070 ms`
- PR: `time: 55.349883 sec, sd: 7.576505 ms` (-1.8%)

#### Changes
- Read and make chunks which have all pairs
- Correct `quals-at-ref` of an overlapping pair by looking up `quals` in the corrected quals map
- Pileup alignments in each chunk
- Add pileup tests

#### Tests
- `lein check` 🆗
- `lein test :all` 🆗
- `lein eastwood` 🆗

codecov[bot]: # [Codecov](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=h1) Report
> Merging [#160](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=desc) into [master](https://codecov.io/gh/chrovis/cljam/commit/773eb527a09b2d206ee567c8a93eeee27162f065?src=pr&el=desc) will **increase** coverage by `0.06%`.
> The diff coverage is `93.27%`.

[![Impacted file tree graph](https://codecov.io/gh/chrovis/cljam/pull/160/graphs/tree.svg?width=650&token=hrCuiz6NWa&height=150&src=pr)](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #160      +/-   ##
==========================================
+ Coverage    85.9%   85.96%   +0.06%     
==========================================
  Files          71       71              
  Lines        5305     5407     +102     
  Branches      485      488       +3     
==========================================
+ Hits         4557     4648      +91     
- Misses        263      271       +8     
- Partials      485      488       +3
```


| [Impacted Files](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [src/cljam/algo/pileup.clj](https://codecov.io/gh/chrovis/cljam/pull/160/diff?src=pr&el=tree#diff-c3JjL2NsamFtL2FsZ28vcGlsZXVwLmNsag==) | `93.11% <93.27%> (-2.75%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=footer). Last update [773eb52...c4f0b63](https://codecov.io/gh/chrovis/cljam/pull/160?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

yito88: @athos @alumi Thank you for trying and reviews!
I'll investigate them.
yito88: The cause of performance degradation for the longer region was duplicating chunks.
When concatenating chunks for overlapping pairs, it repeatedly piled up for multiple regions including these pairs.

I've tried benchmarking after fixing that.

sample: WES `SRP006438.sorted.bam`

- gene: NRAS `chr1:114703000-114723000`

  |branch|time|
  |:-:|:-:|
  |`master`|`time: 2.424547 sec, sd: 14.857820 µs`|
  |`feature/pileup-chunks`|`time: 2.217902 sec, sd: 129.108257 µs`|

- gene: EGFR `chr7:55015000-55215000`

  |branch|time|
  |:-:|:-:|
  |`master`|`time: 5.902383 sec, sd: 30.402538 µs`|
  |`feature/pileup-chunks`|`time: 5.497433 sec, sd: 109.532839 µs`|
alumi: @yito88 Would you squash your commits down before merging this PR? 🙏 

Patch:
diff --git a/src/cljam/algo/pileup.clj b/src/cljam/algo/pileup.clj
index a54dffd..2b880fc 100644
--- a/src/cljam/algo/pileup.clj
+++ b/src/cljam/algo/pileup.clj
@@ -9,28 +9,31 @@
   (:import [cljam.io.protocols SAMAlignment]
            [cljam.io.pileup PileupBase LocusPile]))
 
-(defn- cover-locus? [^long pos ^SAMAlignment aln]
-  (<= (.pos aln) pos))
+(defn- seq-step
+  [^long start ^long end ^long step-size alns]
+  (letfn [(starts-before? [^long pos ^SAMAlignment aln]
+            (<= (.pos aln) pos))
+          (ends-after? [^long pos ^SAMAlignment aln]
+            (<= pos (.end aln)))
+          (step [^long s carried-over rest-alns]
+            (lazy-seq
+             (when (<= s end)
+               (let [e (Math/min (+ s step-size) end)
+                     [in out] (split-with #(starts-before? e %) rest-alns)
+                     pile (into carried-over in)
+                     carry (filterv #(ends-after? (inc e) %) pile)
+                     out-pos (some-> out ^SAMAlignment (first) .pos)
+                     next-s (if-not (seq carry) (or out-pos (inc end)) (inc e))]
+                 (cons [s pile] (step next-s carry out))))))]
+    (let [[x :as xs] (drop-while #(< (.end ^SAMAlignment %) start) alns)]
+      (when x
+        (step (max start (.pos ^SAMAlignment x)) [] xs)))))
 
 (defn pileup-seq
   "Returns a lazy sequence that each element contains [position (alignments
   piled-up at the locus)]."
   [^long start ^long end alns]
-  (letfn [(step [^long pos prev-buf rest-alns]
-            (lazy-seq
-             (when (< (dec pos) end)
-               (let [[i rests] (split-with (partial cover-locus? pos) rest-alns)
-                     p (some-> rests ^SAMAlignment first .pos)
-                     [b :as buf] (into (filterv #(<= pos (.end ^SAMAlignment %))
-                                                prev-buf) i)
-                     next-pos (if (and p (not b)) p (inc pos))]
-                 (if b
-                   (cons [pos buf] (step next-pos buf rests))
-                   (when (and p (<= p end))
-                     (step next-pos buf rests)))))))]
-    (let [[x :as xs] (drop-while #(< (.end ^SAMAlignment %) start) alns)]
-      (when x
-        (step (max start (.pos ^SAMAlignment x)) [] xs)))))
+  (seq-step start end 0 alns))
 
 (defn- quals-at-ref
   [idx ^String qual]
@@ -104,35 +107,6 @@
   (when (seq pile)
     (LocusPile. chr pos pile)))
 
-(defn- correct-qual
-  "Correct quality of two overlapped mate reads by setting zero quality for one
-  of the base."
-  [^PileupBase r1 ^PileupBase r2]
-  (let [b1 (.base r1)
-        q1 (.qual r1)
-        b2 (.base r2)
-        q2 (.qual r2)]
-    (if (= b1 b2)
-      [(assoc r1 :qual (min 200 (+ q1 q2))) (assoc r2 :qual 0)]
-      (if (<= q2 q1)
-        [(assoc r1 :qual (int (* 0.8 q1))) (assoc r2 :qual 0)]
-        [(assoc r1 :qual 0) (assoc r2 :qual (int (* 0.8 q2)))]))))
-
-(defn correct-overlapped-bases
-  "Find out overlapped bases and tweak their base quality scores."
-  [xs]
-  (if (<= (count xs) 1)
-    xs
-    (->> xs
-         (group-by (fn [x] (.qname ^PileupBase x)))
-         (into [] (mapcat
-                   (fn [[_ xs]]
-                     (if (<= (count xs) 1) xs (apply correct-qual xs))))))))
-
-(defn- correct-overlaps
-  [pile]
-  (update pile 1 correct-overlapped-bases))
-
 (defn filter-by-base-quality
   "Returns a predicate for filtering piled-up reads by base quality at its
   position."
@@ -142,6 +116,74 @@
          (partial filterv)
          (update p 1))))
 
+(defn- merge-corrected-quals
+  "Merge corrected quals with the uncorrected part."
+  [^SAMAlignment aln ^long correct-start corrected-quals]
+  (let [quals (:quals-at-ref aln)
+        start (.pos aln)
+        len   (count corrected-quals)]
+    (if (< start correct-start)
+      (-> quals
+          (subvec 0 (- correct-start start))
+          (into corrected-quals)
+          (into (subvec quals (+ (- correct-start start) len))))
+      (into corrected-quals (subvec quals len)))))
+
+(defn- correct-pair-quals
+  "Correct quals of a pair. Returns a map with corrected quals."
+  [^SAMAlignment a1 ^SAMAlignment a2]
+  (when (and (pos? (.pnext a1))
+             (<= (max (.pos a1) (.pos a2)) (.end a1)))
+    (let [tlen1 (.tlen a1)
+          tlen2 (.tlen a2)
+          quals1 (:quals-at-ref a1)
+          quals2 (:quals-at-ref a2)
+          seqs1  (:seqs-at-ref a1)
+          seqs2  (:seqs-at-ref a2)
+          correct-start (max (.pos a1) (.pos a2))
+          new-quals (for [pos (range correct-start
+                                     (inc (min (.end a1) (.end a2))))]
+                      (let [relative-pos1 (- pos (.pos a1))
+                            relative-pos2 (- pos (.pos a2))
+                            q1 (quals1 relative-pos1)
+                            q2 (quals2 relative-pos2)
+                            [b1] (seqs1 relative-pos1)
+                            [b2] (seqs2 relative-pos2)]
+                        (if (= b1 b2)
+                          [(min 200 (+ q1 q2)) 0]
+                          (if (<= q2 q1)
+                            [(int (* 0.8 q1)) 0]
+                            [0 (int (* 0.8 q2))]))))
+          [new1 new2] (apply map vector new-quals)
+          new-quals1 (merge-corrected-quals a1 correct-start new1)
+          new-quals2 (merge-corrected-quals a2 correct-start new2)]
+      (if (flag/r1? (.flag a1))
+        [new-quals1 new-quals2]
+        [new-quals2 new-quals1]))))
+
+(defn- make-corrected-quals-map
+  "Make a map which has corrected quals of all overlapping pairs."
+  [alns]
+  (->> alns
+       (group-by (fn [x] (.qname ^SAMAlignment x)))
+       (into {} (keep
+                 (fn [[qname xs]]
+                   (when (<= 2 (count xs))
+                     [qname (apply correct-pair-quals xs)]))))))
+
+(defn- correct-quals-at-ref
+  "Returns an alignment with corrected quals by looking up quals
+  in the corrected quals map."
+  [corrected-map ^SAMAlignment aln]
+  (if-not (= (.rnext aln) "=")
+    aln
+    (if-let [quals-pair (corrected-map (.qname aln))]
+      (let [new-quals (if (flag/r1? (.flag aln))
+                        (first quals-pair)
+                        (second quals-pair))]
+        (assoc aln :quals-at-ref new-quals))
+      aln)))
+
 (defn pileup
   "Piles up alignments in given region and returns a lazy sequence of
   `cljam.io.pileup.LocusPile`s.
@@ -149,32 +191,38 @@
   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]"
+  - `ignore-overlaps?` Disable detecting overlapped bases of PE reads [false]
+  - `chunk-size` Size of a chunk to pile up at once [5000]"
   ([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 ignore-overlaps?]
-     :or {min-base-quality 13 min-map-quality 0 ignore-overlaps? false}}]
+    {:keys [min-base-quality min-map-quality ignore-overlaps? chunk-size]
+     :or {min-base-quality 13 min-map-quality 0 ignore-overlaps? false
+          chunk-size 5000}}]
    (when-let [len (:len (refs/ref-by-name (sam/read-refs sam-reader) chr))]
      (let [s (max 1 start)
-           e (min len end)]
-       (->> {:chr chr :start s :end e}
-            (sam/read-alignments sam-reader)
+           e (min len end)
+           region {:chr chr :start s :end e}]
+       (->> (sam/read-alignments sam-reader region)
             (sequence
              (comp
               (filter (basic-mpileup-pred min-map-quality))
               (map index-cigar)))
-            (pileup-seq s e)
+            (seq-step start end chunk-size)
             (sequence
-             (comp (map resolve-bases)
-                   (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)))))))))
+             (comp
+              (mapcat (fn [[pos alns]]
+                        (->> (if ignore-overlaps?
+                               alns
+                               (keep (partial correct-quals-at-ref
+                                              (make-corrected-quals-map alns)) alns))
+                             (pileup-seq pos (min end (+ pos chunk-size))))))
+              (map resolve-bases)
+              (if (pos? min-base-quality)
+                (map (filter-by-base-quality min-base-quality))
+                identity)
+              (keep (partial ->locus-pile chr)))))))))
 
 (defn align-pileup-seqs
   "Align multiple piled-up seqs."


Test patch:
diff --git a/test/cljam/algo/pileup_test.clj b/test/cljam/algo/pileup_test.clj
index 269b31f..25d5086 100644
--- a/test/cljam/algo/pileup_test.clj
+++ b/test/cljam/algo/pileup_test.clj
@@ -198,6 +198,31 @@
                                   [{:pos 3 :pile [{:pos 3 :end 3}]}]
                                   [])))))
 
+(deftest chunk-step-seq
+  (testing "empty"
+    (is (nil? (seq (#'plp/seq-step 1 10 2 []))))
+    (is (nil? (seq (#'plp/seq-step 2 10 2 (map ->aln [{:pos 1 :end 1}])))))
+    (is (nil? (seq (#'plp/seq-step 1 10 2 (map ->aln [{:pos 11 :end 11}]))))))
+  (testing "dense"
+    (are [?in ?out]
+        (= ?out (mapv ->pile (#'plp/seq-step 1 10 2 (mapv ->aln (map #(zipmap [:pos :end] %) ?in)))))
+      [[1 1]] [[1 [[1 1]]]]
+      [[3 5]] [[3 [[3 5]]]]
+      [[3 5] [4 4]] [[3 [[3 5] [4 4]]]]
+      [[3 5] [4 6]] [[3 [[3 5] [4 6]]] [6 [[4 6]]]]))
+  (testing "sparse"
+    (are [?in ?out]
+        (= ?out (mapv ->pile (#'plp/seq-step 1 10 2 (mapv ->aln (map #(zipmap [:pos :end] %) ?in)))))
+      [[3 5] [14 16]] [[3 [[3 5]]]]
+      [[3 5] [6 7]] [[3 [[3 5]]] [6 [[6 7]]]]
+      [[3 5] [4 5] [6 7]] [[3 [[3 5] [4 5]]] [6 [[6 7]]]]
+      [[3 5] [7 8] [10 11]] [[3 [[3 5]]] [7 [[7 8]]] [10 [[10 11]]]]))
+  (testing "wide"
+    (are [?in ?out]
+        (= ?out (mapv ->pile (#'plp/seq-step 1 10 20 (mapv ->aln (map #(zipmap [:pos :end] %) ?in)))))
+      [[1 1] [14 16]] [[1 [[1 1]]]]
+      [[3 5] [4 6] [20 21]] [[3 [[3 5] [4 6]]]])))
+
 (deftest about-pileup
   (testing "dense"
     (with-open [br (sam/bam-reader test-sorted-bam-file)
@@ -249,6 +274,9 @@
       (is (= (filter seq (take 40 test-bam-mpileup-seq-ref2-freq))
              (map #(frequencies (map :base (:pile %))) plp-ref2))))))
 
+;;   --------
+;;   --------
+;; 1234567890
 (def ^:private reads-for-pileup
   (mapv
    p/map->SAMAlignment
@@ -257,6 +285,50 @@
     {: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)}]))
 
+;;       ----
+;;   ---
+;; 1234567890
+(def ^:private reads-gap-for-pileup
+  (mapv
+   p/map->SAMAlignment
+   [{:qname "R001" :flag (int 99) :rname "seq1" :pos (int 3) :end (int 5) :seq "AATTG" :qual "AABz~" :cigar "2S3M"
+     :rnext "=" :pnext (int 7) :tlen (int 5) :mapq (int 60)}
+    {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 7) :end (int 10) :seq "TTGGTT" :qual "CCDDz~" :cigar "4M2S"
+     :rnext "=" :pnext (int 3) :tlen (int -5) :mapq (int 20)}]))
+
+;;       ----
+;;   ----
+;; 1234567890
+(def ^:private reads-continal-for-pileup
+  (mapv
+   p/map->SAMAlignment
+   [{:qname "R001" :flag (int 99) :rname "seq1" :pos (int 3) :end (int 6) :seq "AATTGG" :qual "AABBz~" :cigar "2S4M"
+     :rnext "=" :pnext (int 6) :tlen (int 5) :mapq (int 60)}
+    {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 7) :end (int 10) :seq "TTGGTT" :qual "CCDDz~" :cigar "4M2S"
+     :rnext "=" :pnext (int 3) :tlen (int -5) :mapq (int 20)}]))
+
+;;        --------
+;;   --------
+;; 123456789012345
+(def ^:private reads-overlapping-for-pileup
+  (mapv
+   p/map->SAMAlignment
+   [{:qname "R001" :flag (int 147) :rname "seq1" :pos (int 3) :end (int 10) :seq "AATTGGCCAA" :qual "AABBCCDDz~" :cigar "2S8M"
+     :rnext "=" :pnext (int 7) :tlen (int 6) :mapq (int 60)}
+    {:qname "R001" :flag (int 99) :rname "seq1" :pos (int 8) :end (int 15) :seq "CCCTTGGATT" :qual "z~AABBCCDD" :cigar "8M2S"
+     :rnext "=" :pnext (int 3) :tlen (int -6) :mapq (int 20)}]))
+
+;;     ----
+;;   --------
+;; 1234567890
+(def ^:private reads-including-for-pileup
+  (mapv
+   p/map->SAMAlignment
+   [{:qname "R001" :flag (int 99) :rname "seq1" :pos (int 3) :end (int 10) :seq "AATTGGCCAA" :qual "AABBCCDDz~" :cigar "2S8M"
+     :rnext "=" :pnext (int 5) :tlen (int 8) :mapq (int 60)}
+    {:qname "R001" :flag (int 147) :rname "seq1" :pos (int 5) :end (int 8) :seq "GGACTT" :qual "AABBz~" :cigar "4M2S"
+     :rnext "=" :pnext (int 3) :tlen (int -8) :mapq (int 20)}]))
+
 (defn- pileup* [region options xs]
   (plp/pileup
    (reify p/IAlignmentReader
@@ -269,13 +341,58 @@
 
 (deftest pileup-default
   (let [plps (->> reads-for-pileup
-                  (pileup* {:chr "seq1" :start 1 :end 10} {}))]
+                  (pileup* {:chr "seq1" :start 1 :end 10} {}))
+        plps-gap (->> reads-gap-for-pileup
+                      (pileup* {:chr "seq1" :start 1 :end 10} {}))
+        plps-cont (->> reads-continal-for-pileup
+                       (pileup* {:chr "seq1" :start 1 :end 10} {}))
+        plps-overlap (->> reads-overlapping-for-pileup
+                          (pileup* {:chr "seq1" :start 1 :end 10} {}))
+        plps-include (->> reads-including-for-pileup
+                          (pileup* {:chr "seq1" :start 1 :end 10} {}))
+        plps-chunk (->> reads-overlapping-for-pileup
+                        (pileup* {:chr "seq1" :start 1 :end 15} {:chunk-size 2}))]
     (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)))))
+           (map #(map :qual (:pile %)) plps)))
+
+    (is (= (filter pos? [0 0 1 1 1 0 1 1 1 1])
+           (map (comp count :pile) plps-gap)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [] [\T] [\T] [\G] [\G]])
+           (map #(map :base (:pile %)) plps-gap)))
+    (is (= (filter seq [[] [] [33] [89] [93] [] [34] [34] [35] [35]])
+           (map #(map :qual (:pile %)) plps-gap)))
+
+    (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
+           (map (comp count :pile) plps-cont)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\T] [\T] [\G] [\G]])
+           (map #(map :base (:pile %)) plps-cont)))
+    (is (= (filter seq [[] [] [33] [33] [89] [93] [34] [34] [35] [35]])
+           (map #(map :qual (:pile %)) plps-cont)))
+
+    (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
+           (map (comp count :pile) plps-overlap)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\C] [\A]])
+           (map #(map :base (:pile %)) plps-overlap)))
+    (is (= (filter seq [[] [] [33] [33] [34] [34] [35] [124] [74] [74]])
+           (map #(map :qual (:pile %)) plps-overlap)))
+
+    (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1])
+           (map (comp count :pile) plps-include)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\A] [\A]])
+           (map #(map :base (:pile %)) plps-include)))
+    (is (= (filter seq [[] [] [33] [33] [66] [66] [28] [68] [89] [93]])
+           (map #(map :qual (:pile %)) plps-include)))
+
+    (is (= (filter pos? [0 0 1 1 1 1 1 1 1 1 1 1 1 1 1])
+           (map (comp count :pile) plps-chunk)))
+    (is (= (filter seq [[] [] [\T] [\T] [\G] [\G] [\C] [\C] [\C] [\A] [\T] [\T] [\G] [\G] [\A]])
+           (map #(map :base (:pile %)) plps-chunk)))
+    (is (= (filter seq [[] [] [33] [33] [34] [34] [35] [124] [74] [74] [32] [33] [33] [34] [34]])
+           (map #(map :qual (:pile %)) plps-chunk)))))
 
 (deftest overlap-correction
   (let [plps (->> reads-for-pileup
@@ -287,18 +404,6 @@
     (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 (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 (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
   (let [plps (->> reads-for-pileup
                   (pileup* {:chr "seq1" :start 1 :end 10} {:min-base-quality 0}))]