shell: """
# Find differential alternative
# splicing events such as skipped
# exons, alternative 5'/3' splice 
# sites, and retained introns 
flair diffSplice \\
    --threads {threads} \\
    --test \\
    --isoforms {input.isoforms} \\
    --counts_matrix {input.counts} \\
    --conditionA {params.group1} \\
    --conditionB {params.group2} \\
    --out_dir {params.outdir} \\
    --out_dir_force
"""

shell: """
# Subset counts matrix to contain 
# only samples belonging to the 
# contrast of interest, this will
# ensure the correct group is set
# as the baseline for design matrix.
# By default, the first group it 
# encounters in the sample names 
# will be the intercept. 
# Get the first group's samples
grp1_samples=$(
    awk -F '\\t' -v OFS='_' \\
        '$2=="{params.group1}" {{print $1,$2,$3}}' \\
        {input.manifest}
)
# Get the second group's samples
grp2_samples=$(
    awk -F '\\t' -v OFS='_' \\
        '$2=="{params.group2}" {{print $1,$2,$3}}' \\
        {input.manifest}
)
# Create per-contrast counts matrix,
# extracts samples by their col name
{params.script} -c ids ${{grp1_samples}} ${{grp2_samples}} \\
    -i {input.counts} \\
> {output.counts}

# Find differential gene/isoform 
# expression and find differential
# isoform usage
flair diffExp \\
    --threads {threads} \\
    --counts_matrix {output.counts} \\
    --out_dir {params.outdir} \\
    --exp_thresh 10 \\
    --out_dir_force
"""

shell: 
    """
    # Download primary assembly,
    # reference genome sequence
    wget -O - \\
        {params.uri} \\
    | gzip -d \\
    > {output.genome}
    """

shell: 
    """
    # Download transcript sequences
    wget -O - \\
        {params.uri} \\
    | gzip -d \\
    > {output.transcripts}
    """

shell: 
    """
    # Download PRI gene annotation
    wget -O - \\
        {params.uri} \\
    | gzip -d \\
    > {output.gtf}
    """

shell: 
    """
    # Convert SAM to raw features TSV
    samtools view \\
        -@{threads} \\
        -h \\
        -F 4 \\
        {input.bam} \\
    | java -jar ${{PICARDJARPATH}}/sam2tsv.jar -r {input.ref} \\
    | awk 'BEGIN{{FS=OFS="\\t"}} ($9 != "S") && ($9 != "H") && ($9 != "N")' - \\
    | awk 'BEGIN{{FS=OFS="\\t"}} \\
        ($7=="."){{$7="-99";}} \\
        ($4=="."){{$4="-99"}} \\
        ($5=="."){{$5="na"}} \\
        ($8=="."){{$8="na"}} \\
        ($9=="D"){{$6=" "}} \\
        ($2==16){{$2="n"}} \\
        ($2==0){{$2="p"}} 1' \\
    > {output.tsv}
    """

shell: 
    """
    # Setup for the combine step,
    # see github repo for more info:
    # https://github.com/darelab2014/Dinopore/blob/main/code/s2.Combine_raw_nnpl.sh
    cd "{params.outdir}"
    head -n 1 {input.events} \\
        > {params.header}
    noline=$(wc -l {input.events} | cut -d " " -f 1)
    tmp="{params.tmp}"
    num=2
    chunk=10000000
    num1=$(expr $num + $chunk)
    num2=$(expr $num1 + 1)
    filecount=1

    # Combine event-level signal information
    while [ $num -le $noline ]; do
        sed -n "$num,${{num1}}p; ${{num2}}q" {input.events} > "$tmp"
        chk=$(wc -l "$tmp" | cut -f 1 -d " ")
        echo "Processing reads $num to $num1 in file $filecount"
        if [ $chk -gt 0 ]; then
            Rscript ${{DINOPORE_CODE}}/s2.Combine_raw_nanopolish.R \\
                -f "$tmp" \\
                -t {threads} \\
                -o {output.events}.part${{filecount}} \\
                -s ${{noline}} \\
                -n ${{num}} \\
                -c ${{chunk}}
            if test -f tmp.eli; then
                eli=$(cat tmp.eli)
                rm "tmp.eli"
            else
                eli=0
            fi
        fi
        num=$(( $num1 - $eli + 1 ))
        num1=$(expr $num + $chunk)
        num2=$(expr $num1 + 1)
        filecount=$(expr $filecount + 1)
    done

    # Combine results of all parts
    cat ${{DINOPORE_CODE}}/misc/nnpl.header > {output.events}
    cat {output.events}.part* | grep -v contig >> {output.events}
    rm -f tmp.nnpl {output.events}.part*
    """

