shell:
    """
    cat {input} > initial_ALL.fasta
    {params.MAFFT} --nuc --preservecase --maxiterate 1000 --localpair --thread {threads} initial_ALL.fasta > references/initial_aln.fasta 2> >(tee {log.errfile} >&2)
    rm initial_ALL.fasta

    {params.REMOVE_GAPS} references/initial_aln.fasta -o {output} -p 0.5 > {log.outfile} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    mkdir -p {wildcards.dataset}/references/
    {params.EXTRACT_SEQ} {input} -o {output} -s "${{CONSENSUS_NAME}}"
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    mkdir -p {wildcards.dataset}/references/
    cp {input} {output}
    sed -i -e "s/>.*/>${{CONSENSUS_NAME}}/" {output}
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    mkdir -p {wildcards.dataset}/references/
    cp {input} {output}
    sed -i -e "s/>.*/>${{CONSENSUS_NAME}}/" {output}
    """

shell:
    """
    {params.GUNZIP} -c {input} > {output} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    # 1. clean previous run
    rm -rf   {wildcards.dataset}/alignments
    rm -f    {wildcards.dataset}/references/ref_ambig.fasta
    rm -f    {wildcards.dataset}/references/ref_majority.fasta
    mkdir -p {wildcards.dataset}/alignments
    mkdir -p {wildcards.dataset}/references

    # 2. perform alignment # -l = leave temps
    {params.NGSHMMALIGN} -v {params.EXTRA} -R {input.initial_ref} -o {output.good_aln} -w {output.reject_aln} -t {threads} -N "${{CONSENSUS_NAME}}" {params.LEAVE_TEMP} {input.FASTQ} > {log.outfile} 2> >(tee {log.errfile} >&2)

    # 3. move references into place
    mv {wildcards.dataset}/{{alignments,references}}/ref_ambig.fasta
    mv {wildcards.dataset}/{{alignments,references}}/ref_majority.fasta
    """

shell:
    """
    cat {input} > ALL_{wildcards.kind}.fasta
    {params.MAFFT} --nuc --preservecase --maxiterate 1000 --localpair --thread {threads} ALL_{wildcards.kind}.fasta > {output} 2> >(tee {log.errfile} >&2)
    rm ALL_{wildcards.kind}.fasta
    """

shell:
    """
    {params.CONVERT_REFERENCE} -t {params.REF_NAME} -m {input.REF_ambig} -i {input.BAM} -o {output} > {log.outfile} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    echo "Writing BAM file"
    rm -f '{params.sort_tmp}'.[0-9]*.bam
    {params.SAMTOOLS} sort -T "{params.sort_tmp}" -o "{output.BAM}" "{input}" > {log.outfile} 2> >(tee {log.errfile} >&2)
    {params.SAMTOOLS} index "{output.BAM}" >> {log.outfile} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.BWA} index {input} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    {params.BWA} mem -t {threads} {params.EXTRA} -o "{output.TMP_SAM}" "{input.REF}" {input.FASTQ} 2> >(tee {log.errfile} >&2)
    # Filter alignments: (1) remove unmapped reads (single-end) or keep only reads mapped in proper pairs (paired-end), (2) remove supplementary aligments
    {params.SAMTOOLS} view -h {params.FILTER} -F 2048 -o "{output.REF}" "{output.TMP_SAM}" 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.BOWTIE} {input} {input} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    {params.BOWTIE} -x {input.REF} -1 {input.R1} -2 {input.R2} {params.PHRED} {params.PRESET} -X {params.MAXINS} {params.EXTRA} -p {threads} -S {output.TMP_SAM} 2> >(tee {log.errfile} >&2)
    # Filter alignments: (1) keep only reads mapped in proper pairs, and (2) remove supplementary aligments
    {params.SAMTOOLS} view -h -f 2 -F 2048 -o "{output.REF}" "{output.TMP_SAM}" 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.BOWTIE} -x {input.REF} -U {input.R1} {params.PHRED} {params.PRESET} {params.EXTRA} -p {threads} -S {output.TMP_SAM} 2> >(tee {log.errfile} >&2)
    # Filter alignments: (1) remove unmapped reads, and (2) remove supplementary aligments
    {params.SAMTOOLS} view -h -F 4 -F 2048 -o "{output.REF}" "{output.TMP_SAM} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    rm -rf {params.DIR}/*/*/extracted_data
    """

