shell:
    """
    fastqc -t {threads} {input} -o $(dirname {output[0]}) --nogroup 2> {log[0]}
    """

shell:
    """
    multiqc --interactive -f -n {output[0]} $(dirname {input} | tr "\n" " ") --cl_config "max_table_rows: 100000" 2> {log[0]}
    """

import pandas
import re
import itertools

## Define a function to extract the .fastq name pattern
def get_read_naming_patterns(directory):
    result = []
    extension= {}
    for fname in sorted(os.listdir(directory)):
        if fname.endswith("fastq.gz") or fname.endswith("fq.gz") or fname.endswith("fastq") or fname.endswith("fq"):
            regex_str = '(_L0+[1-9]+)?_(R)?(1|2)(\.|_)' #regex for finding R1 and R2, if L001 is present before, it is also included
            regex = re.compile(regex_str)
            ext = re.search(regex, fname)
            if ext is None:
                ext = re.search(r'f(?:ast)?q(?:\.gz)?', fname)
                samp = re.sub("\.$", "", re.search(r'^([^\.]*)\.*', fname).group(0))
                if samp in extension.keys():
                    if ext.group(0).endswith(".gz"):
                        extension[samp] = [ext.group(0)]
                else:
                    extension[samp] = [ext.group(0)]
            else:
                regex_after = re.compile(regex_str+".*")
                regex_before = re.compile(".*"+regex_str)
                read = re.compile(re.search(regex_after, fname).group(0))
                samp = re.sub(regex, '', re.search(regex_before, fname).group(0))
                extension.setdefault(samp, [])
                extension[samp].append(re.sub("^_", "", read.pattern))
    return(extension)


## Check for the presence of a directory for .fastq in config. If none, will be "links".
if "link_directory" in config.keys():
    link_directory = config["link_directory"]
    if not link_directory.endswith("/"):
        link_directory = link_directory + "/"
else:
    link_directory = "links/"

## Check that "tax_DB_path" ends with "/". Otherwise, correct. 
if not config["tax_DB_path"].endswith("/"):
    config["tax_DB_path"] = config["tax_DB_path"] + "/"

## Check for the existence/accessbility of tax_DB_path   
if os.path.isdir(config['tax_DB_path']) is False:
   raise IOError("Cannot access to '{}' variable.".format('tax_DB_path'))


## Check for the presence of a metadata table in in config, either for local reads ("local_samples") or distant reads ("sra_samples")
if "local_samples" not in config.keys() and "sra_samples" not in config.keys():
    raise ValueError("No samples defined in the config file. local_samples or sra_samples must be defined in config")

## Check the local_samples/sra_samples tables for the presence of following variables in config: 'sample_label', 'grouping_column', 'run_column'.
### read sample dataframe from config file
if "local_samples" in config.keys():
    metadata = config["local_samples"] 
elif "sra_samples" in config.keys():
    metadata = config["sra_samples"]    
df = pandas.read_csv(metadata, sep="\t", index_col=0)
df_colnames = set(df.columns)

### Get values for which presence should be tested in dataframe
to_check = ['sample_label', 'grouping_column', 'run_column'] # a list of config values describing columns that should be in the dataframe
to_check_dict = (dict((k, config[k]) for k in to_check))

### Set a function to compare list
def list_diff(list1, list2): 
	return (list(set(list1) - set(list2))) 

### Loop over the values of config variable to check
for key in to_check_dict.keys():

    ## Extract values from dictionnary
    to_check_values = to_check_dict[key]

    ## If it is a single value, test that it is in the columns of the metadata
    if isinstance(to_check_values, str):

        ## if the column name list is not in the column names of df dataframe then raise an error
        if to_check_values not in df_colnames:                
            message = f"{' and '.join(to_check_values)} not available in the sample TSV file."
            raise IOError(message)

    ## If it is a list, check that it contained in metadata
    elif isinstance(to_check_values, list):

        diff = []
        diff = list_diff(to_check_values, df_colnames)

        ## if the set of column names in to_check list is not in the column names of df dataframe then raise an error
        if diff:
            for item in diff:
                """if the set of column names in to_check list is not a subset of column names of df dataframe then raise an error"""
                message = f"{' and '.join(diff)} not available in the sample TSV file."
                raise IOError(message)

    else:
        message = f"{key} should be a string or a list of string"
        raise IOError(message)


