script:
    "../scripts/filter_taxonomy.py"

shell:
    """
    exec 2> {log}

    export BLASTDB={params.taxdb}

    blastdbcmd -db {params.blast_DB} -tax_info -outfmt %T \
    > {output.taxlist}
    """

script:
    "../scripts/make_blast_mask.py"

script:
    "../scripts/apply_blocklist.py"

shell:
    """
    touch {output.mask} 2> {log}
    """

shell:
    """
    touch {output.block} > {log}
    """

shell:
    """
    exec 2> {log}
    if [ -s {input.report} ]; then
      grep -v -f {params.acc_list} {input.report} > {output.report}
    else
      touch {output.report}
    fi
    """

script:
    "../scripts/filter_blast.py"

script:
    "../scripts/min_consensus_filter.py"

shell:
    """
    exec 2> {log}

    if [ -s {input.blast} ]
    then
        # Get list of all OTUs
        OTUs=$(grep "^>" {input.otus} | cut -d";" -f1 | tr -d '>' | sort -u)

        for otu in $OTUs
        do
            size=$(grep -E "^>${{otu}}\>" {input.otus}  | cut -d"=" -f2)
            bhits=$(grep -c -E "^${{otu}};" {input.blast} || true)
            if [ $bhits -eq 0 ]
            then
                # When there is no blast hit
                echo "{wildcards.sample}\t$otu\t$size\t0\t0\t0\t0\t0\t-\t-\t-\t- (1.0)\t../{input.blast}\t../{input.filtered}" >> {output}
            else
                # Otherwise collect and print stats to file
                bit_best=$(grep -E "^${{otu}};" {input.blast} | cut -f5 | cut -d. -f1 | sort -rn | head -n1)
                bit_low=$(grep -E "^${{otu}};" {input.blast} | cut -f5 | cut -d. -f1 | sort -n | head -n1)
                bit_thr=$(($bit_best - {params.bit_diff}))
                shits=$(grep -c -E "^${{otu}}\>" {input.filtered})
                cons=$(grep -E "^${{otu}}\>" {input.lca} | cut -d'\t' -f2-5)

                echo "{wildcards.sample}\t$otu\t$size\t$bhits\t$bit_best\t$bit_low\t$bit_thr\t$shits\t$cons\t../{input.blast}\t../{input.filtered}" >> {output}
            fi
        done
        # Sort by size and add header (just to get hits on top)
        sort -k3,3nr -o {output} {output}
        sed -i '1 i\Sample\tQuery\tCount\tBlast hits\tBest bit-score\tLowest bit-score\tBit-score threshold\tSaved Blast hits\tConsensus\tRank\tTaxid\tDisambiguation\tlink_report\tlink_filtered' {output}

    else
        echo "{wildcards.sample}\t-\t-\t0\t0\t0\t0\t0\t-\t-\t-\t-\t../{input.blast}\t../{input.filtered}" > {output}
        sed -i '1 i\Sample\tQuery\tCount\tBlast hits\tBest bit-score\tLowest bit-score\tBit-score threshold\tSaved Blast hits\tConsensus\tRank\tTaxid\tDisambiguation\tlink_report\tlink_filtered' {output}
    fi
    """

shell:
    """
    head -n 1 {input.report[0]} > {output.agg}
    for i in {input.report}; do 
      cat ${{i}} | tail -n +2 >> {output.agg}
    done
    """

shell:
    """
    exec 2> {log}
    cat {input.report[0]} | head -n 1 > {output.agg}
    for i in {input.report}; do 
        cat ${{i}} | tail -n +2 >> {output.agg}
    done
    """

script:
    "../scripts/summarize_results.py"

shell:
    """
    exec 2> {log}

    cat {input.report[0]} | head -n 1 > {output.agg}
    for i in {input.report}; do 
        cat ${{i}} | tail -n +2 >> {output.agg}
    done
    """

script:
    "../scripts/krona_table.py"

shell:
    """
    exec 2> {log}
    i=0
    for file in {input.report}
    do
        file_list[$i]="${{file}},$(echo ${{file}} | cut -d'/' -f1)"
        ((i+=1))
    done

    ktImportText -o {output.agg} ${{file_list[@]}}
    """

