feat: also write out mapping quality bigWig file (#476) (#480)

bihealth · Jan 22, 2024 · dc1d92b · dc1d92b
1 parent c54bba1
commit dc1d92b
Show file tree

Hide file tree

Showing 4 changed files with 57 additions and 12 deletions.
diff --git a/snappy_pipeline/workflows/ngs_mapping/__init__.py b/snappy_pipeline/workflows/ngs_mapping/__init__.py
@@ -1314,8 +1314,10 @@ def _get_output_files_run_work(self):
         yield "vcf_md5", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.md5"
         yield "vcf_tbi", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.tbi"
         yield "vcf_tbi_md5", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.tbi.md5"
-        yield "bw", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw"
-        yield "bw_md5", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw.md5"
+        yield "cov_bw", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw"
+        yield "cov_bw_md5", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw.md5"
+        yield "mq_bw", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw"
+        yield "mq_bw_md5", "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw.md5"
 
     @dictify
     def get_log_file(self, action):

diff --git a/snappy_wrappers/wrappers/maelstrom/bam_collect_doc/wrapper.py b/snappy_wrappers/wrappers/maelstrom/bam_collect_doc/wrapper.py
@@ -50,7 +50,7 @@
 md5sum $(basename {snakemake.output.vcf}) >$(basename {snakemake.output.vcf_md5})
 md5sum $(basename {snakemake.output.vcf_tbi}) >$(basename {snakemake.output.vcf_tbi_md5})
 
-# Convert to bigWig file
+# Convert coverage to bigWig file
 
 bcftools query -f '%CHROM\t%POS[\t%CV]\n' $(basename {snakemake.output.vcf}) \
 | awk -v span=$WINDOW -F $'\t' 'BEGIN {{ OFS=FS; prev=0; }}
@@ -63,12 +63,33 @@
         old3=$3;
         prev=$1;
     }}' \
-> $TMPDIR/out.wig
+> $TMPDIR/out_cov.wig
 cut -f 1-2 {snakemake.config[static_data_config][reference][path]}.fai \
 > $TMPDIR/chrom.sizes
 
-wigToBigWig $TMPDIR/out.wig $TMPDIR/chrom.sizes $(basename {snakemake.output.bw})
-md5sum $(basename {snakemake.output.bw}) >$(basename {snakemake.output.bw_md5})
+wigToBigWig $TMPDIR/out_cov.wig $TMPDIR/chrom.sizes $(basename {snakemake.output.cov_bw})
+md5sum $(basename {snakemake.output.cov_bw}) >$(basename {snakemake.output.cov_bw_md5})
+
+# Convert mapping quality to bigWig file
+
+bcftools query -f '%CHROM\t%POS[\t%MQ]\n' $(basename {snakemake.output.vcf}) \
+| awk -v span=$WINDOW -F $'\t' 'BEGIN {{ OFS=FS; prev=0; }}
+        {{ if (prev != $1) {{
+            printf("variableStep chrom=%s span=%d\n", $1, span);
+        }} else {{
+            printf("%s\t%f\n", old2, old3);
+        }}
+        old2=$2;
+        old3=$3;
+        prev=$1;
+    }}' \
+> $TMPDIR/out_mq.wig
+cut -f 1-2 {snakemake.config[static_data_config][reference][path]}.fai \
+> $TMPDIR/chrom.sizes
+
+wigToBigWig $TMPDIR/out_mq.wig $TMPDIR/chrom.sizes $(basename {snakemake.output.mq_bw})
+md5sum $(basename {snakemake.output.mq_bw}) >$(basename {snakemake.output.mq_bw_md5})
+
 popd
 
 # Create output links -----------------------------------------------------------------------------

diff --git a/tests/snappy_pipeline/workflows/test_workflows_ngs_mapping.py b/tests/snappy_pipeline/workflows/test_workflows_ngs_mapping.py
@@ -664,8 +664,10 @@ def test_generate_doc_files_step_part_run_get_input_files(ngs_mapping_workflow):
 def test_generate_doc_files_step_part_get_output_files(ngs_mapping_workflow):
     """Tests BamCollectDocStepPart.get_output_files() - run"""
     expected = {
-        "bw": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw",
-        "bw_md5": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw.md5",
+        "cov_bw": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw",
+        "cov_bw_md5": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw.md5",
+        "mq_bw": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw",
+        "mq_bw_md5": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw.md5",
         "vcf": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz",
         "vcf_md5": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.md5",
         "vcf_tbi": "work/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.tbi",
@@ -677,6 +679,8 @@ def test_generate_doc_files_step_part_get_output_files(ngs_mapping_workflow):
             "output/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.vcf.gz.tbi.md5",
             "output/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw",
             "output/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.cov.bw.md5",
+            "output/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw",
+            "output/{mapper}.{library_name}/report/cov/{mapper}.{library_name}.mq.bw.md5",
             "output/{mapper}.{library_name}/log/{mapper}.{library_name}.bam_collect_doc.log",
             "output/{mapper}.{library_name}/log/{mapper}.{library_name}.bam_collect_doc.log.md5",
             "output/{mapper}.{library_name}/log/{mapper}.{library_name}.bam_collect_doc.conda_info.txt",
@@ -779,10 +783,19 @@ def test_ngs_mapping_workflow_files(ngs_mapping_workflow):
         for stats in ("bamstats", "flagstats", "idxstats")
     ]
     expected += [
-        "output/bwa.P00{i}-N1-DNA1-WGS1/report/cov/bwa.P00{i}-N1-DNA1-WGS1.cov.{ext}".format(
+        "output/bwa.P00{i}-N1-DNA1-WGS1/report/cov/bwa.P00{i}-N1-DNA1-WGS1.{ext}".format(
             i=i, ext=ext
         )
-        for ext in ("bw", "bw.md5", "vcf.gz", "vcf.gz.md5", "vcf.gz.tbi", "vcf.gz.tbi.md5")
+        for ext in (
+            "cov.bw",
+            "cov.bw.md5",
+            "cov.vcf.gz",
+            "cov.vcf.gz.md5",
+            "cov.vcf.gz.tbi",
+            "cov.vcf.gz.tbi.md5",
+            "mq.bw",
+            "mq.bw.md5",
+        )
         for i in range(1, 7)
     ]
     expected += [

diff --git a/tests/snappy_pipeline/workflows/test_workflows_ngs_mapping_processed_fastq.py b/tests/snappy_pipeline/workflows/test_workflows_ngs_mapping_processed_fastq.py
@@ -714,10 +714,19 @@ def test_ngs_mapping_workflow_files(ngs_mapping_workflow):
         for stats in ("bamstats", "flagstats", "idxstats")
     ]
     expected += [
-        "output/bwa.{library_name}/report/cov/bwa.{library_name}.cov.{ext}".format(
+        "output/bwa.{library_name}/report/cov/bwa.{library_name}.{ext}".format(
             library_name=library_name, ext=ext
         )
-        for ext in ("bw", "bw.md5", "vcf.gz", "vcf.gz.md5", "vcf.gz.tbi", "vcf.gz.tbi.md5")
+        for ext in (
+            "cov.bw",
+            "cov.bw.md5",
+            "cov.vcf.gz",
+            "cov.vcf.gz.md5",
+            "cov.vcf.gz.tbi",
+            "cov.vcf.gz.tbi.md5",
+            "mq.bw",
+            "mq.bw.md5",
+        )
         for library_name in dna
     ]
     expected += [