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

bihealth · Jan 22, 2024 · 1e7c5f7 · 1e7c5f7
1 parent 3e3cead
commit 1e7c5f7
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 6 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})
+popd
+
+# 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 -----------------------------------------------------------------------------