rnaseq.rules

#
#  Copyright (c) 2016-2021 Sequana Dev Team (https://sequana.readthedocs.io)
#
#  The full license is in the LICENSE file, distributed with this software.
#
#  Website:       https://github.com/sequana/sequana
#  Documentation: http://sequana.readthedocs.io
#  Contributors:  https://github.com/sequana/sequana/graphs/contributors
##############################################################################
# standard modules
import glob
import os
import shutil
import subprocess

import sequana
from sequana_pipetools import snaketools as sm
import sequana.featurecounts as fc

# ========================================================= The main config file
#
configfile: "config.yaml"


# ================================================== The sequana pipeline manager
#
manager = sm.PipelineManager("rnaseq", config)

expected_output = []
# ========================================= Define output of the pipeline
#
manager.globals = {}

if manager.config['general']['aligner'] == 'salmon':
    manager.globals['strand_summary'] = None
    rule rnaseq:
        input:
            "multiqc/multiqc_report.html",
            ".sequana/rulegraph.svg",
            "post_analysis/rnadiff.sh"
else:
    manager.globals['strand_summary'] = "outputs/strand_summary.csv"
    rule rnaseq:
        input:
            "multiqc/multiqc_report.html",
            ".sequana/rulegraph.svg",
            "post_analysis/rnadiff.sh",


# ========================================= Define genome directory and inputs
# Make sure it is absolute
#
genome_directory = os.path.abspath(manager.config["general"]["genome_directory"])
genome_name = genome_directory.rsplit("/", 1)[1]

__prefix_name__ = f"{genome_directory}/{genome_name}"
__fasta_file__ = f"{__prefix_name__}.fa"
__gff_file__   = f"{__prefix_name__}.gff"


# ================================================ Build custom GFF if required
# If we have several features, we need to build a custom GFF file
#
if manager.config['general']['custom_gff']:
    __gff_file__   = manager.config['general']['custom_gff']
assert os.path.exists(__gff_file__), f"GFF file {__gff_file__} does not exist"

# check existence of fasta and gff before starting;
for this in [__fasta_file__, __gff_file__]:
    if os.path.exists(this) is False:
        raise IOError("File {} not found".format(__fasta_file__))


if manager.config.general.contaminant_file and manager.config.general.rRNA_feature:
    logger.error("Either set contaminant_file or rRNA_feature in the config file, not both.")
    sys.exit(1)


# ==================================== search for specific sequences as contaminants
#
if manager.config.general.contaminant_file:
    # Could be a local file of in the genome directory file
    __bowtie1_index_rna__fasta = f"{manager.config.general.contaminant_file}"

    # if not found locally, try to find it in the genome_directory path
    if os.path.exists(__bowtie1_index_rna__fasta) is False:
        __bowtie1_index_rna__fasta = f"{genome_directory}/{manager.config.general.contaminant_file}"
        if os.path.exists(__bowtie1_index_rna__fasta) is False:
            logger.error("File {} does not exists. Check your config file".format(__bowtie1_index_rna__fasta))
            sys.exit(1)

    # we will copy the file to keep the information
    os.makedirs("inputs/contamination_file", exist_ok=True)
    shutil.copy(__bowtie1_index_rna__fasta, "inputs/contamination_file")

    # so we need to rename the input
    __bowtie1_index_rna__fasta = "inputs/contamination_file/" + os.path.basename(__bowtie1_index_rna__fasta)

    bowtie1_index_conta__input_reference = __bowtie1_index_rna__fasta
    bowtie1_index_conta__output = f"{__bowtie1_index_rna__fasta}.1.ebwt"
    rule samtools_faidx:
        input:
            __bowtie1_index_rna__fasta
        output:
            __bowtie1_index_rna__fasta + ".fai"
        container:
            config['apptainers']['sequana_tools']
        shell:
            """
            samtools faidx {input[0]}
            """

elif manager.config.general.rRNA_feature:
    # extract the rRNA feature from the GFF file. Build the corresponding FastA
    # file. if not found, a dummy FastA file with  AAAAAAAAAAAAAA is built

    bowtie1_index_conta__input_reference = f"{__prefix_name__}_{manager.config.general.rRNA_feature}.fa"
    bowtie1_index_conta__input_gff = f"{__prefix_name__}_{manager.config.general.rRNA_feature}.gff"
    bowtie1_index_conta__output = f"{__prefix_name__}_{manager.config.general.rRNA_feature}.1.ebwt"
    rule extract_fasta:
        input:
            fasta = __fasta_file__,
            gff = __gff_file__
        params:
            feature = config['general']['rRNA_feature']
        output:
            fasta = bowtie1_index_conta__input_reference,
            fai = bowtie1_index_conta__input_reference + ".fai",
            gff = bowtie1_index_conta__input_gff
        log:
            "logs/indexing/get_rRNA.log"
        container:
            config['apptainers']['sequana_tools']
        shell:
            """
            # used to be gawk but awk is more generic.
            awk '{{ if ($3=="{params.feature}") print }}' {input.gff} > {output.gff}
            if [ -s {output.gff} ]
            then
                bedtools getfasta -fi {input.fasta} -bed {output.gff}  -fo {output.fasta}
            else :
                echo -e ">empty\\nAAAAAAAAAAAAAA" > {output.fasta}
            fi
            samtools faidx {output.fasta}
            """

# ========================================================= Indexing for rRNA and contmination
#
# redo the indexing whatsover since it is pretty fast
if manager.config.general.contaminant_file or manager.config.general.rRNA_feature:

    # identify ribosomal contamination or contamination
    rule ribosomal_contamination:
        input:
            reference = bowtie1_index_conta__input_reference,
            fai = bowtie1_index_conta__input_reference + ".fai"
        output:
            bowtie1_index_conta__output
        log:
            "logs/indexing/bowtie1_index_conta.log"
        params:
            options=""
        threads: 2
        container:
            config['apptainers']['sequana_tools']
        wrapper:
            f"{manager.wrappers}/wrappers/bowtie1/build"


