Repository: chrovis/cljam
Pull request #165 (2019-05-20 07:47:32)

Improve pileup more
#### Summary
Improve the performance of pile up

#### Changes
- Add `alignment` to `PileupBase` to remove `assoc`
- Improve `correct-pair-quals` by using ~~transient vectors~~ a transducer
- ~~Improve `index-cigar` by making `seqs-at-ref` and `quals-at-ref` at once~~
  -　(Tips) I've tried to make new `index-cigar` for pileup to make them without index, however, it's not effective. That's because these index made by `(memoize to-index*)` makes this process fast.
- Call `clojure.lang.Numbers#bitTest` for `long` directly to avoid many casts.
- Add some type hints

#### Tests
- `lein check` 🆗
- `lein test :all` 🆗
- `lein cloverage` 🆗  

#### Benchmark
- sample: WES SRP006438.sorted.bam
  
    - many alignments chr1:152210000-152214999  (improved by 15.7%)

    |branch|time|
    |----|----|
    |master|time: 7.490895 sec, sd: 46.194343 ms|
    |feature/improve-pileup (new)|time: 6.313420 sec, sd: 351.617645 µs|


    - gene: NRAS chr1:114703000-114723000  (improved by 23.7%)
    
    |branch|time|
    |----|----|
    |master|time: 1.706070 sec, sd: 19.496387 µs|
    |feature/improve-pileup (new)|time: 1.302210 sec, sd: 8.227728 µs|

    - gene: EGFR chr7:55015000-55215000  (improved by 23.2%)
    
    |branch|time|
    |----|----|
    |master|time: 4.215258 sec, sd: 47.407546 µs|
    |feature/improve-pileup (new)|time: 3.235970 sec, sd: 25.867269 µs|


Patch:
diff --git a/src/cljam/algo/pileup.clj b/src/cljam/algo/pileup.clj
index 2b880fc..d6af593 100644
--- a/src/cljam/algo/pileup.clj
+++ b/src/cljam/algo/pileup.clj
@@ -84,18 +84,19 @@
   [^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)
-         (.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))))
+        [base indel] ((:seqs-at-ref aln) relative-pos)
+        deletion? (number? indel)]
+    (plpio/->PileupBase
+     (zero? relative-pos)
+     (.mapq aln)
+     base
+     qual
+     (flag/reversed? (.flag aln))
+     (= ref-pos (.end aln))
+     (when-not deletion? indel)
+     (when deletion? indel)
+     (.qname aln)
+     aln)))
 
 (defn- resolve-bases
   [[ref-pos alns]]
@@ -110,12 +111,30 @@
 (defn filter-by-base-quality
   "Returns a predicate for filtering piled-up reads by base quality at its
   position."
-  [min-base-quality]
+  [^long min-base-quality]
   (fn [p]
     (->> #(<= min-base-quality (.qual ^PileupBase %))
          (partial filterv)
          (update p 1))))
 
+(defn- unzip-2
+  [transform-fn]
+  (fn [rf]
+    (let [va (volatile! (transient []))
+          vb (volatile! (transient []))]
+      (fn
+        ([] (rf))
+        ([acc]
+         (-> acc
+             (rf (persistent! @va))
+             (rf (persistent! @vb))
+             rf))
+        ([acc x]
+         (let [[a b] (transform-fn x)]
+           (vswap! va conj! a)
+           (vswap! vb conj! b)
+           acc))))))
+
 (defn- merge-corrected-quals
   "Merge corrected quals with the uncorrected part."
   [^SAMAlignment aln ^long correct-start corrected-quals]
@@ -129,37 +148,46 @@
           (into (subvec quals (+ (- correct-start start) len))))
       (into corrected-quals (subvec quals len)))))
 
