Snakefile

rule alfa_qc:
    '''
        Run ALFA from stranded bedgraph files
    '''
    input:
        plus = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "{sample}.{unique}.plus.bg"),
        minus = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "{sample}.{unique}.minus.bg"),
        gtf = lambda wildcards:
            os.path.join(
                config["alfa_indexes"],
                get_sample(
                    'organism',
                    search_id='index',
                    search_value=wildcards.sample),
                get_sample(
                    'index_size',
                    search_id='index',
                    search_value=wildcards.sample),
                "ALFA",
                "sorted_genes.stranded.ALFA_index")

    output:
        biotypes = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "ALFA_plots.Biotypes.pdf"),
        categories = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "ALFA_plots.Categories.pdf"),
        table = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "{sample}.ALFA_feature_counts.tsv")

    params:
        out_dir = lambda wildcards, output:
            os.path.dirname(output.biotypes),
        plus = lambda wildcards, input:
            os.path.basename(input.plus),
        minus = lambda wildcards, input:
            os.path.basename(input.minus),
        alfa_orientation = lambda wildcards:
            get_sample(
                'alfa_directionality',
                search_id='index',
                search_value=wildcards.sample),
        genome_index = lambda wildcards, input:
            os.path.abspath(
                os.path.join(
                    os.path.dirname(input.gtf),
                    "sorted_genes")),
        name = "{sample}"

    singularity:
        "docker://zavolab/alfa:1.1.1-slim"

    log:
        os.path.join(
            config["log_dir"],
            "samples",
            "{sample}",
            "alfa_qc.{unique}.log")

    shell:
        "(cd {params.out_dir}; \
        alfa \
        -g {params.genome_index} \
        --bedgraph {params.plus} {params.minus} {params.name} \
        -s {params.alfa_orientation}) &> {log}"


# cd {params.out_dir};
rule alfa_qc_all_samples:
    '''
        Run ALFA from stranded bedgraph files on all samples
    '''
    input:
        tables = lambda wildcards:
            expand(
                os.path.join(
                    config["output_dir"],
                    "samples",
                    "{sample}",
                    "ALFA",
                    "{unique}",
                    "{sample}.ALFA_feature_counts.tsv"),
                sample=pd.unique(samples_table.index.values),
                unique=wildcards.unique)
    output:
        biotypes = os.path.join(
            config["output_dir"],
            "ALFA",
            "{unique}",
            "ALFA_plots.Biotypes.pdf"),
        categories = os.path.join(
            config["output_dir"],
            "ALFA",
            "{unique}",
            "ALFA_plots.Categories.pdf")

    params:
        out_dir = lambda wildcards, output:
            os.path.dirname(output.biotypes)

    log:
        os.path.join(
            config["log_dir"],
            "alfa_qc_all_samples.{unique}.log")

    singularity:
        "docker://zavolab/alfa:1.1.1-slim"

    shell:
        "(alfa -c {input.tables} -o {params.out_dir}) &> {log}"


rule alfa_concat_results:
    input:
        expand(
            os.path.join(
                config["output_dir"],
                "ALFA",
                "{unique}",
                "ALFA_plots.{annotation}.pdf"),
            unique=["Unique", "UniqueMultiple"],
            annotation=["Categories", "Biotypes"])

    output:
        os.path.join(
            config["output_dir"],
            "ALFA",
            "ALFA_plots_mqc.png")

    params:
        density = 300

    log:
        os.path.join(
            config["log_dir"],
            "alfa_qc_all_samples.concat.log")

    singularity:
        "docker://zavolab/imagemagick:7.0.8"

    shell:
        "(convert -append -density {params.density} \
            {input} {output}) &> {log}"


rule prepare_multiqc_config:
    '''
        Prepare config for the MultiQC
    '''
    input:
        script = os.path.join(
            workflow.basedir,
            "workflow",
            "scripts",
            "zarp_multiqc_config.py")

    output:
        multiqc_config = os.path.join(
            config["output_dir"],
            "multiqc_config.yaml")

    params:
        logo_path = config['report_logo'],
        multiqc_intro_text = config['report_description'],
        url = config['report_url']

    log:
        stderr = os.path.join(
            config["log_dir"],
            "prepare_multiqc_config.stderr.log"),
        stdout = os.path.join(
            config["log_dir"],
            "prepare_multiqc_config.stdout.log")

    shell:
        "(python {input.script} \
        --config {output.multiqc_config} \
        --intro-text '{params.multiqc_intro_text}' \
        --custom-logo {params.logo_path} \
        --url '{params.url}') \
        1> {log.stdout} 2> {log.stderr}"


