Guppy

rule minimap2:
    input:
        fa="data/ecoli_k12_mg1655.fasta",
        fq="guppy_results/{sample}.fq"
    output:
        "guppy_results/mapped/{sample}.bam"
    shell:
        "minimap2 -a -x map-ont {input.fa} {input.fq} | samtools sort -T tmp -o {output}"

rule samtools_index_and_extract_methylation_from_fast5:
    input:
        bam="guppy_results/mapped/{sample}.bam",
        script="guppy_results/gcf52ref/scripts/extract_methylation_fast5.py"
    output:
        "guppy_results/mapped/{sample}.bam.bai"
    shell:
        "samtools index {input.bam} && python {input.script} -p 10 guppy_results/workspace/*.fast5"

rule extract_methylation_from_rocks:
    input:
        #rocks="base_mods.rocksdb",
        script="guppy_results/gcf52ref/scripts/extract_methylation_from_rocks.py",
        bam="guppy_results/mapped/{sample}.bam",
        bai="guppy_results/mapped/{sample}.bam.bai",
        fa="data/ecoli_k12_mg1655.fasta"
    output:
        "guppy_results/{sample}_guppy-log-perCG.tsv"
    shell:
        """
        set +e
        python {input.script} -d base_mods.rocksdb/ -a {input.bam} -r {input.fa} -o {output}
        exitcode=$?
        if [ $exitcode -eq 1 ]
        then
            exit 0
        else
            exit 0
        fi
        """

rule calculate_frequency:
    input:
        script="script_in_snakemake/run_guppy.R",
        log="guppy_results/{sample}_guppy-log-perCG.tsv"
    output:
      file1="guppy_results/{sample}_guppy-freq-perCG.tsv",
      file2="guppy_results/{sample}_guppy-freq-perCG-combStrand.tsv"
    shell:
        "Rscript {input.script} {input.log} {output.file1} {output.file2}"

rule create_input_for_comb_model:
    input:
        log="guppy_results/{sample}_guppy-log-perCG.tsv",
        script="script_in_snakemake/format_guppy.R"
    output:
        "guppy_results/{sample}_guppy-perRead-score.tsv"
    shell:
        "Rscript {input.script} {input.log} {output}"