## Generate a list of input files, either local (when "local_samples" is in config, or SRA-hosted (when "sra_samples" is in config))
### Create the final pandas dataframe that will be completed
all_samples=pandas.DataFrame()

### Create a set of dictionaries to store sample characteristics
sras_ext = {}
reads_sra = {}
reads_local = {}
original_names = {}
reads_ext = {}
paths  = {}
layout = {}


### In case of local samples, work our way through the local_samples table to extract read paths (if indicated in the R1/R2 columns) or extract match it with .fastq found in the "links" directory.
if "local_samples" in config.keys():
    ## Read the metadata table
    local_data = pandas.read_csv(config["local_samples"], sep="\t", index_col=0)
    local_data.index = [str(x) for x in local_data.index]
    all_local_sample_names =  "".join(list(local_data.index))
    ## Check for forbidden characters in sample names
    if "(" in all_local_sample_names or ")" in all_local_sample_names or "_-_" in all_local_sample_names:
        raise ValueError("Forbidden character in sample name in sample name file")
    ## Extract path if indicated in "R1" column
    if "R1" in list(local_data):
        for sample, sample_data in local_data.iterrows():
            if sample in paths:
                ## Check for sample names to be unique
                raise IOError("Identical sample name used multiple times: %s" % sample_name)
            paths[sample] =[sample_data.loc["R1"]] 
            reads_ext[sample] = "single"
            layout[sample] = "single"
            if 'R2' in local_data.columns.values:
                if "R1" in str(sample_data.loc["R2"]):
                    raise IOError("ATTENTION! R1 flag within R2 filename: %s", sample_data.loc["R2"])
                if (str(sample_data.loc["R2"]).strip()) != "nan":
                    paths[sample].append(sample_data.loc["R2"])
                    reads_ext[sample] = ["R1", "R2"]
                    layout[sample] = "paired"
        all_samples = local_data   
        paths = {**paths}


    ## In absence of a "R1" column indicating file paths, try to match the sample names with the content of the "links" folder
    else:

        reads_local = get_read_naming_patterns(link_directory)            
        original_names = { x : x for x in reads_local.keys() }
        read_correct = {}
        original_correct = {}

        ## In absence of a "OldSampleName", the leftmost column is used to match with "links" content
        if "OldSampleName" not in list(local_data):

            for sample in list(local_data.index):
                regex = re.compile(r'%s([^a-zA-Z0-9]|$)' % sample) # this regex ensures that the matching of the sample names end at the end of the str, to prevent S1 matching S10 for instance
                match = [bool(re.match(regex, x)) for x in sorted(list(original_names.keys()))]
                if sum(match) != 1: #there must be one and only one entry matching one sample name
                    raise ValueError("Problem matching SampleName to read file names: " +  sample)
                read_name = str(sorted(list(original_names.keys()))[match.index(True)]) # Sample name with _SX
                original_correct[sample] = original_names[read_name]
                read_correct[sample] = reads_local[read_name]
                paths[sample] = expand(link_directory + read_name + "_{reads}" ,  reads = read_correct[sample]) 

                ## In presence of a "LibraryLayout" column, samples can be specifid to be "single" or "paired". 
                if "LibraryLayout" in list(local_data):
                    if local_data.loc[sample, "LibraryLayout"].lower()=="paired":
                        reads_ext[sample]=["R1", "R2"]
                        layout[sample] = "paired"
                    elif local_data.loc[sample, "LibraryLayout"].lower()=="single":
                        reads_ext[sample]=["single"]
                        layout[sample] = "single"
                    else:
                        raise ValueError("Problem in the Local_sample file, LibraryLayout badly defined")
                ## Otherwise, "paired" is assumed
                else:
                    reads_ext[sample]=["R1", "R2"]
                    layout[sample] = "paired"

        ## In presence of a "OldSampleName", the this column is used to match with "links" content
        else:
            for Old in list(local_data["OldSampleName"]):
                regex = re.compile(r'%s([^a-zA-Z0-9]|$)' % Old)
                match = [bool(re.match(regex, x)) for x in sorted(list(original_names.keys()))]
                if sum(match) != 1:
                    raise ValueError("Problem matching OldSampleName to read file names : " + Old)
                read_name=str(sorted(list(original_names.keys()))[match.index(True)]) # Sample with have SX to unify with above
                sample = str(local_data.index[local_data['OldSampleName'] == Old][0])
                original_correct[sample] = original_names[read_name]
                read_correct[sample] = reads_local[read_name]
                paths[sample] = expand(link_directory + read_name + "_{reads}" ,  reads = read_correct[sample])

                ## In presence of a "LibraryLayout" column, samples can be specifid to be "single" or "paired". 
                if "LibraryLayout" in list(local_data):
                    if local_data.loc[sample, "LibraryLayout"].lower()=="paired":
                        reads_ext[sample]=["R1", "R2"]
                        layout[sample] = "paired"
                    elif local_data.loc[sample, "LibraryLayout"].lower()=="single":
                        reads_ext[sample]=["single"]
                        layout[sample] = "single"
                    else:
                        raise ValueError("Problem in the Local_sample file, LibraryLayout badly defined")
                ## Otherwise, "paired" is assumed
                else:
                    reads_ext[sample]=["R1", "R2"]
                    layout[sample] = "paired"

        original_names = original_correct
        reads_local = read_correct
        reads_ext = reads_ext


        all_samples=local_data



### In case of sra-hosted samples, work our way through the "sra_table" to extract LibraryLayout and clean sample names.
if "sra_samples" in config.keys():
    sra_data = pandas.read_csv(config["sra_samples"], sep="\t", index_col=0).drop_duplicates()
    all_sra_sample_names = "".join(list(sra_data.index))
    if "(" in all_sra_sample_names or ")" in all_sra_sample_names or "_-_" in all_sra_sample_names:
        raise ValueError("Forbidden character in sample name in sra file")
    for sra_sample in list(sra_data.index):
        sample_name = str(sra_sample).replace(" ", "_").replace("&", "and").replace(":", "-")
        if sample_name in reads_sra.keys(): # if the sample name is already used, add _(n+1) at the end
            sample_name = sample_name+"_"+str(list(reads_sra.keys()).count(sample_name))
        reads_sra[sample_name]=str(sra_sample)
        if sra_data.loc[sra_sample, "LibraryLayout"].lower()=="paired":
            sras_ext[sample_name]=["1.fastq.gz", "2.fastq.gz"]
            reads_ext[sample_name]=["R1", "R2"]
        elif sra_data.loc[sra_sample, "LibraryLayout"].lower()=="single":
            sras_ext[sample_name] = ["fastq.gz"]
            reads_ext[sample_name]=["single"]
        else:
            raise ValueError("Problem in the sra file, LibraryLayout badly defined")
    all_samples=sra_data
    config["local_samples"] =  config["sra_samples"]

        # all_samples.loc[sample_name, "Replicate"]=sra_data.loc[i, "Replicate"]
read_naming = {**reads_local, **sras_ext}
original_names = {**original_names, **reads_sra}

def get_grouping_key(column_of_interest):
    file_list = []
    ## Convert to list if provided as a string in config
    if isinstance(column_of_interest,str):
        column_of_interest=[column_of_interest]
    unique_col=list(set(column_of_interest))
    for i in unique_col:
        combined_values = expand("{column_of_interest}/{column_values}", column_of_interest = i, column_values = list(set(all_samples[i])))
        file_list.extend(combined_values)
    return(file_list)

### Define rarefaction levels (values in config + no rarefaction)
def get_rarefaction_key(rarefaction_values):
    file_list = []
    for i in set(rarefaction_values):
        combined_values = expand("rarefaction_{rarefaction_values}", rarefaction_values = i)
        file_list = file_list + combined_values
    file_list = file_list + ["norarefaction"]

    return(file_list)


# Transform the collapse levels from the plotting taxonomic rank
def get_taxa_collapse_level_key(collapse_level):

    file_list = []
    for i in set(collapse_level):
        if i == "OTU":
            value = 'no_collapse'
        elif i == "Species" :
            value = 'collap_7'
        elif i == "Genus" :
            value = 'collap_6'
        elif i == "Family" :
            value = 'collap_5'
        elif i == "Order" :
            value = 'collap_4'
        elif i == "Class" :
            value = 'collap_3'
        elif i == "Phylum" :
            value = 'collap_2'
        elif i == "Kingdom" :
            value = 'collap_1'
        else :
            raise ValueError("Forbidden value in taxa level for barplots")
        file_list.append(value)
    file_list = file_list + ["no_collapse"]
    return(file_list)



## Set of function to generate list of output from config ##############################################################

### Light output, including QC, count table, consensus sequences and taxonomic assignement
def light_output_list():
    output = []
    output = MultiQC
    output.append(light_output)
    if "DADA2" in config["denoiser"]:
        output.append("DADA2/2_denoised/DADA2_denoising_stats.tsv")
    return(output)

### Basic output, diagnostic plots, KRONA plots and rarefaction curves
def basic_plots_list():
    output = []
    output = light_output_list()
    output.append(basic_plots)
    return(output)

### Complete set of phyloseq, in option including transposed count table and metadata (wide to long)
def phyloseq_output_list():
    output = []
    output = basic_plots_list()
    output.append(phyloseq)
    if config["melted_phyloseq"] == True:
        output.append(phyloseq_melted)
    if config["transposed_tables"] == True:
        output.append(transposed_output)
    return(output)


### Qiime2 outputs, offering interactive visualization of the data
def Qiime2_output_list():
    output = []
    output = basic_plots_list()
    if config["Qiime2_basic_output_visualization"] == True:
        output.append(Qiime2_vis_qzv)
    return(output)

### PICRUSt2 outputs, to predict genomic content based on taxonomic profiles
def PICRUSt2_list():
    output = []
    output = basic_plots_list()
    output.append(picrust2)
    return(output)

### Rule all, the default rule including all default output (not PICRUSt2, since long to compute)
def rule_all_list():
    output = []
    output.append(basic_plots_list())
    output.append(Qiime2_vis_qzv)
    output.append(phyloseq_output_list())
    return(output)

include: "make_output_fcts.py"

### Basic output
#### MuliQC
MultiQC = expand("QC/{RUN}_multiqc_raw_reads_report.html", RUN = set(all_samples[config["run_column"]]))
MultiQC.append("QC/multiqc_raw_reads_report.html")


#### Light pipeline output
light_output = expand("{denoiser}/2_denoised/dna-sequences.fasta", denoiser = config["denoiser"])
light_output.append(expand("{denoiser}/2_denoised/count_table.tsv", denoiser = config["denoiser"]))
light_output.append(expand("{denoiser}/3_classified/{classifier}_{tax_DB}/dna-sequences_tax_assignments.qzv", classifier = config["classifier"], denoiser = config["denoiser"], tax_DB = config["tax_DB_name"]))

#### DADA2 stat table
DADA2_stats_tables = "DADA2/2_denoised/DADA2_denoising_stats.tsv"


### Basic evaluation plots
basic_plots = expand("{denoiser}/5_visualization/{classifier}_{tax_DB}/reads/reads_plot_with_filtered.png",
                                 classifier = config["classifier"],
                                 denoiser = config["denoiser"],
                                 tax_DB = config["tax_DB_name"])



basic_plots.append(expand("{denoiser}/5_visualization/{classifier}_{tax_DB}/reads/reads_plot_with_{filter_lineage}_in_{filter_tax_rank}_filtered.png",
                                denoiser = config["denoiser"],
                                classifier = config["classifier"],
                                tax_DB = config["tax_DB_name"],
                                filter_tax_rank = config["filter_tax_rank"],
                                filter_lineage = config["filter_lineage"]))


basic_plots.append(expand("{denoiser}/5_visualization/{classifier}_{tax_DB}/rarefaction_curve/nofiltering/{rarefaction_value}_rarefaction_curve_{grouping_column}.png",
                                 denoiser = config["denoiser"],
                                 classifier = config["classifier"],
                                 tax_DB = config["tax_DB_name"],
                                 rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                                 grouping_column = config["grouping_column"][0]))


basic_plots.append(expand("{denoiser}/5_visualization/{classifier}_{tax_DB}/rarefaction_curve/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}_rarefaction_curve_{grouping_column}.png",
                                 denoiser = config["denoiser"],
                                 classifier = config["classifier"],
                                 tax_DB = config["tax_DB_name"],
                                 rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                                 filter_tax_rank = config["filter_tax_rank"],
                                 filter_lineage = config["filter_lineage"],
                                 grouping_column = config["grouping_column"][0]))



