TheiaMeta: Viral Metagenomics workflow

.dockstore.yml

@@ -205,6 +205,11 @@ workflows:
    primaryDescriptorPath: /workflows/phylogenetics/wf_lyveset.wdl
    testParameterFiles:
     - empty.json
+ - name: TheiaMeta_Illumina_PE_PHB
+   subclass: WDL
+   primaryDescriptorPath: /workflows/metagenomics/wf_theiameta_illumina_pe.wdl
+   testParameterFiles:
+    - empty.json
  - name: Snippy_Streamline_PHB
    subclass: WDL
    primaryDescriptorPath: /workflows/phylogenetics/wf_snippy_streamline.wdl

.github/workflows/pytest-workflows.yml

@@ -38,6 +38,11 @@ jobs:
         # For every workflow, test it with MiniWDL and Cromwell
         tag: ["${{ fromJson(needs.changes.outputs.workflows) }}"]
         engine: ["miniwdl", "cromwell"]
+        exclude:
+          - tag: "wf_theiacov_illumina_pe"
+            engine: "miniwdl"
+          - tag: "wf_theiacov_illumina_se"
+            engine: "miniwdl"
     defaults:
       run:
         # Play nicely with miniconda

tasks/alignment/task_minimap2.wdl

@@ -0,0 +1,53 @@
+version 1.0
+
+task minimap2 {
+  meta {
+    description: "Align a query genome to a reference with minimap2"
+  }
+  input {
+    File query1
+    File? query2
+    File reference
+    String samplename
+    String docker = "us-docker.pkg.dev/general-theiagen/staphb/minimap2:2.22" # newer versions seem to be bugged (infinite loop)
+    String mode = "asm20"
+    Boolean output_sam = false
+    Int disk_size = 100
+    Int cpu = 2
+    Int mem = 8
+  }
+  command <<<
+    # Preset options - https://lh3.github.io/minimap2/minimap2.html
+    # Version capture
+    minimap2 --version | tee VERSION
+
+    if [ -z "~{query2}" ] ; then
+      INPUT_QUERY="~{query1}"
+    else
+      INPUT_QUERY="~{query1} ~{query2}"
+    fi
+
+    # Run minimap2 - output can be sam or paf file depending on ~{output_sam}
+    minimap2 \
+      ~{true="-a" false="" output_sam} \
+      -x "~{mode}" \
+      -t "~{cpu}" \
+      "~{reference}" \
+      ${INPUT_QUERY} > "~{samplename}"_minimap2.out 
+
+  >>>
+  output {
+    File minimap2_out = "~{samplename}_minimap2.out"
+    String minimap2_version = read_string("VERSION")
+    String minimap2_docker = "~{docker}"
+  }
+  runtime {
+    docker: "~{docker}"
+    memory: mem + " GB"
+    cpu: cpu
+    disks: "local-disk " + disk_size + " SSD"
+    disk: disk_size + " GB"
+    maxRetries: 3
+    preemptible: 0
+  }
+}

tasks/assembly/task_ivar_consensus.wdl

@@ -13,6 +13,7 @@ task consensus {
     Float? consensus_min_freq 
     Int? consensus_min_depth
     String char_unknown = "N"
+    Boolean skip_N = false
     Int disk_size = 100
   }
   command <<<