# ============================================================================ bowtie2 index
#
if manager.config.general.aligner == "bowtie2":
    if manager.config['bowtie2_mapping']['genome_size_larger_than_4gb']:
        bt2_ext = "bt2l"
    else:
        bt2_ext = "bt2"

    # These two variables are used elsewhere and in the rule below
    # Index creatin may differ from one version to another and one may want to
    # keep track of the index and its version. We try to retrieve the version
    # and if succesful, we will add the version as a suffix, otherwise just the
    # name of bowtie2

    # tested on version 2.4.2
    try:
        p = subprocess.Popen(['bowtie2'], stderr=subprocess.PIPE)
        p.wait()
        stderr = p.stderr.read().decode().split("\n")
        hits = [line for line in stderr if "version" in line and 'Bowtie' in line]
        bowtie2_version = "_" + hits[0].split("version")[1].split()[0].strip()
    except Exception:  # various type of exception may occur here
         logger.warning(f"Could not determine bowtie2 version. Index will be stored in {genome_directory}/bowtie2/")
         bowtie2_version = ""

    bowtie2_index = f"{genome_directory}/bowtie2{bowtie2_version}/{genome_name}"

    if os.path.exists(f"{bowtie2_index}.1.{bt2_ext}"):
        pass # index exists, no need to do it, everything should be fine
    else:
        rule bowtie2_index:
            input:
                reference=__fasta_file__
            output:
                multiext(
                    bowtie2_index,
                    ".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2",
                ),

            log:
                "logs/indexing/bowtie2_genome.log"
            params:
                options=config["bowtie2_index"]["options"]
            threads:
                config["bowtie2_index"]["threads"]
            container:
                config['apptainers']['sequana_tools']
            resources:
                **config['bowtie2_index']['resources']
            wrapper:
                f"{manager.wrappers}/wrappers/bowtie2/build"


# ============================================================================ star index
#
elif manager.config.general.aligner  == "star":
    # tested on version 2.7.8a
    try:
        p = subprocess.Popen(['STAR', '--version'], stdout=subprocess.PIPE)
        p.wait()
        star_version = p.stdout.read().decode().strip()
    except Exception:  # various type of exception may occur here
         logger.warning(f"Could not determine STAR version. Index will be stored in {genome_directory}/star/")
         star_version = ""


    __star_index__dir__ = genome_directory + f"/star{star_version}"
    __star_index__done =  f"{__star_index__dir__}/star.done"

    if not os.path.exists(__star_index__done):

        rule star_index:
            input:
                fasta =  __fasta_file__
            output:
                done = __star_index__done
            params:
                options= config['star_index']['options'],
                wkdir= __star_index__dir__
            threads:
                config["star_index"]['threads']
            log:
                "logs/indexing/star_genome.log"
            container:
                config['apptainers']['sequana_tools']
            resources:
                **config['star_index']['resources']
            wrapper:
                f"{manager.wrappers}/wrappers/star/index"


# ========================================================================== salmon
#
elif manager.config.general.aligner == "salmon":
    #tested on salmon 1.4.0
    try:
        p = subprocess.Popen(['salmon', '--version'], stdout=subprocess.PIPE)
        p.wait()
        salmon_version = p.stdout.read().decode().split()[-1]
    except Exception:  # various type of exception may occur here
         logger.warning(f"Could not determine salmon version. Index will be stored in {genome_directory}/salmon/")
         salmon_version = ""

    if os.path.exists(genome_directory + f"/salmon{salmon_version}/salmon.done"):
        pass # index exists, no need to do it, everything should be fine
    else:
        rule salmon_index:
            input:
                fasta=__fasta_file__,
                gff=__gff_file__
            output:
                done=genome_directory + f"/salmon{salmon_version}/salmon.done"
            threads:
                config['salmon_index']['threads']
            resources:
                **config["salmon_mapping"]['resources']
            params:
                options=config['salmon_index']['options']
            container:
                config['apptainers']['salmon']
            log:
                "logs/salmon_indexing.log"
            wrapper:
                f"{manager.wrappers}/wrappers/salmon/index"



# ===================================================================== FASTQC on input data set
#
if not manager.config['fastqc']['skip_fastqc_raw']:
    rule fastqc_raw:
        input:
            manager.getrawdata()
        output:
            done = "{sample}/fastqc_raw/fastqc.done"
        params:
            options= config["fastqc"]["options"],
            working_directory= "{sample}/fastqc_raw/"
        threads: config["fastqc"]["threads"]
        container:
            config['apptainers']['fastqc']
        log:
            "{sample}/fastqc_raw/fastqc.log"
        wrapper:
            f"{manager.wrappers}/wrappers/fastqc"

    expected_output.extend(expand("{sample}/fastqc_raw/fastqc.done", sample=manager.samples))


# ================================================================== trimming
valid_trimmer = ['cutadapt', 'fastp', 'atropos']
if manager.config.trimming.software_choice not in valid_trimmer:
    print(f"Invalid choice for trimming tool. Choose one in {valid_trimmer}")
    sys.exit(1)

if manager.config.trimming.do is False:
    __clean_fastq__output = manager.getrawdata()
