Add Samtools depth task

Author: yueyaog

CoverageProfiler/CoverageProfiler.wdl

@@ -3,31 +3,53 @@ version 1.0
 workflow coverageProfile {
     input {
         String sampleName
+        String coverageTool = "Samtools"
         File alignedBam
         File alignedBamIndex
         File referenceFasta
         File referenceDict
         File referenceFai
         File intervals
+        Int MinBaseQuality = 20
+        Int MinMappingQuality = 20
     }
-    call DepthOfCoverage {
-        input:
-            sampleName = sampleName,
-            alignedBam = alignedBam,
-            alignedBamIndex = alignedBamIndex,
-            referenceFasta = referenceFasta,
-            referenceDict = referenceDict,
-            referenceFai = referenceFai,
-            intervals = intervals
+    if (coverageTool =="Samtools") {
+        call SamtoolsDepth {
+            input:
+                sampleName = sampleName,
+                alignedBam = alignedBam,
+                alignedBamIndex = alignedBamIndex,
+                referenceFasta = referenceFasta,
+                referenceDict = referenceDict,
+                referenceFai = referenceFai,
+                intervals = intervals,
+                minBaseQuality = MinBaseQuality,
+                minMappingQuality = MinMappingQuality
+        }
+    }
+    if (coverageTool == "DepthOfCoverage") {
+        call DepthOfCoverage {
+            input:
+                sampleName = sampleName,
+                alignedBam = alignedBam,
+                alignedBamIndex = alignedBamIndex,
+                referenceFasta = referenceFasta,
+                referenceDict = referenceDict,
+                referenceFai = referenceFai,
+                intervals = intervals,
+                minBaseQuality = MinBaseQuality,
+                minMappingQuality = MinMappingQuality
+        }
     }
     output {
-        File coveragebyInterval = DepthOfCoverage.sample_interval_summary
-        Float meanCoverage = DepthOfCoverage.mean_coverage
+        File? DepthOfCoverageIntervalCov = DepthOfCoverage.sample_interval_summary
+        Float? DepthOfCoverageMeanCoverage = DepthOfCoverage.mean_coverage
+        File? SamtoolsDepthProfile = SamtoolsDepth.depth_profile
     }
     meta {
         author: "Yueyao Gao"
         email: "[email protected]"
-        description: "Calculates the depth of coverage of an input sample using GATK's DepthOfCoverage tool."
+        description: "Calculates the depth of coverage of an input sample"
 
     }
 }
@@ -40,6 +62,8 @@ workflow coverageProfile {
         File referenceDict
         File referenceFai
         File intervals
+        Int minBaseQuality
+        Int minMappingQuality
         Int? mem_gb
         Int? cpu
         String gatk_docker = "broadinstitute/gatk:4.5.0.0"
@@ -56,6 +80,9 @@ workflow coverageProfile {
             --input ~{alignedBam} \
             --read-index ~{alignedBamIndex} \
             --reference ~{referenceFasta} \
+            --minimum-mapping-quality ~{minMappingQuality} \
+            --min-base-quality ~{minBaseQuality} \
+            --count-type COUNT_READS \ # Count all reads independently (even if from the same fragment). The only option supported by GATK 4.5.0.0.
             --output output/~{sampleName}
 
         cat output/~{sampleName}.sample_interval_summary | awk 'BEGIN {FS = ","}{print $3}' | tail -n 1 > output/mean_coverage.txt
@@ -68,6 +95,48 @@ workflow coverageProfile {
         memory: machine_mem_mb + " MB"
         cpu: select_first([cpu, 1])
         docker: gatk_docker
-        disks: "local-disk 500 HDD"
+        disks: "local-disk 500 SSD"
+    }
+}
+
+
+    task SamtoolsDepth {
+        input {
+            String sampleName
+            File alignedBam
+            File alignedBamIndex
+            File referenceFasta
+            File referenceDict
+            File referenceFai
+            File intervals
+            Int minBaseQuality
+            Int minMappingQuality
+            Int? mem_gb
+            Int? cpu
+            String samtools_docker = "euformatics/samtools:1.20"
+    }
+    command <<<
+        # Create directories for output
+        mkdir input
+        readlink -f ~{alignedBam} > input/bam_path.txt
+
+        # Run samtools depth
+        samtools depth \
+        -b ~{intervals} \
+        -f input/bam_path.txt \
+        --min-BQ ~{minBaseQuality} \
+        --min-MQ ~{minMappingQuality} \
+        -s \ # Remove Overlapping Reads, Count fragment
+        -o output/~{sampleName}_samtools.depth
+
+    >>>
+    output {
+        File depth_profile = "output/~{sampleName}_samtools.depth"
+    }
+    runtime {
+        memory: select_first([mem_gb, 7]) * 1000 + " MB"
+        cpu: select_first([cpu, 1])
+        docker: samtools_docker
+        disks: "local-disk 500 SSD"
     }
 }