Snakefile

from os import path,getcwd

# Run the pipeline line so on your local computer
# snakemake -npr --cores 10 --keep-going --rerun-incomplete --restart-times 3

configfile: "config/config.yaml"


RULES=["Download_silva_database", "Extract_primer_silva_reads", "Train_silva_classifier",
       "Download_unite_database", "Unzip_unite_DB", "modify_and_rename_unite_files","Import_unite_sequences",
       "Import_unite_taxonomy", "Import_unite_taxonomy", "Train_unite_classifier"]



rule all:
    input:
        "logs/Download_silva_database/",
        "logs/Download_unite_database/",
        "logs/modify_and_rename_unite_files/",
        "logs/Train_unite_classifier/",
        "logs/Train_silva_classifier",
        "databases/silva-138-99-nb-classifier.qza",
        "databases/silva-138-99-nb-341-926-classifier.qza",
        "databases/unite-classifier.qza"


# This rule will make rule specific log directories
# # in order to easily store the standard input and stand error
# # generated when submiting jobs to the cluster
rule make_logs_directories:
    output:
        directory("logs/Download_silva_database/"),
        directory("logs/Download_unite_database/"),
        directory("logs/modify_and_rename_unite_files/"),
        directory("logs/Train_unite_classifier/"),
        directory("logs/Train_silva_classifier")
    threads: 1
    shell:
        """
         [ -d logs/ ] || mkdir -p logs/
         cd logs/
         for RULE in {RULES}; do
          [ -d ${{RULE}}/ ] || mkdir -p ${{RULE}}/
         done
        """


# -------------- Get SILVA database for QIIME2 ------------------#

rule Download_silva_database:
    input:
        log_dirs=rules.make_logs_directories.output
    output:
        classifier="databases/silva-138-99-nb-classifier.qza",
        sequences="databases/silva-138-99-seqs.qza",
        taxonomy="databases/silva-138-99-tax.qza"
    threads: 1
    log: "logs/Download_silva_database/Download_silva_database.log"
    params:
        classifier=config["SILVA_CLASSIFIER"],
        sequences=config["SILVA_SEQUENCES"],
        taxonomy=config["SILVA_TAXONOMY"]
    shell:
        """

        # Classifier
        wget -O {output.classifier}  {params.classifier} > {log} 2>&1
        
        # Sequences
        wget -O {output.sequences} {params.sequences} > {log} 2>&1

        # Taxonomy
        wget -O {output.taxonomy} {params.taxonomy} > {log} 2>&1

        """

# ----- A primer specific silva database ---------#

rule Extract_primer_silva_reads:
    input: rules.Download_silva_database.output.sequences # OR  rules.download_silva_database.output.classifier
    output: "databases/ref-seqs-341-926.qza" 
    threads: 1
    log: "logs/Extract_primer_silva_reads/Extract_primer_silva_reads.log"
    params:
        conda_activate=config["QIIME2_ENV"],
        forward_primer=config["parameters"]["extract_sequence"]["forward_primer"],
        reverse_primer=config["parameters"]["extract_sequence"]["reverse_primer"],
        trunc_len=config["parameters"]["extract_sequence"]["trunc_length"],
        min_length=config["parameters"]["extract_sequence"]["min_length"],
        max_length=config["parameters"]["extract_sequence"]["max_length"]
    shell:
        """
        set +u

        {params.conda_activate}

        set -u
        
        qiime feature-classifier extract-reads \
             --i-sequences {input} \
             --p-f-primer {params.forward_primer} \
             --p-r-primer {params.reverse_primer} \
             --p-trunc-len {params.trunc_len} \
             --p-min-length {params.min_length} \
             --p-max-length {params.max_length} \
             --o-reads {output}  > {log} 2>&1
        """


#   Train the classifier
#  Will fail with the error message "Killed" if you are out
# of memory. It is a good idea to provide this step with
# a lot of memory say 50GB and above.
rule Train_silva_classifier:
    input: 
        sequences=rules.Extract_primer_silva_reads.output,
        taxonomy=rules.Download_silva_database.output.taxonomy
    output: "databases/silva-138-99-nb-341-926-classifier.qza"
    threads: 1
    log: "logs/Train_silva_classifier/Train_silva_classifier.log"
    params:
        conda_activate=config["QIIME2_ENV"]
    resources:
        mem_mb=50024 # 50GB of memeory required
    shell:
        """
        set +u

        {params.conda_activate}

        set -u
        
        qiime feature-classifier fit-classifier-naive-bayes \
              --i-reference-reads {input.sequences} \
              --i-reference-taxonomy {input.taxonomy} \
              --o-classifier {output}  > {log} 2>&1
             
        """