shell:
    """
    rm -rf {params.DIR}/*/*/preprocessed_data
    """

shell:
    """
    rm -rf references/ALL_aln_*.fasta
    rm -rf references/MAFFT_*_cohort.*
    """

shell:
    """
    rm -rf {params.DIR}/*/*/alignments
    rm -rf {params.DIR}/*/*/QA_alignments
    rm -rf {params.DIR}/*/*/references/ref_ambig.fasta
    rm -rf {params.DIR}/*/*/references/ref_majority.fasta
    rm -rf {params.DIR}/*/*/references/initial_consensus.fasta
    """

shell:
    """
    rm -f {input}
    rm -rf {params.DIR}/*/*/alignments
    rm -rf {params.DIR}/*/*/references/ref_ambig*.fasta
    rm -rf {params.DIR}/*/*/references/ref_majority*.fasta
    """

shell:
    """
    rm -f {input}
    rm -rf {params.DIR}/*/*/alignments
    rm -rf {params.DIR}/*/*/references/ref_ambig*.fasta
    rm -rf {params.DIR}/*/*/references/ref_majority*.fasta
    """

shell:
    """
    rm -rf {params.DIR}/*/*/variants/SNVs
    """

shell:
    """
    rm -rf {params.DIR}/*/*/variants/global/contigs_stage_?.fasta
    rm -rf {params.DIR}/*/*/variants/global/stage_?
    """

shell:
    """
    rm {params.DIR}/*/*/variants/global/quasispecies.*
    """

shell:
    """
    rm -rf {params.DIR}/*/*/variants/global/predicthaplo/
    """

shell:
    """
    rm -rf {params.DIR}/*/*/visualization
    """

shell:
    """
    {params.bcftools} mpileup \
        --threads {threads} \
        -Ou \
        -f {input.fname_ref} \
        --max-depth {params.max_coverage} \
        --max-idepth {params.max_coverage} \
        --annotate FORMAT/AD,FORMAT/DP,INFO/AD \
        {input.fname_bam} \
    | {params.bcftools} call \
        --threads {threads} \
        -Ou \
        -mv \
        --keep-alts \
    | {params.bcftools} norm \
        --threads {threads} \
        -Ou \
        -f {input.fname_ref} \
    | {params.bcftools} filter \
        --threads {threads} \
        -e 'TYPE="INDEL" & INFO/AD[1]<INFO/AD[0]' \
        -Ob \
        --output {output.fname_temp_bcf} \
        2> >(tee {log.errfile} >&2)
    #bcftools csq -f {input.fname_ref} -g wheretogetthis.gff3.gz in.vcf -Ob -o out.bcf

    {params.gunzip} -c {input.fname_cov} | tail -n +2 \
    | awk -v base={params.tsvbased} \
        '$3 < {params.mask_coverage_threshold} {{printf "%s\\t%d\\t%d\\n", $1, $2 - base, $2 - base + 1}}' \
    > {output.fname_mask_lowcoverage} \
        2> >(tee -a {log.errfile} >&2)

    # preparations
    {params.enhance_bcf} \
       {output.fname_temp_bcf} \
       {output.fname_bcf} \
       {params.ambiguous_base_coverage_threshold} \
       2> >(tee -a {log.errfile} >&2)

    {params.bcftools} index {output.fname_bcf} 2> >(tee -a {log.errfile} >&2)

    common_consensus_params="--fasta-ref {input.fname_ref} --mark-del - --mask {output.fname_mask_lowcoverage} --mask-with n"

    # majority bases
    {params.bcftools} consensus \
        --output {output.fname_fasta} \
        --chain {output.fname_chain} \
        $common_consensus_params \
        -H A \
        -i "INFO/AD[0]<INFO/AD[*]" \
        {output.fname_bcf} \
        2> >(tee -a {log.errfile} >&2)

    # ambiguous bases
    {params.bcftools} consensus \
        --output {output.fname_fasta_ambig} \
        --chain {output.fname_chain_ambig} \
        $common_consensus_params \
        -H I --iupac-codes \
        {output.fname_bcf} \
        2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    {params.EXTRACT_CONSENSUS} -i {input.BAM} -f {input.REF} -c {params.MIN_COVERAGE} -n {params.N_COVERAGE} -q {params.QUAL_THRD} -a {params.MIN_FREQ} -N "${{CONSENSUS_NAME}}" -o {params.OUTDIR} > {log.outfile} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    if tail -n +2 {input.REF_majority_dels} | grep -qE '[^n]'; then
        {params.MATCHER} -asequence {input.REF} -bsequence {input.REF_majority_dels} -outfile {output.REF_matcher} 2> >(tee {log.errfile} >&2)
    else
        touch {output.REF_matcher}
        echo "pure 'nnnn...' consensus, no possible alignement" | tee {log.outfile}
    fi
    """

