Repository: chrovis/cljam
Pull request #307 (2024-04-23 08:49:23)

Support reading CRAM files with embedded references
~~**Note:** This PR will be marked as ready for review after #306 has been merged.~~

This PR allows the CRAM reader to handle files with embedded references. With this update, we no longer need to specify any reference file to read CRAM files with embedded references.

## Specification

According to the CRAM specification, if the encoder embeds a portion of the reference sequence into a slice, it allocates a dedicated block to store the raw (ASCII character) base sequence data for the reference sequence covering that slice. The content ID of this allocated block is pointed to by the "embedded reference bases block content ID" field in the slice header. If there is no embedded reference for the slice, this field's value will be `-1`. (See [§8.5](https://github.com/samtools/hts-specs/blob/401e2ec3e8ecf19fa89def74e9c27c1e40df341e/CRAMv3.tex#L944-L945) for more details)

## Implementation

In this change, the CRAM reader first checks the "embedded reference bases block content ID" field in the slice header to identify the block allocated for the embedded reference. If this value is positive, it creates a dedicated seq resolver based on the content of the block associated with the specified content ID. This new seq resolver is responsible for resolving reference sequences covered by the embedded reference and helps with the decoding of CRAM records in the slice.

## Notes

This PR also makes some slight changes to the signature of the seq resolver protocol method to make the new seq resolver for embedded reference a little bit more efficient.