# --------------- Create Unite database for QIIME2 -----------------#

rule Download_unite_database:
    input:
        log_dirs=rules.make_logs_directories.output
    output:
        temp("databases/unite.gz")
    threads: 1
    log: "logs/Download_unite_database/Download_unite_database.log"
    params:
        url=config["UNITE_URL"]
    shell:
        "wget -O {output} {params.url} > {log} 2>&1 " 

rule Unzip_unite_DB:
    input: rules.Download_unite_database.output
    output: 
        sequences="databases/sh_qiime_release_04.02.2020/sh_refs_qiime_ver8_dynamic_04.02.2020.fasta",
        taxonomy="databases/sh_qiime_release_04.02.2020/sh_taxonomy_qiime_ver8_dynamic_04.02.2020.txt"
    threads: 1
    #log: "logs/Unzip_unite_DB/Unzip_unite_DB.log"
    params:
        out_dir=lambda w, input: path.dirname(input[0]),
        basename=lambda w, input: path.basename(input[0])
    shell:
        """
        cd {params.out_dir}
        #[ -f {log} ] || touch {log}
        tar -xvzf {params.basename}
        """

# modify the taxonomy header such that the first two lines 
# are "Feature ID\tTaxon"
# and rename the sequence.fasta file
rule modify_and_rename_unite_files:
    input:
        sequences=rules.Unzip_unite_DB.output.sequences,
        taxonomy=rules.Unzip_unite_DB.output.taxonomy
    output:
        sequences="databases/unite-sequences.fasta",
        taxonomy="databases/unite-taxonomy.txt"
    threads: 2
    log: "logs/modify_and_rename_unite_files/modify_and_rename_unite_files.log"
    shell:
        """
        # Modify and rename the unite taxonomy file
         (echo -e "Feature ID\tTaxon"; cat {input.taxonomy}) > {output.taxonomy} 2> {log}

        # Copy and rename the unite sequences file
        cat {input.sequences} > {output.sequences} 2> {log}
        """


# Setting up the already trimmed database
rule Import_unite_sequences:
    input: rules.modify_and_rename_unite_files.output.sequences
    output: "databases/unite-sequences.qza"
    threads: 2
    log: "logs/Import_unite_sequences/Import_unite_sequences.log"
    params:
        conda_activate=config["QIIME2_ENV"]
    shell:
        """
        set +u
 
        {params.conda_activate}
 
        set -u

        qiime tools import \
             --type 'FeatureData[Sequence]' \
             --input-path {input} \
             --output-path {output} > {log} 2>&1
        """

        
        
# Import Taxonomy
rule Import_unite_taxonomy:
    input: rules.modify_and_rename_unite_files.output.taxonomy
    output: "databases/unite-taxonomy.qza"
    threads: 2
    log: "logs/Import_unite_taxonomy/Import_unite_taxonomy.log"
    params:
        conda_activate=config["QIIME2_ENV"]
    shell:
        """
        set +u

        {params.conda_activate}

        set -u

        qiime tools import \
             --type 'FeatureData[Taxonomy]' \
             --input-path {input} \
             --output-path {output} > {log} 2>&1
        """


# Train the classifier
rule Train_unite_classifier:
    input: 
        sequences=rules.Import_unite_sequences.output,
        taxonomy=rules.Import_unite_taxonomy.output
    output: "databases/unite-classifier.qza"
    threads: 10
    log: "logs/Train_unite_classifier/Train_unite_classifier.log"
    params:
        conda_activate=config["QIIME2_ENV"]
    shell:
        """
        set +u

        {params.conda_activate}

        set -u
        
        qiime feature-classifier fit-classifier-naive-bayes \
              --i-reference-reads {input.sequences} \
              --i-reference-taxonomy {input.taxonomy} \
              --o-classifier {output} > {log} 2>&1
        """