shell:
    """
    exec 2> {log}

    if [[ {params.method} == "otu" ]] 
    then
        echo "Sample\tQ30 rate\tInsert size peak\tRead number\tPseudo-reads\tReads in OTU\tOTU number\tAssigned reads\t(Sub-)Species consensus\tGenus consensus\tHigher rank consensus\tNo match" > {output.report}

        Q30=$(tail -n +2 {input.fastp} | cut -d'\t' -f9)
        size=$(tail -n +2 {input.fastp} | cut -d'\t' -f11)
        reads=$(tail -n +2 {input.merging} | cut -d'\t' -f2)
        pseudo=$(tail -n +2 {input.merging} | cut -d'\t' -f5)
        clustered=$(tail -n +2 {input.clustering} | cut -d'\t' -f10)
        otu=$(tail -n +2 {input.tax} | cut -d'\t' -f2)
        assigned=$(tail -n +2 {input.compo} | awk '$2 != "No match"' | cut -d'\t' -f5 | awk '{{s+=$1}}END{{print s}}')
        spec=$(tail -n +2 {input.tax} | cut -d'\t' -f5)
        gen=$(tail -n +2 {input.tax} | cut -d'\t' -f7)
        high=$(($(tail -n +2 {input.tax} | cut -d'\t' -f9) + $(tail -n +2 {input.tax} | cut -d'\t' -f11)))
        none=$(tail -n +2 {input.tax} | cut -d'\t' -f3)

        echo "{wildcards.sample}\t$Q30\t$size\t$reads\t$pseudo\t$clustered\t$otu\t$assigned\t$spec\t$gen\t$high\t$none" >> {output.report}
    else
        echo "Sample\tQ30 rate\tInsert size peak\tRead number\tPseudo-reads\tReads in ASV\tASV number\tAssigned reads\t(Sub-)Species consensus\tGenus consensus\tHigher rank consensus\tNo match" > {output.report}

        Q30=$(tail -n +2 {input.fastp} | cut -d'\t' -f9)
        size=$(tail -n +2 {input.fastp} | cut -d'\t' -f11)
        reads=$(tail -n +2 {input.clustering} | cut -d'\t' -f2)
        pseudo=$(tail -n +2 {input.clustering} | cut -d'\t' -f6)
        clustered=$(tail -n +2 {input.clustering} | cut -d'\t' -f16)
        otu=$(tail -n +2 {input.tax} | cut -d'\t' -f2)
        assigned=$(tail -n +2 {input.compo} | awk '$2 != "No match"' | cut -d'\t' -f5 | awk '{{s+=$1}}END{{print s}}')
        spec=$(tail -n +2 {input.tax} | cut -d'\t' -f5)
        gen=$(tail -n +2 {input.tax} | cut -d'\t' -f7)
        high=$(($(tail -n +2 {input.tax} | cut -d'\t' -f9) + $(tail -n +2 {input.tax} | cut -d'\t' -f11)))
        none=$(tail -n +2 {input.tax} | cut -d'\t' -f3)

        echo "{wildcards.sample}\t$Q30\t$size\t$reads\t$pseudo\t$clustered\t$otu\t$assigned\t$spec\t$gen\t$high\t$none" >> {output.report}
    fi
    """

shell:
    """
    exec 2> {log}
    cat {input.report[0]} | head -n 1 > {output.agg}
    for i in {input.report}; do 
        cat ${{i}} | tail -n +2 >> {output.agg}
    done
    """

script:
    "../scripts/write_report.Rmd"

script:
    "../scripts/write_report.Rmd"

script:
    "../scripts/conda_collector.py"

shell:
    """
    exec 2> {log}

    echo "Database\tLast modified\tFull path" \
        > {output.report}

    paste \
        <(echo "BLAST") \
        <(date +%F -r {params.blast}.nto) \
        <(echo {params.blast}) \
        >> {output.report}

    paste \
        <(echo "taxdb.bti") \
        <(date +%F -r {params.taxdb}/taxdb.bti) \
        <(echo {params.taxdb}/taxdb.bti) \
        >> {output.report}

    paste \
        <(echo "taxdb.btd") \
        <(date +%F -r {params.taxdb}/taxdb.btd) \
        <(echo {params.taxdb}/taxdb.btd) \
        >> {output.report}

    paste \
        <(echo "taxdump lineages") \
        <(date +%F -r {params.taxdump_lin}) \
        <(echo {params.taxdump_lin}) \
        >> {output.report}

    paste \
        <(echo "taxdump nodes") \
        <(date +%F -r {params.taxdump_nodes}) \
        <(echo {params.taxdump_nodes}) \
        >> {output.report}
    """

script:
    "../scripts/primer_disambiguation.py"

shell:
    """
    exec 2> {log}

    before_r1=$(zcat {input.before_r1} | echo $((`wc -l`/4)))
    after_r1=$(zcat {input.after_r1} | echo $((`wc -l`/4)))
    before_r2=$(zcat {input.before_r2} | echo $((`wc -l`/4)))
    after_r2=$(zcat {input.after_r2} | echo $((`wc -l`/4)))

    before=$(( before_r1 + before_r2 ))
    after=$(( after_r1 + after_r2 ))

    if [ $after -ne 0 ] 
    then
        perc_discarded=$( python -c "print(f'{{round(100*(1-${{after}}/${{before}}),2)}}')" )
    else
        perc_discarded=0.00
    fi

    echo "Sample\tTotal raw reads\tTotal reads after primer trim\tNo primer found [%]" > {output.report}
    echo "{wildcards.sample}\t$before\t$after\t$perc_discarded" >> {output.report}
    """

script:
    "../scripts/parse_fastp.py"

shell:
    """
    paste {input.cutadapt} {input.fastp} 1> {output.report} 2> {log}
    """

shell:
    """
    exec 2> {log}

    cat {input.report[0]} | head -n 1 > {output.agg}
    for i in {input.report}; do 
        cat ${{i}} | tail -n +2 >> {output.agg}
    done
    """

import sys


sys.stderr = open(snakemake.log[0], "w")


def main(taxids, blocklist, output):
    with open(taxids, 'r') as fi:
        taxs = set([line.strip() for line in fi.readlines()])

    with open(blocklist, 'r') as bl:
        blocks = set([line.split('#')[0].strip() for line in bl.readlines()])

    listout = taxs.difference(blocks)

    with open(output, 'w') as fo:
        for tax in listout:
            fo.write(f"{tax}\n")


if __name__ == '__main__':
    main(snakemake.input["taxids"],
         snakemake.input["blocklist"],
         snakemake.output['mask'])

import sys


sys.stderr = open(snakemake.log[0], "w")


from os import stat
import pandas as pd


def main(report, filtered, bit_diff):
    if stat(report).st_size == 0:
        with open(filtered, "w") as fout:
            fout.write(
                "query\tsubject\tevalue\tidentity\tbitscore\tsubject_acc\t"
                "subject_taxid\talignment_length\tmismatch\tgaps\tsubject_name"
            )
    else:
        df = pd.read_csv(report, sep="\t", header=0)
        if df.empty:
            df.to_csv(filtered, sep="\t", header=True, index=False)
        else:
            sd = dict(tuple(df.groupby("query")))
            dfout = pd.DataFrame()
            for key, val in sd.items():
                dfout = pd.concat(
                    [dfout, val[val["bitscore"] >= max(val["bitscore"]) - bit_diff]]
                )
            dfout["query"] = dfout["query"].str.split(";").str[0]
            dfout.to_csv(filtered, sep="\t", header=True, index=False)


if __name__ == '__main__':
    main(snakemake.input['report'],
         snakemake.output['filtered'],
         snakemake.params['bit_diff'])

import sys


sys.stderr = open(snakemake.log[0], "w")


import taxidTools as txd
import pandas as pd


def get_lineage(taxid, tax):
    if taxid == "-":
        return ["Unassigned"]
    elif taxid == "Undetermined":
        return ["Undetermined"]
    else:
        return [node.name for node in tax.getAncestry(taxid)][::-1]
        # inverting list to have the lineage descending for Krona


def main(input, output, taxonomy):
    tax = txd.load(taxonomy)
    df = pd.read_csv(input, sep='\t', header=0)
    with open(output, "w") as out:
        for index, row in df.iterrows():
            out.write(
                "\t".join(
                    [str(row["Count"])] + get_lineage(row['Taxid'], tax)
                ) + "\n")


if __name__ == '__main__':
    main(snakemake.input['compo'],
         snakemake.output['krt'],
         snakemake.input['tax'])

import sys


sys.stderr = open(snakemake.log[0], "w")


import taxidTools as txd


def main(taxid_file, parent, output, taxonomy):

    tax = txd.load(taxonomy)

    with open(taxid_file, "r") as fin:
        db_entries = set(fin.read().splitlines()[1:])

    with open(output, "w") as fout:
        for taxid in db_entries:
            try:
                if tax.isDescendantOf(str(taxid).strip(), str(parent).strip()):
                    fout.write(taxid + "\n")
                else:
                    pass
            except KeyError:
                pass  # Ignoring missing taxids as they are either not in the
                # taxdumps or actively filtered by the user.


if __name__ == '__main__':
    main(snakemake.input['taxlist'],
         snakemake.params["taxid"],
         snakemake.output['mask'],
         snakemake.input['tax'])

import sys


sys.stderr = open(snakemake.log[0], "w")


import taxidTools as txd
from collections import Counter, defaultdict


def parse_blast(blast_file):
    """
    Parse a BLAST report and returns a dictionnary where Keys are query
    sequence names and values list of taxids for each hit.
    BLAST report must have the following formatting:
        '6 qseqid sseqid evalue pident bitscore sacc
        staxids sscinames scomnames stitle'
    """
    dictout = defaultdict()
    with open(blast_file, 'r') as fi:
        next(fi)  # Skip header
        for line in fi:
            ls = line.split()
            taxids = ls[6].split(";")  # split taxids if nescessary
            # extend taxids list for this OTU
            if ls[0] in dictout.keys():
                dictout[ls[0]].extend(taxids)
            # or inititate the list
            else:
                dictout[ls[0]] = taxids

    # Make sure everything is str formated
    dictout = {k: [str(e) for e in v] for k, v in dictout.items()}

    return dictout


def main(blast_report, output, min_consensus, taxonomy):
    if min_consensus <= 0.5 or min_consensus > 1:
        raise ValueError("'min_consensus' must be in the interval (0.5 , 1]")

    tax = txd.load(taxonomy)
    otu_dict = parse_blast(blast_report)
    with open(output, 'w') as out:
        out.write("queryID\tConsensus\tRank\tTaxid\tDisambiguation\n")

        for queryID, taxid_list in otu_dict.items():
            try:
                consensus = tax.consensus(taxid_list, min_consensus)

            except KeyError:
                # Taxid not present in the Taxdump version
                # used raises a KeyError
                # Filter out missing sequences (verbose)
                taxid_list_new = []
                for taxid in taxid_list:
                    if taxid not in tax.keys():
                        pass  # This is most likely the result of active filtering by the user
                        # No need ot be over verbose with this
                        # print(f"WARNING: taxid {taxid} missing from Taxonomy "
                        #      f"reference, it will be ignored")
                    else:
                        taxid_list_new.append(taxid)

                # Update list
                taxid_list = taxid_list_new

                # Empty list case:
                if not taxid_list:
                    consensus = "Undetermined"
                else:
                    # Get the consensus with the filtered taxids
                    consensus = tax.consensus(taxid_list, min_consensus)

            finally:
                if consensus != "Undetermined":
                    rank = consensus.rank
                    name = consensus.name
                    taxid = consensus.taxid
                else:
                    taxid = "Undetermined"
                    rank = "Undetermined"
                    name = "Undetermined"

                # (freq, name) tuple to sort
                freqs = [((v/len(taxid_list)), tax.getName(k))
                         for k, v in Counter(taxid_list).items()]
                sorted_freqs = sorted(freqs, reverse=True)

                names = "; ".join([f"{f} ({round(n,2)})"
                                   for (n, f) in sorted_freqs])
                out.write(f"{queryID}\t{name}\t{rank}\t{taxid}\t{names}\n")


if __name__ == '__main__':
    main(snakemake.input['blast'],
         snakemake.output['consensus'],
         snakemake.params["min_consensus"],
         snakemake.input['tax'])

import sys


sys.stderr = open(snakemake.log[0], "w")


import pandas as pd


def concatenate_uniq(entries):
    s = "; ".join(entries.to_list())
    df = pd.DataFrame(
        [e.rsplit(" (", 1) for e in s.split("; ")], columns=["name", "freq"]
        )  # parenthesis in names
    df.loc[:, "freq"] = df["freq"].str.replace(")", "", regex=False).astype(float)
    # Aggreagte, normalize, and sort
    tot = df["freq"].sum()
    df = df.groupby("name").apply(lambda x: x.sum() / tot)
    df = df.sort_values(by=["freq"], ascending=False)
    # Format as string
    uniq = df.to_dict()["freq"]
    uniq = [f"{name} ({round(freq, 2)})" for name, freq in uniq.items()]
    return "; ".join(uniq)