shell:
    """
    {params.FRAMESHIFT_DEL_CHECKS} -i {input.BAM} -c {input.REF_majority_dels} -f {input.REF_NAME} -g {input.GENES_GFF} --english=true -o {output.FRAMESHIFT_DEL_CHECK_TSV} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    # 1. cleanup old run
    rm -f {output.SAM} {output.MSA}

    # 2. concatenate references
    mkdir -p {wildcards.dataset}/QA_alignments
    cat {input.patient_ref} {input.virusmix_ref} > {wildcards.dataset}/QA_alignments/bwa_refs_{wildcards.kind}.fasta

    # 3. indexing
    {params.BWA} index {wildcards.dataset}/QA_alignments/bwa_refs_{wildcards.kind}.fasta 2> >(tee {log.errfile} >&2)

    # 4. align
    {params.BWA} mem -t {threads} {wildcards.dataset}/QA_alignments/bwa_refs_{wildcards.kind}.fasta {input.FASTQ} > {output.SAM} 2> >(tee -a {log.errfile} >&2)

    # 5. MSA
    {params.MAFFT} --nuc --preservecase --maxiterate 1000 --localpair --thread {threads} {wildcards.dataset}/QA_alignments/bwa_refs_{wildcards.kind}.fasta > {output.MSA} 2> >(tee -a {log.errfile} >&2)

    # 6. cleanup BWA indices
    rm -f {wildcards.dataset}/QA_alignments/bwa_refs_{wildcards.kind}.fasta.*
    """

shell:
    """
    CONSENSUS_NAME={wildcards.dataset}
    CONSENSUS_NAME="${{CONSENSUS_NAME#*/}}"
    CONSENSUS_NAME="${{CONSENSUS_NAME//\//-}}"

    # 1. clean previous run
    rm -f {output}
    mkdir -p {wildcards.dataset}/QA_alignments

    # 2. collect coverage stats
    # we only collect statistics in the loop regions
    # of HIV-1 in order
    {params.COV_STATS} -t {params.TARGET} -i {input.BAM} -o {output} -m {input.MSA} --select "${{CONSENSUS_NAME}}" > {log.outfile} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    echo "Keep reject  -----------------------------------------------------"
    echo

    {params.SAMTOOLS} bam2fq -@ {threads} \
                             -F 2 \
                             -1 {output.reject_1} \
                             -2 {output.reject_2} \
                             {input.reject_aln} \
                             2> >(tee {log.errfile} >&2)
    echo
    """

shell:
    """
    echo "Filter out virus' reads  -----------------------------------------"
    echo

    {params.BWA} mem -t {threads} \
                     -o {output.tmp_aln} \
                     {input.global_ref} {input.fastq} \
                     > {log.outfile} 2> >(tee {log.errfile} >&2)

    echo
    echo "Keep reject  -----------------------------------------------------"
    echo

    {params.SAMTOOLS} bam2fq -@ {threads} \
                             -F 2 \
                             -1 {output.reject_1} \
                             -2 {output.reject_2} \
                             {output.tmp_aln} \
                             2> >(tee -a {log.errfile} >&2)
    echo
    """

shell:
    """
    curl --output "{output.host_ref}" "{params.host_ref_url}"
    """

