Repository: chrovis/cljam
Pull request #140 (2018-08-10 10:03:21)

Rewrite pileup module
#### Summary
**BREAKING CHANGES**
Completely rewrote `pileup` module in order to improve its performance and usability.

#### Changes
- Don't pile up uncovered regions
    - Sparsely piled-up alignments
- Extract I/O from pileup algorithm
    - Reader/Writer of pileup format
    - Improved serialization with start/end (^/$) and mapq support
- Move everything in `cljam.pileup.mpileup` to `cljam.pileup`
    - Remove `cljam.pileup.common` which is no longer referred
- Add a benchmarking script

##### Performance comparison
Piling up `chr1:1-247249719` of `large.bam` with a single thread.

###### master
```time: 2301.451698 sec, sd: 40390.179732 sec```

###### feature/sparse-mpileup
```time: 48.248668 sec, sd: 648.517263 µs```

#### Affects
- `cljam.algo.pileup`
- `cljam.io.sam.util.cigar`
- `cljam.io.sam.util.quality`

#### Tests

- `lein check` 🆗
- `lein test :all` 🆗

#### Notes
- It could be even faster if we implement the eager version of pileup like `cljam.algo.depth/depth`
- Upgraded cloverage to `1.0.11`


Patch:
diff --git a/.travis.yml b/.travis.yml
index 3379700..ef4a772 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -27,7 +27,7 @@ jobs:
     - stage: coverage
       jdk: oraclejdk8
       script:
