feat: aligners & variant callers

data_structures/read_group.wdl

@@ -416,3 +416,50 @@ task inner_read_group_to_string {
         maxRetries: 1
     }
 }
+
+task read_group_to_array {
+    meta {
+        description: "Converts a `ReadGroup` struct to a `Array[String]` **without any validation**."
+        outputs: {
+            converted_read_group: "Input `ReadGroup` as a `Array[String]`"
+        }
+    }
+
+    parameter_meta {
+        read_group: "`ReadGroup` struct to convert to array"
+    }
+
+    input {
+        ReadGroup read_group
+    }
+
+    String delimiter = "\n"
+
+    command <<<
+        {
+            echo -n "~{"ID:" + read_group.ID}"  # required field. All others optional
+            echo -n "~{delimiter + "BC:" + read_group.BC}"
+            echo -n "~{delimiter + "CN:" + read_group.CN}"
+            echo -n "~{delimiter + "DS:" + read_group.DS}"
+            echo -n "~{delimiter + "DT:" + read_group.DT}"
+            echo -n "~{delimiter + "FO:" + read_group.FO}"
+            echo -n "~{delimiter + "KS:" + read_group.KS}"
+            echo -n "~{delimiter + "LB:" + read_group.LB}"
+            echo -n "~{delimiter + "PG:" + read_group.PG}"
+            echo -n "~{delimiter + "PI:" + read_group.PI}"
+            echo -n "~{delimiter + "PL:" + read_group.PL}"
+            echo -n "~{delimiter + "PM:" + read_group.PM}"
+            echo -n "~{delimiter + "PU:" + read_group.PU}"
+            echo "~{delimiter + "SM:" + read_group.SM}"
+        } >> out.txt
+    >>>
+
+    output {
+        Array[String] converted_read_group = read_lines("out.txt")
+    }
+
+    runtime {
+        container: "ghcr.io/stjudecloud/util:3.0.3"
+        maxRetries: 1
+    }
+}

docker/bwamem2/Dockerfile

@@ -0,0 +1,8 @@
+FROM quay.io/biocontainers/samtools:1.17--h00cdaf9_0 AS samtools
+FROM quay.io/biocontainers/bwa-mem2:2.3--he70b90d_0
+
+COPY --from=samtools /usr/local/bin/ /usr/local/bin/
+COPY --from=samtools /usr/local/lib/ /usr/local/lib/
+COPY --from=samtools /usr/local/libexec/ /usr/local/libexec/
+
+ENTRYPOINT [ "bwa-mem2" ]

docker/bwamem2/package.json

@@ -0,0 +1,5 @@
+{
+    "name": "bwamem2",
+    "version": "2.3",
+    "revision": "0"
+}

docker/hisat2/Dockerfile

@@ -0,0 +1,8 @@
+FROM quay.io/biocontainers/samtools:1.17--h00cdaf9_0 AS samtools
+FROM quay.io/biocontainers/hisat2:2.2.1--hdbdd923_7
+
+COPY --from=samtools /usr/local/bin/ /usr/local/bin/
+COPY --from=samtools /usr/local/lib/ /usr/local/lib/
+COPY --from=samtools /usr/local/libexec/ /usr/local/libexec/
+
+ENTRYPOINT [ "hisat2" ]

docker/hisat2/package.json

@@ -0,0 +1,5 @@
+{
+    "name": "hisat2",
+    "version": "2.2.1",
+    "revision": "0"
+}

docker/minimap2/Dockerfile

@@ -0,0 +1,8 @@
+FROM quay.io/biocontainers/samtools:1.17--h00cdaf9_0 AS samtools
+FROM quay.io/biocontainers/minimap2:2.30--h577a1d6_0
+
+COPY --from=samtools /usr/local/bin/ /usr/local/bin/
+COPY --from=samtools /usr/local/lib/ /usr/local/lib/
+COPY --from=samtools /usr/local/libexec/ /usr/local/libexec/
+
+ENTRYPOINT [ "minimap2" ]

docker/minimap2/package.json

@@ -0,0 +1,5 @@
+{
+    "name": "minimap2",
+    "version": "2.30",
+    "revision": "0"
+}

docker/ngsep/Dockerfile

@@ -0,0 +1,5 @@
+FROM eclipse-temurin:21
+
+RUN wget https://github.com/NGSEP/NGSEPcore/releases/download/v5.1.0/NGSEPcore_5.1.0.jar -O /usr/local/bin/NGSEPcore.jar
+
+ENTRYPOINT [ "java", "-jar", "/usr/local/bin/NGSEPcore.jar" ]

docker/ngsep/package.json

@@ -0,0 +1,5 @@
+{
+    "name": "ngsep",
+    "version": "5.1.0",
+    "revision": "0"
+}

tools/bwa.wdl

@@ -62,7 +62,7 @@ task bwa_aln {
             n_cores=$(nproc)
         fi
         # -1 because samtools uses one more core than `--threads` specifies
-        (( samtools_cores = n_cores - 1 ))
+        (( samtools_cores = n_cores - 1 || 1 ))
 
         mkdir bwa_db
         tar -C bwa_db -xzf "~{bwa_db_tar_gz}" --no-same-owner
@@ -160,7 +160,7 @@ task bwa_aln_pe {
             n_cores=$(nproc)
         fi
         # -1 because samtools uses one more core than `--threads` specifies
-        (( samtools_cores = n_cores - 1 ))
+        (( samtools_cores = n_cores - 1 || 1 ))
 
         mkdir bwa_db
         tar -C bwa_db -xzf "~{bwa_db_tar_gz}" --no-same-owner
@@ -250,14 +250,14 @@ task bwa_mem {
     )
 
     command <<<
-        set -euo pipefail
+        set -xeuo pipefail
 
         n_cores=~{ncpu}
         if ~{use_all_cores}; then
             n_cores=$(nproc)
         fi
         # -1 because samtools uses one more core than `--threads` specifies
-        (( samtools_cores = n_cores - 1 ))
+        (( samtools_cores = n_cores - 1 || 1 ))
 
         mkdir bwa_db
         tar -C bwa_db -xzf "~{bwa_db_tar_gz}" --no-same-owner
@@ -294,7 +294,7 @@ task bwa_mem {
 
     runtime {
         cpu: ncpu
-        memory: "25 GB"
+        memory: "120 GB"
         disks: "~{disk_size_gb} GB"
         container: "ghcr.io/stjudecloud/bwa:0.7.17-2"
         maxRetries: 1

tools/bwamem2.wdl

@@ -0,0 +1,134 @@
+version 1.2
+
+task align {
+    meta {
+        description: "Align DNA sequences against a large reference database using BWA-MEM2"
+        outputs: {
+            alignments: "The output alignment file in SAM format"
+        }
+    }
+
+    parameter_meta {
+        read_one_fastq_gz: "Input gzipped FASTQ read one file to align with BWA-MEM2"
+        reference_index: "The BWA-MEM2 index file for the reference genome"
+        read_group: "The read group string to be included in the SAM header. Format: '@RG\\tID:foo\\tSM:bar'"
+        read_two_fastq_gz: "Input gzipped FASTQ read two file to align with BWA-MEM2"
+        prefix: "Prefix for the BAM file. The extension `.bam` will be added."
+        smart_pairing: "If true, enable smart pairing mode for paired-end reads"
+        skip_mate_rescue: "If true, skip mate rescue for paired-end reads"
+        threads: "Number of threads to use for alignment"
+        modify_disk_size_gb: "Additional disk space to allocate (in GB)"
+        seed_length: "Seed value for the BWA-MEM2 aligner"
+        min_score: "Minimum score threshold for reporting alignments"
+    }
+
+    input {
+        File read_one_fastq_gz
+        File reference_index
+        String read_group
+        File? read_two_fastq_gz
+        String prefix = sub(
+            basename(read_one_fastq_gz),
+            "([_\\.][rR][12])?(\\.subsampled)?\\.(fastq|fq)(\\.gz)?$",
+            ""
+        )
+        Boolean smart_pairing = false
+        Boolean skip_mate_rescue = false
+        Int threads = 4
+        Int modify_disk_size_gb = 0
+        Int seed_length = 19
+        Int min_score = 30
+    }
+
+    String output_name = prefix + ".bam"
+    Int disk_size_gb = ceil((
+            size(read_one_fastq_gz, "GB") + size(read_two_fastq_gz, "GB")
+        ) * 2)
+        + ceil(size(reference_index, "GB"))
+        + 30
+        + modify_disk_size_gb
+
+    command <<<
+        set -euo pipefail
+
+        mkdir bwa_db
+        tar -C bwa_db -xzf "~{reference_index}" --no-same-owner
+        PREFIX=$(basename bwa_db/*.ann ".ann")
+
+        bwa-mem2 mem \
+            -t ~{threads} \
+            -R "~{read_group}" \
+            -k ~{seed_length} \
+            -T ~{min_score} \
+            ~{if smart_pairing then "-p" else ""} \
+            ~{if skip_mate_rescue then "-S" else ""} \
+            bwa_db/"$PREFIX" \
+            "~{read_one_fastq_gz}" \
+            ~{if defined(read_two_fastq_gz) then "\"~{read_two_fastq_gz}\"" else ""} |
+        samtools view -b -o "~{output_name}" -
+
+        rm -r bwa_db
+    >>>
+
+    output {
+        File alignments = output_name
+    }
+
+    requirements {
+        container: "ghcr.io/stjudecloud/bwamem2:branch-minimap2-2.3-0"
+        cpu: threads
+        memory: "~{4 * threads} GB"
+        disks: "~{disk_size_gb} GB"
+    }
+}
+
+task index {
+    meta {
+        description: "Index a reference genome for alignment with minimap2"
+        outputs: {
+            reference_index: "The minimap2 index file for the reference genome"
+        }
+    }
+
+    parameter_meta {
+        reference_fasta: "The reference genome in FASTA format to be indexed"
+        db_name: "The base name for the output index files"
+        modify_disk_size_gb: "Additional disk space to allocate (in GB)"
+    }
+
+    input {
+        File reference_fasta
+        String db_name = "reference"
+        Int modify_disk_size_gb = 0
+    }
+
+    Float input_fasta_size = size(reference_fasta, "GB")
+    Int disk_size_gb = ceil(input_fasta_size * 2) + 10 + modify_disk_size_gb
+    String bwa_db_out_name = db_name + ".tar.gz"
+
+    command <<<
+        set -euo pipefail
+
+        ref_fasta=~{basename(reference_fasta, ".gz")}
+        gunzip -c "~{reference_fasta}" > "$ref_fasta" \
+            || ln -sf "~{reference_fasta}" "$ref_fasta"
+
+        bwa-mem2 index \
+            "$ref_fasta"
+
+        tar -czf "~{bwa_db_out_name}" "$ref_fasta"*
+
+        rm -r "$ref_fasta"
+    >>>
+
+    output {
+        File reference_index = bwa_db_out_name
+    }
+
+    requirements {
+        container: "ghcr.io/stjudecloud/bwamem2:branch-minimap2-2.3-0"
+        cpu: 1
+        memory: "120 GB"
+        disks: "~{disk_size_gb} GB"
+    }
+}