Note that this PR does not include additional test code, although you can easily verify the change using a simple snippet with a test file `0600_mapped.cram` from [hts-specs](https://github.com/samtools/hts-specs/tree/master/test/cram/3.0/passed) (which is [documented](https://github.com/samtools/hts-specs/tree/master/test/cram/3.0#embedded-reference) as being encoded with embedded references) as follows:

```clojure
(require '[cljam.io.cram :as cram])

(with-open [cram-rdr (cram/reader "0600_mapped.cram")]
  (= ["ATTTTTCGGGTTTTTTGAAAATGTAGCTACAGAAAGTTTGTTTAAATATCTTGTTTTCTTGCACTTTGTGCAGAATT"
      "CCCTTTCAGAAAAATTATTTTTAAGAATTTTTCATRYTAGGAATATTGTTATTTCAGAAAATAGCTAAATGTGATTTCTGTAATTTTGCCTGCTAAAGGG"]
     (map :seq (cram/read-alignments cram-rdr))))
```

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/chrovis/cljam/pull/307?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) Report
Attention: Patch coverage is `45.00000%` with `11 lines` in your changes are missing coverage. Please review.
> Project coverage is 88.37%. Comparing base [(`4795708`)](https://app.codecov.io/gh/chrovis/cljam/commit/4795708b09754510df7b1f5ba5dda434e7961fd3?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) to head [(`db98dac`)](https://app.codecov.io/gh/chrovis/cljam/pull/307?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).

| [Files](https://app.codecov.io/gh/chrovis/cljam/pull/307?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) | Patch % | Lines |
|---|---|---|
| [src/cljam/io/cram/reader.clj](https://app.codecov.io/gh/chrovis/cljam/pull/307?src=pr&el=tree&filepath=src%2Fcljam%2Fio%2Fcram%2Freader.clj&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis#diff-c3JjL2NsamFtL2lvL2NyYW0vcmVhZGVyLmNsag==) | 31.25% | [11 Missing :warning: ](https://app.codecov.io/gh/chrovis/cljam/pull/307?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) |

<details><summary>Additional details and impacted files</summary>


```diff
@@                      Coverage Diff                      @@
##           feature/cram-index-access     #307      +/-   ##
=============================================================
- Coverage                      88.48%   88.37%   -0.12%     
=============================================================
  Files                             94       94              
  Lines                           8148     8163      +15     
  Branches                         505      505              
=============================================================
+ Hits                            7210     7214       +4     
- Misses                           433      444      +11     
  Partials                         505      505              
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/chrovis/cljam/pull/307?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).

athos: Rebased onto the latest master branch.
athos: Rebased onto the latest master branch.

Patch:
diff --git a/src/cljam/io/cram/decode/record.clj b/src/cljam/io/cram/decode/record.clj
index 38e2b6c..bca7849 100644
--- a/src/cljam/io/cram/decode/record.clj
+++ b/src/cljam/io/cram/decode/record.clj
@@ -98,8 +98,7 @@
 
 (defn- record-seq
   [seq-resolver {:keys [preservation-map]} {:keys [rname pos end] :as record} features]
-  (let [region {:chr rname :start pos :end end}
-        ref-bases (.getBytes ^String (resolver/resolve-sequence seq-resolver region))
+  (let [^bytes ref-bases (resolver/resolve-sequence seq-resolver rname pos end)
         len (long (::len record))
         bs (byte-array len (byte (int \N)))
         subst (:SM preservation-map)]
diff --git a/src/cljam/io/cram/reader.clj b/src/cljam/io/cram/reader.clj
index 6a5734d..bb4e375 100644
--- a/src/cljam/io/cram/reader.clj
+++ b/src/cljam/io/cram/reader.clj
@@ -3,11 +3,14 @@
             [cljam.io.cram.decode.data-series :as ds]
             [cljam.io.cram.decode.record :as record]
             [cljam.io.cram.decode.structure :as struct]
+            [cljam.io.cram.seq-resolver.protocol :as resolver]
             [cljam.io.protocols :as protocols]
+            [cljam.io.util.byte-buffer :as bb]
             [cljam.util.intervals :as intervals])
   (:import [java.io Closeable]
            [java.nio Buffer ByteBuffer ByteOrder]
-           [java.nio.channels FileChannel FileChannel$MapMode]))
+           [java.nio.channels FileChannel FileChannel$MapMode]
+           [java.util Arrays]))
 
 (declare read-alignments read-alignments-in-region)
 
@@ -59,6 +62,25 @@
   (read-to-buffer rdr 26)
   (struct/decode-file-definition (.-buffer rdr)))
 
+(defn- seq-resolver-for-slice [^CRAMReader rdr slice-header blocks]
+  (let [embedded-ref-content-id (long (:embedded-reference slice-header))]
+    (if (>= embedded-ref-content-id 0)
+      (let [ref-bases-block (->> blocks
+                                 (filter #(= (:content-id %) embedded-ref-content-id))
+                                 first)
+            offset (long (:start slice-header))
+            ^bytes bs (bb/read-bytes (:data ref-bases-block)
+                                     (:raw-size ref-bases-block))]
+        (reify resolver/ISeqResolver
+          ;; According to the CRAM specification v3.1 §8.5, a slice with an embedded
+          ;; reference must not be a multiple reference slice, so the temporary
+          ;; sequence resolver here can safely ignore the passed chr
+          (resolve-sequence [_ _chr start end]
+            (Arrays/copyOfRange bs
+                                (- (long start) offset)
+                                (inc (- (long end) offset))))))
+      (.-seq-resolver rdr))))
+
 (defn- read-slice-records [^CRAMReader rdr bb compression-header slice-header]
   (let [blocks (into [] (map (fn [_] (struct/decode-block bb)))
                      (range (:blocks slice-header)))
@@ -67,7 +89,7 @@
                      (bs/make-bit-stream-decoder (:data core-block)))
         ds-decoders (ds/build-data-series-decoders compression-header bs-decoder blocks)
         tag-decoders (ds/build-tag-decoders compression-header bs-decoder blocks)]
-    (record/decode-slice-records (.-seq-resolver rdr)
+    (record/decode-slice-records (seq-resolver-for-slice rdr slice-header blocks)
                                  @(.-header rdr)
                                  compression-header
                                  slice-header
diff --git a/src/cljam/io/cram/seq_resolver.clj b/src/cljam/io/cram/seq_resolver.clj
index 3dee9f9..0fa95ea 100644
--- a/src/cljam/io/cram/seq_resolver.clj
+++ b/src/cljam/io/cram/seq_resolver.clj
@@ -8,11 +8,12 @@
   (close [_]
     (.close ^Closeable seq-reader))
   proto/ISeqResolver
-  (resolve-sequence [_ region]
-    (cseq/read-sequence seq-reader region)))
+  (resolve-sequence [_ chr start end]
+    (when-let [s (cseq/read-sequence seq-reader {:chr chr :start start :end end})]
+      (.getBytes ^String s))))
 
 (defn seq-resolver
-  "Creates e new sequence resolver from the given sequence file."
+  "Creates a new sequence resolver from the given sequence file."
   [seq-file]
   (->SeqResolver (cseq/reader seq-file)))
 
diff --git a/src/cljam/io/cram/seq_resolver/protocol.clj b/src/cljam/io/cram/seq_resolver/protocol.clj
index ddebb14..f9bf482 100644
--- a/src/cljam/io/cram/seq_resolver/protocol.clj
+++ b/src/cljam/io/cram/seq_resolver/protocol.clj
@@ -1,11 +1,11 @@
 (ns cljam.io.cram.seq-resolver.protocol)
 
 (defprotocol ISeqResolver
-  (resolve-sequence [this region]))
+  (resolve-sequence [this chr start end]))
 
 (extend-protocol ISeqResolver
   nil
-  (resolve-sequence [_ region]
+  (resolve-sequence [_ chr start end]
     (throw
      (ex-info "reference was not specified, but tried to resolve sequence"
-              region))))
+              {:chr chr :start start :end end}))))


Test patch:
diff --git a/test/cljam/io/cram/decode/record_test.clj b/test/cljam/io/cram/decode/record_test.clj
index fd214b7..eefbf03 100644
--- a/test/cljam/io/cram/decode/record_test.clj
+++ b/test/cljam/io/cram/decode/record_test.clj
@@ -48,9 +48,9 @@
   (let [seqs (with-open [r (cseq/reader common/test-fa-file)]
                (into {} (map (juxt :name :sequence)) (cseq/read-all-sequences r)))]
     (reify resolver/ISeqResolver
-      (resolve-sequence [_ {:keys [chr start end]}]
+      (resolve-sequence [_ chr start end]
         (let [s (get seqs chr)]
-          (subs s (dec (long start)) end))))))
+          (.getBytes (subs s (dec (long start)) end)))))))
 
 (deftest record-end-test
   (are [?record ?features ?expected]
diff --git a/test/cljam/io/cram/seq_resolver_test.clj b/test/cljam/io/cram/seq_resolver_test.clj
index 25cffb1..984d108 100644
--- a/test/cljam/io/cram/seq_resolver_test.clj
+++ b/test/cljam/io/cram/seq_resolver_test.clj
@@ -2,20 +2,17 @@
   (:require [cljam.io.cram.seq-resolver :as resolver]
             [cljam.io.cram.seq-resolver.protocol :as proto]
             [cljam.test-common :as common]
-            [clojure.test :refer [are deftest is testing]]))
+            [clojure.test :refer [deftest is testing]]))
 
 (deftest resolve-sequence-test
   (testing "seq resolver resolves sequence for specified region"
     (let [resolver (resolver/seq-resolver common/test-fa-file)
           resolver' (resolver/clone-seq-resolver resolver)]
-      (are [?region ?expected]
-           (= ?expected
-              (proto/resolve-sequence resolver ?region)
-              (proto/resolve-sequence resolver' ?region))
-        {:chr "ref" :start 1 :end 5}
-        "AGCAT"
-
-        {:chr "unknown" :start 1 :end 5}
-        nil)))
+      (is (= "AGCAT"
+             (String. (proto/resolve-sequence resolver "ref" 1 5))
+             (String. (proto/resolve-sequence resolver' "ref" 1 5))))
+      (is (= nil
+             (proto/resolve-sequence resolver "unknown" 1 5)
+             (proto/resolve-sequence resolver' "unknown" 1 5)))))
   (testing "resolver can be omitted, but in that case, trying to resolve seq will end up with error"
-    (is (thrown? Exception (proto/resolve-sequence nil {:chr "ref" :start 1 :end 5})))))
+    (is (thrown? Exception (proto/resolve-sequence nil "ref" 1 5)))))