@@ -44,7 +45,8 @@ task consensus {
     -q ~{min_qual} \
     -t ~{consensus_min_freq} \
     -m ~{consensus_min_depth} \
-    -n ~{char_unknown}
+    -n ~{char_unknown} \
+    ~{true = "-k" false = "" skip_N} 
 
     # clean up fasta header
     echo ">~{samplename}" > ~{samplename}.ivar.consensus.fasta

tasks/assembly/task_metaspades.wdl

@@ -0,0 +1,46 @@
+version 1.0
+
+task metaspades_pe {
+  input {
+    File read1_cleaned
+    File read2_cleaned
+    String samplename
+    String docker = "us-docker.pkg.dev/general-theiagen/biocontainers/spades:3.12.0--h9ee0642_3"
+    Int disk_size = 100
+    Int cpu = 4
+    Int memory = 16
+    String? kmers
+    String? metaspades_opts
+    Int phred_offset = 33
+  }
+  command <<<
+    metaspades.py --version | head -1 | cut -d ' ' -f 2 | tee VERSION
+    metaspades.py \
+      -1 ~{read1_cleaned} \
+      -2 ~{read2_cleaned} \
+      ~{'-k ' + kmers} \
+      -m ~{memory} \
+      -t ~{cpu} \
+      -o metaspades \
+      --phred-offset ~{phred_offset} \
+      ~{metaspades_opts}
+
+    mv metaspades/contigs.fasta ~{samplename}_contigs.fasta
+
+  >>>
+  output {
+    File assembly_fasta = "~{samplename}_contigs.fasta"
+    String metaspades_version = read_string("VERSION")
+    String metaspades_docker = '~{docker}'
+  }
+  runtime {
+    docker: "~{docker}"
+    memory: "~{memory} GB"
+    cpu: "~{cpu}"
+    disks:  "local-disk " + disk_size + " SSD"
+    disk: disk_size + " GB"
+    maxRetries: 3
+    preemptible: 0
+  }
+}
+

tasks/quality_control/task_fastq_scan.wdl

@@ -7,6 +7,7 @@ task fastq_scan_pe {
     String read1_name = basename(basename(basename(read1, ".gz"), ".fastq"), ".fq")
     String read2_name = basename(basename(basename(read2, ".gz"), ".fastq"), ".fq")
     Int disk_size = 100
+    String docker = "quay.io/biocontainers/fastq-scan:0.4.4--h7d875b9_1"
   }
   command <<<
     # capture date and version
@@ -21,9 +22,11 @@ task fastq_scan_pe {
     fi
 
     # capture forward read stats
-    eval "${cat_reads} ~{read1}" | fastq-scan | tee ~{read1_name}_fastq-scan.json >(jq .qc_stats.read_total > READ1_SEQS)
+    eval "${cat_reads} ~{read1}" | fastq-scan | tee ~{read1_name}_fastq-scan.json
+    cat ~{read1_name}_fastq-scan.json | jq .qc_stats.read_total | tee READ1_SEQS
     read1_seqs=$(cat READ1_SEQS)
-    eval "${cat_reads} ~{read2}" | fastq-scan | tee ~{read2_name}_fastq-scan.json >(jq .qc_stats.read_total > READ2_SEQS)
+    eval "${cat_reads} ~{read2}" | fastq-scan | tee ~{read2_name}_fastq-scan.json
+    cat ~{read2_name}_fastq-scan.json | jq .qc_stats.read_total | tee READ2_SEQS
     read2_seqs=$(cat READ2_SEQS)
 
     # capture number of read pairs
@@ -38,11 +41,12 @@ task fastq_scan_pe {
   output {
     File read1_fastq_scan_report = "~{read1_name}_fastq-scan.json"
     File read2_fastq_scan_report = "~{read2_name}_fastq-scan.json"
-    Int read1_seq = read_string("READ1_SEQS")
-    Int read2_seq = read_string("READ2_SEQS")
+    Int read1_seq = read_int("READ1_SEQS")
+    Int read2_seq = read_int("READ2_SEQS")
     String read_pairs = read_string("READ_PAIRS")
     String version = read_string("VERSION")
     String pipeline_date = read_string("DATE")
+    String fastq_scan_docker = docker
   }
   runtime {
     docker: "us-docker.pkg.dev/general-theiagen/biocontainers/fastq-scan:0.4.4--h7d875b9_1"
@@ -60,6 +64,7 @@ task fastq_scan_se {
     File read1
     String read1_name = basename(basename(basename(read1, ".gz"), ".fastq"), ".fq")
     Int disk_size = 100
+    String docker = "quay.io/biocontainers/fastq-scan:0.4.4--h7d875b9_1"
   }
   command <<<
     # capture date and version
@@ -74,13 +79,15 @@ task fastq_scan_se {
     fi
 
     # capture forward read stats
-    eval "${cat_reads} ~{read1}" | fastq-scan | tee ~{read1_name}_fastq-scan.json >(jq .qc_stats.read_total > READ1_SEQS)
+    eval "${cat_reads} ~{read1}" | fastq-scan | tee ~{read1_name}_fastq-scan.json
+    cat ~{read1_name}_fastq-scan.json | jq .qc_stats.read_total | tee READ1_SEQS
   >>>
   output {
     File fastq_scan_report = "~{read1_name}_fastq-scan.json"
-    Int read1_seq = read_string("READ1_SEQS")
+    Int read1_seq = read_int("READ1_SEQS")
     String version = read_string("VERSION")
     String pipeline_date = read_string("DATE")
+    String fastq_scan_docker = docker
   }
   runtime {
     docker: "us-docker.pkg.dev/general-theiagen/biocontainers/fastq-scan:0.4.4--h7d875b9_1"
@@ -91,4 +98,4 @@ task fastq_scan_se {
     preemptible: 0
     maxRetries: 3
   }
-}
+}

tasks/quality_control/task_ncbi_scrub.wdl

@@ -24,7 +24,7 @@ task ncbi_scrub_pe {
     fi
 
     # dehost reads
-    /opt/scrubber/scripts/scrub.sh -n ${read1_unzip} |& tail -n1 | awk -F" " '{print $1}' > FWD_SPOTS_REMOVED
+    /opt/scrubber/scripts/scrub.sh -i ${read1_unzip} |& tail -n1 | awk -F" " '{print $1}' > FWD_SPOTS_REMOVED
 
     # gzip dehosted reads
     gzip ${read1_unzip}.clean -c > ~{samplename}_R1_dehosted.fastq.gz
@@ -40,7 +40,7 @@ task ncbi_scrub_pe {
     fi
 
     # dehost reads
-    /opt/scrubber/scripts/scrub.sh -n ${read2_unzip} |& tail -n1 | awk -F" " '{print $1}' > REV_SPOTS_REMOVED
+    /opt/scrubber/scripts/scrub.sh -i ${read2_unzip} |& tail -n1 | awk -F" " '{print $1}' > REV_SPOTS_REMOVED
 
     # gzip dehosted reads
     gzip ${read2_unzip}.clean -c > ~{samplename}_R2_dehosted.fastq.gz
@@ -51,6 +51,7 @@ task ncbi_scrub_pe {
     Int read1_human_spots_removed = read_int("FWD_SPOTS_REMOVED")
     Int read2_human_spots_removed = read_int("REV_SPOTS_REMOVED")
     String ncbi_scrub_docker = docker
+
   }
   runtime {
       docker: "~{docker}"
@@ -67,7 +68,7 @@ task ncbi_scrub_se {
   input {
     File read1
     String samplename
-    String docker = "gcr.io/ncbi-sys-gcr-public-research/sra-human-scrubber@sha256:b7dba71079344daea4ea3363e1a67fa54edb7ec65459d039669c68a66d38b140"
+    String docker = "us-docker.pkg.dev/general-theiagen/ncbi/sra-human-scrubber:2.1.0"
     Int disk_size = 100
   }
   String r1_filename = basename(read1)
@@ -85,7 +86,7 @@ task ncbi_scrub_se {
     fi
 
     # dehost reads
-    /opt/scrubber/scripts/scrub.sh -n ${read1_unzip} |& tail -n1 | awk -F" " '{print $1}' > FWD_SPOTS_REMOVED
+    /opt/scrubber/scripts/scrub.sh -i ${read1_unzip} |& tail -n1 | awk -F" " '{print $1}' > FWD_SPOTS_REMOVED
 
     # gzip dehosted reads
     gzip ${read1_unzip}.clean -c > ~{samplename}_R1_dehosted.fastq.gz

tasks/quality_control/task_pilon.wdl

@@ -0,0 +1,43 @@
+version 1.0
+
+task pilon {
+  input {
+    File assembly
+    File bam
+    File bai
+    String samplename
+    String docker = "us-docker.pkg.dev/general-theiagen/biocontainers/pilon:1.24--hdfd78af_0"
+    Int cpu = 4
+    Int memory = 8
+    Int disk_size = 100
+  }
+  command <<<
+    # version capture
+    pilon --version | cut -d' ' -f3 | tee VERSION
+
+    # run pilon
+    pilon \
+    --genome ~{assembly} \
+    --frags ~{bam} \
+    --output ~{samplename} \
+    --outdir pilon \
+    --changes --vcf
+
+  >>>
+  output {
+    File assembly_fasta = "pilon/~{samplename}.fasta"
+    File changes = "pilon/~{samplename}.changes"
+    File vcf = "pilon/~{samplename}.vcf"
+    String pilon_version = read_string("VERSION")
+    String pilon_docker = "~{docker}"
+  }
+  runtime {
+    docker: "~{docker}"
+    memory: "~{memory} GB"
+    cpu: cpu
+    disks:  "local-disk " + disk_size + " SSD"
+    disk: disk_size + " GB" # TES
+    preemptible: 0
+    maxRetries: 3
+  }
+}

tasks/quality_control/task_quast.wdl

@@ -4,17 +4,18 @@ task quast {
   input {
     File assembly
     String samplename
+    Int min_contig_len = 500
     String docker = "us-docker.pkg.dev/general-theiagen/staphb/quast:5.0.2"
     Int disk_size = 100
-    Int memory = 2 # added default value
+    Int mem = 2 # added default value
     Int cpu = 2 # added default value
   }
   command <<<
     # capture date and version
     date | tee DATE
     quast.py --version | grep QUAST | tee VERSION
 
-    quast.py ~{assembly} -o .
+    quast.py ~{assembly} -o . --min-contig ~{min_contig_len}
     mv report.tsv ~{samplename}_report.tsv
 
     python <<CODE
@@ -34,7 +35,13 @@ task quast {
               n50_value.write(line[1])
           if "GC" in line[0]:
             with open("GC_PERCENT", 'wt') as gc_percent:
-              gc_percent.write(line[1])       
+              gc_percent.write(line[1])
+          if "Largest contig" in line[0]:
+            with open("LARGEST_CONTIG", 'wt') as largest_contig:
+              largest_contig.write(line[1])
+          if "# N's per 100 kbp" in line[0]:
+            with open("UNCALLED_BASES", "wt") as uncalled_bases:
+              uncalled_bases.write(line[1])
 
     CODE
 
@@ -46,11 +53,14 @@ task quast {
     Int genome_length = read_int("GENOME_LENGTH")
     Int number_contigs = read_int("NUMBER_CONTIGS")
     Int n50_value = read_int("N50_VALUE")
-    Float gc_percent = read_float("GC_PERCENT")    
+    Float gc_percent = read_float("GC_PERCENT")
+    Int largest_contig = read_int("LARGEST_CONTIG")
+    Float uncalled_bases = read_float("UNCALLED_BASES")
+    String quast_docker = docker    
   }
   runtime {
     docker:  "~{docker}"
-    memory:  "~{memory} GB"
+    memory:  "~{mem} GB"
     cpu:   "~{cpu}"
     disks: "local-disk " + disk_size + " SSD"
     disk: disk_size + " GB"

tasks/quality_control/task_readlength.wdl

@@ -0,0 +1,36 @@
+version 1.0
+
+task readlength {
+  input {
+    File read1
+    File read2
+    Int memory = 8
+    String docker = "us-docker.pkg.dev/general-theiagen/staphb/bbtools:38.76"
+    Int disk_size = 100
+  }
+  command <<<
+    # date and version control
+    date | tee DATE
+
+    readlength.sh in=~{read1} > STDOUT_FORWARD
+    readlength.sh in=~{read2} > STDOUT_REVERSE
+
+    avg_forward=$(cat STDOUT_FORWARD | grep "#Avg:" | cut -f 2)
+    avg_reverse=$(cat STDOUT_REVERSE | grep "#Avg:" | cut -f 2)
+
+    result=$(awk "BEGIN { printf \"%.2f\", ($avg_forward +  $avg_reverse ) / 2 }")
+    echo $result | tee AVERAGE_READ_LENGTH
+  >>>
+  output {
+    Float average_read_length = read_string("AVERAGE_READ_LENGTH")
+  }
+  runtime {
+    docker: "~{docker}"
+    memory: "~{memory} GB"
+    cpu: 4
+    disks:  "local-disk " + disk_size + " SSD"
+    disk: disk_size + " GB" # TES
+    preemptible: 0
+    maxRetries: 3
+  }
+}