basic_plots.append(expand("{denoiser}/5_visualization/{classifier}_{tax_DB}/KRONA/{grouping_key}.html",
                                 denoiser = config["denoiser"],
                                 classifier = config["classifier"],
                                 tax_DB = config["tax_DB_name"],
                                 grouping_key = get_grouping_key(config["grouping_column"])))



### Advances Phyloseq set of output
#### Not filtered
phyloseq = expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]))


phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"])))

phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   norm_value = config["normalization"], 
                   abund_value = config["min_abundance"], 
                   prev_value = config["min_prevalence"]))

#### Tax filtered
phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"]))


phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"]))

phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}.rds",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"],
                   norm_value = config["normalization"], 
                   abund_value = config["min_abundance"], 
                   prev_value = config["min_prevalence"]))


### Potentially spurious ASVs
phyloseq.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/norarefaction/no_collapse/base_export/tree_treeshrink/output.txt",
            denoiser = config["denoiser"],
            classifier = config["classifier"],
            tax_DB = config["tax_DB_name"],
            collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
            rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
            filter_tax_rank = config["filter_tax_rank"],
            filter_lineage = config["filter_lineage"],
            norm_value = config["normalization"], 
            abund_value = config["min_abundance"], 
            prev_value = config["min_prevalence"]))


#### Melted Phyloseq set of output
#### Not filtered
phyloseq_melted = expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]))


phyloseq_melted.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"])))

phyloseq_melted.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   norm_value = config["normalization"], 
                   abund_value = config["min_abundance"], 
                   prev_value = config["min_prevalence"]))

#### Tax filtered
phyloseq_melted.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"]))


phyloseq_melted.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"]))

phyloseq_melted.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}_melted.tsv",
                   denoiser = config["denoiser"],
                   classifier = config["classifier"],
                   tax_DB = config["tax_DB_name"],
                   collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
                   rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
                   filter_tax_rank = config["filter_tax_rank"],
                   filter_lineage = config["filter_lineage"],
                   norm_value = config["normalization"], 
                   abund_value = config["min_abundance"], 
                   prev_value = config["min_prevalence"]))


#### Transposed output tables
#### Count-table transposed format - filtered
transposed_output = expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}_export/count_table_transposed.txt",
            denoiser = config["denoiser"],
            classifier = config["classifier"],
            tax_DB = config["tax_DB_name"],
            collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
            rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
            filter_tax_rank = config["filter_tax_rank"],
            filter_lineage = config["filter_lineage"],
            norm_value = config["normalization"], 
            abund_value = config["min_abundance"], 
            prev_value = config["min_prevalence"])


#### Count-table transposed format - not filtered
transposed_output.append(expand("{denoiser}/4_physeq/{classifier}_{tax_DB}/nofiltering/{rarefaction_value}/{collapse_key}/base_with_tree_norm{norm_value}_abund{abund_value}_prev{prev_value}_export/count_table_transposed.txt",
            denoiser = config["denoiser"],
            classifier = config["classifier"],
            tax_DB = config["tax_DB_name"],
            collapse_key = get_taxa_collapse_level_key(config["phyloseq_tax_ranks"]),
            rarefaction_value = get_rarefaction_key(config["rarefaction_value"]),
            filter_tax_rank = config["filter_tax_rank"],
            filter_lineage = config["filter_lineage"],
            norm_value = config["normalization"], 
            abund_value = config["min_abundance"], 
            prev_value = config["min_prevalence"]))


### Qiime2 output
#### Qiime2 interactive visualization
Qiime2_vis_qzv = expand("{denoiser}/2_denoised/dna-sequences.fasta",
                        denoiser = config["denoiser"])
Qiime2_vis_qzv.append(expand("{denoiser}/2_denoised/count-table.qzv",
                             denoiser = config["denoiser"]))
Qiime2_vis_qzv.append(expand("{denoiser}/2_denoised/rep-seqs.qzv",
                             denoiser = config["denoiser"]))