def main(compo, report, sample):
    df = pd.read_csv(compo, sep="\t", header=0).fillna(0)

    # Empty input case
    if len(df["Query"]) == 1 and df["Query"].head(1).item() == "-":
        with open(report, "w") as fout:
            fout.write(
                "Sample\tConsensus\tRank\tTaxid\tCount\tDisambiguation\tPercent of total\tPercent of assigned"
            )

    else:
        groups = df.groupby(["Consensus", "Rank", "Taxid"]).agg(
            {"Count": "sum", "Disambiguation": concatenate_uniq}
        )
        groups = groups.sort_values("Count", ascending=False).reset_index()

        # Get percs of total
        groups["perc"] = round(groups["Count"] / groups["Count"].sum() * 100, 2)

        # Get percs of assigned
        assigned, notassigned = (
            groups[groups["Consensus"] != "-"],
            groups[groups["Consensus"] == "-"],
        )
        assigned["perc_ass"] = round(assigned["Count"] / assigned["Count"].sum() * 100, 2)
        notassigned["perc_ass"] = "-"
        groups = pd.concat([assigned, notassigned])

        # Formatting
        groups.insert(0, "Sample", sample)
        groups.rename(columns={"perc": "Percent of total",
                               "perc_ass": "Percent of assigned"},
                      inplace=True)
        groups["Consensus"].replace({"-": "No match"}, inplace=True)
        groups["Taxid"].replace({0: "-"}, inplace=True)
        groups.to_csv(report, sep="\t", index=False)


if __name__ == '__main__':
    main(snakemake.input['compo'],
         snakemake.output['report'],
         snakemake.params['sample_name'])

htmltools::a(
    href="https://cvua-rrw.github.io/FooDMe/",
    htmltools::img(
        src = knitr::image_uri(snakemake@params[['logo']]), 
        alt = 'FooDMe documentation', 
        style = 'position:absolute; top:0; right:0; padding:10px;',
        width=200
    )
)

workdir <- snakemake@params[["workdir"]]

overview <- snakemake@input[["summary"]]
fastp <- snakemake@input[["fastp"]]
qc_filtering <- snakemake@input[["qc_filtering"]]
clustering <- snakemake@input[["clustering"]]
blast_rep <- snakemake@input[["blast_rep"]]
taxonomy <- snakemake@input[["taxonomy"]]
result <- snakemake@input[["result"]]
db <- snakemake@input[["db"]]
soft <- snakemake@input[["soft"]]

OTU_bool <- snakemake@params[["method"]] == "otu" # store True if using OTU 

# infer run name from workdir
run <- basename(workdir)
#head(tail(strsplit(workdir,"/")[[1]],2),1)

reportAll <- snakemake@params[["sample"]] == "all"

# Number of samples
nsamples <- nrow(read.csv(file = overview, sep = "\t", check.names=FALSE))

data_table <- read.csv(file = overview, sep = "\t", check.names=FALSE)
datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						)
						))

data_table <- read.csv(file = fastp, sep = "\t", check.names=FALSE)

# Create hyperlinks
data_table$links <- paste0("<a href=", data_table$link_to_report, ">file</a>")
data_table$link_to_report = NULL

datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						)
						))

cat("## Read filtering statistics\n")

data_table <- read.csv(file = qc_filtering, sep = "\t", check.names=FALSE)
datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						)
						))

cat("## Clustering statistics\n")

data_table <- read.csv(file = clustering, sep = "\t", check.names=FALSE)
datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						)
						))

cat("## Denoising statistics\n")

data_table <- read.csv(file = clustering, sep = "\t", check.names=FALSE)
datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						)
						))

data_table <- read.csv(file = blast_rep, sep = "\t", check.names=FALSE)
#Process links
data_table$blast_report <- paste0("<a href=", data_table$link_report, ">file</a>")
data_table$link_report = NULL
data_table$filtered_report <- paste0("<a href=", data_table$link_filtered, ">file</a>")
data_table$link_filtered = NULL

datatable(data_table, filter = 'top', rownames = FALSE, escape = FALSE,
		extensions = list("ColReorder" = NULL, "Buttons" = NULL),
		options = list(	
					dom = 'BRrltpi',
					autoWidth=FALSE,
					scrollX = TRUE,
					lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
					ColReorder = TRUE,
					buttons =
					list(
						'copy',
						'print',
						list(
						extend = 'collection',
						buttons = c('csv', 'excel', 'pdf'),
						text = 'Download'
						),
						I('colvis')
						),
						deferRender = TRUE,
						scroller = TRUE
						))