6 add consensus fasta generation

config.yml.template

@@ -10,3 +10,4 @@ min_depth: 4
 log_dir: "logs"
 repeat_bed: "resources/repeat.bed"
 final_bam_extension: ".rmdup"
+consensus_fasta_dir: "results/consensus"

workflow/Snakefile

@@ -3,7 +3,7 @@ import re
 configfile: "config.yml"
 localrules: filtered, index_bam, index_reference, merge_trimmed, retain_list
 
-rule filtered:
+rule filtered: 
     input:
         expand('{vcf_dir}/genome.IF-GF-MM2-RM.vcf.gz', vcf_dir = config["vcf_dir"]),
         expand('{vcf_dir}/genome.IF-GF-MM2-RM.vcf.gz.csi', vcf_dir = config["vcf_dir"])
@@ -14,4 +14,6 @@ include: "rules/align.smk"
 
 include: "rules/call.smk"
 
-include: "rules/filter.smk"
+include: "rules/filter.smk"
+
+include: "rules/consensus.smk"

workflow/rules/consensus.smk

@@ -0,0 +1,22 @@
+
+
+
+def individual_fasta(wildcards):
+    individuals = get_individuals()
+    fasta_list = expand("{fasta_dir}/{individual}.fasta", fasta_dir = config["consensus_fasta_dir"], individual = individuals.keys())
+    return fasta_list 
+
+rule consensus:
+    input:
+        individual_fasta
+
+rule consensus_individual:
+    input:
+        expand("{bam_dir}/{{individual}}{extension}.bam", bam_dir = config["bam_dir"], extension = config["final_bam_extension"]),
+        expand("{bam_dir}/{{individual}}{extension}.bam.bai", bam_dir = config["bam_dir"], extension = config["final_bam_extension"])
+    output:
+        expand("{fasta_dir}/{{individual}}.fasta", fasta_dir = config["consensus_fasta_dir"])
+    log: expand("{logs}/{{individual}}/consensus.log", logs=config["log_dir"])
+    threads: 2
+    shell:
+        "samtools consensus -@ {threads} -o {output} {input[0]} > {log} 2>&1"