Carry aligned reads as CRAM, with --aligned_format cram|bam

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Changed`
 
+- [#207](https://github.com/IntGenomicsLab/lrsomatic/pull/207) - Aligned and haplotagged reads are now CRAM 3.0 with the reference embedded instead of BAM, from `MINIMAP2_ALIGN` onwards; `bamfiles/` holds `<sample>_<type>.cram` + `.crai`. The new `--aligned_format` parameter (`cram` by default, or `bam`) keeps BAM throughout for tools that cannot read CRAM. On chr20 of three real samples CRAM is 33 % (ONT), 52 % (older PacBio) and 69 % (Revio) smaller, with identical records. Version 3.0 because Clair3 and ClairS-TO's Verdict bundle an htslib that cannot decode 3.1 and then call nothing without failing; the embedded reference (1-7 % of the file) because ClairS, ClairS-TO, Severus and SAVANA open the reads without the FASTA. `CRAMINO_POST` and `SAMTOOLS_MERGE` now get the FASTA (`ASCAT` does not: its image has no `ref.fasta` argument, and it reads the embedded reference); `CLAIRSTO` deletes its intermediate haplotagged copy of the tumour reads; the uBAM's own `RG` tag is dropped before alignment, since CRAM keeps one read group per record and kept the basecaller's. That change alters `MINIMAP2_ALIGN` in both formats, so a run resumed from an earlier version realigns every sample, also with `--aligned_format bam` (@ljwharbers).
 - [#201](https://github.com/IntGenomicsLab/lrsomatic/pull/201) - `CLAIRSTO` and `CLAIRSTO_VERDICT_TAG` now pull the fork image from Docker Hub: `oras://docker.io/ljwharbers/clairs-to-sif:0.5.1-verdict-chm13-c0687e8-flat` under Singularity/Apptainer and `docker.io/ljwharbers/clairs-to:0.5.1-verdict-chm13-c0687e8-flat` otherwise, instead of `ghcr.io/ljwharbers/clairs-to`. The `-cpu` SIF on ghcr failed with `PROTOCOL_ERROR` on slow links: ghcr redirects every blob download to an Azure URL that expires at the next 5-minute mark and resets a stream still open then, and Apptainer resumes neither an `oras://` nor a `docker://` download. Docker Hub's download URLs are valid for 50 minutes and only checked when the request starts. `-flat` is the same software copied into an empty image in a few layers (3.3 GB instead of 7 GB); the software and its outputs are unchanged. `docs/usage.md` describes `pullTimeout`, Docker Hub's anonymous pull limit, pre-pulling, and how to recover the remaining `oras://ghcr.io` SIFs resumably (@ljwharbers).
 - [#199](https://github.com/IntGenomicsLab/lrsomatic/pull/199) - `CLAIRSTO` and `CLAIRSTO_VERDICT_TAG` now run the `-cpu` rebuild of the fork image (`0.5.1-verdict-chm13-c0687e8-cpu`), which swaps PyTorch's CUDA build for the CPU build of the same version. The software is otherwise unchanged, but the Apptainer SIF drops from 6.53 GB to 3.46 GB. The old image could not be pulled on a normal VSC link: Apptainer fetches an `oras://` SIF as a single unresumable stream, and the signed blob URL ghcr redirects to expires on a 15-minute wall-clock boundary, so 6.53 GB needed 7.3 MB/s sustained and was otherwise cut mid-transfer with `PROTOCOL_ERROR` (@ljwharbers).
 - [#197](https://github.com/IntGenomicsLab/lrsomatic/pull/197) - `CLAIRSTO` now runs `ghcr.io/ljwharbers/clairs-to:0.5.1-verdict-chm13-c0687e8` (ClairS-TO 0.5.1) instead of `docker.io/hkubal/clairs-to:v0.4.2`: a fork that lets Verdict read its CNA resources from `--cna_resource_dir`, fixes four places where Verdict's Python port of ASCAT departed from R, and disables Verdict with a warning when its resources cannot be read. **GRCh38 results move as well as CHM13 ones.** Revert to the upstream image once HKU-BAL/ClairS-TO carries these changes. The module also selects the SIF under `-profile apptainer` and sets explicit output prefixes (@ljwharbers).

conf/modules.config

@@ -306,6 +306,13 @@ process {
                 "-R '${rg}'"
             ].join(' ').trim()
         }
+        // --aligned_format cram: CRAM 3.0, not 3.1, because the htslib bundled in Clair3 and in ClairS-TO's
+        // Verdict (1.11) cannot decode 3.1 slices and then calls nothing without failing; embed_ref=1 because
+        // ClairS, ClairS-TO, Severus and SAVANA's read counting open the reads without the reference (1-7 % larger)
+        ext.args2  = { params.aligned_format == 'cram' ? '--output-fmt-option version=3.0 --output-fmt-option embed_ref=1' : '' }
+        // Drop the uBAM's own RG tag: -T '*' would carry it through -y next to the @RG from -R, and a record
+        // with two RG tags keeps only the last one in CRAM -- the basecaller's, which has no @RG line here
+        ext.args3  = "-x RG"
         ext.args4  = "-T '*'"
         publishDir = [
             enabled: false
@@ -314,6 +321,7 @@ process {
 
     withName: '.*:SAMTOOLS_MERGE' {
         ext.prefix = { "${meta.id}_${meta.type}_mapped" }
+        ext.args   = { params.aligned_format == 'cram' ? '--output-fmt-option version=3.0 --output-fmt-option embed_ref=1' : '' }
         publishDir = [
             enabled: false
         ]
@@ -416,6 +424,9 @@ process {
                     params.longphase_tag_supplementary ? "--tagSupplementary" : ''
                 ].join(' ').trim()
             }
+        // --aligned_format cram: CRAM through a FIFO, so the reference can be embedded (longphase's own --cram cannot)
+        ext.cram   = { params.aligned_format == 'cram' }
+        ext.args2  = '--output-fmt-option version=3.0 --output-fmt-option embed_ref=1'
         publishDir = [
             path: { "${params.outdir}/${meta.id}/bamfiles" },
             mode: params.publish_dir_mode,

docs/output.md

@@ -152,18 +152,22 @@ The pipeline produces per-sample output directories. Two modes exist depending o
 
 ```
 ├── bamfiles
-│   ├── sample_normal.bam
-│   ├── sample_normal.bam.bai
-│   ├── sample_tumor.bam
-│   ├── sample_tumor.bam.bai
-```
-
-| File                    | Description                                                                                          |
-| ----------------------- | ---------------------------------------------------------------------------------------------------- |
-| `sample_normal.bam`     | Aligned and haplotagged bam file (with methylation and nucleosome predictions) for the normal sample |
-| `sample_normal.bam.bai` | index file for the normal bam file                                                                   |
-| ` sample_tumor.bam`     | Aligned and haplotagged bam file (with methylation and nucleosome predictions) for the tumor sample  |
-| `sample_tumor.bam.bai`  | index file for the tumor bam file                                                                    |
+│   ├── sample_normal.cram
+│   ├── sample_normal.cram.crai
+│   ├── sample_tumor.cram
+│   ├── sample_tumor.cram.crai
+```
+
+| File                      | Description                                                                                           |
+| ------------------------- | ----------------------------------------------------------------------------------------------------- |
+| `sample_normal.cram`      | Aligned and haplotagged CRAM file (with methylation and nucleosome predictions) for the normal sample |
+| `sample_normal.cram.crai` | index file for the normal CRAM file                                                                   |
+| `sample_tumor.cram`       | Aligned and haplotagged CRAM file (with methylation and nucleosome predictions) for the tumor sample  |
+| `sample_tumor.cram.crai`  | index file for the tumor CRAM file                                                                    |
+
+The CRAM files are version 3.0 with the reference embedded, so they can be read without the reference FASTA
+(`samtools view sample_tumor.cram`); `samtools view -b` turns one back into a BAM. With `--aligned_format bam`
+this directory holds `sample_{normal,tumor}.bam` and `.bam.bai` instead.
 
 </details>
 

docs/usage.md

@@ -242,6 +242,13 @@ opt-in. See [VEP plugins](#vep-plugins) for sizes, licence terms and per-assembl
 | `--minimap2_ont_model`       | specifies which model to use minimap2 with for ONT samples. Default = `null`                |
 | `--minimap2_pb_model`        | specifies which model to use minimap2 with for PacBio samples. Default = `null`             |
 | `--save_secondary_alignment` | A boolean to specify if secondary alignments are kept in aligned bam file. Default = `true` |
+| `--aligned_format`           | `cram` or `bam`: format of the aligned and haplotagged reads. Default = `cram` (see below)  |
+
+`--aligned_format cram` writes CRAM 3.0 with the reference embedded, from alignment onwards and in `bamfiles/`. On
+long reads this is a third (ONT) to two thirds (PacBio Revio) smaller than BAM, with identical records and the same
+calls from every tool in the pipeline. Because the reference is embedded, the files can be read without the FASTA.
+Use `--aligned_format bam` if a tool you run on the results cannot read CRAM, or convert afterwards with
+`samtools view -b -o sample_tumor.bam sample_tumor.cram`.
 
 #### ASCAT Options
 

modules.json

@@ -87,7 +87,8 @@
                     "longphase/haplotag": {
                         "branch": "master",
                         "git_sha": "b8d30a43f33aee3148b0e9e9f00587984a4ac195",
-                        "installed_by": ["modules"]
+                        "installed_by": ["modules"],
+                        "patch": "modules/nf-core/longphase/haplotag/longphase-haplotag.diff"
                     },
                     "longphase/phase": {
                         "branch": "master",

modules/local/clairsto/main.nf

@@ -71,6 +71,11 @@ process CLAIRSTO {
             cp -- "\$src" "${prefix}_Tumor_\${table}.txt"
         fi
     done
+
+    # The intermediate haplotagged copy of the tumor reads (phased_bam_output/, always BAM, so larger than a
+    # CRAM input) is unused once ClairS-TO exits; --remove_intermediate_dir would also take cna_output above.
+    # 0.5.1 names the directory tmp_<sample_name>.
+    rm -rf tmp*/phasing_output/phased_bam_output
     """
 
     stub:

modules/local/cramino/main.nf

@@ -9,6 +9,7 @@ process CRAMINO {
 
     input:
     tuple val(meta), path(bam)
+    tuple val(meta2), path(fasta)   // [[:], []] for BAM/uBAM input
 
     output:
     tuple val(meta), path("*.txt"), emit: txt
@@ -21,9 +22,10 @@ process CRAMINO {
     script:
     def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
+    def reference = fasta ? "--reference ${fasta}" : ''
 
     """
-    cramino $args $bam --arrow ${prefix}.arrow > ${prefix}_cramino.txt
+    cramino $args $reference $bam --arrow ${prefix}.arrow > ${prefix}_cramino.txt
     """
 
     stub:

modules/local/cramino/meta.yml

@@ -33,6 +33,15 @@ input:
           - edam: "http://edamontology.org/format_25722"
           - edam: "http://edamontology.org/format_2573"
           - edam: "http://edamontology.org/format_3462"
+  - - meta2:
+        type: map
+        description: |
+          Groovy Map containing reference information
+          e.g. `[ id:'genome' ]`
+    - fasta:
+        type: file
+        description: Reference FASTA to decode CRAM input; empty for BAM/uBAM
+        pattern: "*.{fa,fasta,fna}"
 
 output:
   txt: