"""
Amplicon Metagenomics Workflow for NIOZ MMBL.
Version: 4.4
Author: Julia Engelmann and Alejandro Abdala
Last update: 15/02/2018
"""
run=config["RUN"]
rule all:
    input:
        expand("{PROJECT}/runs/{run}/{sample}_data/report_f_"+config["assignTaxonomy"]["method"]+".pdf" if config["pdfReport"] == "T" else "{PROJECT}/runs/{run}/{sample}_data/report_f_"+config["assignTaxonomy"]["method"]+".html", PROJECT=config["PROJECT"],sample=config["LIBRARY"], run=run)

if len(config["LIBRARY"])==1:
    rule init_structure:
        input:
            fw = config["fw_reads"],
            rv = config["rv_reads"],
            metadata = config["metadata"]
        output:
            r1="{PROJECT}/samples/{sample}/rawdata/fw.fastq",
            r2="{PROJECT}/samples/{sample}/rawdata/rv.fastq",
            metadata="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt"
        shell:
            "./init_sample.sh "+config["PROJECT"]+" "+config["LIBRARY"][0]+" {input.metadata} {input.fw} {input.rv}"
#First we run fastQC over the rawdata
rule fast_qc:
    input:
        r1="{PROJECT}/samples/{sample}/rawdata/fw.fastq",
        r2="{PROJECT}/samples/{sample}/rawdata/rv.fastq"
    output:
        o1="{PROJECT}/samples/{sample}/qc/fw_fastqc.html",
        o2="{PROJECT}/samples/{sample}/qc/rv_fastqc.html",
        s1="{PROJECT}/samples/{sample}/qc/fw_fastqc/summary.txt",
        s2="{PROJECT}/samples/{sample}/qc/rv_fastqc/summary.txt"
    benchmark:
        "{PROJECT}/samples/{sample}/qc/fq.benchmark"
    shell:
        "fastqc {input.r1} {input.r2} --extract -o {wildcards.PROJECT}/samples/{wildcards.sample}/qc/"
#validate qc if to many fails on qc report
rule validateQC:
    input:
        "{PROJECT}/samples/{sample}/qc/fw_fastqc/summary.txt",
        "{PROJECT}/samples/{sample}/qc/rv_fastqc/summary.txt",
        "{PROJECT}/samples/{sample}/qc/fw_fastqc.html",
        "{PROJECT}/samples/{sample}/qc/rv_fastqc.html",
        "{PROJECT}/samples/{sample}/rawdata/fw.fastq",
        "{PROJECT}/samples/{sample}/rawdata/rv.fastq"
    output:
        "{PROJECT}/samples/{sample}/qc/fq_fw_internal_validation.txt",
        "{PROJECT}/samples/{sample}/qc/fq_rv_internal_validation.txt"
    script:
        "Scripts/validateQC.py"
#Run pear to extend fragments
rule pear:
     input:
         r1="{PROJECT}/samples/{sample}/rawdata/fw.fastq",
         r2="{PROJECT}/samples/{sample}/rawdata/rv.fastq",
         tmp1="{PROJECT}/samples/{sample}/qc/fq_fw_internal_validation.txt",
         tmp2="{PROJECT}/samples/{sample}/qc/fq_rv_internal_validation.txt"
     output:
        "{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq",
        "{PROJECT}/runs/{run}/{sample}_data/peared/pear.log"
     benchmark:
        "{PROJECT}/runs/{run}/{sample}_data/peared/pear.benchmark"
     params:
        "{PROJECT}/runs/{run}/{sample}_data/peared/seqs"
     shell:
        "pear -f {input.r1} -r {input.r2} -o {params[0]} "
        "-t {config[pear][t]} -v {config[pear][v]} -j {config[pear][j]} -p {config[pear][p]} {config[pear][extra_params]} > "
        "{output[1]}"
#Validate % of peared reds
rule validatePear:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/peared/pear.log"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/peared/pear.log.validation"
    script:
        "Scripts/validatePear.py"

if config["fastQCPear"] == "T":
    rule fastQCPear:
      input:
          r1="{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq"
      output:
          o1="{PROJECT}/runs/{run}/{sample}_data/peared/qc/seqs.assembled_fastqc.html",
          s2="{PROJECT}/runs/{run}/{sample}_data/peared/qc/seqs.assembled_fastqc/summary.txt"
      params:
          "{PROJECT}/runs/{run}/{sample}_data/peared/qc/"
      benchmark:
          "{PROJECT}/runs/{run}/{sample}_data/peared/qc/fq.benchmark"
      shell:
          "fastqc {input.r1} --extract -o {params}"

    rule validateFastQCPear:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/peared/qc/seqs.assembled_fastqc/summary.txt",
            "{PROJECT}/runs/{run}/{sample}_data/peared/qc/seqs.assembled_fastqc.html",
            "{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/peared/qc/fq_fw_internal_validation.txt"
        script:
            "Scripts/validatePearedQC.py"
else:
    rule skipFastQCPear:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/peared/qc/fq_fw_internal_validation.txt"
        shell:
            "touch {output}"
## check if the mapping file is ok. Header line needs to start with #
rule bc_mapping_validation:
     input:
        mapp="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt"
     output:
        "{PROJECT}/metadata/bc_validation/{sample}/sampleList_mergedBarcodes_{sample}.log" # antes .log pero creo carpeta .log??
     benchmark:
        "{PROJECT}/metadata/bc_validation/{sample}/validation.benchmark"
     params:
        "{PROJECT}/metadata/bc_validation/{sample}/"
     shell:
        "validate_mapping_file.py -o {params} -m {input.mapp}"
#validate bc validation log file stop WF in failure
rule validateBCV:
        input:
            "{PROJECT}/metadata/bc_validation/{sample}/sampleList_mergedBarcodes_{sample}.log"
        output:
            "{PROJECT}/metadata/bc_validation/{sample}/validation.log"
        params:
            "{PROJECT}/metadata/bc_validation/{sample}/"
        script:
            "Scripts/validateBCV.py"

#Extract bc from reads
rule extract_barcodes:
     input:
         tmp2="{PROJECT}/runs/{run}/{sample}_data/peared/pear.log.validation",
         assembly="{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq",
         #tmpinput="{PROJECT}/runs/{run}/{sample}_data/peared/qc/seqs.assembled.html" if config["fastQCExtended"] == "T" else "{PROJECT}/runs/{run}/{sample}_data/peared/pear.log.validation"
         tmpinput="{PROJECT}/runs/{run}/{sample}_data/peared/qc/fq_fw_internal_validation.txt"
     output:
         "{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.fastq",
         "{PROJECT}/runs/{run}/{sample}_data/barcodes/reads.fastq"
     params:
         "{PROJECT}/runs/{run}/{sample}_data/barcodes/"
     benchmark:
        "{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.benchmark"
     shell:
         "extract_barcodes.py -f {input.assembly} -c {config[ext_bc][c]} "
         " {config[ext_bc][bc_length]} {config[ext_bc][extra_params]} -o {params}"
#If allow missmatch correct bar code
if config["bc_missmatch"]:
    rule correct_barcodes:
        input:
            bc="{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.fastq",
            mapp="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.fastq_corrected"
        benchmark:
            "{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes_corrected.benchmark"
        shell:
            "Rscript Scripts/errorCorrectBarcodes.R $PWD {input.mapp} {input.bc} "  + str(config["bc_missmatch"])
#split libraries - demultiplex
rule split_libraries:
     input:
         rFile="{PROJECT}/runs/{run}/{sample}_data/barcodes/reads.fastq",
         mapFile="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt",
         bcFile="{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.fastq_corrected" if config["bc_missmatch"] else "{PROJECT}/runs/{run}/{sample}_data/barcodes/barcodes.fastq",
         tmp3="{PROJECT}/metadata/bc_validation/{sample}/validation.log"
     output:
         seqs="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna", #marc as tmp
         spliLog="{PROJECT}/runs/{run}/{sample}_data/splitLibs/split_library_log.txt",
     params:
         outDir="{PROJECT}/runs/{run}/{sample}_data/splitLibs",
     benchmark:
         "{PROJECT}/runs/{run}/{sample}_data/splitLibs/splitLibs.benchmark"
     shell:
         "split_libraries_fastq.py -m {input.mapFile} -i {input.rFile} "
         "-o {params.outDir} -b {input.bcFile} -q {config[split][q]} -r {config[split][r]} "
         "--retain_unassigned_reads --barcode_type {config[split][barcode_type]} {config[split][extra_params]}"
#split libraries reverse complement
rule get_unassigned:
    input:
        split="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna",
        assembly="{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq"
    output:
        temp("{PROJECT}/runs/{run}/{sample}_data/splitLibs/unassigned.fastq")
    shell:
        "cat {input.split} | grep \"^>Unassigned\" |  sed 's/>Unassigned_[0-9]* /@/g' | "
        "sed 's/ .*//' | grep -F -w -A3  -f - {input.assembly} |  sed '/^--$/d' > {output}"
