Add snakemake logic for coverage bed files

Snakefile

@@ -78,6 +78,7 @@ containers = {
     "debian": "docker://debian:buster-slim",
     "fastqc": "docker://quay.io/biocontainers/fastqc:0.11.7--4",
     "gatk": "docker://broadinstitute/gatk3:3.7-0",
+    "gvcf2coverage": "docker://quay.io/biocontainers/gvcf2coverage:0.1--hfb13731_0",
     "multiqc": "docker://quay.io/biocontainers/multiqc:1.8--py_2",
     "picard": "docker://quay.io/biocontainers/picard:2.22.8--0",
     "python3": "docker://python:3.6-slim",
@@ -128,6 +129,10 @@ rule all:
         cutadapt = (f"{sample}/pre_process/{sample}-{read_group}.txt"
                     for read_group, sample in get_readgroup_per_sample()),
         coverage_stats = coverage_stats,
+        coverage_files = (f"{sample}/vcf/{sample}_{threshold}.bed"
+                          for sample, threshold in itertools.product(
+                              config['samples'], config['coverage_threshold'])
+                          ) if 'coverage_threshold' in config else []
 
 rule create_markdup_tmp:
     """Create tmp directory for mark duplicates"""
@@ -565,3 +570,10 @@ rule stats_tsv:
     container: containers["python3"]
     shell: "python {input.sc} -i {input.stats} > {output}"
 
+
+rule gvcf2coverage:
+    """ Determine coverage from gvcf files """
+    input: rules.gvcf_gather.output.gvcf
+    output: "{sample}/vcf/{sample}_{threshold}.bed"
+    container: containers["gvcf2coverage"]
+    shell: "gvcf2coverage -t {wildcards.threshold} < {input} > {output}"

config/schema.json

@@ -13,6 +13,7 @@
         "scatter_size",
         "female_threshold",
         "bedfile",
+        "coverage_threshold",
         "baitsfile"
     ],
     "properties": {
@@ -67,6 +68,11 @@
         "description": "Bed file of the baits of the capture kit. Used to calculate picard HsMetrics",
         "type": "string"
     },
+    "coverage_threshold": {
+        "description": "One or more thresholds to calculate coverage for, one bedfile per value per sample",
+        "type": "array",
+        "minItems": 1
+    },
     "refflat": {
         "description": "RefFlat file with transcripts",
         "type": "string"

tests/data/config/sample_config_bed_coverage.json

+{
+  "samples": {
+    "micro": {
+      "read_groups": {
+        "lib_01": {
+        "R1": "tests/data/fastq/micro_R1.fq.gz",
+        "R2": "tests/data/fastq/micro_R2.fq.gz"
+        }
+      }
+    }
+  },
+  "reference":"tests/data/reference/ref.fa",
+  "dbsnp": "tests/data/reference/database.vcf.gz",
+  "coverage_threshold": [10, 20, 30, 50],
+  "known_sites": ["tests/data/reference/database.vcf.gz"]
+}

tests/test_dry_run.yml

@@ -39,3 +39,16 @@
   stderr:
     must_not_contain:
       - rror
+
+- name: test-bed-coverage-output
+  tags:
+    - dry-run
+  command: >
+    snakemake -s Snakefile -n --configfile
+    tests/data/config/sample_config_bed_coverage.json
+  stdout:
+    contains:
+    - 'output: micro/vcf/micro_10.bed'
+    - 'output: micro/vcf/micro_20.bed'
+    - 'output: micro/vcf/micro_30.bed'
+    - 'output: micro/vcf/micro_50.bed'