elif manager.config.trimming.software_choice in ["cutadapt", "atropos"]:
    adapter_tool = manager.config.trimming.software_choice

    fwd = manager.config.cutadapt.fwd
    rev = manager.config.cutadapt.rev

    if adapter_tool in ["cutadapt", "atropos"]:
        adapter_tool = "cutadapt"
        __cutadapt__input_fastq = manager.getrawdata()
        __cutadapt__wkdir = "{sample}/cutadapt"
        __cutadapt__output = ["{sample}/cutadapt/{sample}_R1_.clean.fastq.gz"]
        if manager.paired:
            __cutadapt__output += ["{sample}/cutadapt/{sample}_R2_.clean.fastq.gz"]

        # Set the fwd and rev adapters
        __cutadapt__fwd = manager.config.cutadapt.fwd
        __cutadapt__rev = manager.config.cutadapt.rev

        __cutadapt__options = manager.config.cutadapt.options
        __cutadapt__mode = manager.config.cutadapt.mode
        __cutadapt__log = "%s/cutadapt/cutadapt.txt" % manager.sample
        __cutadapt__sample = manager.sample
        __clean_fastq__output = __cutadapt__output
        include: sm.modules["cutadapt"]
elif manager.config.trimming.software_choice == "fastp":

    __clean_fastq__output = ["{sample}/fastp/{sample}_R1_.fastp.fastq.gz"]
    if manager.paired:
        __clean_fastq__output += ["{sample}/fastp/{sample}_R2_.fastp.fastq.gz"]

    _quality = config["fastp"].get("quality", 15)
    _minlen = config["fastp"].get("minimum_length", 20)

    options_fastp = config["fastp"].get("options", "")
    options_fastp += f" --qualified_quality_phred {_quality}"
    options_fastp += f" -l {_minlen}"
    if config["fastp"].get("disable_adapter_trimming", False) is True:
        options_fastp += "--disable_adapter_trimming"
    if config["fastp"].get("disable_quality_filtering", False) is True:
        options_fastp += "--disable_quality_filtering"

    rule fastp:
        input: 
            sample=manager.getrawdata()
        output:
            trimmed=__clean_fastq__output,
            html="{sample}/fastp/fastp_{sample}.html",
            json="{sample}/fastp/fastp_{sample}.json", # must be named fastp
        log:
            "logs/fastp/{sample}.log"
        params:
            options=config['fastp']["options"],
            adapters=config["fastp"]["adapters"]
        threads:
            config["fastp"].get("threads", 4)
        resources:
            **config['fastp']['resources']
        container:
            config['apptainers']['fastp']
        wrapper:
            f"{manager.wrappers}/wrappers/fastp"


# ===================================================== FASTQC fastp results
#
rule fastqc_clean:
    input:
        __clean_fastq__output
    output:
        done = "{sample}/fastqc_clean/fastqc.done"
    params:
        options= config["fastqc"]["options"],
        working_directory= "{sample}/fastqc_clean/"
    threads: config["fastqc"]["threads"]
    log:
        "{sample}/fastqc_clean/fastqc.log"
    resources:
        **config["fastqc"]['resources']
    container:
        config['apptainers']['fastqc']
    wrapper:
        f"{manager.wrappers}/wrappers/fastqc"
expected_output.extend(expand("{sample}/fastqc_clean/fastqc.done", sample=manager.samples))


# ================================= Decompress fastq.gz file before running bowtie1
#
if manager.config.trimming.software_choice == 'cutadapt' and manager.config.trimming.do:
    #__unpigz_R1__input = manager.getname("cutadapt", "_R1_.clean.fastq.gz")
    __unpigz_R1__input = "{sample}/cutadapt/{sample}_R1_.clean.fastq.gz"
elif manager.config.trimming.software_choice == 'fastp' and manager.config.trimming.do:
    __unpigz_R1__input = "{sample}/fastp/{sample}_R1_.fastp.fastq.gz"
elif manager.config.trimming.software_choice == 'atropos' and manager.config.trimming.do:
    __unpigz_R1__input = "{sample}/atropos/{sample}_R1_.clean.fastq.gz"
else:
    __unpigz_R1__input = manager.getrawdata()


# ==========  decompress and sanity check
#
if int(config['bowtie1_mapping_rna']['nreads']) != -1:
    extra = int(config['bowtie1_mapping_rna']['nreads']) % 4
    config['bowtie1_mapping_rna']['nreads'] -= extra


rule sample_rRNA:
    input:
        __unpigz_R1__input
    output:
        fastq=temp("{sample}/data_for_bowtie1/{sample}_R1_.fastq")
    threads: 4
    params:
        nreads = int(config['bowtie1_mapping_rna']['nreads'])
    shell:
        """
        set +o pipefail
        if [[ {params.nreads} == -1 ]]; then
            unpigz -p {threads} -fk --stdout {input[0]} > {output[0]}
        else
            unpigz -p {threads} -fk --stdout {input[0]} | head -n {params.nreads} > {output[0]}
        fi
        """

"""With paired data, alignement on rRNA leads to 0% alignment if we use R1 and
R2. If we use R1 only, the percentage is >0. First reason is that reads are not
trimmed properly. In truth, bowtie2 supports local alignments which means it can
soft-clip non-matching (=adapter) content while still align the local part of
the read that matches the reference. With Bowtie1 the read will probably go
unaligned due to the many mismatches. So we do not include R2 from version
v0.9.14.
"""

# ========================================== bowtie1 mapping to detect rRNA

