Repository: chrovis/cljam
Pull request #308 (2024-05-01 05:09:17)

Automatic read name generation for CRAM records
This PR enables the CRAM reader to automatically generate read names (QNAMEs) when they are absent in CRAM records.

## Specification

According to section [§10.3](https://github.com/samtools/hts-specs/blob/401e2ec3e8ecf19fa89def74e9c27c1e40df341e/CRAMv3.tex#L1382) of the CRAM specification v3.1:

- If the `RN` field in the compression header's preservation map of a slice is set to `true`, the decoder will decode read names from the `RN` data series
- If the `RN` field is set to `false` and a record has a detached flag (`0x2`) set in the `CF` field, the decoder will also decode its read name from the `RN` data series
- For all other records, the decoder will automatically generate read names
  - Although the CRAM specification does not define a specific method for generating read names, it suggests that names "are typically based on the file name and a numeric ID of the read, utilizing the record counter field from the slice header block"
  - Also, the decoder must generate consistent read names for mate reads

Due to the last condition, read name generation must occur after resolving mate reads.

## Implementation

This PR adds a "qname generator" to the CRAM reader. A qname generator is responsible for generating read names based on the file name and the global record counter. If the `RN` field is set to `false` in the compression header, the CRAM record decoder generates read names by calling the qname generator after mate read resolution.

The format for generated read names is `<file name>:<counter>`, which is aligned with [htslib](https://github.com/samtools/htslib/blob/deeb9f01376ca9416315e4c9f5fe489e6f03e05f/cram/cram_decode.c#L3032-L3043). The `<file name>` will be tweaked from the actual file name to ensure it complies with [the specification for QNAMEs in SAM](https://github.com/samtools/hts-specs/blob/401e2ec3e8ecf19fa89def74e9c27c1e40df341e/SAMv1.tex#L425), meaning it does not contain characters outside the `[!-?A-~]` range and does not exceed 254 characters in length.

## Note

The PR also includes new test code based on the `1001_name.cram` test file from [hts-specs](https://github.com/samtools/hts-specs/tree/401e2ec3e8ecf19fa89def74e9c27c1e40df341e/test/cram/3.0/passed).

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/chrovis/cljam/pull/308?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 `84.84848%` with `5 lines` in your changes are missing coverage. Please review.
> Project coverage is 88.47%. Comparing base [(`1dd64cd`)](https://app.codecov.io/gh/chrovis/cljam/commit/1dd64cd1a959e32a60bdccf6f09dee471f975b0d?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) to head [(`d0143c8`)](https://app.codecov.io/gh/chrovis/cljam/pull/308?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/308?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/core.clj](https://app.codecov.io/gh/chrovis/cljam/pull/308?src=pr&el=tree&filepath=src%2Fcljam%2Fio%2Fcram%2Fcore.clj&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis#diff-c3JjL2NsamFtL2lvL2NyYW0vY29yZS5jbGo=) | 66.66% | [5 Missing :warning: ](https://app.codecov.io/gh/chrovis/cljam/pull/308?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             @@
##           master     #308      +/-   ##
==========================================
- Coverage   88.48%   88.47%   -0.02%     
==========================================
  Files          94       94              
  Lines        8148     8175      +27     
  Branches      505      505              
==========================================
+ Hits         7210     7233      +23     
- Misses        433      437       +4     
  Partials      505      505              
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/chrovis/cljam/pull/308?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/core.clj b/src/cljam/io/cram/core.clj
index beb739e..6d3cc7d 100644
--- a/src/cljam/io/cram/core.clj
+++ b/src/cljam/io/cram/core.clj
@@ -5,12 +5,26 @@
             [cljam.io.sam.util.refs :as util.refs]
             [cljam.io.util.byte-buffer :as bb]
             [cljam.util :as util]
-            [clojure.java.io :as cio])
+            [clojure.java.io :as cio]
+            [clojure.string :as str])
   (:import [cljam.io.cram.reader CRAMReader]
            [java.io FileNotFoundException]
+           [java.net URL]
            [java.nio.channels FileChannel]
            [java.nio.file OpenOption StandardOpenOption]))
 
+(defn- qname-generator [^URL url]
+  (let [prefix (-> (.getPath url)
+                   (str/replace #"^.*/" "")
+                   (str/replace #"[^!-?A-~]+" "_"))
+        len (inc (count prefix))]
+    (fn [i]
+      (let [digits (str (inc (long i)))
+            len' (+ len (count digits))]
+        (cond-> (str prefix \: digits)
+          (> len' 254)
+          (subs (- len' 254) len'))))))
+
 (defn reader
   "Creates a new CRAM reader that reads a CRAM file f.
 
@@ -23,6 +37,7 @@
                              (into-array OpenOption [StandardOpenOption/READ]))
         bb (bb/allocate-lsb-byte-buffer 256)
         seq-resolver (some-> reference resolver/seq-resolver)
+        qname-gen (qname-generator url)
         header (volatile! nil)
         refs (delay (vec (util.refs/make-refs @header)))
         offset (volatile! nil)
@@ -31,7 +46,7 @@
                 (crai/read-index (str f ".crai") @refs)
                 (catch FileNotFoundException _
                   nil)))
-        rdr (reader/->CRAMReader url ch bb header refs offset idx seq-resolver)]
+        rdr (reader/->CRAMReader url ch bb header refs offset idx seq-resolver qname-gen)]
     (reader/read-file-definition rdr)
     (vreset! header (reader/read-header rdr))
     (vreset! offset (.position ch))
@@ -51,7 +66,8 @@
                                   (delay @(.-refs rdr))
                                   (delay @(.-offset rdr))
                                   (delay @(.-index rdr))
-                                  seq-resolver)]
+                                  seq-resolver
+                                  (.-qname-generator rdr))]
     (reader/read-file-definition rdr')
     (reader/skip-container rdr')
     rdr'))
diff --git a/src/cljam/io/cram/decode/record.clj b/src/cljam/io/cram/decode/record.clj
index bca7849..b673e64 100644
--- a/src/cljam/io/cram/decode/record.clj
+++ b/src/cljam/io/cram/decode/record.clj
@@ -39,10 +39,10 @@
 (defn- build-read-name-decoder [{:keys [preservation-map]} {:keys [RN]}]
   (if (:RN preservation-map)
     #(assoc % :qname (String. ^bytes (RN)))
-    (throw (ex-info "Omitted read names are not supported yet." {}))))
+    identity))
 