shell:
    """
    # using zcat FILENAME.gz causes issues on Mac, see
    # https://serverfault.com/questions/570024/
    # redirection fixes this:
    unpack_rawreads() {{
        for fq in "${{@}}"; do
            zcat -f < "${{fq}}"
        done
    }}

    echo
    echo "Count aligned reads ---------------------------------------------"
    echo

    count=$({params.SAMTOOLS} view -@ {threads} -c -f {params.F} -F 2304 {input.host_aln} | tee {output.filter_count} 2> >(tee {log.errfile} >&2) )

    if (( count > 0 )); then
        echo
        echo "-----------------------------------------------------------------"
        echo "Needs special care: ${{count}} potential human reads found"
        echo "-----------------------------------------------------------------"
        echo
        echo "Removing identified host reads from raw reads -------------------"
        echo

        # get list
        {params.SAMTOOLS} view -@ {threads} \
                               -f {params.F} \
                               {input.host_aln} \
                               | cut -f 1 > {output.filter_list} 2> >(tee -a {log.errfile} >&2)

        unpack_rawreads {input.R1:q} \
               | {params.remove_reads_script} {output.filter_list} \
               | gzip \
               > {output.filtered_1} 2> >(tee -a {log.errfile} >&2) &

        unpack_rawreads {input.R2:q} \
               | {params.remove_reads_script} {output.filter_list} \
               | gzip \
               > {output.filtered_2} 2> >(tee -a {log.errfile} >&2) &

        wait

        if (( {params.keep_host} )); then
            # keep the rejects for further analysis

            echo
            echo "Keeping host-aligned virus' rejects ------------------------------"
            echo


            # (we compress reference-less, because the reference size is larger
            # than the contaminant reads)
            rm -f '{params.sort_tmp}'.[0-9]*.bam
            FMT=cram,no_ref,use_bzip2,use_lzma,level=9,seqs_per_slice=1000000
            {params.SAMTOOLS} view -@ {threads} \
                                  -h -f 2 \
                                  {input.host_aln} \
                | {params.SAMTOOLS} sort -@ {threads} \
                                -T {params.sort_tmp} \
                                --output-fmt ${{FMT}} \
                                -o {params.host_aln_cram} \
                                2> >(tee -a {log.errfile} >&2)

            echo
            echo "Compressing host-depleted raw reads ------------------------------"
            echo

            {params.SAMTOOLS} index -@ {threads} {params.host_aln_cram} 2> >(tee -a {log.errfile} >&2)
        fi
    else
        echo
        echo "No potential human reads found -----------------------------------"
        echo "Copy raw reads file"
        echo
        unpack_rawreads {input.R1} | gzip > {output.filtered_1} 2> >(tee -a {log.errfile} >&2) &
        unpack_rawreads {input.R2} | gzip > {output.filtered_2} 2> >(tee -a {log.errfile} >&2) &
        wait
        touch {output.filter_list}
    fi
    echo
    """

script:
    "../scripts/compute_diversity_measures.py"

script:
    "../scripts/aggregate_diversity.py"

shell:
    """
    {params.HAPLOCLIQUE} {params.EXTRA_PARAMETERS} {params.RELAX} {params.NO_SINGLETONS} {params.NO_PROB0} --limit_clique_size={params.CLIQUE_SIZE_LIMIT} --max_cliques={params.MAX_NUM_CLIQUES} --log={log.outfile} --bam {input} {params.OUTPREFIX} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    {params.COMPUTE_MDS} -q {params.INPREFIX} -s {params.REGION_START} -e {params.REGION_END} {params.USE_MSA} {params.MSA} -p {output.PDF} -o {params.TSV} > {log.output} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    # Convert BAM to FASTQ without re-reversing reads - SAVAGE expect all reads in the same direction
    source {params.FUNCTIONS}
    SamToFastq {params.PICARD} I={input} FASTQ={output.R1} SECOND_END_FASTQ={output.R2} RC=false 2> >(tee {log.errfile} >&2)
    # Remove /1 and /2 from the read names
    sed -i -e "s:/1$::" {output.R1}
    sed -i -e "s:/2$::" {output.R2}

    R1=${{PWD}}/{output.R1}
    R2=${{PWD}}/{output.R2}
    {params.SAVAGE} -t {threads} --split {params.SPLIT} -p1 ${{R1}} -p2 ${{R2}} -o {params.OUTDIR} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    # Convert BAM to FASTQ without re-reversing reads - SAVAGE expect all reads in the same direction
    source {params.FUNCTIONS}
    SamToFastq {params.PICARD} I={input} FASTQ={output.R1} 2> >(tee {log.errfile} >&2)

    R1=${{PWD}}/{output.R1}
    {params.SAVAGE} -t {threads} --split {params.SPLIT} -s ${{R1}} -o {params.OUTDIR} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.SAMTOOLS} sort -n {input.fname_bam} -o {output.fname_sam} 2> >(tee {log.errfile} >&2)

    mkdir -p {output.OUTPREFIX}
    {params.PREDICTHAPLO} \
        --sam {output.fname_sam} \
        --reference {input.fname_ref} \
        --prefix {output.OUTPREFIX}/ \
        --have_true_haplotypes 0 \
        --min_length {params.read_min_length} \
        2> >(tee -a {log.errfile} >&2)

    # TODO: copy over actual haplotypes
    touch {output.fname_out}
    """

shell:
    """
    {params.ALN2BASECNT} --first "{params.ARRAYBASED}" --basecnt "{output.BASECNT}" --coverage "{output.COVERAGE}" --name "{params.NAME}" --stats "{output.STATS}" "{input.BAM}" > {log.outfile} 2> >(tee {log.errfile} >&2)
    """

run:
    with open(output[0], "w") as out:
        out.write("\n".join(input))

shell:
    """
    {params.GATHER_COVERAGE} --output {output.COVERAGE} --stats {output.COVSTATS} --threads {threads} @{input.COVLIST} > >(tee {log.outfile}) 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    {params.MINORITY_CALLER} -r {input.REF} -c {params.MIN_COVERAGE} -N {params.NAMES} -t {threads} -o {params.OUTDIR} {params.FREQUENCIES} {input.BAM} > >(tee {log.outfile}) 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    echo "Trimming BAM with samtools"

    # samtools complains without that:
    rm -f '{params.sort_tmp}'.[0-9]*.bam

    {params.SAMTOOLS} ampliconclip -b {params.BED_PRIMERS} -f {output.stats} {input.BAM} |
        {params.SAMTOOLS}  sort -o {output.BAM} -T {params.sort_tmp} 2> >(tee {log.errfile} >&2)
    {params.SAMTOOLS}  index {output.BAM} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.script} {params.options} "{output.upload_prepared_touch}" "{params.sample_id}" "{wildcards.dataset}" {input:q}
    """

shell:
    """
    # using zcat FILENAME.gz causes issues on Mac, see
    # https://serverfault.com/questions/570024/
    # redirection fixes this:
    unpack_rawreads() {{
        for fq in "${{@}}"; do
            zcat -f < "${{fq}}"
        done
    }}

    echo "Compress un-filtered sequences -----------------------------------"
    echo

    {params.BWA} mem -t {threads} \
                     -C \
                     -o {output.cram_sam} \
                     {input.global_ref} \
                     <(unpack_rawreads {input.R1:q}) \
                     <(unpack_rawreads {input.R2:q}) \
                     > {log.outfile} 2> >(tee {log.errfile} >&2)

    # HACK handle incompatibilities between:
    #  - Illumina's 'bcl2fastq', which write arbitrary strings
    #  - 'bwa mem' which keep comments in the SAM file verbatim as in the FASTQ file
    #  - 'samtools' which expects comment to be properly marked as 'BC:Z:'
    #    as per SAM format specs
    REGEXP=\'s{{(?<=\\t)([[:digit:]]:[[:upper:]]:[[:digit:]]:([ATCGN]+(\+[ATCGN]+)?|[[:digit:]]+))$}}{{BC:Z:\\1}}\'
    FMT=cram,embed_ref,use_bzip2,use_lzma,level=9,seqs_per_slice=1000000

    rm -f '{params.sort_tmp}'.[0-9]*.bam
    perl -p -e ${{REGEXP}} {output.cram_sam} \
          | {params.SAMTOOLS} sort -@ {threads} \
                                   -T {params.sort_tmp} \
                                   -M \
                                   --reference {input.global_ref} \
                                   --output-fmt ${{FMT}} \
                                   -o {output.final_cram} \
                                   2> >(tee -a {log.errfile} >&2)

    {params.checksum_type}sum {output.final_cram} > {output.checksum} 2> >(tee -a {log.errfile} >&2)

    echo
    echo DONE -------------------------------------------------------------
    echo
    """

