Repository: chrovis/cljam
Pull request #117 (2017-10-25 11:21:23)

Add FASTA-TwoBit converter
Adds FASTA-TwoBit conversion, and test cases also.

As known issue, it is hard to convert a large reference file such as `hg38.fa` because of memory consumption. The reason is that `cljam.io.fasta.reader/sequential-read-string` is non-lazy and `cljam.io.twobit.reader/read-all-sequences` realizes some sequences by chunking. I plan on fixing them on another branch.


Patch:
diff --git a/src/cljam/algo/convert.clj b/src/cljam/algo/convert.clj
index 304c992..428a54a 100644
--- a/src/cljam/algo/convert.clj
+++ b/src/cljam/algo/convert.clj
@@ -1,15 +1,18 @@
 (ns cljam.algo.convert
-  "Format converter from SAM to BAM, and vice versa."
+  "Converters between equivalent formats: SAM/BAM and FASTA/TwoBit."
   (:require [clojure.tools.logging :as logging]
             [cljam.common :refer [*n-threads* get-exec-n-threads]]
             [cljam.io.sam :as sam]
             [cljam.io.bam.encoder :as encoder]
             [cljam.io.sam.util :as sam-util]
+            [cljam.io.sequence :as cseq]
             [cljam.io.util :as io-util]
             [com.climate.claypoole :as cp])
   (:import [java.nio ByteBuffer]
            [cljam.io.bam.writer BAMWriter]))
 
+;;; SAM <-> BAM
+
 (def ^:private default-num-block 100000)
 
 (defn- sam-write-alignments [rdr wtr hdr num-block]
@@ -30,22 +33,45 @@
                             (partition-all num-block (sam/read-alignments rdr {})))]
       (sam/write-blocks wtr blocks))))
 
