import glob
import pandas as pd

"""
Author: Christiane Hassenrueck
Affiliation: Leibniz Institute for Baltic Sea Research Warnemuende (IOW)
Aim: 
"""

# specify config file
# change this path to the location of the config file with your project specific settings
configfile: "/bio/Common_repositories/workflow_templates/Amplicon_dada2_MiSeq/config/config.yaml"

# snakemake rules

rule all:
	input:
		tax_dada2 = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/tax_dada2_" + config["project"] + ".txt",
		tax_blast = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/tax_blast_" + config["project"] + ".txt"

rule blast_classification:
	input:
		fasta = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/seqs_" + config["project"] + ".fasta"
	output:
		blastout = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/blastout_" + config["project"] + ".txt",
		blastout_filt = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/blastout_filt_" + config["project"] + ".txt"
	params:
		db = config["silva_blast"],
		min_blast = config["min_blast"]
	conda:
		config["wdir"] + "/envs/blast.yaml"
	threads:
		config["parallel_threads"]
	shell:
		"""
		blastn -query {input.fasta} -task blastn -db {params.db} -out {output.blastout} -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs qlen gaps" -max_target_seqs 100 -num_threads {threads}
		# filter blast output
		awk -v threshold={params.min_blast} '($3 + $13) / 2 >= threshold' {output.blastout} > {output.blastout_filt}
		"""

rule parse_blastout:
	input:
		blastout_filt = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/blastout_filt_" + config["project"] + ".txt",
		fasta = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/seqs_" + config["project"] + ".fasta"
	output:
		tax_blast = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/tax_blast_" + config["project"] + ".txt"
	params:
		script_I = config["wdir"] + "/scripts/parse_blast_I.R",
		script_II = config["wdir"] + "/scripts/parse_blast_II.R",
		wrapper = config["wdir"] + "/scripts/parse_blast_wrapper.sh",
		project = config["project"],
		adir = config["adir"],
		analysis = config["analysis"],
		lineage = config["lineage"],
		silva_acc2tax = config["silva_acc2tax"],
		silva_taxrank = config["silva_taxrank"]
	conda:
		config["wdir"] + "/envs/r.yaml"
	threads:
		config["parallel_threads"]
	log:
		config["adir"] + "/Intermediate_results/Logfiles/" + config["analysis"] + "/blast_parsing.log"
	shell:
		"""
		TMPDIR="{params.adir}/Intermediate_results/{params.analysis}/tmp"
		mkdir -p $TMPDIR
		cut -f1 {input.blastout_filt} | uniq | parallel -j{threads} '{params.wrapper} {{}} {input.blastout_filt} {params.script_I} {params.silva_acc2tax} {params.silva_taxrank} $TMPDIR' >> {params.adir}/Intermediate_results/{params.analysis}/tax_blast_tmp.txt
		{params.script_II} -n {params.project} -l {params.lineage} -o "{params.adir}/Intermediate_results/{params.analysis}/tax_blast_tmp.txt" -f {input.fasta} -d {params.adir} -v {params.analysis} &>> {log}
		rm {params.adir}/Intermediate_results/{params.analysis}/tax_blast_tmp.txt
		"""

rule dada2_classification:
	input:
		asvtab = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/asvtab_" + config["project"] + ".txt"
	output:
		tax_dada2 = config["adir"] + "/Intermediate_results/" + config["analysis"] + "/tax_dada2_" + config["project"] + ".txt"
	params:
		script = config["wdir"] + "/scripts/classification.R",
		min_boot = config["min_boot"],
		project = config["project"],
		adir = config["adir"],
		analysis = config["analysis"],
		silva_dada2 = config["silva_dada2"],
		lineage = config["lineage"]
	conda:
		config["wdir"] + "/envs/r.yaml"
	threads:
		config["parallel_threads"]
	log:
		config["adir"] + "/Intermediate_results/Logfiles/" + config["analysis"] + "/classification_R.log"
	shell:
		"""
		{params.script} -n {params.project} -a {input.asvtab} -b {params.min_boot} -t {params.silva_dada2} -l "{params.lineage}" -d {params.adir} -v {params.analysis} -c {threads} &>> {log}
		"""