if manager.config.general.rRNA_feature or manager.config.general.contaminant_file:
    # rRNA. Note the list here below because the rule expects a list (in case it
    # is paired

    rule bowtie1_mapping_rna:
        input:
            fastq= rules.sample_rRNA.output.fastq,
            index=bowtie1_index_conta__output
        output:
            bam = "{sample}/bowtie1_mapping_rna/{sample}_rRNA.bam",
            sorted = "{sample}/bowtie1_mapping_rna/{sample}_rRNA.sorted.bam",
        log:
            "{sample}/bowtie1_mapping_rna/{sample}_bowtie1.log"
        params:
            options=""
        threads:
            config['bowtie1_mapping_rna']['threads']
        container:
            config['apptainers']['sequana_tools']
        wrapper:
            f"{manager.wrappers}/wrappers/bowtie1/align"

    rule fix_bowtie1_log:
        input:
            expand("{sample}/bowtie1_mapping_rna/{sample}_bowtie1.log", sample=manager.samples)
        output:
            "logs/fix_bowtie1/fix_bowtie1.log"
        run:

            for filename in input:
                # we read the file
                with open(filename) as fin:
                    data = fin.readlines()
                # we update the file
                with open(filename, "w") as fout:
                    for line in data:
                        if "least one alignment" in line:
                            fout.write(line)
                            fout.write(line.replace("least one alignment", "least one reported alignment"))
                        else:
                            fout.write(line)
            with open(output[0], "w") as fout:
                fout.write("")
    expected_output += ["logs/fix_bowtie1/fix_bowtie1.log"]


# ========================================================== bowtie2 mapping
if manager.config.general.aligner == "bowtie2":

    rule bowtie2_mapping:
        input:
            fastq=__clean_fastq__output,
            idx=multiext(
            bowtie2_index,
            ".1.bt2",
            ".2.bt2",
            ".3.bt2",
            ".4.bt2",
            ".rev.1.bt2",
            ".rev.2.bt2",
        ),

        output:
            bam="{sample}/bowtie2/{sample}.sorted.bam",
        log:
            "{sample}/bowtie2/{sample}.log"
        params:
            options=config["bowtie2_mapping"]["options"],
        threads:
            config["bowtie2_mapping"]["threads"]
        container:
            config['apptainers']['sequana_tools']
        resources:
             **config['bowtie2_mapping']['resources']
        wrapper:
            f"{manager.wrappers}/wrappers/bowtie2/align"

    __mapping_output = "{sample}/bowtie2/{sample}.sorted.bam"


# ========================================================== star mapping
elif manager.config.general.aligner == "star":
    # Mapper rna-star

    rule star_mapping:
        input:
            fastq= __clean_fastq__output,
            reference=__fasta_file__,
            index= __star_index__done
        output:
            bam = "{sample}/star_mapping/{sample}_Aligned.sortedByCoord.out.bam"
        params:
            options=config['star_mapping']['options'],
            # for legacy mapping, set the first and second pass options
            options_first_pass=config['star_mapping']['options'],
            options_second_pass=config['star_mapping']['options'],
            prefix = "{sample}/star_mapping/{sample}",
            legacy=True,
        threads:
            config['star_mapping']['threads']
        log:
            "{sample}/star_mapping/{sample}.log"
        container:
            config['apptainers']['sequana_tools']
        resources:
            **config['star_mapping']['resources']
        wrapper:
            f"{manager.wrappers}/wrappers/star/align"


    expected_output.extend(
        expand(
            "{sample}/star_mapping/{sample}_Aligned.sortedByCoord.out.bam",
            sample=manager.samples
            )
    )

    __mapping_output = "{sample}/star_mapping/{sample}_Aligned.sortedByCoord.out.bam"

# ========================================================== salmon mapping

elif manager.config.general.aligner == "salmon":
    # to be used later by salmon_to_features
    __salmon_mapping__output_counts = "{sample}/salmon_mapping/{sample}_quant.sf"

    rule salmon_mapping:
        input:
            fastq=__clean_fastq__output,
            index=genome_directory + f"/salmon{salmon_version}/salmon.done"
        output:
            quant="{sample}/salmon_mapping/{sample}_quant.sf"
        params:
            options=config['salmon_mapping']['options']
        threads:
            config['salmon_mapping']['threads']
        resources:
            **config["salmon_mapping"]['resources']
        container:
            config['apptainers']['salmon']
        log:
            "{sample}/salmon_mapping/salmon.log"
        wrapper:
            "add_salmon/wrappers/salmon/align"


    expected_output.extend(expand("{sample}/salmon_mapping/{sample}_quant.sf", sample=manager.samples))
    # There is no BAM created
    __mapping_output = None


# ========================================================== add_read_group
# The input is the output of the mapping
# Add Read group on BAM files
if manager.config.general.aligner not in ['salmon']:
    rule add_read_group:
        input:
            __mapping_output
        output:
            "{sample}/add_read_group/{sample}.sorted.bam"
        log:
            "{sample}/add_read_group/{sample}.log"
        params:
            options=config["add_read_group"]["options"],
            SM="{sample}"
        container:
            config['apptainers']['sequana_tools']
        wrapper:
            f"{manager.wrappers}/wrappers/add_read_group"



# we always add read group so input is the read group output
# output is stored in __final_bam__
# duplicates can be from PCR or if SR, by pure chance.
# if Paired, most likely a PCR origin.
# Mark duplicates
if config["mark_duplicates"]["do"]:
    rule mark_duplicates:
        input:
            "{sample}/add_read_group/{sample}.sorted.bam"
        output:
            bam = "{sample}/mark_duplicates/{sample}.sorted.markdup.bam",
            metrics = "{sample}/mark_duplicates/{sample}.sorted.markdup.metrics",
        log:
            out = "{sample}/mark_duplicates/log.out",
            err = "{sample}/mark_duplicates/log.err"
        params:
            remove_dup = "false",
            tmpdir = "{sample}/mark_duplicates/tmp"
        container:
            config['apptainers']['sequana_tools']
        resources:
            **config['mark_duplicates']['resources']
        wrapper:
            f"{manager.wrappers}/wrappers/mark_duplicates"
    __final_bam__ = "{sample}/mark_duplicates/{sample}.sorted.markdup.bam"
elif manager.config.general.aligner not in ['salmon']:
    __final_bam__   = "{sample}/add_read_group/{sample}.sorted.bam"
else:
    __final_bam__ = []


# ====================================================================== generating bigwig files
if manager.config.bam_coverage.do is True and config['general']['aligner'] not in ['salmon']:

    rule bam_coverage:
        input: __final_bam__
        output:
            "{sample}/bam_coverage/{sample}.norm.bw"
        params:
            options = config['bam_coverage']["options"]
        log:
            "{sample}/bam_coverage/{sample}.log"
        threads:
            config['bam_coverage']['threads']
        container:
            config['apptainers']['sequana_tools']
        resources:
            **config['bam_coverage']['resources']
        wrapper:
            f"{manager.wrappers}/wrappers/deeptools/bam_coverage"

    expected_output.extend(
        expand(
            "{sample}/bam_coverage/{sample}.norm.bw",
            sample=manager.samples
            )
    )


# ============================================================= generating IGV plots
if manager.config.igvtools.do and config['general']['aligner'] not in ['salmon']:
    # if nothing provided, it must be an empty string
    if manager.config.igvtools.chrom_sizes_file.strip():
        pass
    else:
        config["igvtools"]["chrom_sizes_file"] = __fasta_file__

    rule igvtools:
        input: __final_bam__
        output:
            "{sample}/igvtools/{sample}.tdf"
        log:
            "{sample}/igvtools/{sample}.log"
        params:
            chromSize=config['igvtools']['chrom_sizes_file']
        container:
            config['apptainers']['igvtools']
        threads: 4
        shell:
            """
            igvtools count -z 5 -w 25 -f mean,max --includeDuplicates {input} {output} {params.chromSize}
            """
    expected_output.extend(expand("{sample}/igvtools/{sample}.tdf", sample=manager.samples))


# ===================================================================== Feature counts from subread suite
if manager.config.general.aligner == "salmon":
    __feature_counts__input = __salmon_mapping__output_counts
else :
    __feature_counts__input = __final_bam__

fc_outdir = "post_analysis/feature_counts/"

if manager.config.feature_counts.do and manager.config.general.aligner not in ['salmon']:
    # Guessing strandness is not always straightfoward; Even when we set it; 
    # collaborators may want to look at the other options. So, we compute
    # everything with the 3 different options of strandness.
    # We will copy one of them based on our criteria, but all 3 will be
    # available

    feature_type = config['feature_counts']['feature']
    if "," in feature_type: # assume this is a custom GFF file
        feature_type = "custom"
    else:
        feature_type = config['feature_counts']["feature"]

    if config['feature_counts']['extra_attributes']:
        fc_options = f" {config['feature_counts']['options']} "
        fc_options += " --extraAttributes {} ".format(config['feature_counts']['extra_attributes'])
    else:
        fc_options = f" {config['feature_counts']['options']} "

    if manager.paired:
        fc_options += " -p "


    # ======================= calls feature counts 3 times Nsamples here below
    strand = [0,1,2]
    rule feature_counts:
        input:
            bam=__feature_counts__input,
            gff=__gff_file__
        output:
            counts="{sample}/feature_counts/{strand}/{sample}_feature.out",
            summary="{sample}/feature_counts/{strand}/{sample}_feature.out.summary"
        params:
            options=fc_options,
            feature=feature_type,
            attribute=config['feature_counts']["attribute"],
            strandness="{strand}"
        threads:
            config["feature_counts"]['threads']
        container:
            config['apptainers']['sequana_tools']
        log:
            "{sample}/feature_counts/{strand}/feature_counts.log"
        wrapper:
            f"{manager.wrappers}/wrappers/feature_counts"


    # ===================== guessing the strand 
    #
    __guess_strandness__output = expand(fc_outdir + "{sample}_feature.out", sample=manager.samples)
    rule guess_strandness:
        """Guessing strandnes"""
        input:
            counts = expand("{sample}/feature_counts/{strand}/{sample}_feature.out", sample=manager.samples, strand=[0,1,2])
        output:
            data=__guess_strandness__output,
            summary=manager.globals['strand_summary']
        run:
            # We compute all strandness
            import sequana.featurecounts as fc

            mfc = fc.MultiFeatureCount(rnaseq_folder=".", 
                        tolerance=manager.config.feature_counts.tolerance)
            mfc.df.to_csv(output.summary)

            try:
                mfc.plot_strandness(savefig=True, output_filename="outputs/strand_summary.png")
            except Exception as err:
                logger.warning("Could not create plot_strandness")

            logger.info(f"strandness inference: {mfc.probable_strand}")
            msg = f"This is {mfc.probable_strand} data (check in the multiqc report)"
            if mfc.probable_strand in [0, 1, 2]:
                choice = mfc.probable_strand
                logger.info(msg)
            else:
                logger.warning("Strandness is apparently neither of 0, 1, 2")
                logger.warning("you will need to copy the feature counts files yourself in ./feature_counts")
                choice = -1

            # If user knowns what he/she wants we overwrite the choice
            if "strandness" in config['feature_counts'] and config["feature_counts"]["strandness"]:
                user_choice =  int(config["feature_counts"]["strandness"])
                if user_choice in [0,1,2]:
                    choice = user_choice
                else:
                    logger.error(f"strandness in the config file must be 0,1,2. You gave {user_choice}")
                    sys.exit(1)

            if choice in {0, 1, 2}:
                for filename in input:
                    if f"feature_counts/{choice}/" in filename:
                        shell(f"cp {filename} {fc_outdir}")
                        shell(f"cp {filename}.summary {fc_outdir}")
            else:
                # if not clear, we copy everything and users should clean up the directory
                for filename in input.fc0:
                    shell("cp {} {}".format(filename, fc_outdir))
                    shell("cp {}.summary {}".format(filename, fc_outdir))
                for filename in input.fc1:
                    shell("cp {} {}".format(filename, fc_outdir))
                    shell("cp {}.summary {}".format(filename, fc_outdir))
                for filename in input.fc2:
                    shell("cp {} {}".format(filename, fc_outdir))
                    shell("cp {}.summary {}".format(filename, fc_outdir))