shell:
    """
    {params.checksum_type}sum {input} > {output}
    """

shell:
    """
    {params.GUNZIP} -c {input} > {output}
    """

shell:
    """
    echo "The length cutoff is: {params.LEN_CUTOFF}" > {log.outfile}

    {params.PRINSEQ} -fastq {input.R1} -fastq2 {input.R2} {params.EXTRA} -out_format 3 -out_good {wildcards.dataset}/preprocessed_data/R -out_bad null -min_len {params.LEN_CUTOFF} -log {log.outfile} 2> >(tee {log.errfile} >&2)

    # make sure that the lock held prinseq has been effectively released and propagated
    # on some networked shares this could otherwise lead to confusion or corruption
    if [[ "$OSTYPE" =~ ^linux ]]; then
        echo "Waiting for unlocks" >&2
        for U in {wildcards.dataset}/preprocessed_data/R_{{1,2}}.fastq; do
            flock -x -o ${{U}} -c "echo ${{U}} unlocked >&2"
        done
    fi

    mv {wildcards.dataset}/preprocessed_data/R{{_,}}1.fastq
    mv {wildcards.dataset}/preprocessed_data/R{{_,}}2.fastq
    rm -f {wildcards.dataset}/preprocessed_data/R_?_singletons.fastq

    gzip {wildcards.dataset}/preprocessed_data/R1.fastq
    gzip {wildcards.dataset}/preprocessed_data/R2.fastq
    """

shell:
    """
    echo "The length cutoff is: {params.LEN_CUTOFF}" > {log.outfile}

    {params.PRINSEQ} -fastq {input.R1} {params.EXTRA} -out_format 3 -out_good {wildcards.dataset}/preprocessed_data/R -out_bad null -min_len {params.LEN_CUTOFF} {params.EXTRA} -log {log.outfile} 2> >(tee {log.errfile} >&2)

    # make sure that the lock held prinseq has been effectively released and propagated
    # on some network shares this could otherwise lead to confusion or corruption
    if [[ "$OSTYPE" =~ ^linux ]]; then
        echo "Waiting for unlocks" >&2
        for U in {wildcards.dataset}/preprocessed_data/R_{{1,2}}.fastq; do
            flock -x -o ${{U}} -c "echo ${{U}} unlocked >&2"
        done
    fi

    mv {wildcards.dataset}/preprocessed_data/R{{,1}}.fastq

    gzip {wildcards.dataset}/preprocessed_data/R1.fastq
    """

shell:
    """
    {params.FASTQC} -o {params.OUTDIR} -t {threads} {params.NOGROUP} {input} 2> >(tee {log.errfile} >&2)
    """

shell:
    """
    TMPTSV=$(mktemp -t XXXXXXXX_cov.tsv)
    TMPINT=$(mktemp -t XXXXXXXX_int.tsv)
    {params.GUNZIP} -c {input.TSV} | cut -f'2-' > $TMPTSV
    mkdir -p "$(dirname "{output}")"
    {params.EXTRACT_COVERAGE_INTERVALS} -c "{params.COVERAGE}" -w "{params.WINDOW_LEN}" -s "{params.SHIFT}" -N "{params.NAME}" {params.LIBERAL} -b {params.ARRAYBASED} {params.OVERLAP} -t "{threads}" -o $TMPINT "{input.BAM}" > >(tee {log.outfile}) 2>&1
    cat $TMPINT >> "{log.outfile}"
    read name intervals < $TMPINT
    IFS=',' read -r -a interarray <<< "$intervals"
    printf "%s\n" "${{interarray[@]}}" > "{output}"
    rm $TMPTSV $TMPINT
    """