+(defn- correct-pair-qual
+  [^SAMAlignment a1 ^SAMAlignment a2]
+  (let [pos1 (.pos a1)
+        pos2 (.pos a2)
+        quals1 (:quals-at-ref a1)
+        quals2 (:quals-at-ref a2)
+        seqs1  (:seqs-at-ref a1)
+        seqs2  (:seqs-at-ref a2)]
+    (fn [^long pos]
+      (let [relative-pos1 (- pos pos1)
+            relative-pos2 (- pos pos2)
+            ^int q1 (quals1 relative-pos1)
+            ^int q2 (quals2 relative-pos2)
+            [b1] (seqs1 relative-pos1)
+            [b2] (seqs2 relative-pos2)]
+        (if (= b1 b2)
+          [(min 200 (+ q1 q2)) 0]
+          (if (<= q2 q1)
+            [(int (* 0.8 q1)) 0]
+            [0 (int (* 0.8 q2))]))))))
+
 (defn- correct-pair-quals
   "Correct quals of a pair. Returns a map with corrected quals."
   [^SAMAlignment a1 ^SAMAlignment a2]
-  (when (and (pos? (.pnext a1))
-             (<= (max (.pos a1) (.pos a2)) (.end a1)))
-    (let [tlen1 (.tlen a1)
-          tlen2 (.tlen a2)
-          quals1 (:quals-at-ref a1)
-          quals2 (:quals-at-ref a2)
-          seqs1  (:seqs-at-ref a1)
-          seqs2  (:seqs-at-ref a2)
-          correct-start (max (.pos a1) (.pos a2))
-          new-quals (for [pos (range correct-start
-                                     (inc (min (.end a1) (.end a2))))]
-                      (let [relative-pos1 (- pos (.pos a1))
-                            relative-pos2 (- pos (.pos a2))
-                            q1 (quals1 relative-pos1)
-                            q2 (quals2 relative-pos2)
-                            [b1] (seqs1 relative-pos1)
-                            [b2] (seqs2 relative-pos2)]
-                        (if (= b1 b2)
-                          [(min 200 (+ q1 q2)) 0]
-                          (if (<= q2 q1)
-                            [(int (* 0.8 q1)) 0]
-                            [0 (int (* 0.8 q2))]))))
-          [new1 new2] (apply map vector new-quals)
-          new-quals1 (merge-corrected-quals a1 correct-start new1)
-          new-quals2 (merge-corrected-quals a2 correct-start new2)]
-      (if (flag/r1? (.flag a1))
-        [new-quals1 new-quals2]
-        [new-quals2 new-quals1]))))
+  (let [correct-start (max (.pos a1) (.pos a2))
+        correct-end (min (.end a1) (.end a2))]
+    (when (and (pos? (.pnext a1))
+               (<= correct-start (.end a1)))
+      (let [[quals1 quals2] (into []
+                                  (unzip-2 (correct-pair-qual a1 a2))
+                                  (range correct-start (inc correct-end)))
+            new-quals1 (merge-corrected-quals a1
+                                              correct-start
+                                              quals1)
+            new-quals2 (merge-corrected-quals a2
+                                              correct-start
+                                              quals2)]
+        (if (flag/r1? (.flag a1))
+          [new-quals1 new-quals2]
+          [new-quals2 new-quals1])))))
 
 (defn- make-corrected-quals-map
   "Make a map which has corrected quals of all overlapping pairs."
@@ -196,11 +224,11 @@
   ([sam-reader region]
    (pileup sam-reader region {}))
   ([sam-reader
-    {:keys [chr start end] :or {start 1 end Integer/MAX_VALUE}}
-    {:keys [min-base-quality min-map-quality ignore-overlaps? chunk-size]
+    {:keys [chr ^long start ^long end] :or {start 1 end Integer/MAX_VALUE}}
+    {:keys [^long min-base-quality ^long min-map-quality ignore-overlaps? ^long chunk-size]
      :or {min-base-quality 13 min-map-quality 0 ignore-overlaps? false
           chunk-size 5000}}]
-   (when-let [len (:len (refs/ref-by-name (sam/read-refs sam-reader) chr))]
+   (when-let [^long len (:len (refs/ref-by-name (sam/read-refs sam-reader) chr))]
      (let [s (max 1 start)
            e (min len end)
            region {:chr chr :start s :end e}]
@@ -212,7 +240,7 @@
             (seq-step start end chunk-size)
             (sequence
              (comp
-              (mapcat (fn [[pos alns]]
+              (mapcat (fn [[^long pos alns]]
                         (->> (if ignore-overlaps?
                                alns
                                (keep (partial correct-quals-at-ref
diff --git a/src/cljam/io/pileup.clj b/src/cljam/io/pileup.clj
index d8e6410..f93317b 100644
--- a/src/cljam/io/pileup.clj
+++ b/src/cljam/io/pileup.clj
@@ -17,7 +17,8 @@
                        ^boolean end?
                        insertion ;; String?
                        deletion ;; int?
-                       qname]) ;; String?
+                       qname ;; String?
+                       alignment]) ;; SAMAlignment?
 
 (defrecord LocusPile [rname
                       ^int pos
@@ -67,7 +68,7 @@
           (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)))
+            (persistent! (conj! results (PileupBase. start? mapq upper-base -1 reverse? false nil nil nil nil)))
             (let [x (.charAt column (+ i (if mapq 3 1)))
                   ins? (= x \+)
                   del? (= x \-)
@@ -90,7 +91,7 @@
                                     (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))))))
+                     (conj! results (PileupBase. start? mapq upper-base -1 reverse? end? (when ins? indel-seq) (when del? indel-num) nil nil))))))
         (persistent! results)))))
 
 (defn- parse-pileup-line
diff --git a/src/cljam/io/sam/util/flag.clj b/src/cljam/io/sam/util/flag.clj
index 106a39e..971eb39 100644
--- a/src/cljam/io/sam/util/flag.clj
+++ b/src/cljam/io/sam/util/flag.clj
@@ -19,10 +19,14 @@
 (def ^:const flag-keywords
   (vec (map vector (map key (sort-by val flags)) (range))))
 
+(defn- long-bit-test
+  [^long x ^long n]
+  (. clojure.lang.Numbers testBit x n))
+
 (defn decode
   "Returns a set of keywords for a given flag integer."
   [^long f]
-  (into #{} (for [[k i] flag-keywords :when (bit-test f i)] k)))
+  (into #{} (for [[k i] flag-keywords :when (long-bit-test f i)] k)))
 
 (defn encode
   "Returns a flag integer encoding set of keywords."
@@ -43,17 +47,17 @@
 (defn multiple?
   "Tests if the template has multiple segments."
   [^long f]
-  (bit-test f 0))
+  (long-bit-test f 0))
 
 (defn properly-aligned?
   "Tests if the paired-end segments are properly aligned."
   [^long f]
-  (bit-test f 1))
+  (long-bit-test f 1))
 
 (defn unmapped?
   "Tests if the segment is unmapped."
   [^long f]
-  (bit-test f 2))
+  (long-bit-test f 2))
 
 (defn both-unmapped?
   "Tests if both the segment and its mate segment are unmapped."
@@ -63,17 +67,17 @@
 (defn reversed?
   "Tests if the segment is reversed."
   [^long f]
-  (bit-test f 4))
+  (long-bit-test f 4))
 
 (defn r1?
   "Tests if the segment is the first in the template."
   [^long f]
-  (bit-test f 6))
+  (long-bit-test f 6))
 
 (defn r2?
   "Tests if the segment is the last in the template."
   [^long f]
-  (bit-test f 7))
+  (long-bit-test f 7))
 
 (defn r1r2
   "Returns 0 for single-end, 1 for R1 and 2 for R2."
@@ -85,4 +89,4 @@
 (defn secondary?
   "Tests if the alignment is secondary."
   [^long f]
-  (bit-test f 8))
+  (long-bit-test f 8))
diff --git a/src/cljam/io/vcf/reader.clj b/src/cljam/io/vcf/reader.clj
index 19a931d..62c58cf 100644
--- a/src/cljam/io/vcf/reader.clj
+++ b/src/cljam/io/vcf/reader.clj
@@ -27,7 +27,7 @@
   protocols/IRegionReader
   (read-in-region [this region]
     (protocols/read-in-region this region {}))
-  (read-in-region [this {:keys [chr start end]} option]
+  (read-in-region [this {:keys [chr start end]} {:keys [depth] :as option}]
     (logging/warn "May cause degradation of performance.")
     (filter
      (fn [v] (and (if chr (= (:chr v) chr) true)
@@ -141,7 +141,7 @@
 ;; ------------------
 
 (defn- parse-data-line
-  [line kws]
+  [line header kws]
   ;; When splitting a string with single-character delimiter,
   ;; `java.lang.String#split` is slightly faster than `clojure.string/split`.
   ;; For more details, please refer to https://github.com/chrovis/cljam/pull/29.
@@ -165,7 +165,7 @@
   [^java.io.BufferedReader rdr header kws]
   (when-let [line (.readLine rdr)]
     (if-not (or (meta-line? line) (header-line? line))
-      (cons (parse-data-line line kws)
+      (cons (parse-data-line line header kws)
             (lazy-seq (read-data-lines rdr header kws)))
       (read-data-lines rdr header kws))))
 
diff --git a/src/cljam/io/vcf/util.clj b/src/cljam/io/vcf/util.clj
index 1af41bd..ba0af93 100644
--- a/src/cljam/io/vcf/util.clj
+++ b/src/cljam/io/vcf/util.clj
@@ -1,5 +1,6 @@
 (ns cljam.io.vcf.util
-  (:require [clojure.string :as cstr]))
+  (:require [clojure.string :as cstr]
+            [cljam.util :as util]))
 
 (definline dot-or-nil?
   "Checks if given string is equal to \".\" or nil."
@@ -89,9 +90,9 @@
   [^String gt]
   (when-not (dot-or-nil? gt)
     (->> gt
-         (re-seq #"([\||/])?(.|\d+)")
-         (map (fn [[_ phase ^String allele]]
-                [(when-not (dot-or-nil? allele) (Integer/parseInt allele))
+         (re-seq #"([\||/])?(\d+)")
+         (map (fn [[_ phase allele]]
+                [(Integer/parseInt allele)
                  (if phase
                    (= phase "|")
                    (neg? (.indexOf ^String gt "/")))])))))
@@ -100,64 +101,7 @@
   "Stringifies genotype map into VCF-style GT string."
   [gt-seq]
   (when-not (or (nil? gt-seq) (empty? gt-seq))
-    (->> gt-seq
-         (mapcat
-          (fn [[allele phase]]
-            [(if phase "|" "/") (or allele \.)]))
-         rest
-         (apply str))))
-
-(defn genotype-seq
-  "Returns a sequence of genotypes represented as sequences of integers."
-  ([^long ploidy ^long n-alt-alleles]
-   (genotype-seq ploidy n-alt-alleles []))
-  ([^long ploidy ^long n-alt-alleles s]
-   (mapcat
-    (fn [i]
-      (if (= 1 ploidy)
-        [(cons i s)]
-        (genotype-seq (dec ploidy) i (cons i s))))
-    (range (inc n-alt-alleles)))))
-
-(defn genotype-index
-  "Returns an index for given genotype."
-  ^long [genotype]
-  {:pre [(seq genotype)
-         (every? (complement neg?) genotype)]}
-  (letfn [(combination ^long [^long n ^long k]
-            (cond
-              (< n k) 0
-              (= k n) 1
-              (or (= k 1) (= k (dec n))) n
-              (< (quot n 2) k) (recur n (- n k))
-              :else (loop [i (inc (- n k)), j 1, r 1]
-                      (if (<= j k)
-                        (recur (inc i) (inc j) (quot (* r i) j))
-                        r))))]
-    (->> genotype
-         sort
-         (map-indexed (fn [^long m ^long k] (combination (+ k m) (inc m))))
-         (apply +))))
-
-(defn biallelic-genotype
-  "Converts a multiallelic `genotype` string into a biallelic one. Ignores all
-  alleles other than the reference allele and the `target-allele`."
-  [genotype ^long target-allele]
-  (->> genotype
-       parse-genotype
-       (map (fn [[allele phased?]]
-              [(when allele (if (= target-allele allele) 1 0)) phased?]))
-       stringify-genotype))
-
-(defn biallelic-coll
-  "Picks up elements in multiallelic `coll` and make a biallelic one. Ignores
-  all alleles other than the reference allele and the `target-allele`."
-  [^long ploidy ^long n-alt-alleles ^long target-allele coll]
-  (keep-indexed
-   (fn [i gt]
-     (when (every? (fn [^long a] (or (= a 0) (= a target-allele))) gt)
-       (nth coll i)))
-   (genotype-seq ploidy n-alt-alleles)))
+    (apply str (rest (mapcat (fn [[allele phase]] [(if phase "|" "/") allele]) gt-seq)))))
 
 (defn genotype->ints
   "Convert genotype to a sequence of integers."
@@ -196,8 +140,8 @@
   (->> formats
        (map (fn [k] [k (sample-map k)]))
        reverse
-       (drop-while (fn [[_ v]] (or (nil? v) (= [nil] v))))
-       (map (fn [[_ v]]
+       (drop-while (fn [[k v]] (or (nil? v) (= [nil] v))))
+       (map (fn [[k v]]
               (cond
                 (sequential? v) (cstr/join "," (map (fn [i] (if (nil? i) "." i)) v))
                 (nil? v) "."


Test patch:
diff --git a/test/cljam/io/pileup_test.clj b/test/cljam/io/pileup_test.clj
index 7ef0cf4..4d195a0 100644
--- a/test/cljam/io/pileup_test.clj
+++ b/test/cljam/io/pileup_test.clj
@@ -34,7 +34,8 @@
      :end? true,
      :insertion nil,
      :deletion nil,
-     :qname nil}]})
+     :qname nil,
+     :alignment nil}]})
 
 ;; Reader
 ;; ------
@@ -42,8 +43,8 @@
 (deftest parse-bases-col
   (testing "without-ref"
     (are [?in ?out]
-        (= (mapv #(into {:qname nil} %) ?out)
-           (mapv rec->map (#'plpio/parse-bases-col nil ?in)))
+         (= (mapv #(into {:qname nil :alignment 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}]
@@ -70,8 +71,8 @@
                    {: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)))
+         (= (mapv #(into {:qname nil :alignment 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}]
@@ -86,7 +87,7 @@
                    {: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))
+         (thrown? Exception (#'plpio/parse-bases-col nil ?in))
       "A+"
       "^"
       "^A"
@@ -105,30 +106,30 @@
 (deftest parse-pileup-line
   (testing "without-ref"
     (are [?in ?out]
-        (= ?out (rec->map (#'plpio/parse-pileup-line ?in)))
+         (= ?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}]}
+                              :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment 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}]}
+                                :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment nil}
+                                       {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment 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}]}))
+                                     :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment nil}
+                                            {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? true, :insertion nil, :deletion 2, :qname nil, :alignment nil}]}))
   (testing "with-ref"
     (are [?in ?out]
-        (= ?out (rec->map (#'plpio/parse-pileup-line ?in)))
+         (= ?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}]}
+                              :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment 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}]}
+                                :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment nil}
+                                       {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment 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}]})))
+                                     :pile [{:start? false, :mapq nil, :base \A, :qual 40, :reverse? true, :end? false, :insertion nil, :deletion nil, :qname nil, :alignment nil}
+                                            {:start? false, :mapq nil, :base \A, :qual 33, :reverse? false, :end? true, :insertion nil, :deletion 2, :qname nil, :alignment nil}]})))
 
 (deftest reader
   (with-open [r (plpio/reader test-pileup-file)]
@@ -152,8 +153,8 @@
 (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)))