elif manager.config.feature_counts.do and manager.config.general.aligner in ['salmon']:

    __salmon_to_features__output = fc_outdir + "{sample}_feature.out"
    rule salmon_to_features:
        input: __salmon_mapping__output_counts
        output: __salmon_to_features__output
        params:
            gff=__gff_file__
        shell:
            """sequana salmon --input {input} --output {output} --gff {params.gff} --attribute ID   """
    expected_output += expand(__salmon_to_features__output, sample=manager.samples)


# ==================================================================== Guess strandness

if manager.config.general.aligner in ['salmon']:
    __guess_strandness__output = expand(__salmon_to_features__output, sample=manager.samples)
rule merge_feature_counts:
    input: __guess_strandness__output
    output: "post_analysis/all_features.out"
    run:
        from sequana.featurecounts import FeatureCountMerger
        fcm = FeatureCountMerger(fof=input)
        fcm.to_tsv(output[0])
expected_output.append("post_analysis/all_features.out")

# ================================================================== rseqc diag tool


if config['rseqc']['do']:

    rule gff2bed:
        input:
            gff=__gff_file__
        output:
            bed="tmp/temp.bed" # config['rseqc'].get('bed_file', "tmp/temp.bed")
        message: "Build BED file from GFF using Sequana"
        run:
            from sequana import GFF3
            g = GFF3(input[0])
            g.to_bed(output[0], 'Name')

    rule rseqc:
        input:
            bam=__final_bam__,
            bed=rules.gff2bed.output.bed
        # no need to put all outputs
        output:
            bam_stat= "{sample}/rseqc/{sample}_bam_stat.txt",
            #read_gc="{sample}/rseqc/{sample}.GC.xls",
            geneBody_coverage="{sample}/rseqc/{sample}.geneBodyCoverage.txt"
        params:
            paired="PE" if manager.paired else "SE"
        log:
            "{sample}/rseqc/{sample}.log"
        container:
            config['apptainers']['sequana_tools']
        shell:
            """
            # for paired data only
            inner_distance.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample} -r {input.bed} &>{log}

            # For now GC not very useful in the output so commented
            # read_GC.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample} &>{log}

            # genebody coverage
            geneBody_coverage.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample} -r {input.bed} &>{log}

            # uses bigwig redundant with geneBody_coverage
            # geneBody_coverage2.py -i {wildcards.sample}/bamCoverage/{wildcards.sample}.norm.bw  -o {wildcards}/rseq/{wildcards.sample} -r test.bed &>{log}

            # Not included in the multiqc module so commented for now
            #clipping_profile.py -i {input.bam} -s {params.paired} -o {wildcards.sample}/rseqc/{wildcards.sample}
            read_duplication.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample} &>{log}
            junction_annotation.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample} -r {input.bed}  &>{log}
            junction_saturation.py -i {input.bam} -o {wildcards.sample}/rseqc/{wildcards.sample}  -r {input.bed} &>{log}
            infer_experiment.py -i {input.bam} -r {input.bed} > {wildcards.sample}/rseqc/{wildcards.sample}.infer.txt &>{log}

            # bam stats  KEEP last since that is the expected output to make sure
            # previous files (not listed in output:) are computed first.
            bam_stat.py -i {input.bam} > {output.bam_stat} &>{log}
            """

    # one series is enough. if bam_stats is created, others are also created
    expected_output.extend(
        expand("{sample}/rseqc/{sample}_bam_stat.txt", sample=manager.samples)
    )

# ========================================================== RNAseqc diag tool
# No need for mark_duplicates for RNASEQC . Just use the BAM file
if config["rnaseqc"]["do"] and config['general']['aligner'] != 'salmon':


    # for multiqc, important that output directories are called rnaseqc
    __gtf_file__ = config['rnaseqc']['gtf_file'].strip()

    # Could be a local file. If provided,
    if __gtf_file__:
        if os.path.exists(__gtf_file__) is False:
            logger.error(f"{__gtf_file__} not found")
            sys.exit(1)
    else: # if gtf not provided, maybe available in the genome directory ?
        __gtf_file__   = __prefix_name__ + ".gtf"
        if os.path.exists(__gtf_file__) is False:
            gdir = config["general"]["genome_directory"]
            from sequana import logger as logs  # to not interfere with snakemake
            logs.critical(f"{__gtf_file__} not found in {gdir}. Trying the GFF file")
            __gtf_file__ = __prefix_name__ + ".gff"


    rule rnaseqc_fixup:
        input:
            gtf = __gtf_file__
        output:
            gtf = temp("tmp/test.gtf")
        run:
            # If input GTF has no exon or genes, an error message is printed and 
            # no files are created. This seems to be an issue in rnaseqc. 
            # So, we create dummy gene and exon
            with open(output.gtf, "w") as ff:
                ff.write(open(input['gtf'], "r").read())
                ff.write('myCHR\tSGD\tgene\t0\t0\t.\t+\t0\tgene_id "dummy"\n')
                ff.write('myCHR\tSGD\texon\t0\t0\t.\t+\t0\texon_id "dummy"\n')
                ff.close()

    rule rnaseqc:
        input:
            bam = __final_bam__,
            gtf = rules.rnaseqc_fixup.output.gtf
        output:
            metrics = "{sample}/rnaseqc/{sample}.metrics.tsv"
        log:
            "{sample}/rnaseqc/{sample}.log",
        params:
            directory = "{sample}/rnaseqc",
            options= config['rnaseqc']['options']
        resources:
            **config["rnaseqc"]["resources"]
        container:
            config["apptainers"]["rnaseqc"]
        shell:
            """
            rnaseqc {input.gtf} {input.bam} {params.directory} -s {wildcards.sample} {params.options} &>{log}
            """

    expected_output.extend(expand("{sample}/rnaseqc/{sample}.metrics.tsv", sample=manager.samples))


