add snp-sites to iqtree env and begin qualimap

workflow/Snakefile

@@ -39,6 +39,7 @@ configfile: os.path.join(project_dir, "config", "snakemake.yaml")
 include: rules_dir + "/download.smk"
 include: rules_dir + "/alignment.smk"
 include: rules_dir + "/phylogeny.smk"
+include: rules_dir + "/qc.smk"
 include: rules_dir + "/functions.smk"
 include: rules_dir + "/targets.smk"
 

workflow/envs/iqtree.yaml

@@ -3,3 +3,4 @@ channels:
   - defaults
 dependencies:
   - iqtree=2.0.3
+  - snp-sites=2.5.1

workflow/envs/qc.yaml

@@ -0,0 +1,9 @@
+channels:
+  - defaults
+  - bioconda
+  - conda-forge
+  - anaconda
+dependencies:
+  - samtools=1.9
+  - qualimap=2.2.2d
+  - multiqc=1.9

workflow/rules/qc.smk

@@ -0,0 +1,35 @@
+# -----------------------------------------------------------------------------#
+# Quality Control                                                              #
+# -----------------------------------------------------------------------------#
+
+rule qualimap:
+    """
+    Run qualimap metrics on the output of snippy pairwise.
+    """
+    input:
+        snippy_dir = results_dir + "/snippy_pairwise_assembly/{sample}",
+    output:
+        bamq = results_dir + "/qualimap/{sample}/{sample}.bam",
+        html = results_dir + "/qualimap/{sample}/qualimapReport.html"
+    conda:
+        os.path.join(envs_dir,"qc.yaml")
+    log:
+        os.path.join(logs_dir, "qualimap/{sample}.log")
+    threads:
+        workflow.cores
+    shell:
+        "samtools view -b -q {config[snippy_map_qual]} {input.snippy_dir}/{wildcards.sample}_snippy.bam > {output.bamq}; "
+        "qualimap bamqc -bam {output.bamq} --skip-duplicated -c -outformat 'HTML' -outdir {results_dir}/qualimap/{wildcards.sample} -nt {threads} 1> {log}; "
+
+
+rule multiqc:
+    """
+    Run multiqc on miscellaneous data files.
+    """
+    input:
+        qualimap_asm_dir = expand(results_dir + "/qualimap/{sample}", sample=identify_assembly_sample()),
+    conda:
+        os.path.join(envs_dir,"qc.yaml")
+    shell:
+        "echo {input}; "
+        "multiqc {input.qualimap_asm_dir}"

workflow/rules/targets.smk

@@ -94,3 +94,15 @@ rule test_iqtree:
     input:
         expand(results_dir + "/iqtree/iqtree.core-filter{missing_data}.treefile",
         missing_data = config["snippy_missing_data"])
+
+#------------------------------------------------------------------------------#
+# QC
+#------------------------------------------------------------------------------#
+rule test_qualimap:
+    input:
+        expand(results_dir + "/qualimap/{sample}/qualimapReport.html",
+        sample=identify_assembly_sample())
+
+rule test_multiqc:
+    input:
+        results_dir + "/multiqc/multiqc.html"