shell:
    """bowtie2 --large-index -x {params.BOWTIE2_DB} --end-to-end --threads {threads} --very-sensitive -U {input.fastq} > $(dirname {output.bam})/AlignedToBowtie2DB.sam 2> {log};"""
    """grep @ $(dirname {output.bam})/AlignedToBowtie2DB.sam | awk '!seen[$2]++' > $(dirname {output.bam})/header_nodups.txt;"""
    """grep -v '^@' $(dirname {output.bam})/AlignedToBowtie2DB.sam > $(dirname {output.bam})/AlignedToBowtie2DB.noheader.sam;"""
    """cat $(dirname {output.bam})/header_nodups.txt $(dirname {output.bam})/AlignedToBowtie2DB.noheader.sam > $(dirname {output.bam})/AlignedToBowtie2DB.nodups.sam;"""
    """samtools view -bS -q 1 -h -@ {threads} $(dirname {output.bam})/AlignedToBowtie2DB.nodups.sam | samtools sort - -@ {threads} -o {output.bam};"""
    """samtools index {output.bam};"""
    """rm $(dirname {output.bam})/header_nodups.txt;"""
    """rm $(dirname {output.bam})/AlignedToBowtie2DB.noheader.sam;"""
    """rm $(dirname {output.bam})/AlignedToBowtie2DB.nodups.sam;"""
    """rm $(dirname {output.bam})/AlignedToBowtie2DB.sam;"""

shell:
    "postprocessing.AMPS.r -m def_anc -r {params.extract} -t {threads} -n {input.node_list} || echo 'postprocessing failed for {wildcards.sample}_{wildcards.taxid}' > {output.analysis}  2> {log}"

shell:
    "(samtools view -h {input.bam} || true) | pmdtools --printDS > {output.scores}"

shell:
    "mkdir {output.dir}; "
    "if [ -s {input.species_tax_id} ]; then "
    'cat {input.species_tax_id} | parallel -j {threads} "grep -w {{}} {params.bowtie2_seqid2taxid} | cut -f1 > {output.dir}/{{}}.seq_ids" ; '
    "for i in $(cat {input.species_tax_id}); do xargs --arg-file={output.dir}/${{i}}.seq_ids samtools view -bh {input.bam} --write-index -@ {threads} -o {output.dir}/${{i}}.tax.bam; done >> {log} 2>&1; "
    "find {output.dir} -name '*.tax.bam' | parallel -j {threads} \"mapDamage {params.options} -i {{}} -r {params.BOWTIE2_DB} --merge-reference-sequences -d {output.dir}/mapDamage_{{}}\" >> {log} 2>&1 || true; "
    "for filename in {output.dir}/*.tax.bam; do newname=`echo $filename | sed 's/tax\.//g'`; mv $filename $newname; done >> {log} 2>&1; "
    "mv {output.dir}/mapDamage_{output.dir}/* {output.dir} >> {log} 2>&1; "
    "rm -r {output.dir}/mapDamage_results >> {log} 2>&1; "
    "else echo NO MICROBES TO AUTHENTICATE > {output.dir}/README.txt; fi"

shell:
    "krakenuniq --preload --db {params.DB} --fastq-input {input.fastq} --threads {threads} --output {output.seqs} --report-file {output.report} --gzip-compressed --only-classified-out &> {log}"

shell:
    "{params.exe} {input.report} {input.seqs} &> {log}; "
    "cat {output.seqs} | cut -f 2,3 > {output.krona}; "
    "ktImportTaxonomy {output.krona} -o {output.html} {params.DB} &>> {log}"

shell:
    "unset DISPLAY; malt-run -at SemiGlobal -m BlastN -i {input.fastq} -o {output.rma6} -a {params.gunzipped_sam} -t {threads} -d {input.db} -sup 1 -mq 100 -top 1 -mpi 85.0 -id 85.0 -v &> {log}"

shell:
    "fastqc {input.fastq} --threads {threads} --nogroup --outdir results/FASTQC_BEFORE_TRIMMING &> {log}"

shell:
    "cutadapt {params.adapters} --minimum-length 31 -o {output.fastq} {input.fastq} &> {log}"

shell:
    'echo {input} | tr " " "\n" > {output.html}.fof;'
    "multiqc -c {params.config} -l {output.html}.fof --verbose --force --outdir results/MULTIQC &> {log}"

__author__ = "Per Unneberg"
__copyright__ = "Copyright 2022, Per Unneberg"
__email__ = "per.unneberg@scilifelab.se"
__license__ = "MIT"

import re
from snakemake.shell import shell
import subprocess as sp

out = sp.run(["malt-build", "--help"], capture_output=True)
regex = re.compile("version (?P<major>\d+)\.(?P<minor>\d+)")
m = regex.search(out.stderr.decode())
if m is None:
    # Assume minor version 4
    minor = 4
else:
    minor = int(m.groupdict()["minor"])

a2t_option = "-a2taxonomy" if minor <= 4 else "-a2t"

log = snakemake.log_fmt_shell(stdout=False, stderr=True, append=True)
options = snakemake.params.get("extra", "")
unique_taxids = snakemake.input.unique_taxids

seqid2taxid = snakemake.params.seqid2taxid
nt_fasta = snakemake.params.nt_fasta
accession2taxid = snakemake.params.accession2taxid

output = snakemake

shell(
    "grep -wFf {unique_taxids} {seqid2taxid} > {snakemake.output.seqid2taxid_project}; "
    "cut -f1 {snakemake.output.seqid2taxid_project} > {snakemake.output.seqids_project}; "
    "grep -Ff {snakemake.output.seqids_project} {nt_fasta} | sed 's/>//g' > {snakemake.output.project_headers}; "
    "seqtk subseq {nt_fasta} {snakemake.output.project_headers} > {snakemake.output.project_fasta} {log}; "
    "unset DISPLAY; "
    "malt-build -i {snakemake.output.project_fasta} {a2t_option} {accession2taxid} -s DNA -t {snakemake.threads} -d {snakemake.output.db} {log}"
)