rule multiqc_report:
    '''
        Create report with MultiQC
    '''
    input:
        fastqc_se = expand(
            os.path.join(
                config['output_dir'],
                "samples",
                "{sample}",
                "fastqc",
                "{mate}"),
            sample=pd.unique(samples_table.index.values),
            mate="fq1"),

        fastqc_pe = expand(
            os.path.join(
                config['output_dir'],
                "samples",
                "{sample}",
                "fastqc",
                "{mate}"),
            sample=[i for i in pd.unique(
                samples_table[samples_table['seqmode'] == 'pe'].index.values)],
            mate="fq2"),

        pseudoalignment = expand(
            os.path.join(
                config['output_dir'],
                "samples",
                "{sample}",
                "quant_kallisto",
                "{sample}.{seqmode}.kallisto.pseudo.sam"),
            zip,
            sample=[i for i in pd.unique(samples_table.index.values)],
            seqmode=[get_sample('seqmode', search_id='index', search_value=i) 
                for i in pd.unique(samples_table.index.values)]),

        TIN_boxplot_PNG = os.path.join(
            config['output_dir'],
            "TIN_scores_boxplot_mqc.png"),

        TIN_boxplot_PDF = os.path.join(
            config['output_dir'],
            "TIN_scores_boxplot_mqc.pdf"),

        alfa_concat_out = os.path.join(
            config["output_dir"],
            "ALFA",
            "ALFA_plots_mqc.png"),

        multiqc_config = os.path.join(
            config["output_dir"],
            "multiqc_config.yaml")

    output:
        multiqc_report = directory(
            os.path.join(
                config["output_dir"],
                "multiqc_summary"))

    params:
        results_dir = os.path.join(
            config["output_dir"]),
        log_dir = config["log_dir"]

    log:
        stderr = os.path.join(
            config["log_dir"],
            "multiqc_report.stderr.log"),
        stdout = os.path.join(
            config["log_dir"],
            "multiqc_report.stdout.log")

    singularity:
        "docker://ewels/multiqc:1.7"

    shell:
        "(multiqc \
        --outdir {output.multiqc_report} \
        --config {input.multiqc_config} \
        {params.results_dir} \
        {params.log_dir};) \
        1> {log.stdout} 2> {log.stderr}"


rule sort_bed_4_big:
    '''
        sort bedGraphs in order to work with bedGraphtobigWig
    '''
    input:
        bg = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "ALFA",
            "{unique}",
            "{sample}.{unique}.{strand}.bg")

    output:
        sorted_bg = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "bigWig",
            "{unique}",
            "{sample}_{unique}_{strand}.sorted.bg")

    singularity:
        "docker://cjh4zavolab/bedtools:2.27"

    log:
        stderr = os.path.join(
            config["log_dir"],
            "samples",
            "{sample}",
            "sort_bg_{unique}_{strand}.stderr.log")

    shell:
        "(sortBed \
         -i {input.bg} \
         > {output.sorted_bg};) 2> {log.stderr}"

rule prepare_bigWig:
    '''
        bedGraphtobigWig, for viewing in genome browsers
    '''
    input:
        sorted_bg = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "bigWig",
            "{unique}",
            "{sample}_{unique}_{strand}.sorted.bg"),
        chr_sizes = lambda wildcards:
            os.path.join(
                config['star_indexes'],
                get_sample(
                    'organism',
                    search_id='index',
                    search_value=wildcards.sample),
                get_sample(
                    'index_size',
                    search_id='index',
                    search_value=wildcards.sample),
                "STAR_index",
                "chrNameLength.txt")

    output:
        bigWig = os.path.join(
            config["output_dir"],
            "samples",
            "{sample}",
            "bigWig",
            "{unique}",
            "{sample}_{unique}_{strand}.bw")

    singularity:
        "docker://zavolab/bedgraphtobigwig:4-slim"

    log:
        stderr = os.path.join(
            config["log_dir"],
            "samples",
            "{sample}",
            "bigwig_{unique}_{strand}.stderr.log"),

        stdout = os.path.join(
            config["log_dir"],
            "samples",
            "{sample}",
            "bigwig_{unique}_{strand}.stdout.log")

    shell:
        "(bedGraphToBigWig \
         {input.sorted_bg} \
         {input.chr_sizes} \
         {output.bigWig};) \
         1> {log.stdout} 2> {log.stderr}"