# ========================================================== multiqc
multiqc_params_options = config['multiqc']['options']
if manager.config.multiqc.config_file:
    multiqc_params_options += f" -c {manager.config.multiqc.config_file}"



rule multiqc:
    input:
        expected_output
    output:
       "multiqc/multiqc_report.html"
    params:
        options=multiqc_params_options,
        input_directory=config['multiqc']['input_directory'],
        config_file=config['multiqc']['config_file'],
        modules=config['multiqc']['modules']
    log:
        "multiqc/multiqc.log"
    resources:
        **config['multiqc']['resources']
    container:
        config["apptainers"]["multiqc"]
    wrapper:
       f"{manager.wrappers}/wrappers/multiqc"

# ========================================================== rulegraph

rule rulegraph:
    input:
        workflow.snakefile,
    output:
        "rulegraph/rulegraph.dot",
    params:
        configname="config.yaml",
        mapper = {"multiqc": "../multiqc/multiqc_report.html"},
    wrapper:
        f"{manager.wrappers}/wrappers/rulegraph"


rule dot2svg:
    input:
        "rulegraph/rulegraph.dot"
    output:
        ".sequana/rulegraph.svg"
    container:
        config['apptainers']['graphviz']
    shell:
        """dot -Tsvg {input} -o {output}"""



rule prepare_DGE_analysis:
    input:
        features="post_analysis/all_features.out",
    output:
        rnadiff="post_analysis/rnadiff.sh",
        design="post_analysis/design.csv"
    run:

	    # ------------------------------------------- RNADIFF
	    # 1. save data for the RNADiff analysis
	    from sequana.featurecounts import FeatureCount
	    try:
	        fc = FeatureCount(input[0], guess_design=True)
	        fc.design_df.to_csv(output[1], index=False)
	    except:
	        msg = "Could not build the design.csv file in rnadiff. You will need to create it manually."
	        logger.warning(msg)
	        with open("post_analysis/README.rst", "w") as fout:
	            fout.write(f"""{msg}
	The design.csv file must be formatted as follows (for 2 conditions with 3 replicates each):

	label,condition
	samplename_1,condition_name_1
	samplename_2,condition_name_1
	samplename_3,condition_name_1
	samplename_4,condition_name_2
	samplename_5,condition_name_2
	samplename_6,condition_name_2
	""")

	    # 2. save the script
	    with open(output.rnadiff, "w") as fout:
	        attribute = config['feature_counts']['attribute']
	        feature = config['feature_counts']['feature']
	        fout.write("#/bin/sh\nsequana rnadiff --features all_features.out " +
	            f" --annotation-file {__gff_file__} --design design.csv --feature-name {feature} --attribute-name {attribute}")
	    shell(f"chmod 755 {output.rnadiff}")



# Those rules takes a couple of seconds so no need for a cluster
localrules:  rulegraph, prepare_DGE_analysis