Qiime2_vis_qzv.append(expand("{denoiser}/3_classified/{classifier}_{tax_DB}/dna-sequences_tax_assignments.qzv",
                      denoiser = config["denoiser"],
                      classifier = config["classifier"],
                      tax_DB = config["tax_DB_name"]))


### Picrust2
picrust2 = expand("{denoiser}/6_picrust2/{classifier}_{tax_DB}/{filter_lineage}_in_{filter_tax_rank}/{rarefaction_value}/picrust/",
                    denoiser = config["denoiser"],
                    classifier = config["classifier"],
                    tax_DB = config["tax_DB_name"],
                    filter_tax_rank = config["filter_tax_rank"],
                    filter_lineage = config["filter_lineage"],
                    rarefaction_value = get_rarefaction_key(config["rarefaction_value"])
                    )

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
from snakemake import shell

if snakemake.params['amplicon_type'] == "ITS":
    print("ITS Trimming")
    forward_primer_compl = Seq.reverse_complement(Seq(snakemake.params['forward_primer'], IUPAC.ambiguous_dna))
    reverse_primer_compl = Seq.reverse_complement(Seq(snakemake.params['reverse_primer'], IUPAC.ambiguous_dna))
    shell("""cutadapt \
    --cores {snakemake.threads} \
    --error-rate 0.1 \
    --times 2 \
    --overlap 3 \
    -o {snakemake.output[R1_trimmed_reads]} \
    -p {snakemake.output[R2_trimmed_reads]} \
    -g '{snakemake.params[forward_primer]}' \
    -a '{reverse_primer_compl}' \
    -G '{snakemake.params[reverse_primer]}' \
    -A '{forward_primer_compl}' \
    --match-read-wildcards \
    --discard-untrimmed \
    {snakemake.input[R1_raw_reads]} \
    {snakemake.input[R2_raw_reads]} >> {snakemake.log[0]}""")

elif snakemake.params['amplicon_type'] == "16S":
    print("16S Trimming")
    shell("""cutadapt \
    --cores {snakemake.threads} \
    --error-rate 0.1 \
    --times 1 \
    --overlap 3 \
    -o {snakemake.output[R1_trimmed_reads]} \
    -p {snakemake.output[R2_trimmed_reads]} \
    -g '{snakemake.params[forward_primer]}' \
    -G '{snakemake.params[reverse_primer]}' \
    --match-read-wildcards \
    --discard-untrimmed \
    {snakemake.input[R1_raw_reads]} \
    {snakemake.input[R2_raw_reads]} >> {snakemake.log[0]}""")

    log <- file(snakemake@log[[1]], open="wt")
    sink(log)
    sink(log, type="message")

## Input
    fnFs <- snakemake@input[[1]]
    fnRs <- snakemake@input[[2]]

## Output
    q_score_filtered_F <- snakemake@output[["q_score_filtered_F"]]
    q_score_filtered_R <- snakemake@output[["q_score_filtered_R"]]
    q_filtering_stats_path <- snakemake@output[["q_filtering_stats"]]

## Parameters
    F_length <- snakemake@params[["F_reads_length_trim"]]
    R_length <- snakemake@params[["R_reads_length_trim"]]
    F_extected_error <- snakemake@params[["F_reads_extected_error"]]
    R_extected_error <- snakemake@params[["R_reads_extected_error"]]
    sample_name <- snakemake@params[["sample_name"]]


## Load needed libraries
    library(dada2);packageVersion("dada2")

## Filter and trim
### Reads are filtered based on the number of errors expected for the read (integration of the qscore and the length of the read). All reads with uncalled nucleotide (N) are removed too. Remaining phiX reads will be removed too. Finally, reads are cut at a length set in config.

### Filter and trim. The filtered reads are directly written while the filtering stats are being save for later compilation.
    filtering_stats <- filterAndTrim(fwd = fnFs, filt = q_score_filtered_F, rev = fnRs, filt.rev = q_score_filtered_R, truncLen=c(F_length,R_length), maxN=0, maxEE=c(F_extected_error,R_extected_error), truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=snakemake@threads, verbose = TRUE)
    filtering_stats <- as.data.frame(filtering_stats)
    filtering_stats$Sample <- sample_name

### Save the stats for this sample in a R object
    saveRDS(filtering_stats, file = q_filtering_stats_path)