-(defn- _convert!
+(defn- convert-sam*
   [rdr wtr num-block write-alignments-fn]
   (let [hdr (sam/read-header rdr)]
     (sam/write-header wtr hdr)
     (sam/write-refs wtr hdr)
     (write-alignments-fn rdr wtr hdr num-block)))
 
-(defn convert
-  "Converts file format from input file to output file by the file extension."
+(defn convert-sam
+  "Converts file format between SAM and BAM based on the file extension."
   [in out & {:keys [n-threads num-block]
              :or {n-threads 0, num-block default-num-block}}]
   (with-open [rdr (sam/reader in)
               wtr (sam/writer out)]
     (binding [*n-threads* n-threads]
       (cond
-        (io-util/sam-writer? wtr) (_convert! rdr wtr num-block sam-write-alignments)
-        (io-util/bam-writer? wtr) (_convert! rdr wtr num-block bam-write-alignments)
+        (io-util/sam-writer? wtr) (convert-sam* rdr wtr num-block sam-write-alignments)
+        (io-util/bam-writer? wtr) (convert-sam* rdr wtr num-block bam-write-alignments)
         :else (throw (ex-info (str "Unsupported output file format " out) {})))))
   nil)
+
+;;; FASTA <-> TwoBit
+
+(defn convert-sequence
+  "Converts file format between FASTA and TwoBit based on the file extension."
+  [in out]
+  (with-open [rdr (cseq/reader in)
+              wtr (cseq/writer out)]
+    (cseq/write-sequences wtr (cseq/read-all-sequences rdr {:mask? true}))))
+
+;;; General converter
+
+(def ^:private alignment-io? (every-pred (comp #{:sam :bam} io-util/file-type)))
+
+(def ^:private sequence-io? (every-pred (comp #{:fasta :2bit} io-util/file-type)))
+
+(defn convert
+  "Converts file format from input file to output file by the file extension."
+  [in out & opts]
+  (cond
+    (alignment-io? in out) (apply convert-sam in out opts)
+    (sequence-io? in out) (convert-sequence in out)
+    :else (throw (ex-info (str "Unsupported I/O pair: " in " and " out) {}))))
diff --git a/src/cljam/tools/cli.clj b/src/cljam/tools/cli.clj
index 5695264..f41ac05 100644
--- a/src/cljam/tools/cli.clj
+++ b/src/cljam/tools/cli.clj
@@ -82,19 +82,15 @@
 ;; ### convert command
 
 (def ^:private convert-cli-options
-  [["-if" "--input-format FORMAT" "Input file format <auto|sam|bam>"
-    :default "auto"]
-   ["-of" "--output-format FORMAT" "Output file format <auto|sam|bam>"
-    :default "auto"]
-   ["-t" "--thread THREAD" "Number of threads (0 is auto)"
+  [["-t" "--thread THREAD" "Number of threads (0 is auto)"
     :default 0
     :parse-fn #(Integer/parseInt %)]
    ["-h" "--help" "Print help"]])
 
 (defn- convert-usage [options-summary]
-  (->> ["Convert SAM to BAM or BAM to SAM."
+  (->> ["Convert file format based on the file extension."
         ""
-        "Usage: cljam convert [-if FORMAT] [-of FORMAT] [-t THREAD] <in.bam|sam> <out.bam|sam>"
+        "Usage: cljam convert [-t THREAD] <in-file> <out-file>"
         ""
         "Options:"
         options-summary]
@@ -381,7 +377,7 @@
                      "Usage: cljam {view,convert,normalize,sort,index,pileup,faidx,dict,level,version} ..."
                      :options  [["-h" "--help" "Show help" :default false :flag true]]
                      :commands [["view"    "Extract/print all or sub alignments in SAM or BAM format."]
-                                ["convert" "Convert SAM to BAM or BAM to SAM."]
+                                ["convert" "Convert file format based on the file extension."]
                                 ["normalize" "Normalize references of alignments."]
                                 ["sort"    "Sort alignments by leftmost coordinates."]
                                 ["index"   "Index sorted alignment for fast random access."]


Test patch:
diff --git a/test/cljam/algo/convert_test.clj b/test/cljam/algo/convert_test.clj
new file mode 100644
index 0000000..de0a328
--- /dev/null
+++ b/test/cljam/algo/convert_test.clj
@@ -0,0 +1,32 @@
+(ns cljam.algo.convert-test
+  (:require [clojure.test :refer :all]
+            [cljam.test-common :refer :all]
+            [cljam.algo.convert :as convert]))
+
+(deftest convert-test
+  (with-before-after {:before (prepare-cache!)
+                      :after (clean-cache!)}
+    (testing "SAM -> BAM"
+      (let [in test-sam-file
+            out (str temp-dir "/test.bam")]
+        (is (not-throw? (convert/convert in out)))
+        (is (same-sam-contents? in out))))
+    (testing "BAM -> SAM"
+      (let [in test-bam-file
+            out (str temp-dir "/test.sam")]
+        (is (not-throw? (convert/convert in out)))
+        (is (same-sam-contents? in out))))
+    (testing "FASTA -> TwoBit"
+      (let [in test-fa-file
+            out (str temp-dir "/test.2bit")]
+        (is (not-throw? (convert/convert in out)))
+        (is (same-sequence-contents? in out))))
+    (testing "TwoBit -> FASTA"
+      (let [in test-twobit-file
+            out (str temp-dir "/test.fa")]
+        (is (not-throw? (convert/convert in out)))
+        (is (same-sequence-contents? in out))))
+    (testing "error"
+      (are [in out] (thrown? Exception (convert/convert in out))
+        test-bam-file (str temp-dir "/test.unknown")
+        test-bam-file (str temp-dir "/test.2bit")))))
diff --git a/test/cljam/test_common.clj b/test/cljam/test_common.clj
index 3ab9b6d..3d98eef 100644
--- a/test/cljam/test_common.clj
+++ b/test/cljam/test_common.clj
@@ -5,6 +5,7 @@
             [clojure.tools.logging.impl :refer [disabled-logger-factory]]
             [cljam.io.protocols :as protocols]
             [cljam.io.sam :as sam]
+            [cljam.io.sequence :as cseq]
             [cavia.core :as cavia :refer [defprofile with-profile]]))
 
 (defn- _in-cloverage? []
@@ -477,7 +478,8 @@
   (= (digest/sha1 (file f1)) (digest/sha1 (file f2))))
 
 (defn same-sam-contents?
-  "Returns true if the contents of two SAM/BAM files are equivalent, false if not."
+  "Returns true if the contents of two SAM/BAM files are equivalent, false if
+  not."
   [f1 f2]
   (with-open [r1 (sam/reader f1)
               r2 (sam/reader f2)]
@@ -486,6 +488,15 @@
          (= (seq (sam/read-alignments r1))
             (seq (sam/read-alignments r2))))))
 
+(defn same-sequence-contents?
+  "Returns true if the contents of two FASTA/TwoBit files are equivalent, false
+  if not."
+  [f1 f2]
+  (with-open [r1 (cseq/reader f1)
+              r2 (cseq/reader f2)]
+    (= (cseq/read-all-sequences r1)
+       (cseq/read-all-sequences r2))))
+
 ;;;; FASTA
 
 (def test-fa-header [{:desc "", :name "ref", :offset 5} {:desc "", :name "ref2", :offset 57}])
diff --git a/test/cljam/tools/cli_test.clj b/test/cljam/tools/cli_test.clj
index b7c9911..19f5074 100644
--- a/test/cljam/tools/cli_test.clj
+++ b/test/cljam/tools/cli_test.clj
@@ -24,6 +24,8 @@
 (def temp-out (str temp-dir "/out"))
 (def temp-bam (str temp-dir "/out.bam"))
 (def temp-sam (str temp-dir "/out.sam"))
+(def temp-fasta (str temp-dir "/out.fa"))
+(def temp-twobit (str temp-dir "/out.2bit"))
 
 (deftest about-view
   (with-before-after {:before (prepare-cache!)
@@ -41,16 +43,24 @@
 (deftest about-convert
   (with-before-after {:before (prepare-cache!)
                       :after (clean-cache!)}
-    ;; sam => bam
-    (is (not-throw? (cli/convert [test-sam-file temp-bam])))
-    (is (= (slurp-bam-for-test temp-bam) (slurp-sam-for-test test-sam-file)))
-    (is (= (slurp-bam-for-test temp-bam) (slurp-bam-for-test test-bam-file)))
-    ;; bam => sam
-    (is (not-throw? (cli/convert [test-bam-file temp-sam])))
-    (is (= (slurp-sam-for-test temp-sam) (slurp-bam-for-test test-bam-file)))
-    (is (= (slurp-sam-for-test temp-sam) (slurp-sam-for-test test-sam-file)))
-    ;; error
-    (is (thrown? IllegalArgumentException (cli/convert [test-bam-file (str temp-dir "/test.unknown")])))))
+    (testing "SAM -> BAM"
+      (is (not-throw? (cli/convert [test-sam-file temp-bam])))
+      (is (= (slurp-bam-for-test temp-bam) (slurp-sam-for-test test-sam-file)))
+      (is (= (slurp-bam-for-test temp-bam) (slurp-bam-for-test test-bam-file))))
+    (testing "BAM -> SAM"
+      (is (not-throw? (cli/convert [test-bam-file temp-sam])))
+      (is (= (slurp-sam-for-test temp-sam) (slurp-bam-for-test test-bam-file)))
+      (is (= (slurp-sam-for-test temp-sam) (slurp-sam-for-test test-sam-file))))
+    (testing "FASTA -> TwoBit"
+      (is (not-throw? (cli/convert [test-fa-file temp-twobit])))
+      (is (same-sequence-contents? test-fa-file temp-twobit)))
+    (testing "TwoBit -> FASTA"
+      (is (not-throw? (cli/convert [test-twobit-file temp-fasta])))
+      (is (same-sequence-contents? test-twobit-file temp-fasta)))
+    (testing "error"
+      (are [in out] (thrown? Exception (cli/convert [in out]))
+        test-bam-file (str temp-dir "/test.unknown")
+        test-bam-file temp-twobit))))
 
 (deftest about-normalize
   (with-before-after {:before (prepare-cache!)