rule rc_unassigned:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs/unassigned.fastq"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs/unassigned.reversed.fastq"
    shell:
        "fastx_reverse_complement -i {input} -o {output}"
rule extract_barcodes_unassigned:
     input:
         assembly="{PROJECT}/runs/{run}/{sample}_data/splitLibs/unassigned.reversed.fastq"
     output:
         "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.fastq",
         "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/reads.fastq"
     params:
         "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/"
     benchmark:
        "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.benchmark"
     shell:
         "extract_barcodes.py -f {input.assembly} -c {config[ext_bc][c]} "
         " {config[ext_bc][bc_length]} {config[ext_bc][extra_params]} -o {params}"
#If allow missmatch correct bar code
if config["bc_missmatch"]:
    rule correct_barcodes_unassigned:
        input:
            bc="{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.fastq",
            mapp="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.fastq_corrected"
        benchmark:
            "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes_corrected.benchmark"
        shell:
            "Rscript Scripts/errorCorrectBarcodes.R $PWD {input.mapp} {input.bc} "  + str(config["bc_missmatch"])

rule remove_unassigned:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.assigned.fna"
    shell:
        "cat {input} | grep -P -A1 \"(?!>Unass)^>\" | sed '/^--$/d' > {output}"

#This rule will call a script in order to execute the librarie splitting for the
#RC sequences. It still needs the seqs.fna file bz this file contains the exact
#number of sequences assigned during the first split, then the script takes taht
#number and start to assign reads for the new splitting starting at the previous
#number...
rule split_libraries_rc:
    input:
        spplited="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna",
        rFile="{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/reads.fastq",
        mapFile="{PROJECT}/metadata/sampleList_mergedBarcodes_{sample}.txt",
        bcFile="{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.fastq_corrected" if config["bc_missmatch"]
        else "{PROJECT}/runs/{run}/{sample}_data/barcodes_unassigned/barcodes.fastq"
    output:
         seqsRC="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/seqs.fna", #marc as tmp
         spliLog="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/split_library_log.txt"
    params:
         outDirRC="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC"
    benchmark:
         "{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/splitLibs.benchmark"
    script:
         "Scripts/splitRC.py"
         #"split_libraries_fastq.py -m {input.mapFile} -i {input.rFile} "
         #"-o {params.outDirRC} -b {input.bcFile} -q {config[split][q]} "
         #"--barcode_type {config[split][barcode_type]} {config[split][extra_params]} --rev_comp_mapping_barcodes --rev_comp"
rule validateDemultiplex:
    input:
        #split="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna",
        split="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.assigned.fna",
        splitRC="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/seqs.fna",
        logSplit="{PROJECT}/runs/{run}/{sample}_data/splitLibs/split_library_log.txt",
        logSplitRC="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/split_library_log.txt",
        allreads="{PROJECT}/runs/{run}/{sample}_data/barcodes/reads.fastq"
    output:
        temp("{PROJECT}/runs/{run}/{sample}_data/splitLibs/split_library_log.txt.validation")
    params:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs",
        "{PROJECT}/runs/{run}/{sample}_data/splitLibsRC"
    script:
        "Scripts/validateSplitNew.py"

rule combine_accepted_reads:
     input:
        seqs="{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.assigned.fna", #marc as tmp
        seqsRC="{PROJECT}/runs/{run}/{sample}_data/splitLibsRC/seqs.fna", #marc as
        tmpFlow="{PROJECT}/runs/{run}/{sample}_data/splitLibs/split_library_log.txt.validation"
     output:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted.fna"
     benchmark:
        "{PROJECT}/runs/{run}/{sample}_data/combine_seqs_fw_rev.benchmark"
     shell:
        "cat {input.seqs} {input.seqsRC}  > {output}"
if config["cutAdapters"] == "T":
    rule cutadapt:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted.fna"
        output:
            out="{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_no_adapters.fna",
            log="{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_no_adapters.log"
        benchmark:
           "{PROJECT}/runs/{run}/{sample}_data/cutadapt.benchmark"
        shell:
           "cutadapt {config[cutadapt][adapters]} "
           "{config[cutadapt][extra_params]} -o {output.out} {input}  > {output.log}"