shell: 
    """
    # Setups temporary directory for
    # intermediate files with built-in 
    # mechanism for deletion on exit
    if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
    tmp=$(mktemp -d -p "{params.tmpdir}")
    trap 'rm -rf "${{tmp}}"' EXIT

    # Setup for the feature table step,
    # see github repo for more info:
    # https://github.com/darelab2014/Dinopore/blob/main/code/S3.Generate_raw_features.sh
    cd "{params.outdir}"
    num=1
    chunk=100000
    num1=$(expr $num + $chunk - 1)
    num2=$(expr $num1 + 1)
    filecount=1

    # Extract common reads between filtered
    # alignment TSV and nanopolish events 
    awk -F '\\t' 'NR!=1 {{print $1}}' {input.tsv} \\
        | LC_ALL=C uniq \\
        | LC_ALL=C sort -S {params.memory} --parallel={threads}  \\
    > {params.tsv}
    awk '{{print $2}}' {input.events} \\
        | LC_ALL=C uniq \\
        | LC_ALL=C sort -S {params.memory} --parallel={threads}  \\
    > {params.events}
    # Prints lines/reads common in both files
    comm -12 \\
        {params.tsv} \\
        {params.events} \\
    > {params.common}
    noline=$(wc -l {params.common} | cut -d " " -f 1)

    # Generate raw features files
    while [ $num -le $noline ]; do
        echo "Processing reads $num to $num1 in file $filecount"
        sed -n "$num,${{num1}}p; ${{num2}}q" {params.common} \\
        > {params.tmp}
        chk=$(wc -l {params.tmp} | cut -f 1 -d " ")
        if [ $chk -gt 0 ]; then
            # The headers can also be found here:
            # https://github.com/darelab2014/Dinopore/tree/main/code/misc
            cat ${{DINOPORE_CODE}}/misc/nnpl.header \\
            > {params.innnpl}
            cat ${{DINOPORE_CODE}}/misc/tsv.header \\
            > {params.intsv}
            LC_ALL=C grep \\
                -h \\
                -F \\
                -w \\
                -f {params.tmp} \\
                {input.tsv} \\
            >> {params.intsv}
            LC_ALL=C grep \\
                -h \\
                -F \\
                -w \\
                -f {params.tmp} \\
                {input.events} \\
            >> {params.innnpl}

            # Create chunked raw features files,
            # Rscript cannot use files with abs paths
            rel_intsv="$(basename {params.intsv})"
            rel_innnpl="$(basename {params.innnpl})"
            Rscript ${{DINOPORE_CODE}}/s3.Generate_raw_feature_table.R \\
                -n ${{rel_innnpl}} \\
                -t ${{rel_intsv}} \\
                -o {params.outname}.part${{filecount}}

            num=$num2
            num1=$(expr $num + $chunk - 1)
            num2=$(expr $num1 + 1)
            filecount=$(expr $filecount + 1)
            rm {params.tmp}
        fi
    done

    # Combine results of all parts
    cat ${{DINOPORE_CODE}}/misc/inae.header \\
    > {output.table}
    cat {params.outname}.part*.positive \\
        | LC_ALL=C sort -T "${{tmp}}" -k1,1 -k3,3n \\
        | grep -v contig \\
    >> {output.table}
    cat {params.outname}.part*.negative \\
        | LC_ALL=C sort -T "${{tmp}}" -k1,1 -k3,3n \\
        | grep -v contig \\
    >> {output.table}

    # Cleanup intermediate files
    rm -f {params.tsv} \\
        {params.events} \\
        {params.common} \\
        {params.intsv} \\
        {params.innnpl} \\
        {params.outname}.part*.positive \\
        {params.outname}.part*.negative
    """