-        - CLOVERAGE_VERSION=1.0.9 lein cloverage -e '^cljam.main' --codecov
+        - lein cloverage --codecov
         - bash <(curl -s https://codecov.io/bash) -f target/coverage/codecov.json
     - stage: deploy
       jdk: oraclejdk8
diff --git a/bench/cljam/algo/pileup_bench.clj b/bench/cljam/algo/pileup_bench.clj
new file mode 100644
index 0000000..675ed93
--- /dev/null
+++ b/bench/cljam/algo/pileup_bench.clj
@@ -0,0 +1,18 @@
+(ns cljam.algo.pileup-bench
+  (:require [libra.bench :refer :all]
+            [libra.criterium :as c]
+            [criterium.core :as criterium]
+            [cljam.test-common :as tcommon]
+            [cljam.algo.pileup :as pileup]
+            [cljam.io.sam :as sam]))
+
+(defbench mpileup-large-bench
+  (tcommon/with-before-after {:before (tcommon/prepare-cavia!)}
+    (let [region {:chr "chr1"}]
+      (with-open [r (sam/reader tcommon/large-bam-file)]
+        (dorun (sam/read-alignments r region)) ;; warm up
+        (criterium/with-progress-reporting
+          (is
+           (c/quick-bench
+            (dorun (pileup/pileup r region)))))))))
+
diff --git a/project.clj b/project.clj
index 1336c58..8695fdc 100644
--- a/project.clj
+++ b/project.clj
@@ -21,7 +21,7 @@
                    :plugins [[lein-binplus "0.6.4" :exclusions [org.clojure/clojure]]
                              [lein-codox "0.10.4"]
                              [lein-marginalia "0.9.1" :exclusions [org.clojure/clojure]]
-                             [lein-cloverage "1.0.10" :exclusions [org.clojure/clojure org.tcrawley/dynapath]]
+                             [lein-cloverage "1.0.11" :exclusions [org.clojure/clojure org.tcrawley/dynapath]]
                              [net.totakke/lein-libra "0.1.2"]]
                    :test-selectors {:default #(not-any? % [:slow :remote])
                                     :slow :slow ; Slow tests with local resources
@@ -32,10 +32,10 @@
              :1.7 {:dependencies [[org.clojure/clojure "1.7.0"]]}
              :1.8 {:dependencies [[org.clojure/clojure "1.8.0"]]}
              :1.9 {:dependencies [[org.clojure/clojure "1.9.0"]]}
-             :1.10 {:dependencies [[org.clojure/clojure "1.10.0-alpha4"]]}
+             :1.10 {:dependencies [[org.clojure/clojure "1.10.0-alpha6"]]}
              :uberjar {:dependencies [[org.clojure/clojure "1.9.0"]
-                                      [org.apache.logging.log4j/log4j-api "2.11.0"]
-                                      [org.apache.logging.log4j/log4j-core "2.11.0"]]
+                                      [org.apache.logging.log4j/log4j-api "2.11.1"]
+                                      [org.apache.logging.log4j/log4j-core "2.11.1"]]
                        :resource-paths ["bin-resources"]
                        :main cljam.tools.main
                        :jvm-opts ["-Dclojure.compiler.direct-linking=true"]
diff --git a/src/cljam/algo/pileup.clj b/src/cljam/algo/pileup.clj
index d101eb2..2e9ecac 100644
--- a/src/cljam/algo/pileup.clj
+++ b/src/cljam/algo/pileup.clj
@@ -1,17 +1,209 @@
 (ns cljam.algo.pileup
   "Functions to calculate pileup from the BAM."
   (:require [cljam.io.sam :as sam]
-            [cljam.algo.pileup.mpileup :as mplp]))
+            [cljam.io.sam.util.cigar :as cigar]
+            [cljam.io.sam.util.flag :as flag]
+            [cljam.io.sam.util.quality :as qual]
+            [cljam.io.sam.util.refs :as refs]
+            [cljam.io.pileup :as plpio])
+  (:import [cljam.io.protocols SAMAlignment]
+           [cljam.io.pileup PileupBase LocusPile]))
+
+(defn- cover-locus? [^long pos ^SAMAlignment aln]
+  (<= (.pos aln) pos))
+
+(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)))))
+
+(defn- quals-at-ref
+  [idx ^String qual]
+  (let [empty-qual? (and (= (.length qual) 1)
+                         (= (.charAt qual 0) \*))]
+    (if empty-qual?
+      (vec (repeat (count idx) 93)) ;; \~
+      (mapv (fn [[_ x]]
+              (if (number? x)
+                (qual/fastq-char->phred-byte (.charAt qual x))
+                93))
+            idx))))
+
+(defn- seqs-at-ref
+  [idx ^String s]
+  (mapv (fn [[op x xs]]
+          (let [c (if (number? x) (.charAt s x) x)]
+            (case op
+              :m [c]
+              :d [c xs]
+              :i [c (subs s (first xs) (last xs))]))) idx))
+
+(defn index-cigar
+  "Align bases and base quality scores with the reference coordinate."
+  [^SAMAlignment aln]
+  (let [idx (cigar/to-index (.cigar aln))]
+    (assoc aln
+           :seqs-at-ref (seqs-at-ref idx (:seq aln))
+           :quals-at-ref (quals-at-ref idx (.qual aln)))))
+
+(defn basic-mpileup-pred
+  "Basic predicate function for filtering alignments for mpileup."
+  [^long min-mapq]
+  (fn [^SAMAlignment aln]
+    (let [flag (.flag aln)]
+      (and flag
+           (zero? (bit-and (flag/encoded #{:unmapped
+                                           :filtered-out
+                                           :duplicated
+                                           :supplementary
+                                           :secondary}) flag))
+           (or (not (flag/multiple? flag)) (flag/properly-aligned? flag))
+           (not-empty (.cigar aln))
+           (<= min-mapq (.mapq aln))))))
+
+(defn resolve-base
+  "Find a piled-up base and an indel from an alignment."
+  [^long ref-pos ^SAMAlignment aln]
+  (let [relative-pos (- ref-pos (.pos aln))
+        qual ((:quals-at-ref aln) relative-pos)
+        [base indel] ((:seqs-at-ref aln) relative-pos)]
+    (-> (PileupBase.
+         (zero? relative-pos)
+         (when (zero? relative-pos) (.mapq aln))
+         base
+         qual
+         (flag/reversed? (.flag aln))
+         (= ref-pos (.end aln))
+         (when-not (number? indel) indel)
+         (when (number? indel) indel)
+         (.qname aln))
+        (assoc :alignment aln))))
+
+(defn- resolve-bases
+  [[ref-pos alns]]
+  [ref-pos (mapv (partial resolve-base ref-pos) alns)])
+
+(defn ->locus-pile
+  "Convert a pile into `cljam.io.pileup.LocusPile`."
+  [chr [pos pile]]
+  (let [c (count pile)]
+    (when (pos? c)
+      (LocusPile. chr pos \N c 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."
+  [min-base-quality]
+  (fn [p]
+    (->> #(<= min-base-quality (.qual ^PileupBase %))
+         (partial filterv)
+         (update p 1))))
+
+(defn pileup
+  "Piles up alignments in given region and returns a lazy sequence of
+  `cljam.io.pileup.LocusPile`s."
+  ([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}}]
+   (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)
+            (sequence
+             (comp
+              (filter (basic-mpileup-pred min-map-quality))
+              (map index-cigar)))
+            (pileup-seq s e)
+            (sequence
+             (comp (map resolve-bases)
+                   (map correct-overlaps)
+                   (map (filter-by-base-quality min-base-quality))
+                   (keep (partial ->locus-pile chr)))))))))
+
+(defn align-pileup-seqs
+  "Align multiple piled-up seqs."
+  [& xs]
+  (if (<= (count xs) 1)
+    (map (fn [{:keys [pos] :as m}] [pos [m]]) (first xs))
+    (letfn [(step [xs]
+              (lazy-seq
+               (when (some seq xs)
+                 (let [min-pos (apply min (keep (comp :pos first) xs))]
+                   (cons [min-pos (mapv
+                                   (fn [[{:keys [pos] :as m}]]
+                                     (when (= pos min-pos) m)) xs)]
+                         (step (mapv
+                                (fn [[{:keys [pos]} :as ys]]
+                                  (if (= pos min-pos) (next ys) ys)) xs)))))))]
+      (step xs))))
 
 (defn mpileup
-  "Returns a lazy sequence of cljam.algo.pileup.mpileup.MPileupElement
-  calculated from FASTA and BAM. If start and end are not supplied, piles whole
-  range up."
-  ([bam-reader region] (mpileup nil bam-reader region))
-  ([ref-reader bam-reader region] (mplp/pileup ref-reader bam-reader region)))
+  "Pile up alignments from multiple sources."
+  [region & sam-readers]
+  (apply align-pileup-seqs (map #(pileup % region) sam-readers)))
 
 (defn create-mpileup
   "Creates a mpileup file from the BAM file."
-  [in-bam out-mplp]
-  (with-open [r (sam/bam-reader in-bam)]
-    (mplp/create-mpileup out-mplp nil r)))
+  ([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]
+   (with-open [s (sam/reader in-sam)
+               w (plpio/writer out-mplp in-ref)]
+     (let [regs (if region
+                  [region]
+                  (map
+                   (fn [{:keys [name len]}]
+                     {:chr name :start 1 :end len})
+                   (sam/read-refs s)))]
+       (doseq [reg regs]
+         (plpio/write-piles w (pileup s reg)))))))
diff --git a/src/cljam/algo/pileup/common.clj b/src/cljam/algo/pileup/common.clj
deleted file mode 100644
index c237352..0000000
--- a/src/cljam/algo/pileup/common.clj
+++ /dev/null
@@ -1,4 +0,0 @@
-(ns cljam.algo.pileup.common)
-
-;; TODO: estiamte from actual data
-(def ^:const step 10000)
diff --git a/src/cljam/algo/pileup/mpileup.clj b/src/cljam/algo/pileup/mpileup.clj
deleted file mode 100644
index 91f8712..0000000
--- a/src/cljam/algo/pileup/mpileup.clj
+++ /dev/null
@@ -1,176 +0,0 @@
-(ns cljam.algo.pileup.mpileup
-  (:require [clojure.java.io :as cio]
-            [clojure.string :as cstr]
-            [clojure.tools.logging :as logging]
-            [cljam.io.sam :as sam]
-            [cljam.io.sam.util.cigar :as cigar]
-            [cljam.io.sam.util.flag :as flag]
-            [cljam.io.sam.util.refs :as refs]
-            [cljam.io.sequence :as cseq]))
-
-(defn to-mpileup
-  "Stringify mpileup sequence."
-  [x]
-  (if (vector? x)
-    (let [[y op xs] x] (apply str y op (count xs) xs))
-    (str x)))
-
-(defn substitute-seq
-  "Substitute sequence with mpileup index."
-  [[^Character r & refs] ^String s [op x xs]]
-  (letfn [(get-char [i] (if (number? i) (let [c (.charAt s i)] (if (= c r) \. c)) i))]
-    (case op
-      :m (get-char x)
-      :d [(get-char x) \- (take (count xs) refs)]
-      :i [(get-char x) \+ (map #(.charAt s %) xs)]
-      x)))
-
-(defn substitute-qual
-  "Substitute base quality with mpileup index."
-  [^String q [op x]]
-  (if (and (number? x) (not= q "*")) (.charAt q x) \~))
-
-(defn pileup-seq
-  "Returns a lazy sequence that each element contains reads piled-up at the locus."
-  [^long start ^long end reads]
-  (letfn [(cover-locus? [pos aln]
-            (<= (:pos aln) pos (:end aln)))
-          (step [pos buf alns]
-            (lazy-seq
-             (when (< pos end)
-               (let [[i o] (split-with (partial cover-locus? (inc pos)) alns)
-                     b (doall (concat (filter #(<= (inc pos) (:end %)) buf) i))]
-                 (cons b (step (inc pos) b o))))))]
-    (->> reads
-         (sequence
-          (comp
-           (remove #(empty? (:cigar %)))
-           (drop-while #(< (:end %) start))))
-         (step (dec start) []))))
-
-(defrecord MPileupElement [^String rname ^long pos ^Character ref pile])
-
-(defn gen-mpileup
-  "Compute mpileup info from piled-up reads and reference."
-  [^String rname ^long locus [^Character ref-base :as refs] reads]
-  (->> reads
-       (map (fn gen-mpileup-seq [{:keys [^long pos seq qual cig-index] :as aln}]
-              {:qual (- (int (substitute-qual qual (nth cig-index (- locus pos)))) 33)
-               :seq (substitute-seq refs seq (nth cig-index (- locus pos)))
-               :read aln}))
-       (MPileupElement. rname locus ref-base)))
-
-(defn pileup*
-  "Internal mpileup function independent from I/O.
-   Can take multiple alignments seqs."
-  [refseq {:keys [chr start end]} & aln-seqs]
-  (->> aln-seqs
-       (map
-        (fn [alns]
-          (->> alns
-               (map (fn [a] (assoc a :cig-index (cigar/to-index (:cigar a)))))
-               (pileup-seq start end))))
-       (apply map
-              (fn [index refs & plps]
-                (map (fn [plp] (gen-mpileup chr index refs plp)) plps))
-              (range start (inc end))
-              (partition 100 1 (concat refseq (repeat \N))))))
-
-(defn basic-mpileup-pred
-  "Basic predicate function for filtering alignments for mpileup."
-  [{:keys [flag]}]
-  (and flag
-       (zero? (bit-and (flag/encoded #{:unmapped :filtered-out :duplicated :supplementary :secondary}) flag))
-       (or (not (flag/multiple? flag)) (flag/properly-aligned? flag))))
-
-(defn- correct-qual
-  "Correct quality of two overlapped mate reads by setting zero quality for one of the base."
-  [{s1 :seq q1 :qual :as r1} {s2 :seq q2 :qual :as r2}]
-  (if (= s1 s2)
-    [(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-reads
-  "Correct quality of two overlapped mate reads in piled-up reads."
-  [pile]
-  (if (<= (count (:pile pile)) 1)
-    pile
-    (->> (:pile pile)
-         (group-by (comp :qname :read))
-         (mapcat
-          (fn [[qname xs]]
-            (if (<= (count xs) 1) xs (apply correct-qual xs))))
-         (assoc pile :pile))))
-
-(defn filter-by-base-quality
-  "Returns a predicate for filtering piled-up reads by base quality at its position."
-  [min-base-quality]
-  (fn [p] (update p :pile (fn [pp] (filter #(<= min-base-quality (:qual %)) pp)))))
-
-(defn transpose-pile
-  "Converts a pile {:pile [{:seq SEQ, :qual QUAL, :read READ}]}
-  into {:seq [], :qual [], :reads []}."
-  [p]
-  (if (zero? (count (:pile p)))
-    (assoc p :seq [] :qual [] :reads [] :count 0)
-    (->> (:pile p)
-         (map (juxt :seq #(char (+ (:qual %) 33)) :read))
-         (apply map vector)
-         (zipmap [:seq :qual :reads])
-         (merge p {:count (count (:pile p))}))))
-
-(defn pileup
-  "Returns a lazy sequence of MPileupElement calculated from FASTA and BAM."
-  ([bam-reader region]
-   (pileup nil bam-reader region))
-  ([ref-reader bam-reader {:keys [chr start end] :or {start -1 end -1}}]
-   (try
-     (if-let [r (refs/ref-by-name (sam/read-refs bam-reader) chr)]
-       (let [s (if (neg? start) 1 start)
-             e (if (neg? end) (:len r) end)
-             refseq (if ref-reader
-                      (cseq/read-sequence ref-reader {:chr chr :start s :end e})
-                      (repeat \N))]
-         (->> (sam/read-alignments bam-reader {:chr chr :start s :end e})
-              (sequence
-               (comp
-                (filter basic-mpileup-pred)
-                (filter (fn [a] (<= 0 (:mapq a))))))
-              (pileup* refseq {:chr chr :start s :end e})
-              (sequence
-               (comp (map first)
-                     (map correct-overlapped-reads)
-                     (map (filter-by-base-quality 13))
-                     (map transpose-pile)
-                     (map (fn [p] (update p :seq (fn [s] (map to-mpileup s))))))))))
-     (catch bgzf4j.BGZFException _
-       (throw (RuntimeException. "Invalid file format"))))))
-
-;; Writing
-;; -------
-
-(defn- write-line!
-  [^java.io.BufferedWriter w line]
-  (.write w (cstr/join \tab [(:rname line)
-                             (:pos line)
-                             (:ref line)
-                             (:count line)
-                             (cstr/join (:seq line))
-                             (cstr/join (:qual line))]))
-  (.newLine w))
-
-(defn create-mpileup
-  "Creates a mpileup file from the BAM file."
-  [f fa-rdr bam-reader]
-  (try
-    (with-open [w (cio/writer f)]
-      (doseq [rname (map :name (sam/read-refs bam-reader))]
-        (doseq [line (pileup fa-rdr bam-reader {:chr rname})]
-          (when-not (zero? (:count line))
-            (write-line! w line)))))
-    (catch Exception e (do
-                         (cio/delete-file f)
-                         (logging/error "Failed to create mpileup")
-                         (throw e)))))
diff --git a/src/cljam/io/pileup.clj b/src/cljam/io/pileup.clj
new file mode 100644
index 0000000..4c17db5
--- /dev/null
+++ b/src/cljam/io/pileup.clj
@@ -0,0 +1,199 @@
+(ns cljam.io.pileup
+  (:require [clojure.java.io :as cio]
+            [clojure.string :as cstr]
+            [cljam.io.sam.util.quality :as qual]
+            [cljam.io.sequence :as cseq]
+            [proton.core :as p])
+  (:import [java.io Closeable BufferedWriter]))
+
+;; Records
+;; -------
+
+(defrecord PileupBase [^boolean start?
+                       mapq ;; byte?
+                       ^char base
+                       ^byte qual
+                       ^boolean reverse?
+                       ^boolean end?
+                       insertion ;; String?
+                       deletion ;; int?
+                       qname]) ;; String?
+
+(defrecord LocusPile [rname
+                      ^int pos
+                      ^char ref
+                      ^long count
+                      pile])
+
+;; Reader
+;; ------
+
+(deftype PileupReader [reader]
+  Closeable
+  (close [this]
+    (.close ^Closeable (.reader this))))
+
+(defn ^PileupReader reader
+  "Returns an open instance of cljam.io.pileup.PileupReader of f. Should be used
+  inside with-open to ensure the reader is properly closed."
+  [f]
+  (PileupReader. (cio/reader f)))
+
+(def ^:private
+  upper-table
+  (let [ba (byte-array 128)]
+    (doseq [c "ATGCN"]
+      (aset-byte ba (byte c) (byte c))
+      (aset-byte ba (+ (byte c) (- (byte \a) (byte \A))) (byte c)))
+    (doseq [[from to] [[\, -1] [\. -1] [\< \>] [\> \>] [\* \*]]]
+      (aset-byte ba (byte from) (byte to)))
+    ba))
+
+(defn- parse-bases-col
+  [ref-base ^String column]
+  ;; TODO: Too slow and unreadable
+  (let [len (.length column)]
+    (loop [i 0
+           results (transient [])]
+      (if (< i len)
+        (let [c (.charAt column i)
+              start? (= \^ c)
+              mapq (when start? (qual/fastq-char->phred-byte (.charAt column (inc i))))
+              base (if mapq (.charAt column (+ i 2)) c)
+              base-int (int base)
+              reverse? (or (= base-int 60)
+                           (= base-int 44)
+                           (<= 97 base-int))
+              upper-base-int (aget ^bytes upper-table base-int)
+              upper-base (if (neg? upper-base-int) ref-base (char upper-base-int))]
+          (when (zero? upper-base-int)
+            (throw (ex-info (format "Invalid character %s in %s" base column) {:column column :base base})))
+          (if (= len (+ i (if mapq 3 1)))
+            (persistent! (conj! results (PileupBase. start? mapq upper-base -1 reverse? false nil nil nil)))
+            (let [x (.charAt column (+ i (if mapq 3 1)))
+                  ins? (= x \+)
+                  del? (= x \-)
+                  [indel-num-chars indel-num] (when (or ins? del?)
+                                                (loop [j (+ i (if mapq 4 2))
+                                                       k 0
+                                                       results 0]
+                                                  (let [c (.charAt column j)]
+                                                    (if (Character/isDigit c)
+                                                      (recur (inc j) (inc k) (+ (* 10 results) (- (int c) 48)))
+                                                      [k results]))))
+                  indel-seq (when indel-num
+                              (let [s (+ i (if mapq 4 2) indel-num-chars)]
+                                (cstr/upper-case (.substring column s (+ s indel-num)))))
+                  end? (boolean (or (= \$ x)
+                                    (and indel-num
+                                         (not= len (+ i (if mapq 4 2) indel-num-chars indel-num))
+                                         (= \$ (.charAt column (+ i (if mapq 4 2) indel-num-chars indel-num))))))
+                  next-pos (long (+ i 1 (if mapq 2 0)
+                                    (if indel-num (+ 1 indel-num-chars indel-num) 0)
+                                    (if end? 1 0)))]
+              (recur next-pos
+                     (conj! results (PileupBase. start? mapq upper-base -1 reverse? end? (when ins? indel-seq) (when del? indel-num) nil))))))
+        (persistent! results)))))
+
+(defn- parse-pileup-line
+  [line]
+  (let [[rname pos [ref-base] cnt bases quals & _] (cstr/split line #"\t")
+        ;; TODO: handle extra columns like mapq, qname...
+        upper-ref-base (Character/toUpperCase ^char ref-base)
+        pile (mapv (fn [b q] (assoc b :qual q))
+                   (some->> bases (parse-bases-col upper-ref-base))
+                   (some-> quals qual/fastq->phred))]
+    (LocusPile. rname (p/as-int pos) ref-base (p/as-long cnt) pile)))
+
+(defn read-piles
+  "Reads piled-up bases of the pileup file, returning them as a lazy sequence of
+  cljam.io.pileup.LocusPile."
+  [^PileupReader reader]
+  (->> reader
+       .reader
+       line-seq
+       (map parse-pileup-line)))
+
+
+;; Writer
+;; ------
+
+(defn- write-mpileup-alignment!
+  [^BufferedWriter w ref-reader rname ref-pos ref
+   {:keys [base reverse? mapq start? end? insertion deletion]}]
+  (let [case-base-fn (if-not reverse?
+                       identity
+                       (fn [c]
+                         (if (= c \>)
+                           \<
+                           (Character/toLowerCase ^char c))))
+        case-fn (if-not reverse? identity cstr/lower-case)]
+    (when start?
+      (.append w \^)
+      (.append w (qual/phred-byte->fastq-char mapq)))
+    (if (= base ref) ;; match / deleted / skipped
+      (.append w (if reverse? \, \.))
+      (.append w (unchecked-char (case-base-fn base))))
+    (when deletion
+      (do (.append w \-)
+          (.append w (String/valueOf ^int deletion))
+          (.append w ^String (case-fn
+                              (if ref-reader
+                                (->> {:chr rname
+                                      :start (inc ref-pos)
+                                      :end (+ ref-pos deletion)}
+                                     (cseq/read-sequence ref-reader))
+                                (apply str (repeat deletion \N)))))))
+    (when insertion
+      (do (.append w \+)
+          (.append w (String/valueOf (count insertion)))
+          (.append w ^String (case-fn insertion))))
+    (when end?
+      (.append w \$))))
+
+(defn- write-mpileup-line!
+  [^BufferedWriter w ref-reader {:keys [rname pos pile]}]
+  (let [ref-base (some-> ref-reader
+                         (cseq/read-sequence
+                          {:chr rname :start pos :end pos}
+                          {:mask? true}))
+        ref-char (some-> ref-base cstr/upper-case first)]
+    (.write w ^String rname)
+    (.append w \tab)
+    (.write w (String/valueOf (long pos)))
+    (.append w \tab)
+    (.write w (or ^String ref-base "N"))
+    (.append w \tab)
+    (.write w (String/valueOf (count pile)))
+    (.append w \tab)
+    (doseq [p pile]
+      (write-mpileup-alignment! w ref-reader rname pos ref-char p))
+    (.append w \tab)
+    (doseq [p pile]
+      (.append w (qual/phred-byte->fastq-char (:qual p))))))
+
+(deftype PileupWriter [writer ref-reader]
+  Closeable
+  (close [this]
+    (.close ^Closeable (.writer this))
+    (some-> ^Closeable (.ref-reader this) .close)))
+
+(defn ^PileupWriter writer
+  "Returns an open instance of cljam.io.pileup.PileupWriter of f. Should be used
+  inside with-open to ensure the writer is properly closed. A reference sequence
+  is used when the second arg `reference-path` is supplied."
+  ([f]
+   (writer f nil))
+  ([f reference-path]
+   (PileupWriter. (cio/writer f) (some-> reference-path cseq/reader))))
+
+(defn write-piles
+  "Writes piled-up bases to a file. `pileup-writer` must be an instance of
+  cljam.io.pileup.PileupWriter."
+  [^PileupWriter pileup-writer piles]
+  (let [^BufferedWriter writer (.writer pileup-writer)
+        ref-reader (.ref-reader pileup-writer)]
+    (doseq [pile piles]
+      (write-mpileup-line! writer ref-reader pile)
+      (.newLine writer))
+    (.flush writer)))
diff --git a/src/cljam/io/sam/util/cigar.clj b/src/cljam/io/sam/util/cigar.clj
index 6a9e4e6..fbd7cb2 100644
--- a/src/cljam/io/sam/util/cigar.clj
+++ b/src/cljam/io/sam/util/cigar.clj
@@ -39,10 +39,10 @@
       (if (and l op)
         (condp get op
           #{\M \= \X} (recur xs (+ l r) (+ l s) (concat! idx (map (fn [x] [:m x]) (range s (+ l s)))))
-          #{\D} (recur xs (+ r l) s (concat! (update-last! idx (fn [x] [:d x (range r (+ l r))])) (repeat l [:m \*])))
+          #{\D} (recur xs (+ r l) s (concat! (update-last! idx (fn [x] [:d x l])) (repeat l [:m \*])))
           #{\N} (recur xs (+ r l) s (concat! idx (repeat l [:m \>])))
           #{\S} (recur xs r (+ s l) idx)
-          #{\I} (recur xs r (+ s l) (update-last! idx (fn [x] [:i x (range s (+ l s))]))))
+          #{\I} (recur xs r (+ s l) (update-last! idx (fn [x] [:i x [s (+ l s)]]))))
         (persistent! idx)))))
 
 (def to-index (memoize to-index*))
diff --git a/src/cljam/io/sam/util/quality.clj b/src/cljam/io/sam/util/quality.clj
index e19cdb4..5a7c0e9 100644
--- a/src/cljam/io/sam/util/quality.clj
+++ b/src/cljam/io/sam/util/quality.clj
@@ -2,11 +2,11 @@
   "Utility functions for phred quality strings."
   (:import [java.nio ByteBuffer CharBuffer]))
 
-(defmacro fastq-char->phred-byte [ch]
+(definline fastq-char->phred-byte [ch]
   `(byte (- (int ~ch) 33)))
 
-(defmacro phred-byte->fastq-char [b]
-  `(unchecked-char (unchecked-add ~b 33)))
+(definline phred-byte->fastq-char [b]
+  `(unchecked-char (unchecked-add (Math/min 93 (Math/max (long ~b) 0)) 33)))
 
 (defn fastq->phred ^bytes [^String fastq]
   (let [b (.getBytes fastq)
diff --git a/src/cljam/tools/cli.clj b/src/cljam/tools/cli.clj
index e7572bd..3bba382 100644
--- a/src/cljam/tools/cli.clj
+++ b/src/cljam/tools/cli.clj
@@ -7,17 +7,17 @@
             [cljam.io.sam :as sam]
             [cljam.io.sam.util :as sam-util]
             [cljam.io.sam.util.header :as header]
-            [cljam.io.sequence :as cseq]
             [cljam.algo.bam-indexer :as bai]
             [cljam.algo.normal :as normal]
             [cljam.algo.sorter :as sorter]
             [cljam.algo.fasta-indexer :as fai]
             [cljam.algo.dict :as dict]
             [cljam.algo.depth :as depth]
-            [cljam.algo.pileup :as plp]
             [cljam.algo.convert :as convert]
             [cljam.algo.level :as level]
-            [cljam.util.region :as region])
+            [cljam.algo.pileup :as plp]
+            [cljam.util.region :as region]
+            [clojure.java.io :as cio])
   (:import [java.io Closeable BufferedWriter OutputStreamWriter]))
 
 ;; CLI functions
@@ -36,6 +36,12 @@
   (str "The following errors occurred while parsing your command:\n\n"
        (cstr/join \newline errors)))
 
+(defn- parse-region [region-str]
+  (when region-str
+    (if-let [reg (region/parse-region region-str)]
+      reg
+      (exit 1 (str "Invalid region format: " region-str)))))
+
 ;; Sub-commands
 ;; ------------
 
@@ -70,12 +76,10 @@
                                  "bam"  (sam/bam-reader f))]
         (when (:header options)
           (println (header/stringify-header (sam/read-header r))))
-        (doseq [aln (if (:region options)
-                      (if-let [region (region/parse-region (:region options))]
-                        (if (sam/indexed? r)
-                          (sam/read-alignments r region)
-                          (exit 1 "Random alignment retrieval only works for indexed BAM."))
-                        (exit 1 "Invalid region format"))
+        (doseq [aln (if-let [region (parse-region (:region options))]
+                      (if (sam/indexed? r)
+                        (sam/read-alignments r region)
+                        (exit 1 "Random alignment retrieval only works for indexed BAM."))
                       (sam/read-alignments r))]
           (println (sam-util/stringify-alignment aln))))))
   nil)
@@ -214,77 +218,34 @@
         options-summary]
        (cstr/join \newline)))
 
-(defn- pileup-simple
-  ([rdr n-threads]
-   (doseq [rname (map :name (sam/read-refs rdr))]
-     (pileup-simple rdr n-threads {:chr rname})))
-  ([rdr n-threads region]
-   (binding [*out* (BufferedWriter. (OutputStreamWriter. System/out))
-             *flush-on-newline* false]
-     (doseq [line (depth/lazy-depth rdr region {:n-threads n-threads})]
-       (println line))
-     (flush))))
-
-(defn- pileup-with-ref
-  ([rdr ref-fa]
-   (with-open [fa-rdr (cseq/reader ref-fa)]
-     (doseq [rname (map :name (sam/read-refs rdr))
-             line  (plp/mpileup fa-rdr rdr {:chr rname})]
-       (if-not (zero? (:count line))
-         (println (cstr/join \tab (map #(case %
-                                          :qual (cstr/join (% line))
-                                          :seq (cstr/join (% line))
-                                          (% line)) [:rname :pos :ref :count :seq :qual])))))))
-  ([rdr ref-fa region]
-   (with-open [fa-rdr (cseq/reader ref-fa)]
-     (doseq [line  (plp/mpileup fa-rdr rdr region)]
-       (if-not (zero? (:count line))
-         (println (cstr/join \tab (map #(case %
-                                          :qual (cstr/join (% line))
-                                          :seq (cstr/join (% line))
-                                          (% line)) [:rname :pos :ref :count :seq :qual]))))))))
-
-(defn- pileup-without-ref
-  ([rdr]
-   (doseq [rname (map :name (sam/read-refs rdr))
-           line  (plp/mpileup rdr {:chr rname})]
-     (if-not (zero? (:count line))
-       (println (cstr/join \tab (map #(case %
-                                        :qual (cstr/join (% line))
-                                        :seq (cstr/join (% line))
-                                        (% line)) [:rname :pos :ref :count :seq :qual]))))))
-  ([rdr region]
-   (doseq [line  (plp/mpileup nil rdr region)]
-     (if-not (zero? (:count line))
-       (println (cstr/join \tab (map #(case %
-                                        :qual (cstr/join (% line))
-                                        :seq (cstr/join (% line))
-                                        (% line)) [:rname :pos :ref :count :seq :qual])))))))
+(defn- depth
+  [f region n-threads]
+  (with-open [r (sam/reader f)]
+    (when-not (sam/indexed? r)
+      (exit 1 "Random alignment retrieval only works for indexed BAM."))
+    (when-not (sorter/sorted? r)
+      (exit 1 "Not sorted"))
+    (let [regs (or (some-> region parse-region vector)
+                   (map (fn [{:keys [name len]}] {:chr name :start 1 :end len}) (sam/read-refs r)))]
+      (binding [*out* (BufferedWriter. (OutputStreamWriter. System/out))
+                *flush-on-newline* false]
+        (doseq [reg regs
+                line (depth/lazy-depth r reg {:n-threads n-threads})]
+          (println line))
+        (flush)))))
 
 (defn pileup [args]
-  (let [{:keys [options arguments errors summary]} (parse-opts args pileup-cli-options)]
+  (let [{:keys [arguments errors summary]
+         {:keys [help region simple ref thread]} :options} (parse-opts args pileup-cli-options)]
     (cond
-     (:help options) (exit 0 (pileup-usage summary))
-     (not= (count arguments) 1) (exit 1 (pileup-usage summary))
-     errors (exit 1 (error-msg errors)))
+      help (exit 0 (pileup-usage summary))
+      (not= (count arguments) 1) (exit 1 (pileup-usage summary))
+      errors (exit 1 (error-msg errors)))
     (let [f (first arguments)]
-      (with-open [r (sam/reader f)]
-        (when (= (type r) cljam.io.sam.reader.SAMReader)
-          (exit 1 "Not support SAM file"))
-        (when-not (sorter/sorted? r)
-          (exit 1 "Not sorted"))
-        (if (:region options)
-          (if-let [region (region/parse-region (:region options))]
-            (cond
-              (:simple options) (pileup-simple r (:thread options) region)
-              (:ref options) (pileup-with-ref r (:ref options) region)
-              :else (pileup-without-ref r region))
-            (exit 1 "Invalid region format"))
-          (cond
-            (:simple options) (pileup-simple r (:thread options))
-            (:ref options) (pileup-with-ref r (:ref options))
-            :else (pileup-without-ref r))))))
-  nil)
+      (if simple
+        (depth f region thread)
+        (with-open [w (cio/writer (cio/output-stream System/out))]
+          (plp/create-mpileup f ref w (parse-region region)))))))
 
 ;; ### faidx command
 


Test patch:
diff --git a/test-resources/pileup/f.pileup b/test-resources/pileup/f.pileup
index 7f0defa..2f9f514 100644
--- a/test-resources/pileup/f.pileup
+++ b/test-resources/pileup/f.pileup
@@ -1,75 +1,75 @@
-ref	7	T	1	.	~
+ref	7	T	1	^?.	~
 ref	8	T	1	.	~
-ref	9	A	3	...	~~~
+ref	9	A	3	.^?.^?.	~~~
 ref	10	G	3	...	~~~
 ref	11	A	3	..C	~~~
 ref	12	T	3	...	~~~
 ref	13	A	3	...	~~~
-ref	14	A	3	.+4AGAG.+2GG.	~~~
+ref	14	A	3	.+4AGAG.+2GG.$	~~~
 ref	15	G	2	..	~~
-ref	16	A	3	...	~~~
+ref	16	A	3	..^?.	~~~
 ref	17	T	3	...	~~~
-ref	18	A	3	.-1G.+2AA.	~~~
+ref	18	A	3	.-1G.+2AA$.	~~~
 ref	19	G	2	*.	~~
 ref	20	C	2	..	~~
 ref	21	T	2	..	~~
-ref	22	G	2	.>	~~
+ref	22	G	2	.$>	~~
 ref	23	T	1	>	~
 ref	24	G	1	>	~
 ref	25	C	1	>	~
 ref	26	T	1	>	~
 ref	27	A	1	>	~
 ref	28	G	1	>	~
-ref	29	T	2	>.	~~
-ref	30	A	2	>.	~~
-ref	31	G	2	>.	~~
-ref	32	G	2	>.	~~
-ref	33	C	2	>.	~~
+ref	29	T	2	>^?,	~~
+ref	30	A	2	>,	~~
+ref	31	G	2	>,	~~
+ref	32	G	2	>,	~~
+ref	33	C	2	>,$	~~
 ref	34	A	1	>	~
 ref	35	G	1	>+1C	~
 ref	36	T	1	.	~
-ref	37	C	2	..	~~
-ref	38	A	2	..	~~
-ref	39	G	2	..	~~
-ref	40	C	2	..	~~
-ref	41	G	1	.	~
-ref	42	C	1	.	~
-ref	43	C	1	.	~
-ref	44	A	1	.	~
-ref	45	T	1	.	~
-ref2	1	A	1	.	?
-ref2	2	G	2	..	??
-ref2	3	G	2	..	??
-ref2	4	T	2	..	??
-ref2	5	T	2	..	??
-ref2	6	T	3	...	???
-ref2	7	T	3	...	???
-ref2	8	A	3	...	???
-ref2	9	T	3	...	???
-ref2	10	A	4	...C	????
-ref2	11	A	4	....	????
-ref2	12	A	5	.....	?????
-ref2	13	A	5	.....	?????
-ref2	14	C	6	...+4AAATTTT	??????
-ref2	15	A	6	......	??????
-ref2	16	A	6	......	??????
-ref2	17	T	6	AA....	??????
-ref2	18	T	6	......	??????
-ref2	19	A	6	......	??????
-ref2	20	A	6	......	??????
-ref2	21	G	5	T....	?????
-ref2	22	T	5	.....	?????
-ref2	23	C	4	....	????
-ref2	24	T	4	....	????
-ref2	25	A	4	....	????
-ref2	26	C	4	....	????
-ref2	27	A	4	....	????
-ref2	28	G	3	...	???
-ref2	29	A	3	...	???
-ref2	30	G	3	...	???
-ref2	31	C	3	...	???
-ref2	32	A	3	...	???
-ref2	33	A	3	...	???
-ref2	34	C	3	...	???
-ref2	35	T	2	..	??
-ref2	36	A	1	.	?
+ref	37	C	2	.^?,	~~
+ref	38	A	2	.,	~~
+ref	39	G	2	.,	~~
+ref	40	C	2	.$,	~~
+ref	41	G	1	,	~
+ref	42	C	1	,	~
+ref	43	C	1	,	~
+ref	44	A	1	,	~
+ref	45	T	1	,$	~
+ref2	1	a	1	^?.	?
+ref2	2	g	2	.^?.	??
+ref2	3	g	2	..	??
+ref2	4	t	2	..	??
+ref2	5	t	2	..	??
+ref2	6	t	3	..^?.	???
+ref2	7	t	3	...	???
+ref2	8	a	3	...	???
+ref2	9	t	3	...	???
+ref2	10	a	4	...^?C	????
+ref2	11	a	4	....	????
+ref2	12	a	5	....^?.	?????
+ref2	13	a	5	.....	?????
+ref2	14	c	6	...+4AAATTT^?T	??????
+ref2	15	a	6	......	??????
+ref2	16	a	6	......	??????
+ref2	17	t	6	AA....	??????
+ref2	18	t	6	......	??????
+ref2	19	a	6	......	??????
+ref2	20	a	6	.$.....	??????
+ref2	21	g	5	T....	?????
+ref2	22	t	5	.$....	?????
+ref2	23	c	4	....	????
+ref2	24	t	4	....	????
+ref2	25	a	4	....	????
+ref2	26	c	4	....	????
+ref2	27	a	4	.$...	????
+ref2	28	g	3	...	???
+ref2	29	a	3	...	???
+ref2	30	g	3	...	???
+ref2	31	c	3	...	???
+ref2	32	a	3	...	???
+ref2	33	a	3	...	???
+ref2	34	c	3	.$..	???
+ref2	35	t	2	.$.	??
+ref2	36	a	1	.$	?
diff --git a/test-resources/pileup/r1.pileup b/test-resources/pileup/r1.pileup
index a3ac12d..a833e1f 100644
--- a/test-resources/pileup/r1.pileup
+++ b/test-resources/pileup/r1.pileup
@@ -1,36 +1,36 @@
-ref2	1	N	1	A	?
-ref2	2	N	2	GG	??
+ref2	1	N	1	^?A	?
+ref2	2	N	2	G^?G	??
 ref2	3	N	2	GG	??
 ref2	4	N	2	TT	??
 ref2	5	N	2	TT	??
-ref2	6	N	3	TTT	???
+ref2	6	N	3	TT^?T	???
 ref2	7	N	3	TTT	???
 ref2	8	N	3	AAA	???
 ref2	9	N	3	TTT	???
-ref2	10	N	4	AAAC	????
+ref2	10	N	4	AAA^?C	????
 ref2	11	N	4	AAAA	????
-ref2	12	N	5	AAAAA	?????
+ref2	12	N	5	AAAA^?A	?????
 ref2	13	N	5	AAAAA	?????
-ref2	14	N	6	CCC+4AAATTTT	??????
+ref2	14	N	6	CCC+4AAATTT^?T	??????
 ref2	15	N	6	AAAAAA	??????
 ref2	16	N	6	AAAAAA	??????
 ref2	17	N	6	AATTTT	??????
 ref2	18	N	6	TTTTTT	??????
 ref2	19	N	6	AAAAAA	??????
-ref2	20	N	6	AAAAAA	??????
+ref2	20	N	6	A$AAAAA	??????
 ref2	21	N	5	TGGGG	?????
-ref2	22	N	5	TTTTT	?????
+ref2	22	N	5	T$TTTT	?????
 ref2	23	N	4	CCCC	????
 ref2	24	N	4	TTTT	????
 ref2	25	N	4	AAAA	????
 ref2	26	N	4	CCCC	????
-ref2	27	N	4	AAAA	????
+ref2	27	N	4	A$AAA	????
 ref2	28	N	3	GGG	???
 ref2	29	N	3	AAA	???
 ref2	30	N	3	GGG	???
 ref2	31	N	3	CCC	???
 ref2	32	N	3	AAA	???
 ref2	33	N	3	AAA	???
-ref2	34	N	3	CCC	???
-ref2	35	N	2	TT	??
-ref2	36	N	1	A	?
+ref2	34	N	3	C$CC	???
+ref2	35	N	2	T$T	??
+ref2	36	N	1	A$	?
diff --git a/test-resources/pileup/r1f.pileup b/test-resources/pileup/r1f.pileup
index eb40366..427ab05 100644
--- a/test-resources/pileup/r1f.pileup
+++ b/test-resources/pileup/r1f.pileup
@@ -1,36 +1,36 @@
-ref2	1	A	1	.	?
-ref2	2	G	2	..	??
-ref2	3	G	2	..	??
-ref2	4	T	2	..	??
-ref2	5	T	2	..	??
-ref2	6	T	3	...	???
-ref2	7	T	3	...	???
-ref2	8	A	3	...	???
-ref2	9	T	3	...	???
-ref2	10	A	4	...C	????
-ref2	11	A	4	....	????
-ref2	12	A	5	.....	?????
-ref2	13	A	5	.....	?????
-ref2	14	C	6	...+4AAATTTT	??????
-ref2	15	A	6	......	??????
-ref2	16	A	6	......	??????
-ref2	17	T	6	AA....	??????
-ref2	18	T	6	......	??????
-ref2	19	A	6	......	??????
-ref2	20	A	6	......	??????
-ref2	21	G	5	T....	?????
-ref2	22	T	5	.....	?????
-ref2	23	C	4	....	????
-ref2	24	T	4	....	????
-ref2	25	A	4	....	????
-ref2	26	C	4	....	????
-ref2	27	A	4	....	????
-ref2	28	G	3	...	???
-ref2	29	A	3	...	???
-ref2	30	G	3	...	???
-ref2	31	C	3	...	???
-ref2	32	A	3	...	???
-ref2	33	A	3	...	???
-ref2	34	C	3	...	???
-ref2	35	T	2	..	??
-ref2	36	A	1	.	?
+ref2	1	a	1	^?.	?
+ref2	2	g	2	.^?.	??
+ref2	3	g	2	..	??
+ref2	4	t	2	..	??
+ref2	5	t	2	..	??
+ref2	6	t	3	..^?.	???
+ref2	7	t	3	...	???
+ref2	8	a	3	...	???
+ref2	9	t	3	...	???
+ref2	10	a	4	...^?C	????
+ref2	11	a	4	....	????
+ref2	12	a	5	....^?.	?????
+ref2	13	a	5	.....	?????
+ref2	14	c	6	...+4AAATTT^?T	??????
+ref2	15	a	6	......	??????
+ref2	16	a	6	......	??????
+ref2	17	t	6	AA....	??????
+ref2	18	t	6	......	??????
+ref2	19	a	6	......	??????
+ref2	20	a	6	.$.....	??????
+ref2	21	g	5	T....	?????
+ref2	22	t	5	.$....	?????
+ref2	23	c	4	....	????
+ref2	24	t	4	....	????
+ref2	25	a	4	....	????
+ref2	26	c	4	....	????
+ref2	27	a	4	.$...	????
+ref2	28	g	3	...	???
+ref2	29	a	3	...	???
+ref2	30	g	3	...	???
+ref2	31	c	3	...	???
+ref2	32	a	3	...	???
+ref2	33	a	3	...	???
+ref2	34	c	3	.$..	???
+ref2	35	t	2	.$.	??
+ref2	36	a	1	.$	?
diff --git a/test-resources/pileup/r2.pileup b/test-resources/pileup/r2.pileup
index 0cb3b30..e589959 100644
--- a/test-resources/pileup/r2.pileup
+++ b/test-resources/pileup/r2.pileup
@@ -1,27 +1,27 @@
-ref2	10	N	4	AAAC	????
+ref2	10	N	4	AAA^?C	????
 ref2	11	N	4	AAAA	????
-ref2	12	N	5	AAAAA	?????
+ref2	12	N	5	AAAA^?A	?????
 ref2	13	N	5	AAAAA	?????
-ref2	14	N	6	CCC+4AAATTTT	??????
+ref2	14	N	6	CCC+4AAATTT^?T	??????
 ref2	15	N	6	AAAAAA	??????
 ref2	16	N	6	AAAAAA	??????
 ref2	17	N	6	AATTTT	??????
 ref2	18	N	6	TTTTTT	??????
 ref2	19	N	6	AAAAAA	??????
-ref2	20	N	6	AAAAAA	??????
+ref2	20	N	6	A$AAAAA	??????
 ref2	21	N	5	TGGGG	?????
-ref2	22	N	5	TTTTT	?????
+ref2	22	N	5	T$TTTT	?????
 ref2	23	N	4	CCCC	????
 ref2	24	N	4	TTTT	????
 ref2	25	N	4	AAAA	????
 ref2	26	N	4	CCCC	????
-ref2	27	N	4	AAAA	????
+ref2	27	N	4	A$AAA	????
 ref2	28	N	3	GGG	???
 ref2	29	N	3	AAA	???
 ref2	30	N	3	GGG	???
 ref2	31	N	3	CCC	???
 ref2	32	N	3	AAA	???
 ref2	33	N	3	AAA	???
-ref2	34	N	3	CCC	???
-ref2	35	N	2	TT	??
-ref2	36	N	1	A	?
+ref2	34	N	3	C$CC	???
+ref2	35	N	2	T$T	??
+ref2	36	N	1	A$	?
diff --git a/test-resources/pileup/r2f.pileup b/test-resources/pileup/r2f.pileup
index 234d180..fbbf360 100644
--- a/test-resources/pileup/r2f.pileup
+++ b/test-resources/pileup/r2f.pileup
@@ -1,27 +1,27 @@
-ref2	10	A	4	...C	????
-ref2	11	A	4	....	????
-ref2	12	A	5	.....	?????
-ref2	13	A	5	.....	?????
-ref2	14	C	6	...+4AAATTTT	??????
-ref2	15	A	6	......	??????
-ref2	16	A	6	......	??????
-ref2	17	T	6	AA....	??????
-ref2	18	T	6	......	??????
-ref2	19	A	6	......	??????
-ref2	20	A	6	......	??????
-ref2	21	G	5	T....	?????
-ref2	22	T	5	.....	?????
-ref2	23	C	4	....	????
-ref2	24	T	4	....	????
-ref2	25	A	4	....	????
-ref2	26	C	4	....	????
-ref2	27	A	4	....	????
-ref2	28	G	3	...	???
-ref2	29	A	3	...	???
-ref2	30	G	3	...	???
-ref2	31	C	3	...	???
-ref2	32	A	3	...	???
-ref2	33	A	3	...	???
-ref2	34	C	3	...	???
-ref2	35	T	2	..	??
-ref2	36	A	1	.	?
+ref2	10	a	4	...^?C	????
+ref2	11	a	4	....	????
+ref2	12	a	5	....^?.	?????
+ref2	13	a	5	.....	?????
+ref2	14	c	6	...+4AAATTT^?T	??????
+ref2	15	a	6	......	??????
+ref2	16	a	6	......	??????
+ref2	17	t	6	AA....	??????
+ref2	18	t	6	......	??????
+ref2	19	a	6	......	??????
+ref2	20	a	6	.$.....	??????
+ref2	21	g	5	T....	?????
+ref2	22	t	5	.$....	?????
+ref2	23	c	4	....	????
+ref2	24	t	4	....	????
+ref2	25	a	4	....	????
+ref2	26	c	4	....	????
+ref2	27	a	4	.$...	????
+ref2	28	g	3	...	???
+ref2	29	a	3	...	???
+ref2	30	g	3	...	???
+ref2	31	c	3	...	???
+ref2	32	a	3	...	???
+ref2	33	a	3	...	???
+ref2	34	c	3	.$..	???
+ref2	35	t	2	.$.	??
+ref2	36	a	1	.$	?
diff --git a/test-resources/pileup/test.pileup b/test-resources/pileup/test.pileup
index d42f2cc..6a4eede 100644
--- a/test-resources/pileup/test.pileup
+++ b/test-resources/pileup/test.pileup
@@ -1,75 +1,75 @@
-ref	7	N	1	T	~
+ref	7	N	1	^?T	~
 ref	8	N	1	T	~
-ref	9	N	3	AAA	~~~
+ref	9	N	3	A^?A^?A	~~~
 ref	10	N	3	GGG	~~~
 ref	11	N	3	AAC	~~~
 ref	12	N	3	TTT	~~~
 ref	13	N	3	AAA	~~~
-ref	14	N	3	A+4AGAGA+2GGA	~~~
+ref	14	N	3	A+4AGAGA+2GGA$	~~~
 ref	15	N	2	GG	~~
-ref	16	N	3	AAA	~~~
+ref	16	N	3	AA^?A	~~~
 ref	17	N	3	TTT	~~~
-ref	18	N	3	A-1NA+2AAA	~~~
+ref	18	N	3	A-1NA+2AA$A	~~~
 ref	19	N	2	*G	~~
 ref	20	N	2	CC	~~
 ref	21	N	2	TT	~~
-ref	22	N	2	G>	~~
+ref	22	N	2	G$>	~~
 ref	23	N	1	>	~
 ref	24	N	1	>	~
 ref	25	N	1	>	~
 ref	26	N	1	>	~
 ref	27	N	1	>	~
 ref	28	N	1	>	~
-ref	29	N	2	>T	~~
-ref	30	N	2	>A	~~
-ref	31	N	2	>G	~~
-ref	32	N	2	>G	~~
-ref	33	N	2	>C	~~
+ref	29	N	2	>^?t	~~
+ref	30	N	2	>a	~~
+ref	31	N	2	>g	~~
+ref	32	N	2	>g	~~
+ref	33	N	2	>c$	~~
 ref	34	N	1	>	~
 ref	35	N	1	>+1C	~
 ref	36	N	1	T	~
-ref	37	N	2	CC	~~
-ref	38	N	2	AA	~~
-ref	39	N	2	GG	~~
-ref	40	N	2	CC	~~
-ref	41	N	1	G	~
-ref	42	N	1	C	~
-ref	43	N	1	C	~
-ref	44	N	1	A	~
-ref	45	N	1	T	~
-ref2	1	N	1	A	?
-ref2	2	N	2	GG	??
+ref	37	N	2	C^?c	~~
+ref	38	N	2	Aa	~~
+ref	39	N	2	Gg	~~
+ref	40	N	2	C$c	~~
+ref	41	N	1	g	~
+ref	42	N	1	c	~
+ref	43	N	1	c	~
+ref	44	N	1	a	~
+ref	45	N	1	t$	~
+ref2	1	N	1	^?A	?
+ref2	2	N	2	G^?G	??
 ref2	3	N	2	GG	??
 ref2	4	N	2	TT	??
 ref2	5	N	2	TT	??
-ref2	6	N	3	TTT	???
+ref2	6	N	3	TT^?T	???
 ref2	7	N	3	TTT	???
 ref2	8	N	3	AAA	???
 ref2	9	N	3	TTT	???
-ref2	10	N	4	AAAC	????
+ref2	10	N	4	AAA^?C	????
 ref2	11	N	4	AAAA	????
-ref2	12	N	5	AAAAA	?????
+ref2	12	N	5	AAAA^?A	?????
 ref2	13	N	5	AAAAA	?????
-ref2	14	N	6	CCC+4AAATTTT	??????
+ref2	14	N	6	CCC+4AAATTT^?T	??????
 ref2	15	N	6	AAAAAA	??????
 ref2	16	N	6	AAAAAA	??????
 ref2	17	N	6	AATTTT	??????
 ref2	18	N	6	TTTTTT	??????
 ref2	19	N	6	AAAAAA	??????
-ref2	20	N	6	AAAAAA	??????
+ref2	20	N	6	A$AAAAA	??????
 ref2	21	N	5	TGGGG	?????
-ref2	22	N	5	TTTTT	?????
+ref2	22	N	5	T$TTTT	?????
 ref2	23	N	4	CCCC	????
 ref2	24	N	4	TTTT	????
 ref2	25	N	4	AAAA	????
 ref2	26	N	4	CCCC	????
-ref2	27	N	4	AAAA	????
+ref2	27	N	4	A$AAA	????
 ref2	28	N	3	GGG	???
 ref2	29	N	3	AAA	???
 ref2	30	N	3	GGG	???
 ref2	31	N	3	CCC	???
 ref2	32	N	3	AAA	???
 ref2	33	N	3	AAA	???
-ref2	34	N	3	CCC	???
-ref2	35	N	2	TT	??
-ref2	36	N	1	A	?
+ref2	34	N	3	C$CC	???
+ref2	35	N	2	T$T	??
+ref2	36	N	1	A$	?
diff --git a/test/cljam/algo/pileup_test.clj b/test/cljam/algo/pileup_test.clj
index 01aeb88..bbc5aa6 100644
--- a/test/cljam/algo/pileup_test.clj
+++ b/test/cljam/algo/pileup_test.clj
@@ -1,19 +1,40 @@
 (ns cljam.algo.pileup-test
   (:require [clojure.test :refer :all]
             [cljam.test-common :refer :all]
+            [cljam.io.protocols :as p]
             [cljam.io.sam :as sam]
+            [cljam.io.sam.util.quality :as qual]
             [cljam.io.sequence :as cseq]
             [cljam.algo.pileup :as plp]
-            [cljam.algo.pileup.mpileup :as mplp]))
+            [clojure.string :as cstr]
+            [cljam.io.pileup :as plpio]))
 
-(def test-bam-pileup-ref '(0 0 0 0 0 0 1 1 3 3 3 3 3 3 2 3 3 3 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 1 1 1 2 2 2 2 1 1 1 1 1))
-(def test-bam-pileup-ref2 '(1 2 2 2 2 3 3 3 3 4 4 5 5 6 6 6 6 6 6 6 5 5 4 4 4 4 4 3 3 3 3 3 3 3 2 1 0 0 0 0))
+(def test-bam-pileup-ref [0 0 0 0 0 0 1 1 3 3 3 3 3 3 2 3 3 3 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 1 1 1 2 2 2 2 1 1 1 1 1])
+(def test-bam-pileup-ref2 [1 2 2 2 2 3 3 3 3 4 4 5 5 6 6 6 6 6 6 6 5 5 4 4 4 4 4 3 3 3 3 3 3 3 2 1 0 0 0 0])
 (def test-bam-mpileup-seq-ref
   '(() () () () () () ("T") ("T") ("A" "A" "A") ("G" "G" "G") ("A" "A" "C") ("T" "T" "T") ("A" "A" "A")
     ("A+4AGAG" "A+2GG" "A") ("G" "G") ("A" "A" "A") ("T" "T" "T") ("A-1N" "A+2AA" "A") ("*" "G") ("C" "C")
     ("T" "T") ("G" ">") (">") (">") (">") (">") (">") (">") (">" "T") (">" "A") (">" "G") (">" "G") (">" "C")
     (">") (">+1C") ;; adjacent indel >+1T
     ("T") ("C" "C") ("A" "A") ("G" "G") ("C" "C") ("G") ("C") ("C") ("A") ("T")))
+(def test-bam-mpileup-seq-ref-freq
+  [{} {} {} {} {} {} {\T 1} {\T 1} {\A 3} {\G 3} {\A 2 \C 1} {\T 3} {\A 3}
+   {\A 3} {\G 2} {\A 3} {\T 3} {\A 3} {\G 1 \* 1} {\C 2}
+   {\T 2} {\G 1 \> 1} {\> 1} {\> 1} {\> 1} {\> 1} {\> 1} {\> 1} {\T 1 \> 1} {\A 1 \> 1} {\G 1 \> 1} {\G 1 \> 1} {\C 1 \> 1}
+   {\> 1} {\> 1}
+   {\T 1} {\C 2} {\A 2} {\G 2} {\C 2} {\G 1} {\C 1} {\C 1} {\A 1} {\T 1}])
+(def test-bam-mpileup-seq-ref-ins
+  [{} {} {} {} {} {} {} {} {} {} {} {} {}
+   {"AGAG" 1 "GG" 1} {} {} {} {"AA" 1} {} {}
+   {} {} {} {} {} {} {} {} {} {} {} {} {}
+   {} {"C" 1}
+   {} {} {} {} {} {} {} {} {} {}])
+(def test-bam-mpileup-seq-ref-del
+  [{} {} {} {} {} {} {} {} {} {} {} {} {}
+   {} {} {} {} {1 1} {} {}
+   {} {} {} {} {} {} {} {} {} {} {} {} {}
+   {} {}
+   {} {} {} {} {} {} {} {} {} {}])
 (def test-bam-mpileup-seq-ref2
   '(("A") ("G" "G") ("G" "G") ("T" "T") ("T" "T") ("T" "T" "T") ("T" "T" "T") ("A" "A" "A") ("T" "T" "T")
     ("A" "A" "A" "C") ("A" "A" "A" "A") ("A" "A" "A" "A" "A") ("A" "A" "A" "A" "A") ("C" "C" "C+4AAAT" "T" "T" "T")
@@ -21,6 +42,13 @@
     ("A" "A" "A" "A" "A" "A") ("A" "A" "A" "A" "A" "A") ("T" "G" "G" "G" "G") ("T" "T" "T" "T" "T")
     ("C" "C" "C" "C") ("T" "T" "T" "T") ("A" "A" "A" "A") ("C" "C" "C" "C") ("A" "A" "A" "A") ("G" "G" "G")
     ("A" "A" "A") ("G" "G" "G") ("C" "C" "C") ("A" "A" "A") ("A" "A" "A") ("C" "C" "C") ("T" "T") ("A") () () () ()))
+(def test-bam-mpileup-seq-ref2-freq
+  [{\A 1} {\G 2} {\G 2} {\T 2} {\T 2} {\T 3} {\T 3} {\A 3} {\T 3}
+   {\A 3 \C 1} {\A 4} {\A 5} {\A 5} {\C 3 \T 3}
+   {\A 6} {\A 6} {\A 2 \T 4} {\T 6}
+   {\A 6} {\A 6} {\T 1 \G 4} {\T 5}
+   {\C 4} {\T 4} {\A 4} {\C 4} {\A 4} {\G 3}
+   {\A 3} {\G 3} {\C 3} {\A 3} {\A 3} {\C 3} {\T 2} {\A 1} {} {} {} {}])
 (def test-bam-mpileup-seq-ref-with-ref
   '(() () () () () () (".") (".") ("." "." ".") ("." "." ".") ("." "." "C") ("." "." ".") ("." "." ".")
     (".+4AGAG" ".+2GG" ".") ("." ".") ("." "." ".") ("." "." ".") (".-1G" ".+2AA" ".") ("*" ".") ("." ".")
@@ -33,27 +61,44 @@
     [\~] [\~] [\~ \~] [\~ \~] [\~ \~] [\~ \~] [\~ \~] [\~] [\~] [\~] [\~ \~] [\~ \~] [\~ \~] [\~ \~]
     [\~] [\~] [\~] [\~] [\~]))
 
+(defn- ->aln [m]
+  (-> {:qname "DUMMY", :flag 0, :rname "DUMMY", :pos -1, :end -1, :rnext "DUMMY", :pnext -1, :tlen -1, :mapq -1, :cigar "10M", :seq "AAAAAAAAAA", :qual "IIIIIIIIII"}
+      (merge m)
+      (update :pos int)
+      (update :end int)
+      (update :flag int)
+      (update :tlen int)
+      (update :pnext int)
+      (update :mapq int)
+      p/map->SAMAlignment))
+
+(defn- ->pbase [m]
+  (plpio/map->PileupBase (merge {:start? false :reverse? false :end? false} m)))
+
+(defn- ->pile [[pos pile]]
+  [pos (map (juxt :pos :end) pile)])
+
 (deftest about-pileup-seq
 
   ;; ----------
   ;; 1234567890...
-  (is (= (map (fn [xs] (map #(dissoc % :end) xs))
-              (mplp/pileup-seq 1 2 [{:pos 1 :cigar "10M" :end 10}]))
-         [[{:pos 1 :cigar "10M"}] [{:pos 1 :cigar "10M"}]]))
+  (is (= (map ->pile (plp/pileup-seq 1 2 (mapv ->aln [{:pos 1 :cigar "10M" :end 10}])))
+         [[1 [[1 10]]]
+          [2 [[1 10]]]]))
 
   ;;    ----------
   ;;   ----------
   ;;  ----------
   ;; ----------
   ;; 1234567890123...
-  (is (= (map count
-              (mplp/pileup-seq 1 20 (map #(hash-map :pos (inc %) :cigar "10M" :end (+ % 10)) (range))))
+  (is (= (map (comp count second)
+              (plp/pileup-seq 1 20 (map #(->aln (hash-map :pos (inc %) :cigar "10M" :end (+ % 10))) (range))))
          [1 2 3 4 5 6 7 8 9 10 10 10 10 10 10 10 10 10 10 10]))
-  (is (= (map count
-              (mplp/pileup-seq 101 120 (map #(hash-map :pos (inc %) :cigar "10M" :end (+ % 10)) (range))))
+  (is (= (map (comp count second)
+              (plp/pileup-seq 101 120 (map #(->aln (hash-map :pos (inc %) :cigar "10M" :end (+ % 10))) (range))))
          [10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10]))
   (is (= (map :pos
-              (last (mplp/pileup-seq 1 100000 (map #(hash-map :pos (inc %) :cigar "10M" :end (+ % 10)) (range)))))
+              (second (last (plp/pileup-seq 1 100000 (map #(->aln (hash-map :pos (inc %) :cigar "10M" :end (+ % 10))) (range))))))
          [99991 99992 99993 99994 99995 99996 99997 99998 99999 100000]))
 
   ;;     -----
@@ -62,8 +107,8 @@
   ;;  --
   ;; -
   ;; 1234567890...
-  (is (= (map count
-              (mplp/pileup-seq 1 10 (map #(hash-map :pos (inc %) :cigar (str (inc %) "M") :end (+ % (inc %))) (range))))
+  (is (= (map (comp count second)
+              (plp/pileup-seq 1 10 (map #(->aln (hash-map :pos (inc %) :cigar (str (inc %) "M") :end (+ % (inc %)))) (range))))
          [1 1 2 2 3 3 4 4 5 5]))
 
   ;;       --------
@@ -74,114 +119,204 @@
   ;;  --------
   ;; ----------
   ;; 1234567890...
-  (is (= (map count
-              (mplp/pileup-seq 1 10 (map #(hash-map :pos (inc %) :cigar (str (- 10 (* (mod % 5) 2)) "M") :end (+ % (- 10 (* (mod % 5) 2)))) (range))))
+  (is (= (map (comp count second)
+              (plp/pileup-seq 1 10 (map #(->aln (hash-map :pos (inc %) :cigar (str (- 10 (* (mod % 5) 2)) "M") :end (+ % (- 10 (* (mod % 5) 2))))) (range))))
          [1 2 3 4 5 6 6 6 6 6])))
 
-(deftest about-mpileup
-  (with-open [br (sam/bam-reader test-sorted-bam-file)
-              fr (cseq/fasta-reader test-fa-file)]
-    (let [mplp-ref (doall (plp/mpileup br {:chr "ref"}))
-          mplp-ref2 (doall (plp/mpileup br {:chr "ref2"}))]
-      (is (= (map :rname mplp-ref) (repeat 45 "ref")))
-      (is (= (map :ref mplp-ref) (repeat 45 \N)))
-      (is (= (map :count mplp-ref) test-bam-pileup-ref))
-      (is (= (map :pos mplp-ref) (range 1 46)))
-      (is (= (map :seq mplp-ref) test-bam-mpileup-seq-ref))
-      (is (= (map :qual mplp-ref) test-bam-mpileup-qual-ref))
-      (is (= (map :count mplp-ref2) test-bam-pileup-ref2))
-      (is (= (map :seq mplp-ref2) test-bam-mpileup-seq-ref2)))
-
-    (let [mplp-ref (doall (plp/mpileup fr br {:chr "ref"}))]
-      (is (= (map :rname mplp-ref) (repeat 45 "ref")))
-      ;;                                  123456789012345678901234567890123456789012345
-      (is (= (apply str (map :ref mplp-ref))
-             "AGCATGTTAGATAAGATAGCTGTGCTAGTAGGCAGTCAGCGCCAT"))
-      (is (= (map :count mplp-ref) test-bam-pileup-ref))
-      (is (= (map :pos mplp-ref) (range 1 46)))
-      (is (= (map :seq mplp-ref) test-bam-mpileup-seq-ref-with-ref))
-      (is (= (map :qual mplp-ref) test-bam-mpileup-qual-ref)))))
-
-(deftest mpileup-region
-  (with-open [br (sam/bam-reader test-sorted-bam-file)
-              fr (cseq/fasta-reader test-fa-file)]
-    ;; 1234567890123456789012345678901234567890
-    ;; aggttttataaaacaattaagtctacagagcaactacgcg
-    (let [mplp-ref1 (doall (plp/mpileup fr br {:chr "ref" :start 1 :end 40}))
-          mplp-ref2 (doall (plp/mpileup fr br {:chr "ref2" :start 1 :end 40}))]
-      (is (= (apply str (map :ref mplp-ref1))
-             "AGCATGTTAGATAAGATAGCTGTGCTAGTAGGCAGTCAGC"))
-      (is (= (map :count mplp-ref1) (take 40 test-bam-pileup-ref)))
-      (is (= (map :pos mplp-ref1) (range 1 41)))
-      (is (= (map :seq mplp-ref1) (take 40 test-bam-mpileup-seq-ref-with-ref)))
-      (is (= (apply str (map :ref mplp-ref2))
-             "AGGTTTTATAAAACAATTAAGTCTACAGAGCAACTACGCG")))))
-
-(def reads-for-pileup
-  [{:qname "R001" :flag 99 :rname "seq1" :pos 3 :end 10 :seq "AATTGGCCAA" :qual "AABBCCDDEE" :cigar "2S8M"
-    :rnext "=" :pnext 1 :tlen 8 :mapq 60}
-   {:qname "R001" :flag 147 :rname "seq1" :pos 3 :end 10 :seq "TTGGCCAATT" :qual "AABBCCDDEE" :cigar "8M2S"
-    :rnext "=" :pnext 3 :tlen -8 :mapq 20}])
+(deftest pileup-seq
+  (testing "empty"
+    (is (nil? (seq (plp/pileup-seq 1 10 []))))
+    (is (nil? (seq (plp/pileup-seq 2 10 (map ->aln [{:pos 1 :end 1}])))))
+    (is (nil? (seq (plp/pileup-seq 1 10 (map ->aln [{:pos 11 :end 11}]))))))
+  (testing "dense"
+    (are [?in ?out]
+        (= ?out (mapv ->pile (plp/pileup-seq 1 10 (mapv ->aln (map #(zipmap [:pos :end] %) ?in)))))
+      [[1 1]] [[1 [[1 1]]]]
+      [[3 5]] [[3 [[3 5]]] [4 [[3 5]]] [5 [[3 5]]]]
+      [[3 5] [4 4]] [[3 [[3 5]]] [4 [[3 5] [4 4]]] [5 [[3 5]]]]
+      [[3 5] [4 6]] [[3 [[3 5]]] [4 [[3 5] [4 6]]] [5 [[3 5] [4 6]]] [6 [[4 6]]]]))
+  (testing "sparse"
+    (are [?in ?out]
+        (= ?out (mapv ->pile (plp/pileup-seq 1 10 (mapv ->aln (map #(zipmap [:pos :end] %) ?in)))))
+      [[3 5] [14 16]] [[3 [[3 5]]] [4 [[3 5]]] [5 [[3 5]]]]
+      [[3 5] [6 7]] [[3 [[3 5]]] [4 [[3 5]]] [5 [[3 5]]] [6 [[6 7]]] [7 [[6 7]]]]
+      [[3 5] [4 5] [6 7]] [[3 [[3 5]]] [4 [[3 5] [4 5]]] [5 [[3 5] [4 5]]] [6 [[6 7]]] [7 [[6 7]]]]
+      [[3 5] [7 8] [10 11]] [[3 [[3 5]]] [4 [[3 5]]] [5 [[3 5]]] [7 [[7 8]]] [8 [[7 8]]] [10 [[10 11]]]])))
+
+(deftest align-pileup-seqs
+  (testing "empty"
+    (is (empty? (plp/align-pileup-seqs)))
+    (is (empty? (plp/align-pileup-seqs [])))
+    (is (empty? (plp/align-pileup-seqs [] []))))
+  (testing "single"
+    (are [?in ?out]
+        (= ?out (plp/align-pileup-seqs ?in))
+      [{:pos 3 :pile [{:pos 3 :end 5}]}]
+      [[3 [{:pos 3 :pile [{:pos 3 :end 5}]}]]]
+
+      [{:pos 3 :pile [[{:pos 3 :end 5}]]}
+       {:pos 4 :pile [[{:pos 3 :end 5}]]}
+       {:pos 5 :pile [[{:pos 3 :end 5}]]}]
+      [[3 [{:pos 3 :pile [[{:pos 3 :end 5}]]}]]
+       [4 [{:pos 4 :pile [[{:pos 3 :end 5}]]}]]
+       [5 [{:pos 5 :pile [[{:pos 3 :end 5}]]}]]]
+
+      [{:pos 3 :pile [{:pos 3 :end 5}]}
+       {:pos 4 :pile [{:pos 3 :end 5} {:pos 4 :end 4}]}
+       {:pos 5 :pile [{:pos 3 :end 5}]}
+       {:pos 8 :pile [{:pos 8 :end 8}]}]
+      [[3 [{:pos 3 :pile [{:pos 3 :end 5}]}]]
+       [4 [{:pos 4 :pile [{:pos 3 :end 5} {:pos 4, :end 4}]}]]
+       [5 [{:pos 5 :pile [{:pos 3 :end 5}]}]]
+       [8 [{:pos 8 :pile [{:pos 8 :end 8}]}]]]))
+  (testing "multiple"
+    (is (= [[3 [{:pos 3, :pile [{:pos 3, :end 5}]} nil]]
+            [4 [{:pos 4, :pile [{:pos 3, :end 5}]} {:pos 4, :pile [{:pos 4, :end 6}]}]]
+            [5 [{:pos 5, :pile [{:pos 3, :end 5}]} {:pos 5, :pile [{:pos 4, :end 6}]}]]
+            [6 [nil {:pos 6, :pile [{:pos 4, :end 6}]}]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 5}]}
+                                   {:pos 4 :pile [{:pos 3 :end 5}]}
+                                   {:pos 5 :pile [{:pos 3 :end 5}]}]
+                                  [{:pos 4 :pile [{:pos 4 :end 6}]}
+                                   {:pos 5 :pile [{:pos 4 :end 6}]}
+                                   {:pos 6 :pile [{:pos 4 :end 6}]}])))
+    (is (= [[3 [{:pos 3, :pile [{:pos 3, :end 5}]} nil]]
+            [4 [{:pos 4, :pile [{:pos 3, :end 5}]} nil]]
+            [5 [{:pos 5, :pile [{:pos 3, :end 5}]} nil]]
+            [8 [nil {:pos 8, :pile [{:pos 8, :end 9}]}]]
+            [9 [nil {:pos 9, :pile [{:pos 8, :end 9}]}]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 5}]}
+                                   {:pos 4 :pile [{:pos 3 :end 5}]}
+                                   {:pos 5 :pile [{:pos 3 :end 5}]}]
+                                  [{:pos 8 :pile [{:pos 8 :end 9}]}
+                                   {:pos 9 :pile [{:pos 8 :end 9}]}])))
+    (is (= [[3 [{:pos 3, :pile [{:pos 3, :end 4}]} nil {:pos 3, :pile [{:pos 3, :end 3}]} nil]]
+            [4 [{:pos 4, :pile [{:pos 3, :end 4}]} {:pos 4, :pile [{:pos 4, :end 4}]} nil nil]]]
+           (plp/align-pileup-seqs [{:pos 3 :pile [{:pos 3 :end 4}]}
+                                   {:pos 4 :pile [{:pos 3 :end 4}]}]
+                                  [{:pos 4 :pile [{:pos 4 :end 4}]}]
+                                  [{:pos 3 :pile [{:pos 3 :end 3}]}]
+                                  [])))))
+
+(deftest about-pileup
+  (testing "dense"
+    (with-open [br (sam/bam-reader test-sorted-bam-file)
+                fr (cseq/fasta-reader test-fa-file)]
+      (let [plp-ref (doall (plp/pileup br {:chr "ref"}))
+            plp-ref2 (doall (plp/pileup br {:chr "ref2"}))]
+        (is (= (filter pos? test-bam-pileup-ref)
+               (map (comp count :pile) plp-ref)))
+        (is (= (keep-indexed #(when (pos? %2) (inc %1)) test-bam-pileup-ref)
+               (map :pos plp-ref)))
+        (is (= (remove empty? test-bam-mpileup-seq-ref-freq)
+               (map #(frequencies (map :base (:pile %))) plp-ref)))
+        (is (= (keep #(seq (map qual/fastq-char->phred-byte %)) test-bam-mpileup-qual-ref)
+               (map #(map :qual (:pile %)) plp-ref)))
+        (is (= (filter pos? test-bam-pileup-ref2)
+               (map (comp count :pile) plp-ref2)))
+        (is (= (remove empty? test-bam-mpileup-seq-ref2-freq)
+               (map #(frequencies (map :base (:pile %))) plp-ref2)))))))
+
+(deftest mpileup
+  (with-open [r1 (sam/reader test-sorted-bam-file)
+              r2 (sam/reader r1)]
+    (let [mplp (plp/mpileup {:chr "ref"} r1 r2)]
+      (is (= (keep-indexed #(when (pos? %2) (inc %1)) test-bam-pileup-ref)
+             (map first mplp))))))
+
+(deftest pileup-region
+  (with-open [br (sam/bam-reader test-sorted-bam-file)]
+    (let [plp-ref1 (doall (plp/pileup br {:chr "ref" :start 1 :end 40}))
+          plp-ref2 (doall (plp/pileup br {:chr "ref2" :start 1 :end 40}))]
+      (is (= (filter pos? (take 40 test-bam-pileup-ref))
+             (map (comp count :pile) plp-ref1)))
+      (is (= (filter pos? (take 40 test-bam-pileup-ref2))
+             (map (comp count :pile) plp-ref2)))
+      (is (= (->> test-bam-pileup-ref
+                  (take 40)
+                  (map vector (range))
+                  (filter (comp pos? second))
+                  (map (comp inc first)))
+             (map :pos plp-ref1)))
+      (is (= (->> test-bam-pileup-ref2
+                  (take 40)
+                  (map vector (range))
+                  (filter (comp pos? second))
+                  (map (comp inc first)))
+             (map :pos plp-ref2)))
+      (is (= (filter seq (take 40 test-bam-mpileup-seq-ref-freq))
+             (map #(frequencies (map :base (:pile %))) plp-ref1)))
+      (is (= (filter seq (take 40 test-bam-mpileup-seq-ref2-freq))
+             (map #(frequencies (map :base (:pile %))) plp-ref2))))))
+
+(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"
+     :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"
+     :rnext "=" :pnext (int 3) :tlen (int -8) :mapq (int 20)}]))
+
+(defn- pileup* [region xs]
+  (plp/pileup
+   (reify p/IAlignmentReader
+     (p/read-refs [_]
+       [{:name "seq1",  :len 249250621}])
+     (p/read-alignments [_ _]
+       xs))
+   region))
 
 (deftest overlap-correction
   (let [plps (->> reads-for-pileup
-                  (filter mplp/basic-mpileup-pred)
-                  (mplp/pileup* "AATTGGCCAATTGGCC" {:chr "seq1" :start 1 :end 10})
-                  (map (comp mplp/transpose-pile mplp/correct-overlapped-reads first)))]
-    (is (= (map :ref plps) [\A \A \T \T \G \G \C \C \A \A]))
-    (is (= (map :count plps) [0 0 2 2 2 2 2 2 2 2]))
-    (is (= (map :seq plps) [[] [] [\. \.] [\. \.] [\. \.] [\. \.] [\. \.] [\. \.] [\. \.] [\. \.]]))
-    (is (= (map :qual plps) [[] [] [\b \!] [\b \!] [\d \!] [\d \!] [\f \!] [\f \!] [\h \!] [\h \!]]))))
+                  (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)))))
 
 (deftest overlap-qual-correction
-  (let [aplp (mplp/correct-overlapped-reads
-              {:pile
-               [{:seq \A :qual (- (int \I) 33) :read {:qname "R001" :flag 99}}
-                {:seq \T :qual (- (int \A) 33) :read {:qname "R001" :flag 147}}]})]
-    (is (= (map :seq (:pile aplp)) [\A \T]))
-    (is (= (map :qual (:pile aplp)) [32 0])))
-  (let [aplp (mplp/correct-overlapped-reads
-              {:pile
-               [{:seq \T :qual (- (int \A) 33) :read {:qname "R001" :flag 147}}
-                {:seq \A :qual (- (int \I) 33) :read {:qname "R001" :flag 99}}]})]
-    (is (= (map :seq (:pile aplp)) [\T \A]))
-    (is (= (map :qual (:pile aplp)) [0 32]))))
+  (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})})]])]
+    (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})})]])]
+    (is (= (map :base aplp) [\T \A]))
+    (is (= (map :qual aplp) [0 32]))))
 
 (deftest filter-by-base-quality
   (let [plps (->> reads-for-pileup
-                  (filter mplp/basic-mpileup-pred)
-                  (mplp/pileup* "AATTGGCCAATTGGCC" {:chr "seq1" :start 1 :end 10})
-                  (sequence
-                   (comp
-                    (map first)
-                    (map mplp/correct-overlapped-reads)
-                    (map (mplp/filter-by-base-quality 1))
-                    (map mplp/transpose-pile))))]
-    (is (= (map :ref plps) [\A \A \T \T \G \G \C \C \A \A]))
-    (is (= (map :count plps) [0 0 1 1 1 1 1 1 1 1]))
-    (is (= (map :seq plps) [[] [] [\.] [\.] [\.] [\.] [\.] [\.] [\.] [\.]]))
-    (is (= (map :qual plps) [[] [] [\b] [\b] [\d] [\d] [\f] [\f] [\h] [\h]]))))
+                  (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)))))
 
 (deftest filter-by-map-quality
   (let [plps (->> reads-for-pileup
-                  (sequence
-                   (comp
-                    (filter mplp/basic-mpileup-pred)
-                    (filter (fn [a] (<= 30 (:mapq a))))))
-                  (mplp/pileup* "AATTGGCCAATTGGCC" {:chr "seq1" :start 1 :end 10})
-                  (sequence
-                   (comp
-                    (map first)
-                    (map mplp/correct-overlapped-reads)
-                    (map mplp/transpose-pile))))]
-    (is (= (map :ref plps) [\A \A \T \T \G \G \C \C \A \A]))
-    (is (= (map :count plps) [0 0 1 1 1 1 1 1 1 1]))
-    (is (= (map :seq plps) [[] [] [\.] [\.] [\.] [\.] [\.] [\.] [\.] [\.]]))
-    (is (= (map :qual plps) [[] [] [\B] [\B] [\C] [\C] [\D] [\D] [\E] [\E]]))))
+                  (filter (fn [a] (<= 30 (:mapq a))))
+                  (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 [[] [] [33] [33] [34] [34] [35] [35] [36] [36]])
+           (map #(map :qual (:pile %)) plps)))))
 
 (deftest about-create-mpileup
   (with-before-after {:before (prepare-cache!)
                       :after (clean-cache!)}
     (let [out-file (str temp-dir "/test.pileup")]
       (is (not-throw? (plp/create-mpileup test-sorted-bam-file out-file)))
-      (is (same-file? out-file test-pileup-file)))))
+      (doseq [[r1 r2] (->> [test-pileup-file out-file]
+                           (map (comp cstr/split-lines slurp))
+                           (apply map vector))
+              :when (not (cstr/starts-with? r1 "ref\t35\t"))] ;; adjacent indel
+        (is (= r1 r2))))))
diff --git a/test/cljam/io/pileup_test.clj b/test/cljam/io/pileup_test.clj
new file mode 100644
index 0000000..7ef0cf4
--- /dev/null
+++ b/test/cljam/io/pileup_test.clj
@@ -0,0 +1,319 @@
+(ns cljam.io.pileup-test
+  (:require [clojure.test :refer :all]
+            [clojure.java.io :as cio]
+            [clojure.string :as cstr]
+            [clojure.walk :as walk]
+            [cljam.test-common :refer :all]
+            [cljam.io.pileup :as plpio]
+            [cljam.io.sequence :as cseq]
+            [cljam.io.fasta-index.core :as fai]
+            [cljam.io.protocols :as p])
+  (:import [java.io File StringWriter StringReader]
+           [cljam.io.pileup PileupReader PileupWriter]))
+
+(defn- rec->map [coll]
+  (walk/postwalk
+   (fn [x]
+     (if (record? x)
+       (into {} x)
+       x))
+   coll))
+
+(def ^:private ^:const the-last-pile
+  ;; ref2	36	N	1	A$	?
+  {:rname "ref2",
+   :pos 36,
+   :ref \N,
+   :count 1,
+   :pile
+   [{:start? false,
+     :mapq nil,
+     :base \A,
+     :qual 30,
+     :reverse? false,
+     :end? true,
+     :insertion nil,
+     :deletion nil,
+     :qname nil}]})
+
+;; Reader
+;; ------
+
+(deftest parse-bases-col
+  (testing "without-ref"
+    (are [?in ?out]
+        (= (mapv #(into {:qname nil} %) ?out)
+           (mapv rec->map (#'plpio/parse-bases-col nil ?in)))
+      "" []
+      "A" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "a" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]
+      "^]A" [{:start? true, :mapq 60, :base \A, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "a$" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? true, :end? true, :insertion nil, :deletion nil}]
+      "A$C" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion nil}
+             {:start? false, :mapq nil, :base \C, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "A$^]C" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion nil}
+               {:start? true, :mapq 60, :base \C, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "^BN" [{:start? true, :mapq 33, :base \N, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      ">" [{:start? false, :mapq nil, :base \>, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "*" [{:start? false, :mapq nil, :base \*, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "<" [{:start? false, :mapq nil, :base \>, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]
+      "A+2TG" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? false, :insertion "TG", :deletion nil}]
+      "a-12nnnnnnnnnnnn$c$" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? true, :end? true, :insertion nil, :deletion 12}
+                             {:start? false, :mapq nil, :base \C, :qual -1, :reverse? true, :end? true, :insertion nil, :deletion nil}]
+      "A+14AAAAATTTTTGGGG" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? false, :insertion "AAAAATTTTTGGGG", :deletion nil}]
+      "^]A+14AAAAATTTTTGGGG" [{:start? true, :mapq 60, :base \A, :qual -1, :reverse? false, :end? false, :insertion "AAAAATTTTTGGGG", :deletion nil}]
+      "A+14AAAAATTTTTGGGG$" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? true, :insertion "AAAAATTTTTGGGG", :deletion nil}]
+      "A+14AAAAATTTTTGGGGG" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? false, :insertion "AAAAATTTTTGGGG", :deletion nil}
+                             {:start? false, :mapq nil, :base \G, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "A-2NN$" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion 2}]
+      "A-2NN$^Ca" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion 2}
+                   {:start? true, :mapq 34, :base \A, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]))
+  (testing "with-ref"
+    (are [?in ?out]
+        (= (mapv #(into {:qname nil} %) ?out)
+           (mapv rec->map (#'plpio/parse-bases-col \T ?in)))
+      "" []
+      "A" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      "a$" [{:start? false, :mapq nil, :base \A, :qual -1, :reverse? true, :end? true, :insertion nil, :deletion nil}]
+      "," [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]
+      ".$" [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion nil}]
+      ".$^]," [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion nil}
+               {:start? true, :mapq 60, :base \T, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]
+      "^B." [{:start? true, :mapq 33, :base \T, :qual -1, :reverse? false, :end? false, :insertion nil, :deletion nil}]
+      ",+2tg" [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? true, :end? false, :insertion "TG", :deletion nil}]
+      ".-2AT$" [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? false, :end? true, :insertion nil, :deletion 2}]
+      ",-2at$^Ca" [{:start? false, :mapq nil, :base \T, :qual -1, :reverse? true, :end? true, :insertion nil, :deletion 2}
+                   {:start? true, :mapq 34, :base \A, :qual -1, :reverse? true, :end? false, :insertion nil, :deletion nil}]))
+  (testing "invalid arguments"
+    (are [?in]
+        (thrown? Exception (#'plpio/parse-bases-col nil ?in))
+      "A+"
+      "^"
+      "^A"
+      "^A+1A"
+      "+"
+      "-"
+      "++"
+      "--"
+      "A$$A"
+      "A+2A^]A"
+      "A+2A"
+      "A+2"
+      "A-2N"
+      "A-2")))
+
+(deftest parse-pileup-line
+  (testing "without-ref"
+    (are [?in ?out]
+        (= ?out (rec->map (#'plpio/parse-pileup-line ?in)))
+      "chr1\t10\tN\t0\t\t" {:rname "chr1", :pos 10, :ref \N, :count 0,
+                            :pile []}
+      "chr1\t10\tN\t1\tA\tI" {:rname "chr1", :pos 10, :ref \N, :count 1,
+                              :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil}]}
+      "chr1\t10\tN\t2\taA\tIB" {:rname "chr1", :pos 10, :ref \N, :count 2,
+                                :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil}
+                                       {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil}]}
+      "chr1\t10\tN\t2\taA-2NN$\tIB" {:rname "chr1", :pos 10, :ref \N, :count 2,
+                                     :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil}
+                                            {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? true, :insertion nil, :deletion 2, :qname nil}]}))
+  (testing "with-ref"
+    (are [?in ?out]
+        (= ?out (rec->map (#'plpio/parse-pileup-line ?in)))
+      "chr1\t10\tA\t0\t\t" {:rname "chr1", :pos 10, :ref \A, :count 0,
+                            :pile []}
+      "chr1\t10\ta\t1\t.\tI" {:rname "chr1", :pos 10, :ref \a, :count 1,
+                              :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil}]}
+      "chr1\t10\tA\t2\t,.\tIB" {:rname "chr1", :pos 10, :ref \A, :count 2,
+                                :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil}
+                                       {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil}]}
+      "chr1\t10\tA\t2\t,.-2CA$\tIB" {:rname "chr1", :pos 10, :ref \A, :count 2,
+                                     :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil}
+                                            {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? true, :insertion nil, :deletion 2, :qname nil}]})))
+
+(deftest reader
+  (with-open [r (plpio/reader test-pileup-file)]
+    (is (instance? PileupReader r))))
+
+(deftest read-piles
+  (with-open [r (plpio/reader test-pileup-file)]
+    (is (= the-last-pile
+           (rec->map (last (plpio/read-piles r)))))))
+
+;; Writer
+;; ------
+
+(defmacro with-string-writer [symbol & exprs]
+  `(with-open [sw# (StringWriter.)
+               ~symbol (cio/writer sw#)]
+     ~@exprs
+     (.flush ~symbol)
+     (str sw#)))
+
+(deftest write-mpileup-alignment!
+  (testing "without-ref"
+    (are [?in ?out]
+        (= ?out (with-string-writer w
+                  (#'plpio/write-mpileup-alignment! w nil "chr1" 10 nil ?in)))
+      {:base \A :reverse? false :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A"
+      {:base \N :reverse? false :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "N"
+      {:base \A :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a"
+      {:base \N :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "n"
+      {:base \* :reverse? false :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "*"
+      {:base \* :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "*"
+      {:base \> :reverse? false :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} ">"
+      {:base \> :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "<"
+      {:base \A :reverse? false :insertion "A" :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A+1A"
+      {:base \A :reverse? false :deletion 2 :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A-2NN"
+      {:base \A :reverse? true :insertion "AT" :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a+2at"
+      {:base \A :reverse? true :insertion "AAAATTTTGGGGCCCC" :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a+16aaaattttggggcccc"
+      {:base \A :reverse? true :deletion 3 :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a-3nnn"
+      {:base \A :reverse? true :deletion 10 :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a-10nnnnnnnnnn"
+      {:base \A :reverse? false :start? true :mapq 60 :alignment {:flag 0 :mapq 60 :pos 10 :end 15}} "^]A"
+      {:base \A :reverse? false :start? true :mapq 93 :alignment {:flag 0 :mapq 93 :pos 10 :end 15}} "^~A"
+      {:base \A :reverse? false :start? true :mapq 94 :alignment {:flag 0 :mapq 94 :pos 10 :end 15}} "^~A"
+      {:base \A :reverse? true :start? true :mapq 40 :insertion "AT" :alignment {:flag 16 :mapq 40 :pos 10 :end 15}} "^Ia+2at"
+      {:base \A :reverse? false :end? true :alignment {:flag 0 :mapq 40 :pos 5 :end 10}} "A$"
+      {:base \A :reverse? true :end? true :deletion 4 :alignment {:flag 16 :mapq 40 :pos 5 :end 10}} "a-4nnnn$"))
+  (testing "with-ref"
+    (let [r (reify p/ISequenceReader
+              (p/read-sequence [this {:keys [start end]}]
+                (subs "ATGCATGCATGCATGCATGCATGCATGC" (dec start) end)))]
+      (are [?in ?out]
+          (= ?out (with-string-writer w
+                    (#'plpio/write-mpileup-alignment! w r "chr1" 10 \T ?in)))
+        {:base \A :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A"
+        {:base \T :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "."
+        {:base \T :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} ","
+        {:base \N :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "N"
+        {:base \A :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a"
+        {:base \N :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "n"
+        {:base \* :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "*"
+        {:base \* :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "*"
+        {:base \> :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} ">"
+        {:base \> :reverse? true :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "<"
+        {:base \A :insertion "A" :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A+1A"
+        {:base \A :deletion 2 :alignment {:flag 0 :mapq 60 :pos 5 :end 15}} "A-2GC"
+        {:base \A :reverse? true :insertion "AT" :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a+2at"
+        {:base \A :reverse? true :insertion "AAAATTTTGGGGCCCC" :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a+16aaaattttggggcccc"
+        {:base \A :reverse? true :deletion 3 :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a-3gca"
+        {:base \A :reverse? true :deletion 10 :alignment {:flag 16 :mapq 60 :pos 5 :end 15}} "a-10gcatgcatgc"
+        {:base \A :start? true :mapq 60 :alignment {:flag 0 :mapq 60 :pos 10 :end 15}} "^]A"
+        {:base \A :start? true :mapq 93 :alignment {:flag 0 :mapq 93 :pos 10 :end 15}} "^~A"
+        {:base \A :start? true :mapq 94 :alignment {:flag 0 :mapq 94 :pos 10 :end 15}} "^~A"
+        {:base \A :reverse? true :start? true :mapq 40 :insertion "AT" :alignment {:flag 16 :mapq 40 :pos 10 :end 15}} "^Ia+2at"
+        {:base \A :end? true :alignment {:flag 0 :mapq 40 :pos 5 :end 10}} "A$"
+        {:base \A :reverse? true :end? true :deletion 4 :alignment {:flag 16 :mapq 40 :pos 5 :end 10}} "a-4gcat$"))))
+
+(defn string-sequence-reader [s]
+  (reify p/ISequenceReader
+    (p/read-sequence [this region]
+      (p/read-sequence this region {}))
+    (p/read-sequence [this {:keys [start end]} {:keys [mask?]}]
+      ((if mask? identity cstr/upper-case)
+       (subs s (dec start) end)))))
+
+(deftest write-mpileup-line!
+  (testing "without-ref"
+    (are [?in ?out]
+        (= ?out (with-string-writer w
+                  (#'plpio/write-mpileup-line! w nil {:rname (first ?in)
+                                                      :pos (second ?in)
+                                                      :pile (last ?in)})))
+      ["chr1" 10 []] "chr1\t10\tN\t0\t\t"
+      ["chr1" 10 [{:base \A :qual 40 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tN\t1\tA\tI"
+      ["chr1" 10 [{:base \A :qual 93 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tN\t1\tA\t~"
+      ["chr1" 10 [{:base \A :qual 94 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tN\t1\tA\t~"
+      ["chr1" 10 [{:base \A :qual 40 :alignment {:flag 0 :pos 5}}
+                  {:base \A :qual 33 :reverse? true :start? true :mapq 60 :alignment {:flag 16 :mapq 60 :pos 10}}]] "chr1\t10\tN\t2\tA^]a\tIB"))
+  (testing "with-ref"
+    (let [r (string-sequence-reader "NNNNNNNNNAtGCATGCAT")]
+      (are [?in ?out]
+          (= ?out (with-string-writer w
+                    (#'plpio/write-mpileup-line! w r {:rname (first ?in)
+                                                      :pos (second ?in)
+                                                      :pile (last ?in)})))
+        ["chr1" 10 []] "chr1\t10\tA\t0\t\t"
+        ["chr1" 10 [{:base \A :qual 40 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tA\t1\t.\tI"
+        ["chr1" 10 [{:base \A :qual 93 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tA\t1\t.\t~"
+        ["chr1" 10 [{:base \A :qual 94 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tA\t1\t.\t~"
+        ["chr1" 10 [{:base \A :insertion "AA" :qual 40 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tA\t1\t.+2AA\tI"
+        ["chr1" 10 [{:base \A :deletion 2 :qual 40 :alignment {:flag 0 :pos 5}}]] "chr1\t10\tA\t1\t.-2TG\tI"
+        ["chr1" 11 [{:base \T :qual 40 :alignment {:flag 0 :pos 5}}]] "chr1\t11\tt\t1\t.\tI"
+        ["chr1" 10 [{:base \A :qual 40 :alignment {:flag 0 :pos 5}}
+                    {:base \A :qual 33 :reverse? true :start? true :mapq 60 :alignment {:flag 16 :mapq 60 :pos 10}}]] "chr1\t10\tA\t2\t.^],\tIB"))))
+
+(deftest writer
+  (let [tmp (File/createTempFile "writer-test" ".mpileup")]
+    (try
+      (with-open [w (plpio/writer tmp)]
+        (is (instance? PileupWriter w)))
+      (finally
+        (when (.isFile (cio/file tmp))
+          (cio/delete-file (cio/file tmp)))))))
+
+(deftest write-piles
+  (is (= "chr1\t10\tN\t1\tA\tI\n"
+         (with-open [sw (StringWriter.)
+                     w (plpio/writer sw)]
+           (plpio/write-piles w [{:rname "chr1"
+                                  :pos 10
+                                  :pile [{:base \A :qual 40 :alignment {:flag 0 :pos 5}}]}])
+           (.toString sw)))))
+
+
+;; Read & Write
+;; ------------
+
+(deftest regression
+  (testing "without-ref"
+    (are [?input]
+        (= ?input
+           (with-open [r (plpio/reader (StringReader. ?input))
+                       sw (StringWriter.)
+                       w (plpio/writer sw)]
+             (plpio/write-piles w (plpio/read-piles r))
+             (str sw)))
+      "chr1\t10\tN\t1\tA\tI\n"
+      "chr1\t10\tN\t4\tAaTt\tIABC\n"
+      "chr1\t10\tN\t4\t^]A+3TTTa-2nn$Tt\tIABC\n"))
+  (testing "with-ref"
+    (are [?input]
+        (= ?input
+           (let [tmp (File/createTempFile "pileup-regression" ".fa")
+                 idx (File/createTempFile "pileup-regression" ".fai")]
+             (try
+               (with-open [w (cseq/writer (.getCanonicalPath tmp))]
+                 (cseq/write-sequences w [{:name "chr1" :sequence "NNNNNNNNNATGC"}]))
+               (fai/create-index tmp idx)
+               (with-open [r (plpio/reader (StringReader. ?input))
+                           sw (StringWriter.)
+                           w (plpio/writer sw tmp)]
+                 (plpio/write-piles w (plpio/read-piles r))
+                 (str sw))
+               (finally
+                 (when (.isFile (cio/file tmp))
+                   (cio/delete-file (cio/file tmp)))))))
+      "chr1\t10\tA\t1\t.\tI\n"
+      "chr1\t10\tA\t4\t.,Tt\tIABC\n"
+      "chr1\t10\tA\t4\t^].+3TTT,-2tg$Tt\tIABC\n")))
+
+(deftest source-type-test
+  (testing "reader"
+    (with-open [server (http-server)]
+      (are [x] (= the-last-pile
+                  (with-open [r (plpio/reader x)]
+                    (rec->map (last (plpio/read-piles r)))))
+        test-pileup-file
+        (cio/file test-pileup-file)
+        (cio/as-url (cio/file test-pileup-file))
+        (cio/as-url (str (:uri server) "/pileup/test.pileup")))))
+
+  (testing "writer"
+    (let [tmp-pileup-file (cio/file temp-dir "pileup-source-type-writer.mpileup")]
+      (are [x] (with-before-after {:before (prepare-cache!)
+                                   :after (clean-cache!)}
+                 (with-open [w (plpio/writer x)]
+                   (not-throw? (plpio/write-piles w [the-last-pile]))))
+        (.getCanonicalPath tmp-pileup-file)
+        tmp-pileup-file
+        (cio/as-url tmp-pileup-file)))))
diff --git a/test/cljam/io/sam/util/cigar_test.clj b/test/cljam/io/sam/util/cigar_test.clj
index 34dd987..d060a59 100644
--- a/test/cljam/io/sam/util/cigar_test.clj
+++ b/test/cljam/io/sam/util/cigar_test.clj
@@ -6,34 +6,34 @@
   (is (= (cigar/parse  "1S2I6M1P11I") '([1 \S] [2 \I] [6 \M] [1 \P] [11 \I]))))
 
 (deftest cigar-to-index
-  (are [?cigar ?index] (= (map (fn [[op x :as xs]] (if (= op :m) x xs)) (cigar/to-index* ?cigar)) ?index)
+  (are [?cigar ?index] (= ?index (map (fn [[op x :as xs]] (if (= op :m) x xs)) (cigar/to-index* ?cigar)))
     "4M"                 [0 1 2 3]
-    "1M3I"               [[:i 0 [1 2 3]]]
-    "1M3D"               [[:d 0 [1 2 3]] \* \* \*]
-    "2M3I"               [0 [:i 1 [2 3 4]]]
+    "1M3I"               [[:i 0 [1 4]]]
+    "1M3D"               [[:d 0 3] \* \* \*]
+    "2M3I"               [0 [:i 1 [2 5]]]
     "4D5M"               [\* \* \* \* 0 1 2 3 4] ;; ^]* * * * A T G C A$
-    "4D5I"               [\* \* \* [:i \* [0 1 2 3 4]]] ;; TODO: ^]* * * *+5TGCA=$
+    "4D5I"               [\* \* \* [:i \* [0 5]]] ;; TODO: ^]* * * *+5TGCA=$
     "4S4M"               [4 5 6 7]
     "4H4M"               [0 1 2 3]
-    "1M3I1M"             [[:i 0 [1 2 3]] 4]
-    "2M3I1M"             [0 [:i 1 [2 3 4]] 5]
+    "1M3I1M"             [[:i 0 [1 4]] 4]
+    "2M3I1M"             [0 [:i 1 [2 5]] 5]
     "4I2D5M"             [\* \* 4 5 6 7 8] ;; ^]* * A T G C A$
-    "1M4D5I"             [[:d 0 [1 2 3 4]] \* \* \* [:i \* [1 2 3 4 5]]] ;; TODO: ^]* * * *+5TGCA=$
-    "1M3D3M"             [[:d 0 [1 2 3]] \* \* \* 1 2 3]
-    "1M2D1I2M"           [[:d 0 [1 2]] \* [:i \* [1]] 2 3] ;; TODO: ^]A-2NN * *+1G G C$
-    "1M2N1I2M"           [0 \> [:i \> [1]] 2 3] ;; TODO: ^]A > >+1G G C$
-    "1M1I2D2M"           [[:i 0 [1]] \* \* 2 3] ;; ^]A+1T * * G C$
-    "1M1I2N2M"           [[:i 0 [1]] \> \> 2 3] ;; ^]A+1T > > G C$
-    "1M4I2D5M"           [[:i 0 [1 2 3 4]] \* \* 5 6 7 8 9] ;; ^]A+4TGCA * * C A T G C$
+    "1M4D5I"             [[:d 0 4] \* \* \* [:i \* [1 6]]] ;; TODO: ^]* * * *+5TGCA=$
+    "1M3D3M"             [[:d 0 3] \* \* \* 1 2 3]
+    "1M2D1I2M"           [[:d 0 2] \* [:i \* [1 2]] 2 3] ;; TODO: ^]A-2NN * *+1G G C$
+    "1M2N1I2M"           [0 \> [:i \> [1 2]] 2 3] ;; TODO: ^]A > >+1G G C$
+    "1M1I2D2M"           [[:i 0 [1 2]] \* \* 2 3] ;; ^]A+1T * * G C$
+    "1M1I2N2M"           [[:i 0 [1 2]] \> \> 2 3] ;; ^]A+1T > > G C$
+    "1M4I2D5M"           [[:i 0 [1 5]] \* \* 5 6 7 8 9] ;; ^]A+4TGCA * * C A T G C$
     "1S4I2D5M"           [\* \* 5 6 7 8 9] ;; ^]* * T G C A T$
-    "3M2P2I3M"           [0 1 [:i 2 [3 4]] 5 6 7] ;; ^]A T G+2CA T G C$
+    "3M2P2I3M"           [0 1 [:i 2 [3 5]] 5 6 7] ;; ^]A T G+2CA T G C$
     "1P4I1D5M"           [\* 4 5 6 7 8] ;; ^]* A T G C A$
-    "5M1D3M4S"           [0 1 2 3 [:d 4 [5]] \* 5 6 7] ;; ^]A T G C A-1N * T G C$
+    "5M1D3M4S"           [0 1 2 3 [:d 4 1] \* 5 6 7] ;; ^]A T G C A-1N * T G C$
     "5M1N3M4S"           [0 1 2 3 4 \> 5 6 7] ;; ^]A T G C A > T G C$
-    "3H2S3M1D2S4H"       [2 3 [:d 4 [3]] \*] ;; ^]G C A-1N \*$
-    "1M1P1I1P1I2M"       [[:i 0 [1 2]] 3 4] ;; ^]A+2TG C A$
-    "2M1I1D2I2D1M"       [0 [:i 1 [2]] [:i \* [3 4]] \* \* 5] ;; TODO: ^]A T+1G *+2AT * * T$
-    "1S2I6M1P1I1P1I4M2I" [3 4 5 6 7 [:i 8 [9 10]] 11 12 13 [:i 14 [15 16]]] ;; ^]C A T G C A+2TG C A T G+2CA$
-    "6M14D1I5M"          [0 1 2 3 4 [:d 5 [6 7 8 9 10 11 12 13 14 15 16 17 18 19]] \* \* \* \* \* \* \* \* \* \* \* \* \* [:i \* [6]] 7 8 9 10 11]  ;; ^]A T G C A T-14NNNNNNNNNNNNNN * * * * * * * * * * * * * *+1C C A T G C$
-    "6M14N1I5M"          [0 1 2 3 4 5 \> \> \> \> \> \> \> \> \> \> \> \> \> [:i \> [6]] 7 8 9 10 11] ;; ^]A T G C A T > > > > > > > > > > > > > >+1C C A T G C$
-    "1S2M2D3M1D4M3D2M3D3M3I4M3D2M3S" [1 [:d 2 [2 3]] \* \* 3 4 [:d 5 [7]] \* 6 7  8 [:d 9 [12 13 14]] \* \* \* 10 [:d 11 [17 18 19]] \* \* \* 12 13 [:i 14 [15 16 17]] 18 19 20 [:d 21 [27 28 29]] \* \* \* 22 23]))
+    "3H2S3M1D2S4H"       [2 3 [:d 4 1] \*] ;; ^]G C A-1N \*$
+    "1M1P1I1P1I2M"       [[:i 0 [1 3]] 3 4] ;; ^]A+2TG C A$
+    "2M1I1D2I2D1M"       [0 [:i 1 [2 3]] [:i \* [3 5]] \* \* 5] ;; TODO: ^]A T+1G *+2AT * * T$
+    "1S2I6M1P1I1P1I4M2I" [3 4 5 6 7 [:i 8 [9 11]] 11 12 13 [:i 14 [15 17]]] ;; ^]C A T G C A+2TG C A T G+2CA$
+    "6M14D1I5M"          [0 1 2 3 4 [:d 5 14] \* \* \* \* \* \* \* \* \* \* \* \* \* [:i \* [6 7]] 7 8 9 10 11]  ;; ^]A T G C A T-14NNNNNNNNNNNNNN * * * * * * * * * * * * * *+1C C A T G C$
+    "6M14N1I5M"          [0 1 2 3 4 5 \> \> \> \> \> \> \> \> \> \> \> \> \> [:i \> [6 7]] 7 8 9 10 11] ;; ^]A T G C A T > > > > > > > > > > > > > >+1C C A T G C$
+    "1S2M2D3M1D4M3D2M3D3M3I4M3D2M3S" [1 [:d 2 2] \* \* 3 4 [:d 5 1] \* 6 7  8 [:d 9 3] \* \* \* 10 [:d 11 3] \* \* \* 12 13 [:i 14 [15 18]] 18 19 20 [:d 21 3] \* \* \* 22 23]))