diff --git a/.github/ISSUE_TEMPLATE/bug_report.yml b/.github/ISSUE_TEMPLATE/bug_report.yml
index 9c53149..bb8336a 100644
--- a/.github/ISSUE_TEMPLATE/bug_report.yml
+++ b/.github/ISSUE_TEMPLATE/bug_report.yml
@@ -45,13 +45,19 @@ body:
       label: Workflow Execution
       description: Where are you running the workflow?
       options:
-        - EPI2ME Desktop application
-        - Command line
-        - EPI2ME cloud agent
+        - EPI2ME Desktop (Local)
+        - EPI2ME Desktop (Cloud)
+        - Command line (Local)
+        - Command line (Cluster)
         - Other (please describe)
     validations:
       required: true
-
+  - type: input
+    id: other-workflow-execution
+    attributes:
+      label: Other workflow execution
+      description: If "Other", please describe
+      placeholder: Tell us where / how you are running the workflow.
 
   - type: markdown
     attributes:
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8288b59..990f235 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,14 +4,14 @@ include:
       file: "wf-containers.yaml"
 
 variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  differential_expression/differential_expression_fastq \
-              --de_analysis --ref_genome differential_expression/hg38_chr20.fa \
-              --transcriptome-source reference-guided \
-              --ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \
-              --direct_rna --minimap2_index_opts '-k 15' --sample_sheet differential_expression/sample_sheet.csv \
-              --jaffal_refBase differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
-              -c demo.nextflow.config"
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/  &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
+              --de_analysis --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
+              --transcriptome-source ${CI_PROJECT_NAME}/data/reference-guided \
+              --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
+              --direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
+              --jaffal_refBase ${CI_PROJECT_NAME}/data/differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
+              -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               CI_FLAVOUR: "new"
 
 macos-run:
@@ -63,120 +63,120 @@ docker-run:
           when: never
         - if: $MATRIX_NAME == "isoforms"
           variables:
-              NF_BEFORE_SCRIPT: wget -O test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  test_data.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  ERR6053095_chr20.fastq --transcriptome-source reference-guided \
-                  --ref_genome chr20/hg38_chr20.fa --ref_annotation chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm \
-                  -c demo.nextflow.config"
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome-source reference-guided \
+                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm \
+                  -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "no_ref_annotation"
           variables:
-              NF_BEFORE_SCRIPT: wget -O test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  test_data.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  ERR6053095_chr20.fastq --transcriptome-source reference-guided \
-                  --ref_genome chr20/hg38_chr20.fa \
-                  -c demo.nextflow.config"
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome-source reference-guided \
+                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa \
+                  -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: run_gffcompare,preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "fusions"
           variables:
-              NF_BEFORE_SCRIPT: wget -O test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  test_data.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB --fastq  ERR6053095_chr20.fastq --transcriptome-source reference-guided \
-                  --ref_genome chr20/hg38_chr20.fa --ref_annotation chr20/gencode.v22.annotation.chr20.gtf \
-                  --jaffal_refBase chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
-                  -c demo.nextflow.config"
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB --fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome-source reference-guided \
+                  --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf \
+                  --jaffal_refBase ${CI_PROJECT_NAME}/data/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
+                  -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "differential_expression"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  differential_expression/differential_expression_fastq \
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
                   --de_analysis \
-                  --ref_genome differential_expression/hg38_chr20.fa --transcriptome-source reference-guided \
-                --ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \
-                --direct_rna --minimap2_index_opts '-k 15'  --sample_sheet differential_expression/sample_sheet.csv \
-                -c demo.nextflow.config"
+                  --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa --transcriptome-source reference-guided \
+                --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
+                --direct_rna --minimap2_index_opts '-k 15'  --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
+                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "only_differential_expression"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  differential_expression/differential_expression_fastq \
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/  &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
                   --transcriptome-source precomputed \
                   --de_analysis \
-                  --ref_genome differential_expression/hg38_chr20.fa \
-                --ref_annotation differential_expression/gencode.v22.annotation.chr20.gff \
+                  --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
+                --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gff \
                 --direct_rna --minimap2_index_opts '-k 15' \
-                --ref_transcriptome differential_expression/ref_transcriptome.fasta \
+                --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/ref_transcriptome.fasta \
                 --sample_sheet test_data/sample_sheet.csv \
-                -c demo.nextflow.config"
+                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "differential_expression_gff3"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
-              NF_WORKFLOW_OPTS: "--fastq  differential_expression/differential_expression_fastq \
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_WORKFLOW_OPTS: "--fastq   ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
                   --transcriptome-source precomputed \
                   --de_analysis \
-                  --ref_genome differential_expression/hg38_chr20.fa \
-                --ref_annotation differential_expression/gencode.v22.annotation.chr20.gff3 \
+                  --ref_genome  ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
+                --ref_annotation  ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gff3 \
                 --direct_rna --minimap2_index_opts '-k 15' \
-                --ref_transcriptome differential_expression/ref_transcriptome.fasta \
+                --ref_transcriptome  ${CI_PROJECT_NAME}/data/differential_expression/ref_transcriptome.fasta \
                 --sample_sheet test_data/sample_sheet.csv \
-                -c demo.nextflow.config"
+                -c  ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "ncbi_gzip"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf differential_expression_ncbi.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O  ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C  ${CI_PROJECT_NAME}/data/ &&  wget -O  ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
               NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \
-                  --fastq  differential_expression_ncbi/differential_expression_fastq \
+                  --fastq   ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
                   --transcriptome-source precomputed \
                   --de_analysis \
-                  --ref_genome differential_expression_ncbi/GRCh38.p14.NCBI_test.fna.gz \
-                --ref_annotation differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \
+                  --ref_genome  ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14.NCBI_test.fna.gz \
+                --ref_annotation  ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \
                 --direct_rna --minimap2_index_opts '-w 25' \
                --sample_sheet test_data/sample_sheet.csv \
-               -c demo.nextflow.config"
+               -c  ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome           
         - if: $MATRIX_NAME == "ncbi_no_gene_id"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf differential_expression_ncbi.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C  ${CI_PROJECT_NAME}/data/ &&  wget -O  ${CI_PROJECT_NAME}/data/demo.nextflow.config  https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
               NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \
-                --fastq differential_expression_ncbi/differential_expression_fastq \
+                --fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
                 --transcriptome-source precomputed --de_analysis \
-                --ref_genome differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
-                --ref_annotation differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
-                --direct_rna --ref_transcriptome differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
+                --ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
+                --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
+                --direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
                 --transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv \
-                -c demo.nextflow.config"
+                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome      
         - if: $MATRIX_NAME == "ensembl_with_versions"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C  ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
               NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \
-                --fastq differential_expression/differential_expression_fastq \
-                --transcriptome-source precomputed --de_analysis \
-                --ref_genome differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
-                --ref_annotation differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \
-                --direct_rna --ref_transcriptome differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.gz \
+                --fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
+                --de_analysis \
+                --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
+                --ref_annotation  ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \
+                --direct_rna --ref_transcriptome  ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.gz \
                 --transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv \
-                -c demo.nextflow.config"
+                -c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "differential_expression_mouse"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression_mouse.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_mouse.tar.gz && tar -xzvf differential_expression_mouse.tar.gz &&  wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
+              NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression_mouse.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_mouse.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/differential_expression_mouse.tar.gz -C ${CI_PROJECT_NAME}/data/&&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
               NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \
-                --fastq differential_expression_mouse/differential_expression_fastq \
+                --fastq  ${CI_PROJECT_NAME}/data/differential_expression_mouse/differential_expression_fastq \
                 --transcriptome-source precomputed --de_analysis \
-                --ref_genome differential_expression_mouse/GRCm39.genome.fa.gz \
-                --ref_annotation differential_expression_mouse/gencode.vM33.annotation.gtf \
-                --direct_rna --ref_transcriptome differential_expression_mouse/gencode.vM33.transcripts.fa.gz \
-                --transcriptome_assembly false --sample_sheet differential_expression_mouse/sample_sheet.csv \
-                -c demo.nextflow.config"
+                --ref_genome  ${CI_PROJECT_NAME}/data/differential_expression_mouse/GRCm39.genome.fa.gz \
+                --ref_annotation  ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.annotation.gtf \
+                --direct_rna --ref_transcriptome  ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.transcripts.fa.gz \
+                --transcriptome_assembly false --sample_sheet  ${CI_PROJECT_NAME}/data/differential_expression_mouse/sample_sheet.csv \
+                -c  ${CI_PROJECT_NAME}/data/demo.nextflow.config"
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,merge_transcriptomes,assemble_transcripts,decompress_annotation,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam    
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 2ca2e89..26625ea 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -8,7 +8,7 @@ repos:
         always_run: true
         pass_filenames: false
         additional_dependencies:
-          - epi2melabs>=0.0.51
+          - epi2melabs>=0.0.53
       - id: build_models
         name: build_models
         entry: datamodel-codegen --strict-nullable --base-class workflow_glue.results_schema_helpers.BaseModel --use-schema-description --disable-timestamp --input results_schema.yml --input-file-type openapi --output bin/workflow_glue/results_schema.py
diff --git a/README.md b/README.md
index 8b1b42e..7a53e55 100644
--- a/README.md
+++ b/README.md
@@ -122,6 +122,13 @@ input_reads.fastq   ─── input_directory  ─── input_directory
 | analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False |
 
 
+### Output Options
+
+| Nextflow parameter name  | Type | Description | Help | Default |
+|--------------------------|------|-------------|------|---------|
+| out_dir | string | Directory for output of all user-facing files. |  | output |
+
+
 ### Sample Options
 
 | Nextflow parameter name  | Type | Description | Help | Default |
@@ -174,13 +181,6 @@ input_reads.fastq   ─── input_directory  ─── input_directory
 | isoform_table_nrows | integer | Maximum rows to dispay in the isoform report table |  | 5000 |
 
 
-### Miscellaneous Options
-
-| Nextflow parameter name  | Type | Description | Help | Default |
-|--------------------------|------|-------------|------|---------|
-| disable_ping | boolean | Enable to prevent sending a workflow ping. |  | False |
-
-
 
 
 
diff --git a/bin/workflow_glue/check_bam_headers_in_dir.py b/bin/workflow_glue/check_bam_headers_in_dir.py
index d582de1..5297658 100755
--- a/bin/workflow_glue/check_bam_headers_in_dir.py
+++ b/bin/workflow_glue/check_bam_headers_in_dir.py
@@ -8,9 +8,10 @@ import pysam
 from .util import get_named_logger, wf_parser  # noqa: ABS101
 
 
-def get_sq_lines(xam_file):
-    """Extract the `@SQ` lines from the header of a XAM file."""
-    return pysam.AlignmentFile(xam_file, check_sq=False).header["SQ"]
+def get_sq_hd_lines(xam_file):
+    """Extract the `@SQ` and `@HD` lines from the header of a XAM file."""
+    alignments = pysam.AlignmentFile(xam_file, check_sq=False)
+    return alignments.header["SQ"], alignments.header["HD"]
 
 
 def main(args):
@@ -27,10 +28,17 @@ def main(args):
     # Set `is_unaligned` accordingly. If there are mixed headers (either with some files
     # containing `@SQ` lines and some not or with different files containing different
     # `@SQ` lines), set `mixed_headers` to `True`.
+    # Also check if there is the SO line, to validate whether the file is (un)sorted.
     first_sq_lines = None
     mixed_headers = False
+    sorted_xam = False
     for xam_file in target_files:
-        sq_lines = get_sq_lines(xam_file)
+        sq_lines, hd_lines = get_sq_hd_lines(xam_file)
+        # Check if it is sorted.
+        # When there is more than one BAM, merging/sorting
+        # will happen regardless of this flag.
+        if hd_lines.get('SO') == 'coordinate':
+            sorted_xam = True
         if first_sq_lines is None:
             # this is the first file
             first_sq_lines = sq_lines
@@ -46,7 +54,9 @@ def main(args):
     # write `is_unaligned` and `mixed_headers` out so that they can be set as env.
     # variables
     sys.stdout.write(
-        f"IS_UNALIGNED={int(is_unaligned)};MIXED_HEADERS={int(mixed_headers)}"
+        f"IS_UNALIGNED={int(is_unaligned)};" +
+        f"MIXED_HEADERS={int(mixed_headers)};" +
+        f"IS_SORTED={int(sorted_xam)}"
     )
     logger.info(f"Checked (u)BAM headers in '{args.input_path}'.")
 
diff --git a/bin/workflow_glue/check_sample_sheet.py b/bin/workflow_glue/check_sample_sheet.py
index 62e3483..7fa1605 100755
--- a/bin/workflow_glue/check_sample_sheet.py
+++ b/bin/workflow_glue/check_sample_sheet.py
@@ -49,6 +49,21 @@ def main(args):
     try:
         encoding = determine_codec(args.sample_sheet)
         with open(args.sample_sheet, "r", encoding=encoding) as f:
+            try:
+                # Excel files don't throw any error until here
+                csv.Sniffer().sniff(f.readline())
+                f.seek(0)  # return to initial position again
+            except Exception as e:
+                # Excel fails with UniCode error
+                sys.stdout.write(
+                    "The sample sheet doesn't seem to be a CSV file.\n"
+                    "The sample sheet has to be a CSV file.\n"
+                    "Please verify that the sample sheet is a CSV file.\n"
+                    f"Parsing error: {e}"
+                 )
+
+                sys.exit()
+
             csv_reader = csv.DictReader(f)
             n_row = 0
             for row in csv_reader:
diff --git a/bin/workflow_glue/check_xam_index.py b/bin/workflow_glue/check_xam_index.py
new file mode 100755
index 0000000..4d027d2
--- /dev/null
+++ b/bin/workflow_glue/check_xam_index.py
@@ -0,0 +1,43 @@
+"""Validate a single (u)BAM file index."""
+
+from pathlib import Path
+import sys
+
+import pysam
+
+from .util import get_named_logger, wf_parser  # noqa: ABS101
+
+
+def validate_xam_index(xam_file):
+    """Use fetch to validate the index.
+
+    Invalid indexes will fail the call with a ValueError:
+    ValueError: fetch called on bamfile without index
+    """
+    alignments = pysam.AlignmentFile(xam_file, check_sq=False)
+    try:
+        alignments.fetch()
+        has_valid_index = True
+    except ValueError:
+        has_valid_index = False
+    return has_valid_index
+
+
+def main(args):
+    """Run the entry point."""
+    logger = get_named_logger("checkBamIdx")
+
+    # Check if a XAM has a valid index
+    has_valid_index = validate_xam_index(args.input_xam)
+    # write `has_valid_index` out so that they can be set as env.
+    sys.stdout.write(
+        f"HAS_VALID_INDEX={int(has_valid_index)}"
+    )
+    logger.info(f"Checked (u)BAM index for: '{args.input_xam}'.")
+
+
+def argparser():
+    """Argument parser for entrypoint."""
+    parser = wf_parser("check_bam_headers")
+    parser.add_argument("input_xam", type=Path, help="Path to target XAM")
+    return parser
diff --git a/docs/06_input_parameters.md b/docs/06_input_parameters.md
index 22b58a2..cc538fd 100644
--- a/docs/06_input_parameters.md
+++ b/docs/06_input_parameters.md
@@ -11,6 +11,13 @@
 | analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False |
 
 
+### Output Options
+
+| Nextflow parameter name  | Type | Description | Help | Default |
+|--------------------------|------|-------------|------|---------|
+| out_dir | string | Directory for output of all user-facing files. |  | output |
+
+
 ### Sample Options
 
 | Nextflow parameter name  | Type | Description | Help | Default |