if config["chimera"]["search"] == "T":
#check for chimeric sequences
    rule search_chimera:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_no_adapters.fna"
            if config["cutAdapters"] == "T" else
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted.fna"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/chimera/chimeras.txt"
        params:
            "{PROJECT}/runs/{run}/{sample}_data/chimera"
        benchmark:
            "{PROJECT}/runs/{run}/{sample}_data/chimera/chimera.benchmark"
        shell:
            "identify_chimeric_seqs.py -m {config[chimera][method]} -i {input}  -o {params} --threads {config[chimera][threads]} {config[chimera][extra_params]}"
        #remove chimeras
    rule remove_chimera:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted.fna",
            "{PROJECT}/runs/{run}/{sample}_data/chimera/chimeras.txt"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_nc.fna",
            "{PROJECT}/runs/{run}/{sample}_data/chimera/chimera.log"
        script:
        #"filter_fasta.py -f {input[1]} -s {input[0]} -n -o {output}"
            "Scripts/remove_chimera.py"
else:
    rule rename_file:
        input:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_no_adapters.fna"
            if config["cutAdapters"] == "T" else
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted.fna"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_nc.fna"
        shell:
            "mv {input} {output}"
#Creates file with sequence length ditribution
rule histogram:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_nc.fna"
    output:
        temp("{PROJECT}/runs/{run}/{sample}_data/seqs_accepted_lengths.txt"),
        temp("{PROJECT}/runs/{run}/{sample}_data/seqs_accepted_hist.txt")
    shell:
        "cat {input} | grep -v '^>' | awk '{{print length}}' > {output[0]} "
        "&&  sort -g {output[0]} | uniq -c > {output[1]}"

rule histogram_chart:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_accepted_hist.txt",
        "{PROJECT}/runs/{run}/{sample}_data/seqs_accepted_lengths.txt"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_dist_hist.png",
        "{PROJECT}/runs/{run}/{sample}_data/seqs_statistics.txt"
    params:
        "{PROJECT}/runs/{run}/{sample}_data/"
    shell:
        "Rscript Scripts/histogram.R $PWD {input[0]} {input[1]} {params}"

#remove to long and to short sequences
rule remove_short_long_reads:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_statistics.txt",
        "{PROJECT}/runs/{run}/{sample}_data/seqs_dist_hist.png",
        "{PROJECT}/runs/{run}/{sample}_data/seqs_accepted_hist.txt",
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_nc.fna"
    params:
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.fasta",
        "{PROJECT}/runs/{run}/{sample}_data/filter.log"
    benchmark:
        temp("{PROJECT}/runs/{run}/{sample}_data/filter.benchmark")
    script:
        "Scripts/rmShortLong.py"

rule count_samples_final:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.fasta"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.dist.txt"
    shell:
        "cat {input} | grep '^>' |  cut -d'_' -f1 | sed 's/>//g' "
        "| sort | uniq -c | sort -nr | awk '{{print $1\"\\t\"$2}}' > {output}"

rule distribution_chart:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.dist.txt"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.dist.png"
    params:
        "{PROJECT}/runs/{run}/{sample}_data/"
    shell:
        "Rscript Scripts/sampleDist.R $PWD {input} {output} {config[sample_chart]}"

rule combine_filtered_samples:
        input:
            allFiltered = expand("{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.fasta",PROJECT=config["PROJECT"],sample=config["LIBRARY"], run=run)
        output:
            "{PROJECT}/runs/{run}/seqs_fw_rev_filtered.fasta"
        benchmark:
            "{PROJECT}/runs/{run}/combine_seqs_fw_rev.benchmark"
        script:
            "Scripts/combineAllReads.py"

#cluster OTUs
rule cluster_OTUs:
    input:
        "{PROJECT}/runs/{run}/seqs_fw_rev_filtered.fasta"
        #"{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.fasta"
    output:
        "{PROJECT}/runs/{run}/otu/seqs_fw_rev_filtered_otus.txt"
    params:
        trieDir="{PROJECT}/runs/{run}/otu/"
    benchmark:
        "{PROJECT}/runs/{run}/otu.benchmark"
    shell:
        "pick_otus.py -m {config[pickOTU][m]} -i {input} "
        "-o {params.trieDir}  -s {config[pickOTU][s]} {config[pickOTU][extra_params]}"