-(defn- build-mate-read-decoder [cram-header {:keys [MF NS NP TS NF]}]
-  (fn [{:keys [rname] :as record}]
+(defn- build-mate-read-decoder [cram-header {:keys [MF NS RN NP TS NF]}]
+  (fn [{:keys [rname qname] :as record}]
     (let [flag (long (::flag record))]
       (if (pos? (bit-and flag 0x02))
         (let [mate-flag (long (MF))
@@ -54,6 +54,7 @@
                          (bit-or (sam.flag/encoded #{:next-unmapped})))
               rnext (ref-name cram-header (NS))]
           (assoc record
+                 :qname (or qname (String. ^bytes (RN)))
                  :flag bam-flag
                  :rnext (if (and (= rname rnext) (not= rname "*")) "=" rnext)
                  :pnext (NP)
@@ -337,9 +338,23 @@
           (aset records i record')
           (aset records j mate'))))))
 
+(defn- assign-record-names [qname-generator slice-header ^objects records]
+  (let [counter-start (long (:counter slice-header))]
+    (dotimes [i (alength records)]
+      (let [record (aget records i)]
+        (when (nil? (:qname record))
+          ;; detached records always have qname, so if the control comes in here,
+          ;; it means this record is not detached and must have ::next-fragment
+          (let [j (inc (+ i (long (::next-fragment record))))
+                mate (aget records j)
+                qname (qname-generator (+ counter-start i))]
+            (aset records i (assoc record :qname qname))
+            (aset records j (assoc mate :qname qname))))))))
+
 (defn decode-slice-records
   "Decodes CRAM records in a slice all at once and returns them as a sequence of maps."
-  [seq-resolver cram-header compression-header slice-header ds-decoders tag-decoders]
+  [seq-resolver qname-generator cram-header compression-header slice-header
+   ds-decoders tag-decoders]
   (let [record-decoder (build-cram-record-decoder seq-resolver
                                                   cram-header
                                                   compression-header
@@ -351,4 +366,6 @@
     (dotimes [i n]
       (aset records i (record-decoder)))
     (resolve-mate-records records)
+    (when-not (get-in compression-header [:preservation-map :RN])
+      (assign-record-names qname-generator slice-header records))
     (map #(dissoc % ::flag ::len ::next-fragment) records)))
diff --git a/src/cljam/io/cram/reader.clj b/src/cljam/io/cram/reader.clj
index bb4e375..e40c613 100644
--- a/src/cljam/io/cram/reader.clj
+++ b/src/cljam/io/cram/reader.clj
@@ -14,7 +14,7 @@
 
 (declare read-alignments read-alignments-in-region)
 
-(deftype CRAMReader [url channel buffer header refs offset index seq-resolver]
+(deftype CRAMReader [url channel buffer header refs offset index seq-resolver qname-generator]
   Closeable
   (close [_]
     (when seq-resolver
@@ -90,6 +90,7 @@
         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-for-slice rdr slice-header blocks)
+                                 (.-qname-generator rdr)
                                  @(.-header rdr)
                                  compression-header
                                  slice-header


Test patch:
diff --git a/test/cljam/io/cram/decode/record_test.clj b/test/cljam/io/cram/decode/record_test.clj
index eefbf03..877b5f1 100644
--- a/test/cljam/io/cram/decode/record_test.clj
+++ b/test/cljam/io/cram/decode/record_test.clj
@@ -426,6 +426,7 @@
                :cigar "5M", :rnext "*", :pnext 0, :tlen 0, :seq "CTGTG", :qual "AEEEE"
                :options []}]
              (record/decode-slice-records test-seq-resolver
+                                          nil
                                           cram-header
                                           compression-header
                                           slice-header
@@ -483,6 +484,73 @@
                :cigar "*", :rnext "*", :pnext 0, :tlen 0, :seq "GCACA", :qual "BCCFD"
                :options []}]
              (record/decode-slice-records test-seq-resolver
+                                          nil
+                                          cram-header
+                                          compression-header
+                                          slice-header
+                                          ds-decoders
+                                          tag-decoders)))))
+  (testing "unnamed reads"
+    (let [cram-header {:SQ
+                       [{:SN "ref"}
+                        {:SN "ref2"}]}
+          compression-header {:preservation-map
+                              {:RN false
+                               :AP true
+                               :RR true
+                               :SM {\A {0 \C, 1 \G, 2 \T, 3 \N}
+                                    \C {0 \A, 1 \G, 2 \T, 3 \N}
+                                    \G {0 \A, 1 \C, 2 \T, 3 \N}
+                                    \T {0 \A, 1 \C, 2 \G, 3 \N}
+                                    \N {0 \A, 1 \C, 2 \G, 3 \T}}
+                               :TD [[]]}}
+          slice-header {:ref-seq-id 0
+                        :start 10
+                        :records 7
+                        :counter 100}
+          ds-decoders (build-stub-decoders
+                       {:BF [67 67 147 147 147 67 0]
+                        :CF [5 5 1 1 3 3 3]
+                        :RL [5 5 5 5 5 5 5]
+                        :AP [0 0 20 0 0 0 0]
+                        :RG [-1 -1 -1 -1 -1 -1 -1]
+                        :RN (->> [nil nil nil nil "q003" "q004" "q005"]
+                                 (keep #(some-> ^String % .getBytes)))
+                        :MF [nil nil nil nil 0 1 2]
+                        :NF [1 1 nil nil nil nil nil]
+                        :NS [nil nil nil nil 0 1 -1]
+                        :NP [nil nil nil nil 50 100 0]
+                        :TS [nil nil nil nil 25 0 0]
+                        :TL [0 0 0 0 0 0 0]
+                        :FN [0 0 0 0 0 0 0]
+                        :QS (->> (repeat 7 "#####")
+                                 (mapcat #(.getBytes ^String %))
+                                 (map #(- (long %) 33)))
+                        :MQ [40 40 40 40 40 40 40]})
+          tag-decoders (build-stub-tag-decoders {})]
+      (is (= [{:qname "gen:101", :flag 99, :rname "ref", :pos 10, :end 14, :mapq 40,
+               :cigar "5M", :rnext "=", :pnext 30, :tlen 25, :seq "GATAA", :qual "#####"
+               :options []}
+              {:qname "gen:102", :flag 99, :rname "ref", :pos 10, :end 14, :mapq 40,
+               :cigar "5M", :rnext "=", :pnext 30, :tlen 25, :seq "GATAA", :qual "#####"
+               :options []}
+              {:qname "gen:101", :flag 147, :rname "ref", :pos 30, :end 34, :mapq 40,
+               :cigar "5M", :rnext "=", :pnext 10, :tlen -25, :seq "AGGCA", :qual "#####"
+               :options []}
+              {:qname "gen:102", :flag 147, :rname "ref", :pos 30, :end 34, :mapq 40,
+               :cigar "5M", :rnext "=", :pnext 10, :tlen -25, :seq "AGGCA", :qual "#####"
+               :options []}
+              {:qname "q003", :flag 147, :rname "ref", :pos 30, :end 34, :mapq 40,
+               :cigar "5M", :rnext "=", :pnext 50, :tlen 25, :seq "AGGCA", :qual "#####"
+               :options []}
+              {:qname "q004", :flag 99, :rname "ref", :pos 30, :end 34, :mapq 40,
+               :cigar "5M", :rnext "ref2", :pnext 100, :tlen 0, :seq "AGGCA", :qual "#####"
+               :options []}
+              {:qname "q005", :flag 8, :rname "ref", :pos 30, :end 34, :mapq 40,
+               :cigar "5M", :rnext "*", :pnext 0, :tlen 0, :seq "AGGCA", :qual "#####"
+               :options []}]
+             (record/decode-slice-records test-seq-resolver
+                                          #(str "gen:" (inc %))
                                           cram-header
                                           compression-header
                                           slice-header