mapping: minor change in targets (ave cov for filtered contigs); new rule/targets (issue #84)

workflow/rules/mapping.smk

@@ -311,6 +311,24 @@ rule mapping_assembly_genomecov_average:
     shell:
         "(date && cat {input} | awk -f {params.script} | tail -n+2 > {output} && date) &> {log}"
 
+rule mapping_assembly_genomecov_bg:
+    input:
+        os.path.join(RESULTS_DIR, "mapping/{mtype}/{rtype}/{tool}/{atype}.sr.bam")
+    output:
+        os.path.join(RESULTS_DIR, "mapping/{mtype}/{rtype}/{tool}/{atype}.sr.cov.bg")
+    wildcard_constraints:
+        mtype="|".join(META_TYPES),
+        rtype="|".join(READ_TYPES),
+        tool="|".join(ASSEMBLERS),
+        atype="ASSEMBLY|ASSEMBLY.FILTERED"
+    threads: 1
+    conda:
+        os.path.join(ENV_DIR, "bedtools.yaml")
+    message:
+        "Mapping: compute assembly coverage"
+    shell:
+        "bedtools genomecov -bg -ibam {input} > {output}"
+
 ##################################################
 # Mappability index
 

workflow/steps/mapping.smk

@@ -18,14 +18,22 @@ rule MAPPING:
 
 rule COVERAGE:
     input:
+        # default: map metaG (LR/SR/hybrid) to filtered contigs
         metag=expand(
             os.path.join(RESULTS_DIR, "mapping/metag/{rtype_tool}.cov.ave.txt"),
-            rtype_tool=["{r}/{t}/ASSEMBLY.{r}".format(r=r, t=t) for r, t in READ_ASSEMBLERS]
+            rtype_tool=["{r}/{t}/ASSEMBLY.FILTERED.{r}".format(r=r, t=t) for r, t in READ_ASSEMBLERS]
         ),
+        # if metaT: map metaT (SR) to filtered contigs
         metat=expand(
             os.path.join(RESULTS_DIR, "mapping/metat/{rtype_tool}.cov.ave.txt"),
             rtype_tool=["{r}/{t}/ASSEMBLY.FILTERED.sr".format(r=r, t=t) for r, t in READ_ASSEMBLERS]
-        ) if "metat" in META_TYPES else []
+        ) if "metat" in META_TYPES else [],
+        # metaG/metaT SR coverage
+        srcov=expand(
+            os.path.join(RESULTS_DIR, "mapping/{mtype}/{rtype_tool}.cov.bg"),
+            mtype=META_TYPES,
+            rtype_tool=["{r}/{t}/ASSEMBLY.FILTERED.sr".format(r=r, t=t) for r, t in READ_ASSEMBLERS]
+        )
     output:
         touch("status/coverage.done")