#pick representative OTUs
rule pick_representatives:
    input:
        otus="{PROJECT}/runs/{run}/otu/seqs_fw_rev_filtered_otus.txt",
        filtered="{PROJECT}/runs/{run}/seqs_fw_rev_filtered.fasta"
        #filtered="{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.fasta"
    output:
        reps="{PROJECT}/runs/{run}/otu/representative_seq_set.fasta",
        log="{PROJECT}/runs/{run}/otu/representative_seq_set.log"
    benchmark:
        "{PROJECT}/runs/{run}/pick_reps.benchmark"
    shell:
        "pick_rep_set.py -m {config[pickRep][m]} -i {input.otus} "
        "-f {input.filtered} -o {output.reps} --log_fp {output.log} {config[pickRep][extra_params]}"
#assign taxonomy to representative OTUs
rule assign_taxonomy:
    input:
        "{PROJECT}/runs/{run}/otu/representative_seq_set.fasta"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_tax_assignments.txt",
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_tax_assignments.log"
    params:
        outdir="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/assign_taxa.benchmark"
    shell:
        "parallel_assign_taxonomy_{config[assignTaxonomy][method]}.py -i {input} --id_to_taxonomy_fp {config[assignTaxonomy][mapFile]} "
        "{config[assignTaxonomy][dbType]} {config[assignTaxonomy][dbFile]} --jobs_to_start {config[assignTaxonomy][jobs]} "
        "--output_dir {params.outdir}  {config[assignTaxonomy][extra_params]}"

#The script make_otu_table.py tabulates the number of times an OTU is found in each sample,
#and adds the taxonomic predictions for each OTU in the last column if a taxonomy file is supplied
rule make_otu_table:
    input:
        tax="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_tax_assignments.txt",
        otus="{PROJECT}/runs/{run}/otu/seqs_fw_rev_filtered_otus.txt"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.biom"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.biom.benchmark"
    shell:
        "make_otu_table.py -i {input.otus} -t {input.tax} -o {output} {config[makeOtu][extra_params]}"

#filter OTU table
rule summarize_taxa:
    input:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.biom"
    output:
        "{PROJECT}/runs/{run}/otu/taxa_"+config["assignTaxonomy"]["method"]+"/otuTable_L6.txt"
    params:
        "{PROJECT}/runs/{run}/otu/taxa_"+config["assignTaxonomy"]["method"]+"/"
    benchmark:
        "{PROJECT}/runs/{run}/otu/taxa_"+config["assignTaxonomy"]["method"]+"/summarize_taxa.benchmark"
    shell:
        "summarize_taxa.py -i {input} -o {params} {config[summTaxa][extra_params]}"

#Converts otu table from biom format to tsv
rule convert_table:
    input:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.biom",
        #"{PROJECT}/runs/{run}/otu/taxa_"+config["assignTaxonomy"]["method"]+"/"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.txt"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.txt.benchmark"
    shell:
        "biom convert -i {input[0]} -o {output} {config[biom][tableType]} "
        "{config[biom][headerKey]} {config[biom][outFormat]} {config[biom][extra_params]}"

#filter OTU table
rule filter_otu:
    input:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.biom"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_noSingletons.biom"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_nosingletons.bio.benchmark"
    shell:
        "filter_otus_from_otu_table.py -i {input} -o {output} -n {config[filterOtu][n]} {config[filterOtu][extra_params]}"
#Convert to/from the BIOM table format.
rule convert_filter_otu:
    input:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_noSingletons.biom"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_noSingletons.txt"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_nosingletons.txt.benchmark"
    shell:
        "biom convert -i {input} -o {output} {config[biom][tableType]} "
        "{config[biom][headerKey]} {config[biom][outFormat]} {config[biom][extra_params]}"
#OTU map-based filtering: Keep all sequences that show up in an OTU map.
rule filter_rep_seqs:
    input:
        fastaRep="{PROJECT}/runs/{run}/otu/representative_seq_set.fasta",
        otuNoSingleton="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_noSingletons.biom"
    output:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_noSingletons.fasta"
    benchmark:
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_noSingletons.benchmark"
    shell:
        "filter_fasta.py -f {input.fastaRep} -o {output} -b {input.otuNoSingleton} {config[filterFasta][extra_params]}"
if config["alignRep"]["align"] == "T":
#Align representative sequences
    rule align_rep_seqs:
        input:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_noSingletons.fasta"
        output:
            aligned="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/representative_seq_set_noSingletons_aligned.fasta",
            log="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/representative_seq_set_noSingletons_log.txt"
        params:
            outdir="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/"
        benchmark:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/align_rep_seqs.benchmark"
        shell:
            "align_seqs.py -m {config[alignRep][m]} -i {input} -o {params.outdir} {config[alignRep][extra_params]}"