shell:
    """
    let "WINDOW_SHIFTS=({params.READ_LEN} * 4/5 + {params.SHIFT}) / {params.SHIFT}"
    let "WINDOW_LEN=WINDOW_SHIFTS * {params.SHIFT}"

    echo "Windows are shifted by: ${{WINDOW_SHIFTS}} bp" > {log.outfile}
    echo "The window length is: ${{WINDOW_LEN}} bp" >> {log.outfile}

    # Get absolute path for input files
    CWD=${{PWD}}
    BAM=${{PWD}}/{input.BAM}
    REF={input.REF}; [[ ${{REF}} =~ ^/ ]] || REF=${{PWD}}/${{REF}}
    OUTFILE=${{PWD}}/{log.outfile}
    ERRFILE=${{PWD}}/{log.errfile}
    WORK_DIR=${{PWD}}/{params.WORK_DIR}

    if [[ -n "{params.LOCALSCRATCH}" ]]; then
        # put input files in localscratch
        rsync -aq "${{BAM}}" "${{REF}}" "{params.LOCALSCRATCH}/"
        BAM="{params.LOCALSCRATCH}/${{BAM##*/}}"
        REF="{params.LOCALSCRATCH}/${{REF##*/}}"
    fi

    # Run ShoRAH in each of the predetermined regions (regions with sufficient coverage)
    LINE_COUNTER=0
    FILES=( )
    FILES_VCF=( )
    while read -r region || [[ -n ${{region}} ]]
    do
        echo "Running ShoRAH on region: ${{region}}" >> $OUTFILE
        (( ++LINE_COUNTER ))
        # Create directory for running ShoRAH in a corresponding region (if doesn't exist)
        DIR=${{WORK_DIR}}/REGION_${{LINE_COUNTER}}
        if [[ ! -d "${{DIR}}" ]]; then
            echo "Creating directory ${{DIR}}" >> $OUTFILE
            mkdir -p ${{DIR}}
        else
            # Results from previous runs
            if [[ {params.KEEP_FILES} == "true" ]]; then
                DIR_DST=${{WORK_DIR}}/old
                echo "Moving results from a previous run to ${{DIR_DST}}" >> $OUTFILE
                rm -rf ${{DIR_DST}}/REGION_${{LINE_COUNTER}}
                mkdir -p ${{DIR_DST}}
                mv -f ${{DIR}} ${{DIR_DST}}
                mkdir -p ${{DIR}}
            fi
        fi
        # Change to the directory where ShoRAH is to be executed
        if [[ -z "{params.LOCALSCRATCH}" ]]; then
            cd ${{DIR}}
        else
            # special case: go in local scratch instead
            rsync -aq "${{DIR}}" "{params.LOCALSCRATCH}/"
            mkdir -p "{params.LOCALSCRATCH}/REGION_${{LINE_COUNTER}}"
            cd "{params.LOCALSCRATCH}/REGION_${{LINE_COUNTER}}"
        fi

        # NOTE: Execution command for ShoRAH2 valid from v1.99.0 and above
        {params.SHORAH} -t {threads} -a {params.ALPHA} -w ${{WINDOW_LEN}} -x 100000 {params.IGNORE_INDELS} -p {params.POSTHRESH} -c {params.COVERAGE} -r ${{region}} -R 42 -b ${{BAM}} -f ${{REF}} >> $OUTFILE 2> >(tee -a $ERRFILE >&2)
        if [[ -n "{params.LOCALSCRATCH}" ]]; then
            # copyback from localscratch
            rsync -auq "{params.LOCALSCRATCH}/REGION_${{LINE_COUNTER}}" "${{WORK_DIR}}"
            cd ${{DIR}}
        fi
        if [[ -s ${{DIR}}/reads.fas && -f ${{DIR}}/snv/SNVs_0.010000_final.csv ]]; then
            # Non empty reads: run had enough data and should have produced SNVs
            {params.BCFTOOLS} view ${{DIR}}/snv/SNVs_0.010000_final.vcf -Oz -o ${{DIR}}/snv/SNVs_0.010000_final.vcf.gz
            {params.BCFTOOLS} index ${{DIR}}/snv/SNVs_0.010000_final.vcf.gz
            FILES+=("${{DIR}}/snv/SNVs_0.010000_final.csv")
            FILES_VCF+=("${{DIR}}/snv/SNVs_0.010000_final.vcf.gz")
        elif (( {params.COVINT} == 0 && LINE_COUNTER == 1 )) && [[ -f ${{DIR}}/reads.fas && ( ! -s ${{DIR}}/reads.fas ) ]]; then
            # if we have disabled coverage intervales entirely, the first and only line might have no reads
            # (e.g.: in negative controls )

            echo "No reads while coverage intervals disabled (possible negative control sample)" 2> >(tee -a $ERRFILE >&2)
            cd ${{CWD}}
            (( --LINE_COUNTER )) || true # Strict mode : (( 0 )) = fail
            break
        else
            echo "ERROR: unsuccesful execution of ShoRAH" 2> >(tee -a $ERRFILE >&2)
            exit 1
        fi

        # Change back to working directory
        cd ${{CWD}}
    done < {input.TSV}

    # Aggregate results from different regions
    if (( ${{#FILES[@]}} )); then
        echo "Intermediate csv files: ${{FILES[*]}}" >> {log.outfile}
        echo "Intermediate vcf files: ${{FILES_VCF[*]}}" >> {log.outfile}
        (head -n 1 "${{FILES[0]}}"; tail -q -n +2 "${{FILES[@]}}" | sort -t, -nk2) > {output.CSV}
        {params.BCFTOOLS} concat -o ${{WORK_DIR}}/snvs_tmp.vcf "${{FILES_VCF[@]}}"
        {params.BCFTOOLS} sort ${{WORK_DIR}}/snvs_tmp.vcf  -o {output.VCF}
        rm -f ${{WORK_DIR}}/snvs_tmp.vcf
    elif (( LINE_COUNTER )); then
        echo "ERROR: unsuccesful execution of ShoRAH" 2> >(tee -a {log.errfile} >&2)
        exit 1
    else
        echo "No alignment region reports sufficient coverage" >> {log.outfile}
        touch {output.CSV}
        touch {output.VCF}
    fi
    """