+         (= ?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"
@@ -179,8 +180,8 @@
               (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)))
+           (= ?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}} ","
@@ -215,10 +216,10 @@
 (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)})))
+         (= ?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~"
@@ -228,10 +229,10 @@
   (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)})))
+           (= ?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~"
@@ -264,35 +265,36 @@
 ;; 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)))
+         (= ?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)))))))
+         (= ?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")))
diff --git a/test/cljam/io/vcf/util_test.clj b/test/cljam/io/vcf/util_test.clj
index b6b9946..35411bf 100644
--- a/test/cljam/io/vcf/util_test.clj
+++ b/test/cljam/io/vcf/util_test.clj
@@ -1,7 +1,8 @@
 (ns cljam.io.vcf.util-test
-  (:require [clojure.test :refer :all]
-            [cljam.test-common :refer :all]
-            [cljam.io.vcf.util :as vcf-util]))
+  (:require  [clojure.test :refer :all]
+             [clojure.string :as cstr]
+             [cljam.test-common :refer :all]
+             [cljam.io.vcf.util :as vcf-util]))
 
 (deftest about-parse-info
   (let [parse-info (vcf-util/info-parser [{:id "NS", :number 1, :type "Integer"}
@@ -73,11 +74,9 @@
     "0/0" [[0 false] [0 false]]
     "0/1" [[0 false] [1 false]]
     "1/1" [[1 false] [1 false]]
-    "./." [[nil false] [nil false]]
     "0|0" [[0 true] [0 true]]
     "0|1" [[0 true] [1 true]]
     "1|1" [[1 true] [1 true]]
-    ".|." [[nil true] [nil true]]
     "0/1/2" [[0 false] [1 false] [2 false]]
     "0/1|2" [[0 false] [1 false] [2 true]]))
 
@@ -92,99 +91,12 @@
     [[0 true] [0 true]] "0|0"
     [[0 true] [1 true]] "0|1"
     [[1 true] [1 true]] "1|1"
-    [[nil true] [nil true]] ".|."
     [[0 false] [0 false]] "0/0"
     [[0 false] [1 false]] "0/1"
     [[1 false] [1 false]] "1/1"
-    [[nil false] [nil false]] "./."
     [[0 false] [1 false] [2 false]] "0/1/2"
     [[0 false] [1 false] [2 true]] "0/1|2"))
 
-(deftest genotype-seq
-  (are [?ploidy ?n-alt-alleles ?expected]
-      (= ?expected (vcf-util/genotype-seq ?ploidy ?n-alt-alleles))
-    1 1 [[0] [1]]
-    1 2 [[0] [1] [2]]
-    2 1 [[0 0] [0 1] [1 1]]
-    2 2 [[0 0] [0 1] [1 1] [0 2] [1 2] [2 2]]
-    2 3 [[0 0] [0 1] [1 1] [0 2] [1 2] [2 2] [0 3] [1 3] [2 3] [3 3]]
-    3 1 [[0 0 0] [0 0 1] [0 1 1] [1 1 1]]
-    3 2 [[0 0 0] [0 0 1] [0 1 1] [1 1 1] [0 0 2]
-         [0 1 2] [1 1 2] [0 2 2] [1 2 2] [2 2 2]]))
-
-(deftest genotype-index
-  (are [?genotype ?expected]
-      (= ?expected (vcf-util/genotype-index ?genotype))
-    [0] 0
-    [1] 1
-    [0 0] 0
-    [0 1] 1
-    [1 1] 2
-    [0 2] 3
-    [1 2] 4
-    [2 2] 5
-    [3 3] 9
-    [0 0 0] 0
-    [1 1 2] 6
-    [1 2 2] 8))
-
-(deftest about-genotypes
-  (are [?ploidy ?n-alt-alleles]
-      (let [x (vcf-util/genotype-seq ?ploidy ?n-alt-alleles)]
-        (= (range (count x)) (map vcf-util/genotype-index x)))
-    1 0
-    1 1
-    1 2
-    1 3
-    2 0
-    2 1
-    2 2
-    2 3
-    2 4
-    3 0
-    3 1
-    3 2
-    3 3
-    3 4
-    4 0
-    4 1
-    4 2
-    4 3
-    4 4))
-
-(deftest biallelic-genotype
-  (are [?genotype ?target-allele ?expected]
-      (= ?expected (vcf-util/biallelic-genotype ?genotype ?target-allele))
-    "0" 1 "0"
-    "1" 1 "1"
-    "2" 1 "0"
-    "0/0" 1 "0/0"
-    "0/1" 1 "0/1"
-    "0/2" 1 "0/0"
-    "1/2" 1 "1/0"
-    "2|2" 1 "0|0"
-    "0/1" 2 "0/0"
-    "0|2" 2 "0|1"
-    "1/2" 2 "0/1"
-    "0/1/2" 2 "0/0/1"
-    "2|3|4" 3 "0|1|0"
-    ;; a call cannot be made
-    nil 0 nil
-    "./." 1 "./."
-    ;; hemizygous
-    "./1" 1 "./1"
-    "./2" 2 "./1"
-    "./1" 2 "./0"))
-
-(deftest biallelic-coll
-  (are [?ploidy ?n-alt-alleles ?target-allele ?coll ?expected]
-      (= ?expected
-         (vcf-util/biallelic-coll ?ploidy ?n-alt-alleles ?target-allele ?coll))
-    2 1 1 [10 20 30] [10 20 30]
-    2 2 1 [10 20 30 40 50 60] [10 20 30]
-    2 2 2 [10 20 30 40 50 60] [10 40 60]
-    3 2 2 [1 2 3 4 5 6 7 8 9 10] [1 5 8 10]))
-
 (deftest about-parse-variant-v4_3
   (let [parse-variant (vcf-util/variant-parser test-vcf-v4_3-meta-info test-vcf-v4_3-header)]
     (are [?variant ?expected]