;; fastq sensitive short read SARS-CoV-2 mutation calling

(use-modules (ccwl ccwl))

(define fastq-1
  (input "fastq_1"))

(define fastq-2
  (input "fastq_2"))

(define bowtie2-index
  (input "bowtie2_index"
         #:other '((secondary-files . #(".1.bt2" ".2.bt2" ".3.bt2" ".4.bt2"
                                        ".rev.1.bt2" ".rev.2.bt2")))))

(define bwa-mem-index
  (input "bwa_mem_index"
         #:other '((secondary-files . #(".amb" ".ann" ".bwt" ".pac" ".sa")))))

(define reference
  (input "reference"
         #:other '((secondary-files . #(".fai")))))

(define (bam2vcf prefix)
  (workflow (string-append prefix "_bam2vcf")
            (list (pipeline "samtools"
                            (list (command "samtools_view"
                                           (list "samtools" "view" "-h" "-b" "-F" "4"
                                                 (input (string-append prefix "_input"))))
                                  (command "samtools_sort"
                                           (list "samtools" "sort" "-n"))
                                  (command "samtools_fixmate"
                                           (list "samtools" "fixmate" "-m" "-" "-"))
                                  (command "samtools_sort2"
                                           (list "samtools" "sort" "-"))
                                  (command "samtools_markdup"
                                           (list "samtools" "markdup" "-" "-")
                                           #:other '((stdout . "dedup.bam")))))
                  (command "samtools_index"
                           (list "samtools" "index" (input "samtools_stdout"))
                           #:outputs (list (output "indexed_dedup_bam"
                                                   #:binding '((glob . "$(inputs.samtools_stdout.basename)"))
                                                   #:other '((secondary-files . #(".bai")))))
                           #:other `((hints (Initial-work-dir-requirement
                                             (listing . #("$(inputs.samtools_stdout)"))))))
                  (pipeline "freebayes"
                            (list (command "freebayes"
                                           (list "freebayes" "--ploidy" "1"
                                                 "--bam" (input "indexed_dedup_bam"
                                                                #:other '((secondary-files . #("bai"))))
                                                 "-f" reference))
                                  (command "bcftools_view"
                                           (list "bcftools" "view"
                                                 "-i" "QUAL > 10"
                                                 "-Oz")
                                           #:other `((stdout . ,(string-append prefix "_output.vcf.gz"))))))
                  (command "tabix"
                           (list "tabix" (input "freebayes_stdout"))
                           #:outputs (list (output "indexed_vcf"
                                                   #:binding '((glob . "$(inputs.freebayes_stdout.basename)"))
                                                   #:other '((secondary-files . #(".tbi")))))
                           #:other `((hints (Initial-work-dir-requirement
                                             (listing . #("$(inputs.freebayes_stdout)")))))))
            (list (output (string-append prefix "_output")
                          #:source "tabix/indexed_vcf"))))

(define variant-calling
  (workflow "variant-calling"
            (list (command "bowtie"
                           (list "bowtie2" "-p" "16"
                                 "-x" bowtie2-index
                                 "-1" fastq-1
                                 "-2" fastq-2)
                           #:outputs (list (output "bowtie_input" #:type 'stdout)))
                  (bam2vcf "bowtie")
                  (command "bwa_mem"
                           (list "bwa" "mem" "-K" "100000000" "-v" "3" "-t" "16"
                                 "-Y" bwa-mem-index fastq-1 fastq-2)
                           #:outputs (list (output "bwa_mem_input" #:type 'stdout)))
                  (bam2vcf "bwa_mem"))
            (list (output "bowtie_output"
                          #:source "bowtie_bam2vcf/bowtie_output")
                  (output "bwa_mem_output"
                          #:source "bwa_mem_bam2vcf/bwa_mem_output"))))

(write-cwl variant-calling
           "workflows/variant-calling-bowtie-bwa-mem.cwl")