Snakefile_org

SAMPLES_DIR = "/home/genomics/mhannaert/snakemake/Illuminapipeline/data/sampels"

# sample names
IDS = glob_wildcards("{id}_1.fq.gz")
CONDITIONS = ["1", "2"]


rule all:
    input:
        "results/00_fastqc/",
        "results/01_multiqc/",
        "results/02_kraken2/",
        "results/03_krona/",
        "results/04_fastp/",
        "results/05_shovill/",
        "results/06_skani/",
        "results/07_quast/",
        "results/08_busco/",


rule fastqc:
    input:
        "{id}_{conditions}.fq.gz"
    output:
        directory("results/00_fastqc/")
    params:
        extra="-t 32"
    log:
        "logs/fastqc.log"
    shell:
        """
        fastqc {params.extra} {input} --extract -o {output} 2>> {log}
        """


rule multiqc:
    input:
        "results/00_fastqc/"
    output:
        directory("results/01_multiqc/")
    log:
        "logs/multiqc.log"
    conda:
        "envs/multiqc.yml"
    shell:
        """
        multiqc {input} 2>> {log}
        rm -rf {input} 2>> {log}
        mv {input}/multiqc_report.html {output}
        """


rule Kraken2:
    input:
        "{id}_{CONDITIONS}.fq.gz"
    output:
        "results/02_kraken2/{id}_kraken2.report"
    params:
        threads=16
    log:
        "logs/Kraken2.log"
    shell:
        """
        kraken2 --gzip-compressed {input} --db /var/db/kraken2/Standard --report {output} --threads {params.threads} --quick --memory-mapping 2>> {log}
        """


rule Krona:
    input:
        "results/02_kraken2/{id}_kraken2.report"
    output:
        "results/03_krona/{id}_krona.html"
    params:
        extra="-t 5 -m 3"
    log:
        "logs/Krona.log"
    conda:
        "envs/krona.yml"
    shell:
        """
        ktImportTaxonomy {params.extra} -o {output.html} {input} 2>> {log}
        rm {input}
        """


rule fastp:
    input:
        first="{id}_1.fq.gz",
        second="{id}_2.fq.gz"
    output:
        first="results/04_fastp/{id}_1.fq.gz",
        second="results/04_fastp/{id}_2.fq.gz",
        html="results/04_fastp/{id}_fastp.html",
        json="results/04_fastp/{id}_fastp.json"
    params:
        extra="-w 32"
    log:
        "logs/fastp.log"
    shell:
        """
        fastp {params.extra} -i {input.first} -I {input.second} -o {output.first} -O {output.second} -h {output.html} -j {output.json} --detect_adapter_for_pe 2>> {log}
        """


rule shovill:
    input:
        first="results/04_fastp/{id}_1.fq.gz",
        second="results/04_fastp/{id}_2.fq.gz"
    output:
        directory("results/05_shovill/")
    params:
        extra="--cpus 16 --ram 16 --minlen 500 --trim"
    log:
        "logs/shovill.log"
    conda:
        "envs/shovill.yml"
    shell:
        """
        shovill --R1 {input.first} --R2 {input.second} {params.extra} -outdir {output}
        """


rule assemblies:
    input:
        copie="results/05_shovill/",
        remove="results/04_fastp/.fq.gz"
    output:
        "results/assemblies/{id}.fna"
    shell:
        """
        for d in $(ls -d {input.copie}*); do cp "$d"/contigs.fa {output}; done
        rm {' '.join(input.remove)}
        """


rule skANI:
    input:
        "results/assemblies/{id}.fna"
    output:
        "results/06_skani/skani_results_file.txt"
    params:
        extra="-t 24 -n 1"
    log:
        "logs/skani.log"
    conda:
        "envs/skani.yml"
    shell:
        """
        skani search {input} -d /home/genomics/bioinf_databases/skani/skani-gtdb-r214-sketch-v0.2 -o {output} {params.extra} 2>> {log}
        """


rule Quast:
    input:
        "results/assemblies/{id}.fna"
    output:
        directory("results/07_quast/")
    conda:
        "envs/quast.yml"
    shell:
        """
        for f in {input}; do quast.py "$f" -o {output}/$(basename "$f" .fna); done
        """


rule quast_summarie:
    input:
        "results/assemblies/{id}.fna"
    output:
        "results/07_quast/quast_summary_table.txt"
    params:
        header="Assembly\tcontigs (>= 0 bp)\tcontigs (>= 1000 bp)\tcontigs (>= 5000 bp)\tcontigs (>= 10000 bp)\tcontigs (>= 25000 bp)\tcontigs (>= 50000 bp)\tTotal length (>= 0 bp)\tTotal length (>= 1000 bp)\tTotal length (>= 5000 bp)\tTotal length (>= 10000 bp)\tTotal length (>= 25000 bp)\tTotal length (>= 50000 bp)\tcontigs\tLargest contig\tTotal length\tGC (%)\tN50\tN90\tauN\tL50\tL90\tN's per 100 kbp",
    shell:
        """
        # Create the output file and add the header
        echo -e "{params.header}" >> {output}

        # Initialize a counter for the number of files processed
        counter=0

        # Find and process all transposed_report.tsv files
        for file in $(find -type f -name "transposed_report.tsv"); do
            # Add the content of each file to the summary table (excluding the header)
            tail -n +2 "$file" >> {output}
            # Increment the counter
            counter=$((counter+1))
        done
    """


rule excel_and_beeswarm:
    input:
        py_excel="skani_quast_to_xlsx.py",
        beeswarm="beeswarm_vis_assemblies.R",
        directory="results/",
        quast="07_quast/quast_summary_table.txt"
    shell:
        """
        {input.py_excel} {input.directory}
        {input.beeswarm} {input.directory}{input.quast}
        """


rule Busco:
    input:
        "results/assemblies/{id}.fna"
    output:
        "results/08_busco/"
    params:
        extra="-m genome --auto-lineage-prok -c 32"
    conda:
        "envs/busco.yml"
    shell:
        """
        for sample in $(ls {input} | awk 'BEGIN{FS=".fna"}{print $1}'); do busco -i "$sample".fna -o {output} {params.extra} ; done
        """


rule Busco_visualisation:
    input:
        sample="results/08_busco/",
        script="generate_plot.py"
    output:
        directory("results/busco_summaries")
    conda:
        "envs/busco.yml"
    shell:
        """
        cp {input.sample}*/*/short_summary.specific.burkholderiales_odb10.*.txt {output}
        cd {output}
        for i in $(seq 1 15 $(ls -1 | wc -l)); do
            echo "Processing files $i to $((i+14))"
            mkdir -p part_"$i-$((i+14))"
            ls -1 | tail -n +$i | head -15 | while read file; do
                echo "Processing file: $file"
                mv "$file" part_"$i-$((i+14))"
            done
            {input.script} -wd part_"$i-$((i+14))"
        done
        """