shell: 
    """
    # Aggregate features of all
    # reads at the same position
    # across multiple samples in 
    # the same group
    cd "{params.outdir}"
    Rscript ${{DINOPORE_CODE}}/s4.Aggregating_reads_pos.R \\
        -t {threads} \\
        -o {output.table} \\
        -r _{params.group} 
    """

shell: 
    """
    # Transform/pre-process features
    # file to create input for CNN,
    # regression model
    cd "{params.outdir}"
    Rscript ${{DINOPORE_CODE}}/s5.Preprocess_data_matrix_inputCNN.R \\
        -t {threads} \\
        -i {input.table} \\
        -o {output.rdata} \\
        -c {params.truth}
    """

shell: 
    """
    # Predict A-to-I RNA-editing
    # Class 0: A (unmodified)
    # Class 1: I (editied)
    # Class 2: G (SNP)
    cd "{params.outdir}"
    DINOPORE_CODE="{params.code}"
    Rscript ${{DINOPORE_CODE}}/s6.Predict_test_data.R \\
        -t {threads} \\
        -i {input.rdata}
    """

shell: 
    """
    # Setups temporary directory for
    # intermediate files with built-in 
    # mechanism for deletion on exit
    if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
    tmp=$(mktemp -d -p "{params.tmpdir}")
    trap 'rm -rf "${{tmp}}"' EXIT

    # Align against reference genome,
    # minimap2 automatically handles
    # conversion of U to T bps, See 
    # this issue for the author of 
    # minimap2's recomendations for 
    # aligning direct RNA reads against 
    # the transcriptome:
    # https://github.com/lh3/minimap2/issues/340  
    minimap2 \\
        -ax map-ont \\
        -N 10 \\
        -k10 \\
        {input.ref} \\
        {input.fq} \\
    | samtools sort -@{threads} \\
        -T "${{tmp}}" \\
        -O bam \\
        --write-index \\
        -o {output.bam}##idx##{output.bai} \\
        -
    """

shell: """
# Build an index to map basecalled
# reads to ONT device's raw signal 
nanopolish index \\
    --directory={params.f5} \\
    {input.fq}
"""

shell: """
# Estimate polyA tail lengths
# for each transcript
nanopolish polya \\
    --threads={threads} \\
    --reads={input.fq} \\
    --bam={input.bam} \\
    --genome={input.ref} \\
> {output.polya}
"""

shell: """
pycoQC \\
    -f {input.tsv} \\
    -o {output.html} \\
    -j {output.json} 
"""

shell: """
fastqc \\
    -t {threads} \\
    -o {params.outdir} \\
    {input.fq}
"""

shell: """
NanoCount \\
    -i {input.bam} \\
    --extra_tx_info \\
    -e {params.em_iter} \\
    -o {output.counts} 
"""

shell: """
# Estimated counts matrix,
# raw counts 
{params.script} \\
    --input {input.counts} \\
    --output {output.est_count} \\
    --join-on transcript_name \\
    --extract est_count \\
    --clean-suffix '.nanocount.transcripts.tsv' \\
    --nan-values 0.0
# TPM counts matrix,
# normalized counts
{params.script} \\
    --input {input.counts} \\
    --output {output.tpm} \\
    --join-on transcript_name \\
    --extract tpm \\
    --clean-suffix '.nanocount.transcripts.tsv' \\
    --nan-values 0.0
"""