#This step should be applied to generate a useful tree when aligning against a template alignment (e.g., with PyNAST)
    rule filter_alignment:
        input:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/representative_seq_set_noSingletons_aligned.fasta"
        output:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/representative_seq_set_noSingletons_aligned_pfiltered.fasta"
        params:
            outdir="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/"
        benchmark:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/align_rep_seqs.benchmark"
        shell:
            "filter_alignment.py -i {input} -o {params.outdir} {config[filterAlignment][extra_params]}"

#Many downstream analyses require that the phylogenetic tree relating the OTUs in a study be present.
#The script make_phylogeny.py produces this tree from a multiple sequence alignment
    rule make_tree:
        input:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/representative_seq_set_noSingletons_aligned_pfiltered.fasta"
        output:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/representative_seq_set_noSingletons_aligned_pfiltered.tre"
        benchmark:
            "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/representative_seq_set_noSingletons_aligned_pfiltered.benchmark"
        shell:
            "make_phylogeny.py -i {input} -o {output} -t {config[makeTree][method]} {config[makeTree][extra_params]}"

rule report_all:
   input:
        a="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable.txt",
        b="{PROJECT}/runs/{run}/otu/taxa_"+config["assignTaxonomy"]["method"]+"/otuTable_L6.txt",
        c="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/otuTable_noSingletons.txt",
        e="{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/aligned/filtered/representative_seq_set_noSingletons_aligned_pfiltered.tre"
        if config["alignRep"]["align"] == "T" else
        "{PROJECT}/runs/{run}/otu/"+config["assignTaxonomy"]["method"]+"/representative_seq_set_noSingletons.fasta"
   output:
        temp("{PROJECT}/runs/{run}/reporttmp_all.html")
   benchmark:
        "{PROJECT}/runs/{run}/report_all.benchmark"
   script:
        "Scripts/report_all.py"

rule tune_report_all:
    input:
        "{PROJECT}/runs/{run}/reporttmp_all.html"
    output:
        "{PROJECT}/runs/{run}/report_all_"+config["assignTaxonomy"]["method"]+".html"
    script:
        "Scripts/tuneReport.py"
if config["pdfReport"] == "T":
    rule translate_to_pdf:
        input:
            "{PROJECT}/runs/{run}/report_all_"+config["assignTaxonomy"]["method"]+".html"
        output:
            "{PROJECT}/runs/{run}/report_all_"+config["assignTaxonomy"]["method"]+".pdf"
        shell:
            "wkhtmltopdf {input} {output}"
rule count_reads:
     params:
        raw=lambda wildcards: ",".join(["{PROJECT}/samples/{sample}/rawdata/fw.fastq", "{PROJECT}/samples/{sample}/rawdata/rv.fastq",
        "{PROJECT}/runs/{run}/{sample}_data/peared/seqs.assembled.fastq","{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_accepted_nc.fna",
        "{PROJECT}/runs/{run}/{sample}_data/splitLibs/seqs.fna"])
     output:
        "{PROJECT}/runs/{run}/{sample}_data/allcounts.txt"
     shell:
        "count_seqs.py -i {params.raw} -o {output}"

rule report:
    input:
        counts="{PROJECT}/runs/{run}/{sample}_data/allcounts.txt",
        report_all="{PROJECT}/runs/{run}/report_all_"+config["assignTaxonomy"]["method"]+".html",
        dist_chart="{PROJECT}/runs/{run}/{sample}_data/seqs_fw_rev_filtered.dist.png"
    output:
        temp("{PROJECT}/runs/{run}/{sample}_data/report.html")
    script:
        "Scripts/report.py"
rule tune_report:
    input:
        "{PROJECT}/runs/{run}/{sample}_data/report.html"
    output:
        "{PROJECT}/runs/{run}/{sample}_data/report_f_"+config["assignTaxonomy"]["method"]+".html"
    script:
        "Scripts/tuneReport.py"
if config["pdfReport"] == "T":
    rule translate_pdf_final_report:
        input:
            toTranslate="{PROJECT}/runs/{run}/{sample}_data/report_f_"+config["assignTaxonomy"]["method"]+".html",
            tmp="{PROJECT}/runs/{run}/report_all_"+config["assignTaxonomy"]["method"]+".pdf"
        output:
            "{PROJECT}/runs/{run}/{sample}_data/report_f_"+config["assignTaxonomy"]["method"]+".pdf"
        shell:
            "wkhtmltopdf {input.toTranslate} {output}"