onsuccess:
    # Create plots about stats
    from sequana import logger as log
    from sequana.modules_report.summary import SequanaReport

    import colorlog
    log = colorlog.getLogger("sequana.rnaseq")
    log.setLevel("INFO")
    manager.teardown(
        extra_files_to_remove=["requirements.txt"],
        extra_dirs_to_remove=[".genomes", "tmp", "logs"])
    manager.clean_multiqc("multiqc/multiqc_report.html")


    try:
        import pandas as pd
        df = pd.read_csv(manager.globals['strand_summary'])
        guess = df['strand'].value_counts().idxmax()
        names = {0: 'stranded', 1: 'unstranded', 2: 'reversely stranded'}
        guess = names[guess]
    except Exception as err:
        if config['general']['aligner'] == "salmon":
            logger.info("Salmon aligner used. No strandness information available")
        else:
            logger.warning(err)
            guess = "?"


    intro = f"""
    <h2>Overview</h2>
    <p>
    The RNA-seq pipeline maps the reads on the provided reference (called <i>{__fasta_file__.split("/")[-1]}</i>). Features counts were extracted and are available in the <a href="./post_analysis/feature_counts">feature counts</a> directory; those files are entry points for differential gene expression analysis. The differential analysis, if performed, should be available in the <a href="post_analysis/rnadiff/">DGE analysis directory</a>. In addition, if enrichment was performed (GO or Kegg pathways), it should be available in the <a href="post_analysis"></a> directory as well.
    </p>

    <p>A <a href="multiqc/multiqc_report.html">multiqc report</a> is available, where various QC and mapping quality plots can be visualised. Some important plots are also in the HTML page here below.
    </p>"""


    rRNA_done = True
    intro += """<h2>ribosomal / contaminant content</h2>"""
    try:
        if os.path.exists("multiqc/multiqc_report_data/multiqc_bowtie1.txt"):
            df = pd.read_csv("multiqc/multiqc_report_data/multiqc_bowtie1.txt", sep='\t')
        elif os.path.exists("multiqc/multiqc_data/multiqc_bowtie1.txt"):
            df = pd.read_csv("multiqc/multiqc_data/multiqc_bowtie1.txt", sep='\t')
        else:
            rRNA_done = False

        if rRNA_done:
            if "reads_aligned_percentage" in df.columns:
                rRNA = int(df.reads_aligned_percentage.mean()*100)/100
            else:
                rRNA = int((100 - df.not_aligned_percentage.mean())*100)/100

            if rRNA < 10:
                intro += f"<p>rRNA content (or contaminant provided) represents {rRNA}%, which is low (as expected). "
            elif rRNA < 20:
                intro += f"<p>rRNA content (or contaminant provided) represents {rRNA}%, which is moderately low. "
            elif rRNA > 20 and rRNA <50:
                rRNA += f"<p>rRNA content (or contaminant provided) represents {rRNA}%, which is relatively high. "
            elif rRNA >= 50:
                rRNA += f"<p>rRNA content (or contaminant provided) represents {rRNA}%, which is very high. "
        else:
            intro += f"<p>rRNA content not computed (no rRNA gene or contaminant provided)</p>"

    except Exception as err:
        print(err)
        pass

    try:
        from sequana.multiqc.plots import Bowtie1Reader
        if rRNA_done:
            filename = "multiqc/multiqc_data/multiqc_bowtie1.txt"
            if not os.path.exists(filename):
                filename = "multiqc/multiqc_report_data/multiqc_bowtie1.txt"
            br = Bowtie1Reader(filename)
            br.df.Sample = [str(x).replace("_bowtie1","") for x in br.df.Sample]
            fig = br.plot_bar(html_code=True)
            from plotly import offline
            intro += offline.plot(fig, output_type="div", include_plotlyjs=True)
    except Exception as err:
        print(err)


    # Include the bowtie plot
    intro += """<h2>Mapping rate</h2>"""
    if config['general']['aligner'] == "bowtie2":
        from sequana.multiqc.plots import Bowtie2
        if not manager.paired:
            filename = "multiqc/multiqc_data/mqc_bowtie2_se_plot_1.txt"
            if not os.path.exists(filename):
                filename = "multiqc/multiqc_report_data/mqc_bowtie2_se_plot_1.txt"
        else:
            filename = "multiqc/multiqc_data/mqc_bowtie2_pe_plot_1.txt"
            if not os.path.exists(filename):
                filename = "multiqc/multiqc_report_data/mqc_bowtie2_pe_plot_1.txt"
        br = Bowtie2(filename)
        fig = br.plot(html_code=True)
        from  plotly import offline
        intro += """<p>The mapping was performed with bowtie2. Here below are the percentage of mapping for each sample. See also the multiqc report.</p>""" + offline.plot(fig, output_type="div", include_plotlyjs=True)

    elif config["general"]["aligner"] == "star":
        from sequana.multiqc.plots import STAR
        filename = "multiqc/multiqc_report_data/multiqc_star.txt"
        if not os.path.exists(filename):
            filename = "multiqc/multiqc_data/multiqc_star.txt"
        br = STAR(filename)
        fig = br.plot(html_code=True)
        from  plotly import offline
        intro += """<p>The mapping was performed with STAR. Here below are the percentage of mapping for each sample. See also the multiqc report.</p>""" + offline.plot(fig, output_type="div", include_plotlyjs=True)

    try:
        from sequana.multiqc.plots import FeatureCounts
        intro += """<h2>Annotation rate</h2>"""
        filename = "multiqc/multiqc_report_data/mqc_featureCounts_assignment_plot_1.txt"
        if not os.path.exists(filename):
            filename = "multiqc/multiqc_data/mqc_featureCounts_assignment_plot_1.txt"

        br = FeatureCounts(filename)
        fig = br.plot(html_code=True)
        from  plotly import offline
        intro += """<p>The annotation was performed with subread/feature counts software. Here below is the percentage of reads assigned to the requested feature (usually gene; see the config file here below). </p>""" + offline.plot(fig, output_type="div", include_plotlyjs=True)
    except Exception as err:
        print(err)

    if manager.globals['strand_summary']:
        intro += """<h2>Strandness</h2>"""
        intro+="""
        <p>Here below is a QC plot related to the strandness found for each samples. The red dotted lines indicate a tolerance. The 0.5 vertical line correspond to an <b>unstranded</b> case. A value close to 0 indicates a <b>reversely stranded</b> case, and a value close to 1 indicates a <b>stranded</b> case.
    """.format(config["general"]['aligner'], guess)

        if "strandness" in config['feature_counts'] and config["feature_counts"]["strandness"]:
            choice = config["feature_counts"]["strandness"]
            intro += "User decided to set strandness to: {} </p>".format(choice)
        else:
            intro += "Strandness was guessed from the data. </p>"

        image = SequanaReport.png_to_embedded_png("strand", "outputs/strand_summary.png",
                     style="width:80%; height:40%")
        intro += image

    intro+=    """<h2>Differential analysis</h2>
    <p>
    Differentially expressed genes analysis is not performed automatically with this pipeline. However, information, feature counts, and other materials can be found in the directory <a href="./post_analysis/.">post_analysis</a> where the standalone 'sequana rnadiff' can be used with DeSEq2. Most probably an analysis is present. If so, please open the directory in  <a href="post_analysis/rnadiff/">rnadiff</a> report.
    </p>
    """

    # Now the final report. add the original command in the HTML report
    data = manager.getmetadata()
    s = SequanaReport(data, intro)

    shell("chmod -R g+w .")
    shell("rm -rf rulegraph")

onerror:
    manager.onerror()