shell: """
# Convert BAM into BED12
bam2Bed12 \\
    -i {input.bam} \\
> {output.bed12}
# Correct misaligned splice
# sites with flair correct 
flair correct \\
    -q {output.bed12} \\
    -g {input.ref} \\
    -f {input.gtf} \\
    -t {threads} \\
    -o {params.prefix} \\
    --nvrna
"""

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Merge the corrected BED files
cat {input.corrected} > {output.merged}

# Find high-confidence isoforms.
# Flair does not like transcript 
# IDs reported by GENCODE  
flair collapse \\
    --temp_dir "${{tmp}}" \\
    --threads {threads} \\
    --gtf {input.gtf} \\
    --genome {input.genome} \\
    --annotation_reliant generate \\
    --check_splice \\
    --stringent \\
    --support 1 \\
    --reads {input.reads} \\
    --query {output.merged} \\
    --output {params.prefix}
"""

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Create sample manifest file,
# TSV containing each samples 
# basename, group, batch, and 
# path to its FastQ file, the
# sample manifest file, cannot
# contain any underscores in the 
# first three columns; however,
# hypens cannot/should not used
# in the groups/condition column
awk -F '\\t' -v OFS='\\t' \\
    '{{print \\
        $1, \\
        $3, \\
        "batch1", \\
        "{params.workdir}/"$1"/fastqs/"$1".filtered.fastq.gz" \\
    }}' {params.groups} \\
    | awk -F '\\t' -v OFS='\\t' \\
        '{{gsub("_","-",$1); \\
            gsub("_","",$2); \\
            gsub("_","-",$3); \\
            print \\
        }}' \\
>  {output.manifest}

# Create isoform-by-sample counts 
# matrix that can be used in the 
# diffExp and diffSplice modules
flair quantify \\
    -r {output.manifest} \\
    -i {input.fa} \\
    --isoform_bed {input.bed} \\
    --check_splice \\
    --stringent \\
    --temp_dir "${{tmp}}" \\
    --threads {threads} \\
    --tpm \\
    --output {params.prefix}
"""

shell: """
# Structural and quality annotation
# of collapsed transcript from flair
sqanti3_qc.py \\
    {input.fa} \\
    {input.gtf} \\
    {input.genome} \\
    --aligner_choice=minimap2 \\
    --cpus {threads} \\
    --output {params.prefix} \\
    --dir {params.outdir} \\
    --CAGE_peak {params.cage_peak} \\
    --polyA_motif_list {params.polya_motif} \\
    --isoAnnotLite \\
    --force_id_ignore \\
    --fasta
"""

shell: """
# Applies ML filter on selected
# SQANTI3 QC classification file
sqanti3_filter.py ML \\
    --output {params.prefix} \\
    --dir {params.outdir} \\
    --gtf {input.gtf} \\
    --isoforms {input.fa} \\
    {input.txt}
"""

shell: 
    """
    # Setups temporary directory for
    # intermediate files with built-in 
    # mechanism for deletion on exit
    if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
    tmp=$(mktemp -d -p "{params.tmpdir}")
    trap 'rm -rf "${{tmp}}"' EXIT

    # Align against reference genome,
    # minimap2 automatically handles
    # conversion of U to T bps, See 
    # this issue for the author of 
    # minimap2's recomendations for 
    # aligning direct RNA reads against 
    # the transcriptome:
    # https://github.com/lh3/minimap2/issues/340  
    minimap2 \\
        -ax map-ont \\
        -N 10 \\
        -k10 \\
        {input.ref} \\
        {input.fq} \\
    | samtools sort -@{threads} \\
        -T "${{tmp}}" \\
        -O bam \\
        --write-index \\
        -o {output.bam}##idx##{output.bai} \\
        -
    """

shell: 
    """
    {params.prefix} {input.fq} {params.suffix} {output.fq}
    """