@@ -63,10 +70,3 @@
 | isoform_table_nrows | integer | Maximum rows to dispay in the isoform report table |  | 5000 |
 
 
-### Miscellaneous Options
-
-| Nextflow parameter name  | Type | Description | Help | Default |
-|--------------------------|------|-------------|------|---------|
-| disable_ping | boolean | Enable to prevent sending a workflow ping. |  | False |
-
-
diff --git a/lib/ingress.nf b/lib/ingress.nf
index 0e492f4..993b811 100644
--- a/lib/ingress.nf
+++ b/lib/ingress.nf
@@ -24,7 +24,8 @@ def is_target_file(Path file, List extensions) {
 
 
 /**
- * Take a channel of the shape `[meta, reads, path-to-stats-dir | null]` and extract the
+ * Take a channel of the shape `[meta, reads, path-to-stats-dir | null]` (or
+ * `[meta, [reads, index], path-to-stats-dir | null]` in the case of XAM) and extract the
  * run IDs from the `run_ids` file in the stats directory into the metamap. If the path
  * to the stats dir is `null`, add an empty list.
  *
@@ -58,6 +59,53 @@ def add_run_IDs_to_meta(ch) {
 }
 
 
+/**
+ * Take a channel of the shape `[meta, reads, path-to-stats-dir | null]` and do the
+ * following:
+ *  - For `fastcat`, extract the number of reads from the `n_seqs` file.
+ *  - For `bamstats`, extract the number of primary alignments and unmapped reads from
+ *    the `bamstats.flagstat.tsv` file.
+ * Then, add add these metrics to the meta map. If the path to the stats dir is `null`,
+ * set the values to 0 when adding them.
+ *
+ * @param ch: input channel of shape `[meta, reads, path-to-stats-dir | null]`
+ * @return: channel with a list of number of reads added to the metamap
+ */
+def add_number_of_reads_to_meta(ch, String input_type_format) {
+    // extract reads from fastcat stats / bamstats results and add to metadata
+    ch = ch | map { meta, reads, stats ->
+        // Check that stats directory is present.
+        if (stats) {
+            if (input_type_format == "fastq") {
+                // Stats from fastcat
+                Integer n_seqs = stats.resolve("n_seqs").splitText()[0] as Integer
+                // `meta + [...]` returns a new map which is handy to avoid any
+                // modifying-maps-in-closures weirdness
+                // See https://github.com/nextflow-io/nextflow/issues/2660
+                [meta + [n_seqs: n_seqs], reads, stats]
+            } else {
+                // or bamstats
+                ArrayList stats_csv = stats.resolve("bamstats.flagstat.tsv").splitCsv(header: true, sep:'\t')
+                // get primary alignments and unmapped and sum them
+                Integer n_primary = stats_csv["primary"].collect{it as Integer}.sum()
+                Integer n_unmapped = stats_csv["unmapped"].collect{it as Integer}.sum()
+                // `meta + [...]` returns a new map which is handy to avoid any
+                // modifying-maps-in-closures weirdness
+                // See https://github.com/nextflow-io/nextflow/issues/2660
+                [meta + [n_primary: n_primary, n_unmapped: n_unmapped], reads, stats]
+            }
+         } else {
+            // return defaults if stats is not there
+            if (input_type_format == "fastq") {
+                [meta + [n_seqs: null], reads, stats]
+            } else {
+                [meta + [n_primary: null, n_unmapped: null], reads, stats]
+            }
+        }
+    }
+    return ch
+}
+
 /**
  * Take a map of input arguments, find valid FASTQ inputs, and return a channel
  * with elements of `[metamap, seqs.fastq.gz | null, path-to-fastcat-stats | null]`.
@@ -86,7 +134,7 @@ def add_run_IDs_to_meta(ch) {
 def fastq_ingress(Map arguments)
 {
     // check arguments
-    Map margs = parse_arguments(arguments, ["fastcat_extra_args": ""])
+    Map margs = parse_arguments("fastq_ingress", arguments, ["fastcat_extra_args": ""])
 
     ArrayList fq_extensions = [".fastq", ".fastq.gz", ".fq", ".fq.gz"]
 
@@ -109,7 +157,9 @@ def fastq_ingress(Map arguments)
     // add sample sheet entries without barcode dirs to the results channel and extract
     // the run IDs into the metamaps before returning
     ch_result = ch_result.mix(input.missing.map { [*it, null] })
-    return add_run_IDs_to_meta(ch_result)
+    ch_result_run_IDs = add_run_IDs_to_meta(ch_result)
+    // add number of reads after potential filtering under the field n_seqs
+    return add_number_of_reads_to_meta(ch_result_run_IDs, "fastq")
 }
 
 
@@ -118,7 +168,8 @@ def fastq_ingress(Map arguments)
  * with elements of `[metamap, reads.bam | null, path-to-bamstats-results | null]`.
  * The second item is `null` for sample sheet entries without a matching barcode
  * directory or samples containing only uBAM files when `keep_unaligned` is `false`.
- * The last item is `null` if `bamstats` was not run (it is only run when `stats: true`).
+ * The last item is `null` if `bamstats` was not run (it is only run when `stats:
+ * true`).
  *
  * @param arguments: map with arguments containing
  *  - "input": path to either: (i) input (u)BAM file, (ii) top-level directory
@@ -129,6 +180,9 @@ def fastq_ingress(Map arguments)
  *  - "analyse_unclassified": boolean whether to keep unclassified reads
  *  - "stats": boolean whether to run `bamstats`
  *  - "keep_unaligned": boolean whether to include uBAM files
+ *  - "return_fastq": boolean whether to convert to FASTQ (this will always run
+ *    `fastcat`)
+ *  - "fastcat_extra_args": string with extra arguments to pass to `fastcat`
  *  - "required_sample_types": list of required sample types in the sample sheet
  *  - "watch_path": boolean whether to use `watchPath` and run in streaming mode
  * @return: channel of `[Map(alias, barcode, type, ...), Path|null, Path|null]`.
@@ -142,7 +196,11 @@ def fastq_ingress(Map arguments)
 def xam_ingress(Map arguments)
 {
     // check arguments
-    Map margs = parse_arguments(arguments, ["keep_unaligned": false])
+    Map margs = parse_arguments(
+        "xam_ingress",
+        arguments,
+        ["keep_unaligned": false, "return_fastq": false, "fastcat_extra_args": ""]
+    )
 
     // we only accept BAM or uBAM for now (i.e. no SAM or CRAM)
     ArrayList xam_extensions = [".bam", ".ubam"]
@@ -153,11 +211,31 @@ def xam_ingress(Map arguments)
     ch_result = input.dirs
     | map { meta, path -> [meta, get_target_files_in_dir(path, xam_extensions)] }
     | mix(input.files)
+    | map{
+        // If there is more than one BAM in each folder we ignore
+        // the indices. For single BAM we add it as a string to the
+        // metadata for later use. If then the BAM returns as position
+        // sorted, the index will be used.
+        meta, paths -> 
+        boolean is_array = paths instanceof ArrayList
+        String xai_fn
+        // Using `.uri` or `.Uri()` leads to S3 paths to be prefixed with `s3:///`
+        // instead of `s3://`, causing the workflow to not find the index file.
+        // `.toUriString()` returns the correct path.
+        if (!is_array){
+            def xai = file(paths.toUriString() + ".bai")
+            if (xai.exists()){
+                xai_fn = xai.toUriString()
+            }
+        }
+        [meta + [xai_fn: xai_fn], paths]
+    }
     | checkBamHeaders
-    | map { meta, paths, is_unaligned_env, mixed_headers_env ->
+    | map { meta, paths, is_unaligned_env, mixed_headers_env, is_sorted_env ->
         // convert the env. variables from strings ('0' or '1') into bools
         boolean is_unaligned = is_unaligned_env as int as boolean
         boolean mixed_headers = mixed_headers_env as int as boolean
+        boolean is_sorted = is_sorted_env as int as boolean
         // throw an error if there was a sample with mixed headers
         if (mixed_headers) {
             error "Found mixed headers in (u)BAM files of sample '${meta.alias}'."
@@ -165,7 +243,7 @@ def xam_ingress(Map arguments)
         // add `is_unaligned` to the metamap (note the use of `+` to create a copy of
         // `meta` to avoid modifying every item in the channel;
         // https://github.com/nextflow-io/nextflow/issues/2660)
-        [meta + [is_unaligned: is_unaligned], paths]
+        [meta + [is_unaligned: is_unaligned, is_sorted: is_sorted], paths]
     }
     | branch { meta, paths ->
         // set `paths` to `null` for uBAM samples if unallowed (they will be added to
@@ -176,13 +254,12 @@ def xam_ingress(Map arguments)
         }
         // get the number of files (`paths` can be a list, a single path, or `null`)
         int n_files = paths instanceof List ? paths.size() : (paths ? 1 : 0)
-        // Preparations finished; we can do the branching now. There will be 3 branches
+        // Preparations finished; we can do the branching now. There will be 5 branches
         // depending on the number of files per sample and whether the reads are already
         // aligned:
-        // * no_op_needed: no need to do anything; just add to the final results channel
-        //   downstream
-        //   - no files
-        //   - a single unaligned file
+        // * no files: no need to do anything
+        // * indexed: a single sorted and indexed BAM file. Index will be validated.
+        // * to_index: a single sorted, but not indexed, BAM file
         // * to_catsort: `samtools cat` into `samtools sort`
         //  - a single aligned file
         //  - more than one unaligned file
@@ -191,12 +268,36 @@ def xam_ingress(Map arguments)
         //    open file descriptors)
         // * to_merge: flatMap > sort > group > merge
         //  - between 1 and `N_OPEN_FILES_LIMIT` aligned files
-        no_op_needed: (n_files == 0) || (n_files == 1 && meta["is_unaligned"])
+        no_files: n_files == 0
+        indexed: \
+            n_files == 1 && (meta["is_unaligned"] || meta["is_sorted"]) && meta["xai_fn"]
+        to_index: 
+            n_files == 1 && (meta["is_unaligned"] || meta["is_sorted"]) && !meta["xai_fn"]
         to_catsort: \
             (n_files == 1) || (n_files > N_OPEN_FILES_LIMIT) || meta["is_unaligned"]
         to_merge: true
     }
 
+    if (margs["return_fastq"]) {
+        // only run samtools fastq on samples with at least one file
+        ch_to_fastq = ch_result.indexed.mix(
+            ch_result.to_index,
+            ch_result.to_merge,
+            ch_result.to_catsort
+        )
+    
+        // input.missing: sample sheet entries without barcode dirs
+        ch_result = input.missing
+        | mix(ch_result.no_files)
+        | map { [*it, null] }
+        | mix(bamToFastq(ch_to_fastq, margs["fastcat_extra_args"]))
+        | map{
+            meta, path, stats ->
+            [meta.findAll { it.key !in ['xai_fn', 'is_sorted'] }, path, stats]
+        }
+        return add_number_of_reads_to_meta(add_run_IDs_to_meta(ch_result), "fastq")
+    }
+
     // deal with samples with few-enough files for `samtools merge` first
     ch_merged = ch_result.to_merge
     | flatMap { meta, paths -> paths.collect { [meta, it] } }
@@ -208,9 +309,59 @@ def xam_ingress(Map arguments)
     ch_catsorted = ch_result.to_catsort
     | catSortBams
 
-    ch_result = input.missing
+    // Validate the index of the input BAM.
+    // If the input BAM index is invalid, regenerate it.
+    // First separate the BAM from the null input channels.
+    ch_to_validate = ch_result.indexed
+    | map{
+        meta, paths ->
+        bai = paths && meta.xai_fn ? file(meta.xai_fn) : null
+        [meta, paths, bai]
+    }
+    | branch {
+        meta, paths, bai ->
+        to_validate: paths && bai
+        no_op_needed: true
+    }
+
+    // Validate non-null files with index
+    ch_validated = ch_to_validate.to_validate
+    | validateIndex
+    | branch {
+        meta, bam, bai, has_valid_index_env -> 
+        boolean has_valid_index = has_valid_index_env as int as boolean
+        // Split if it is a valid index
+        valid_idx: has_valid_index
+            return [meta, bam, bai]
+        invalid_idx: true
+            return [meta, bam]
+    }
+
+    // Create channel for no_op needed (null channels and valid indexes)
+    ch_no_op = ch_validated.valid_idx
+    | mix(ch_to_validate.no_op_needed)
+
+    // Re-index sorted-not-indexed BAM file
+    ch_indexed = ch_result.to_index
+    | mix( ch_validated.invalid_idx )
+    | samtools_index
+
+    // Add extra null for the missing index to input.missing
+    // as well as the missing metadata.
+    // input.missing: sample sheet entries without barcode dirs
+    ch_missing = input.missing
     | mix(
-        ch_result.no_op_needed,
+        ch_result.no_files,
+    )
+    | map{
+        meta, paths ->
+        [meta + [xai_fn: null, is_sorted: false], paths, null]
+    }
+
+    // Combine all possible inputs
+    ch_result = ch_missing | mix(
+        ch_no_op,
+        ch_indexed,
         ch_merged,
         ch_catsorted,
     )
@@ -218,19 +369,97 @@ def xam_ingress(Map arguments)
     // run `bamstats` if requested
     if (margs["stats"]) {
         // branch and run `bamstats` only on the non-`null` paths
-        ch_result = ch_result.branch { meta, path ->
+        ch_result = ch_result.branch { meta, path, index ->
             has_reads: path
             is_null: true
         }
         ch_bamstats = bamstats(ch_result.has_reads)
-        ch_result = add_run_IDs_to_meta(ch_bamstats) | mix(ch_result.is_null)
+
+        // the channel comes from xam_ingress also have the BAM index in it.
+        // Handle this by placing them in a nested array, maintaining the structure 
+        // from fastq_ingress. We do not use variable name as assigning variable
+        // name with a tuple not matching (e.g. meta, bam, bai, stats <- [meta, bam, stats] )
+        // causes the workflow to crash.
+        ch_result = ch_bamstats
+        | map{
+            it[3] ? [it[0], [it[1], it[2]], it[3]] : it
+        }
+        | add_run_IDs_to_meta
+        | map{
+            it.flatten()
+        }
+        | mix(
+            ch_result.is_null.map{it + [null]}
+        )
     } else {
         // add `null` instead of path to `bamstats` results dir
-        ch_result = ch_result | map { meta, bam -> [meta, bam, null] }
+        ch_result = ch_result | map { meta, bam, bai -> [meta, bam, bai, null] }
     }
+
+    // Remove metadata that are unnecessary downstream:
+    // meta.xai_fn: not needed, as it will be part of the channel as a file
+    // meta.is_sorted: if data are aligned, they will also be sorted/indexed
+    //
+    // The output meta can contain the following flags:
+    // [
+    //     barcode: always present
+    //     type: always present
+    //     run_id: always present, but can be empty (i.e. `[]`)
+    //     alias: always present
+    //     n_primary: always present, but can be `null`
+    //     n_unmapped: always present, but can be `null`
+    //     is_unaligned: present if there is a (u)BAM file
+    // ]
+    ch_result = add_number_of_reads_to_meta(
+        ch_result
+            | map{
+                meta, bam, bai, stats ->
+                [meta.findAll { it.key !in ['xai_fn', 'is_sorted'] }, [bam, bai], stats]
+            }, 
+        "xam"
+    )
+    | map{
+        it.flatten()
+    }
+
     return ch_result
 }
 
+process bamToFastq {
+    label "ingress"
+    label "wf_common"
+    cpus 4
+    memory "2 GB"
+    input:
+        tuple val(meta), path(bams, stageAs: "input_dir/reads*.bam")
+        val extra_args
+    output: tuple val(meta), path("seqs.fastq.gz"), path("fastcat_stats")
+    script:
+    """
+    mkdir fastcat_stats
+
+    # Save file as compressed fastq
+    fastcat \
+        -s ${meta["alias"]} \
+        -r >(bgzip -c > fastcat_stats/per-read-stats.tsv.gz) \
+        -f fastcat_stats/per-file-stats.tsv \
+        --histograms histograms \
+        $extra_args \
+        <( 
+            samtools cat -b <(find input_dir -name 'reads*.bam') | \
+            samtools fastq - -n -T '*' -o - -0 - 
+        ) \
+    | bgzip -c > seqs.fastq.gz
+
+    mv histograms/* fastcat_stats
+
+    # extract the run IDs and number of sequences (n_seqs) from the per-read stats
+    csvtk freq -tf runid fastcat_stats/per-read-stats.tsv.gz \
+    | csvtk del-header \
+    | tee >(cut -f 1 | sort > "fastcat_stats/run_ids") \
+    | awk 'BEGIN{n=0}; {n+=\$2}; END{print n}' > "fastcat_stats/n_seqs"
+    """
+}
 
 process checkBamHeaders {
     label "ingress"
@@ -246,6 +475,7 @@ process checkBamHeaders {
             path("input_dir/reads*.bam", includeInputs: true),
             env(IS_UNALIGNED),
             env(MIXED_HEADERS),
+            env(IS_SORTED),
         )
     script:
     """
@@ -255,17 +485,41 @@ process checkBamHeaders {
 }
 
 
+process validateIndex {
+    label "ingress"
+    label "wf_common"
+    cpus 1
+    memory "2 GB"
+    input: tuple val(meta), path("reads.bam"), path("reads.bam.bai")
+    output:
+        // set the two env variables by `eval`-ing the output of the python script
+        // checking the XAM headers
+        tuple(
+            val(meta),
+            path("reads.bam", includeInputs: true),
+            path("reads.bam.bai", includeInputs: true),
+            env(HAS_VALID_INDEX)
+        )
+    script:
+    """
+    workflow-glue check_xam_index reads.bam > env.vars
+    source env.vars
+    """
+}
+
+
 process mergeBams {
     label "ingress"
     label "wf_common"
     cpus 3
     memory "4 GB"
-    input: tuple val(meta), path("input_bams/reads*.bam")
-    output: tuple val(meta), path("reads.bam")
-    shell:
+    input: tuple val(meta), path("input_bams/reads*.bam"), path("input_bams/reads*.bam.bai")
+    output: tuple val(meta), path("reads.bam"), path("reads.bam.bai")
+    script:
+    def merge_threads = Math.max(1, task.cpus - 1)
     """
-    samtools merge -@ ${task.cpus - 1} \
-        -b <(find input_bams -name 'reads*.bam') -o reads.bam
+    samtools merge -@ ${merge_threads} \
+        -b <(find input_bams -name 'reads*.bam') --write-index -o reads.bam##idx##reads.bam.bai
     """
 }
 
@@ -276,11 +530,12 @@ process catSortBams {
     cpus 4
     memory "4 GB"
     input: tuple val(meta), path("input_bams/reads*.bam")
-    output: tuple val(meta), path("reads.bam")
+    output: tuple val(meta), path("reads.bam"), path("reads.bam.bai")
     script:
+    def sort_threads = Math.max(1, task.cpus - 2)
     """
     samtools cat -b <(find input_bams -name 'reads*.bam') \
-    | samtools sort - -@ ${task.cpus - 2} -o reads.bam
+    | samtools sort - -@ ${sort_threads} --write-index -o reads.bam##idx##reads.bam.bai
     """
 }
 
@@ -291,10 +546,11 @@ process sortBam {
     cpus 3
     memory "4 GB"
     input: tuple val(meta), path("reads.bam")
-    output: tuple val(meta), path("reads.sorted.bam")
+    output: tuple val(meta), path("reads.sorted.bam"), path("reads.sorted.bam.bai")
     script:
+    def sort_threads = Math.max(1, task.cpus - 1)
     """
-    samtools sort -@ ${task.cpus - 1} reads.bam -o reads.sorted.bam
+    samtools sort --write-index -@ ${sort_threads} reads.bam -o reads.sorted.bam##idx##reads.sorted.bam.bai
     """
 }
 
@@ -305,10 +561,11 @@ process bamstats {
     cpus 3
     memory "4 GB"
     input:
-        tuple val(meta), path("reads.bam")
+        tuple val(meta), path("reads.bam"), path("reads.bam.bai")
     output:
         tuple val(meta),
               path("reads.bam"),
+              path("reads.bam.bai"),
               path("bamstats_results")
     script:
         def bamstats_threads = Math.max(1, task.cpus - 1)
@@ -325,8 +582,6 @@ process bamstats {
     | csvtk del-header | sort | uniq > bamstats_results/run_ids
     """
 }
-
-
 /**
  * Run `watchPath` on the input directory and return a channel of shape [metamap,
  * path-to-target-file]. The meta data is taken from the sample sheet in case one was
@@ -473,9 +728,11 @@ process fastcat {
             | bgzip > $out
 
         mv histograms/* $fastcat_stats_outdir
-        # extract the run IDs from the per-read stats
-        csvtk cut -tf runid $fastcat_stats_outdir/per-read-stats.tsv.gz \
-        | csvtk del-header | sort | uniq > $fastcat_stats_outdir/run_ids
+        # extract the run IDs and number of sequences (n_seqs) from the per-read stats
+        csvtk freq -tf runid $fastcat_stats_outdir/per-read-stats.tsv.gz \
+        | csvtk del-header \
+        | tee >(cut -f 1 | sort > "$fastcat_stats_outdir/run_ids") \
+        | awk 'BEGIN{n=0}; {n+=\$2}; END{print n}' > "$fastcat_stats_outdir/n_seqs"
         """
 }
 
@@ -489,7 +746,7 @@ process fastcat {
  *  the argument-parsing to be tailored to a particular ingress function)
  * @return: map of parsed arguments
  */
-Map parse_arguments(Map arguments, Map extra_kwargs=[:]) {
+Map parse_arguments(String func_name, Map arguments, Map extra_kwargs=[:]) {
     ArrayList required_args = ["input"]
     Map default_kwargs = [
         "sample": null,
@@ -502,7 +759,7 @@ Map parse_arguments(Map arguments, Map extra_kwargs=[:]) {
     ArgumentParser parser = new ArgumentParser(
         args: required_args,
         kwargs: default_kwargs + extra_kwargs,
-        name: "fastq_ingress")
+        name: func_name)
     return parser.parse_args(arguments)
 }
 
@@ -726,10 +983,10 @@ def get_sample_sheet(Path sample_sheet, ArrayList required_sample_types) {
     // in STDOUT. Thus, we use the somewhat clunky construct with `concat` and `last`
     // below. This lets the CSV channel only start to emit once the error checking is
     // done.
-    ch_err = validate_sample_sheet(sample_sheet, required_sample_types).map {
+    ch_err = validate_sample_sheet(sample_sheet, required_sample_types).map { stdoutput, sample_sheet_file ->
         // check if there was an error message
-        if (it) error "Invalid sample sheet: ${it}."
-        it
+        if (stdoutput) error "Invalid sample sheet: ${stdoutput}."
+        stdoutput
     }
     // concat the channel holding the path to the sample sheet to `ch_err` and call
     // `.last()` to make sure that the error-checking closure above executes before
@@ -743,13 +1000,14 @@ def get_sample_sheet(Path sample_sheet, ArrayList required_sample_types) {
 /**
  * Python script for validating a sample sheet. The script will write messages
  * to STDOUT if the sample sheet is invalid. In case there are no issues, no
- * message is emitted.
+ * message is emitted. The sample sheet will be published to the output dir.
  *
  * @param: path to sample sheet CSV
  * @param: list of required sample types (optional)
  * @return: string (optional)
  */
 process validate_sample_sheet {
+    publishDir params.out_dir, mode: 'copy', overwrite: true
     cpus 1
     label "ingress"
     label "wf_common"
@@ -757,10 +1015,27 @@ process validate_sample_sheet {
     input:
         path "sample_sheet.csv"
         val required_sample_types
-    output: stdout
+    output: 
+        tuple stdout, path("sample_sheet.csv")
     script:
     String req_types_arg = required_sample_types ? "--required_sample_types "+required_sample_types.join(" ") : ""
     """
     workflow-glue check_sample_sheet sample_sheet.csv $req_types_arg
     """
 }
+
+// Generate an index for an input XAM file
+process samtools_index {
+    cpus 4
+    label "ingress"
+    label "wf_common"
+    memory 4.GB
+    input:
+        tuple val(meta), path("reads.bam")
+    output:
+        tuple val(meta), path("reads.bam"), path("reads.bam.bai")
+    script:
+    """
+    samtools index -@ $task.cpus reads.bam
+    """
+}
diff --git a/nextflow.config b/nextflow.config
index 4eb4c73..3aec771 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -105,7 +105,7 @@ params {
       ]
     agent = null
     container_sha = "shae7c9f184996a384e99be68e790f0612f0c732867"
-	common_sha = "sha1c5febff9f75143710826498b093d9769a5edbb9"
+	common_sha = "sha645176f98b8780851f9c476a064d44c2ae76ddf6"
     }
 }
 
@@ -204,4 +204,5 @@ trace {
 
 env {
     PYTHONNOUSERSITE = 1
+    JAVA_TOOL_OPTIONS = "-Xlog:disable -Xlog:all=warning:stderr"
 }
diff --git a/subworkflows/JAFFAL/gene_fusions.nf b/subworkflows/JAFFAL/gene_fusions.nf
index b186740..7c5919f 100644
--- a/subworkflows/JAFFAL/gene_fusions.nf
+++ b/subworkflows/JAFFAL/gene_fusions.nf
@@ -13,6 +13,9 @@ process jaffal{
         tuple val(sample_id), path("jaffal_output_$sample_id/*jaffa_results.csv"), emit: results_csv
     script:
     """
+    # Jaffa requires a virtual environment
+    # so override JAVA_TOOL_OPTIONS env variable.
+    JAVA_TOOL_OPTIONS=""
     JAFFAOUT=jaffal_output_$sample_id
 
     # JAFFAL exists with status code 1 when there's 0 fusion hits. Prevent this with '||:'