shell:
    """
    {params.SAMTOOLS} faidx {input} -o {output} > {log.outfile} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    # Add qualities to indels
    {params.LOFREQ} indelqual --dindel -f {input.REF} -o {output.BAM} --verbose {input.BAM} > {log.outfile} 2> >(tee {log.errfile} >&2)
    # Index bam file
    {params.SAMTOOLS} index {output.BAM} 2> >(tee {log.errfile} >&2)

    # Run Lofreq
    echo "Running LoFreq" >> {log.outfile}
    {params.LOFREQ} call {params.EXTRA} --call-indels -f {input.REF} -o {output.SNVs} --verbose {output.BAM} >> {log.outfile} 2> >(tee -a {log.errfile} >&2)
    """

shell:
    """
    {params.EXTRACT_COVERAGE_INTERVALS} -b {params.ARRAYBASED} -cf {input.TSV} -c {params.COVERAGE} --no-shorah -N {params.NAMES} -o {output} {input.BAM} > >(tee {log.outfile}) 2>&1
    """

shell:
    """
    SAMPLE_ID="{params.ID}"

    # Number of input reads
    LINECOUNT=$( cat {input.R1} | wc -l )
    let "INPUT=LINECOUNT / {params.FACTOR}"

    # Number of reads after QC
    # For portability reason not using zcat
    {params.GUNZIP} -c {input.R1_QC} > {params.R1_temp}
    LINECOUNT=$( cat {params.R1_temp} | wc -l )
    let "READCOUNT=LINECOUNT / {params.FACTOR}"
    rm {params.R1_temp}

    # Number of aligned reads
    ALNCOUNT=$( {params.SAMTOOLS} view {input.BAM} | wc -l )

    echo -e "${{SAMPLE_ID}}\t${{INPUT}}\t${{READCOUNT}}\t${{ALNCOUNT}}" > {output}
    """

shell:
    """
    cat {input} > stats/temp
    echo -e "ID\tInput\tQC\tAlignments" | cat - stats/temp > {output}
    rm stats/temp
    """