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

# Build an index optimized
# for small reference genomes
genomeSize=$(grep -v '^>' {input.fa} | awk '{{sum+=length($0)}}; END {{print sum}}')
# Calculate min(14, log2(GenomeSize)/2 - 1)
genomeSAindexNbases=$(python -c "import math; print(min(14,int((math.log($genomeSize,2)/2)-1)))")
echo "Building index with --genomeSAindexNbases $genomeSAindexNbases"
STAR \\
    --runThreadN {threads} \\
    --runMode genomeGenerate \\
    --genomeSAindexNbases ${{genomeSAindexNbases}} \\
    --genomeDir STAR/2.7.6a/genome \\
    --genomeFastaFiles {input.fa} \\
    --outFileNamePrefix STAR/2.7.6a/build_genome_ \\
    --outTmpDir ${{tmp}}/tmp_genome
"""

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

STAR \\
    --runThreadN {threads} \\
    --runMode genomeGenerate \\
    --genomeDir STAR/2.7.6a/genome \\
    --genomeFastaFiles {input.fa} \\
    --outFileNamePrefix STAR/2.7.6a/build_genome_ \\
    --outTmpDir ${{tmp}}/tmp_genome
"""

shell: """
python3 {params.create_rRNA} \\
    {input.fa} \\
    {input.gtf} \\
    {params.genome} > {params.genome}.rRNA_interval_list
"""

shell: """
python3 {params.get_karyoplot} {input.gtf} > karyoplot_gene_coordinates.txt
"""

shell: """
mkdir -p karyobeds && cd karyobeds
python3 {params.get_karyoplot} {input.gtf}
"""

shell: """
python {params.gtf2protein} {input.gtf} > protein_coding_genes.lst
gtfToGenePred -ignoreGroupsWithoutExons -genePredExt {input.gtf} genes.gtf.genePred
genePredToBed genes.gtf.genePred genes.gtf.genePred.bed

awk -F '\\t' -v OFS='\\t' '{{print $12,$1}}' genes.gtf.genePred \\
    | sort -k1,1n > gene2transcripts

while read gene; do
    grep "${{gene}}" gene2transcripts;
done < protein_coding_genes.lst > gene2transcripts.protein_coding_only

python {params.gene2transcripts} \\
    gene2transcripts.protein_coding_only \\
    genes.gtf.genePred.bed \\
    > gene2transcripts.protein_coding_only.with_len

sort -k1,1 -k3,3nr gene2transcripts.protein_coding_only.with_len | \\
    awk -F '\\t' '{{if (!seen[$1]) {{seen[$1]++; print $2}}}}' > protein_coding_only.txt

while read transcript; do
    grep -m1 "${{transcript}}" genes.gtf.genePred.bed;
done < <(awk -F '.' '{{print $1}}' protein_coding_only.txt) > transcripts.protein_coding_only.bed12

rm -f "protein_coding_genes.lst" "genes.gtf.genePred" "genes.gtf.genePred.bed" \\
    "gene2transcripts" "gene2transcripts.protein_coding_only" "protein_coding_only.txt" \\
    "gene2transcripts.protein_coding_only.with_len" "protein_coding_only.txt"
"""

shell: """
python3 {params.generate_qualimap} \\
    -g {input.gtf} \\
    -f {input.fa} \\
    -o {output} \\
    --ignore-strange-chrom 2> qualimap_error.log
"""

shell: """
wget https://hpc.nih.gov/~OpenOmics/common/{params.tarname} -O {params.tarfile}
tar vxf {params.tarfile} -C {params.outdir} && rm {params.tarfile}
"""

shell: """
wget https://hpc.nih.gov/~OpenOmics/common/{params.tarname} -O {params.tarfile}
tar vxf {params.tarfile} -C {params.outdir} && rm {params.tarfile}
"""

shell: """
mkdir -p "{params.outdir}"
wget https://hpc.nih.gov/~OpenOmics/common/{params.conf1} -O {output.conf1}
wget https://hpc.nih.gov/~OpenOmics/common/{params.conf2} -O {output.conf2}
sed -i 's@/data/OpenOmics/references/common@{params.new}@g' {output.conf1}
sed -i 's@/data/OpenOmics/references/common@{params.new}@g' {output.conf2}
"""

shell: """
wget https://hpc.nih.gov/~OpenOmics/common/{params.tarfile} -O {params.outfh}
tar vxf {params.outfh} -C {params.outdir} && rm {params.outfh}
"""

shell: """
python {params.get_flowcell_lanes} {input.R1} {wildcards.name} > {output.fqinfo}
"""

shell: """
python {params.pythonscript} {params.annotate} {params.inputdir} {params.inputdir}
sed 's/\\t/|/1' {output.gene_counts_matrix} | \
    sed '1 s/^gene_id|GeneName/symbol/' > {output.reformatted}
"""

shell: """
infer_experiment.py -r {params.bedref} -i {input.file1} -s 1000000 > {output.out1}
read_distribution.py -i {input.file1} -r {params.bedref} > {output.out2}
"""

shell: """
# tin.py writes to current working directory
cd {params.outdir}
tin.py -i {input.bam} -r {params.bedref}
"""

shell: """
python {params.create_matrix} {input.tins} > {output.matrix}
"""

run:
    with open(output.out1, "w") as out:
        out.write("sampleName\tfileName\tcondition\tlabel\n")
        i=0
        for f in input.files:
            out.write("{}\t".format(params.allsamples[i]))
            out.write("{}/{}.star.count.txt\t".format(params.pathprefix, params.allsamples[i]))
            out.write("{}\t".format(params.groups[i]))
            out.write("{}\n".format(params.labels[i]))
            i=i+1
        out.close()

shell: """
cd {params.outdir}
module load {params.rver}
Rscript {params.rscript} '{params.outdir}' '{input.files}' '{params.annotate}' '{input.sampletable}'
"""

shell: """
cd {params.outdir}
module load {params.rver}
Rscript {params.rscript1} '{params.outdir}' '{output.outhtml}' \
'{input.file1}' '{input.file2}' '{params.projectId}' '{params.projDesc}' '{params.rscript2}'
"""

shell: """
mkdir -p {params.outdir}
fastQValidator --noeof --file {input.R1} > {output.out1}
fastQValidator --noeof --file {input.R2} > {output.out2}
"""

shell: """
fastqc {input.R1} {input.R2} -t {threads} -o {params.outdir};
"""

shell: """
# Get encoding of Phred Quality Scores
encoding=$(python {params.encoding} {input.R1})
echo "Detected Phred+${{encoding}} ASCII encoding"

bbtools bbmerge-auto in1={input.R1} in2={input.R2} qin=${{encoding}} \
    ihist={output} k=62 extend2=200 rem ecct -Xmx{params.memory}G
"""

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

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield}  \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --twopassMode Basic \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

# SAMtools sort (uses less memory than STAR SortedByCoordinate)
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam {params.prefix}.Aligned.out.bam \
    > {output.out1}

rm {params.prefix}.Aligned.out.bam
mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

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

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield} \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --sjdbGTFfile {params.gtffile} \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}
"""

shell: """
cat {input.files} | \
    sort | \
    uniq | \
    awk -F \"\\t\" '{{if ($5>0 && $6==1) {{print}}}}'| \
    cut -f1-4 | sort | uniq | \
grep \"^chr\" | grep -v \"^chrM\" > {output.out1}
"""

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

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield}  \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --sjdbFileChrStartEnd {input.tab} \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

# SAMtools sort (uses less memory than STAR SortedByCoordinate)
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam {params.prefix}.Aligned.out.bam \
    > {output.out1}

rm {params.prefix}.Aligned.out.bam
mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

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

# Get strandedness to calculate Forward Probability
fp=$(tail -n1 {input.file2} | awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "1.0"; else print "0.5";}}')

echo "Forward Probability Passed to RSEM: $fp"
rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  \
    --bam --paired-end -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time \
    --temporary-folder ${{tmp}} --keep-intermediate-files --forward-prob=${{fp}} --estimate-rspd
"""

shell: """
inner_distance.py -i {input.bam} -r {params.genemodel} \
    -k 10000000 -o {params.prefix}
"""

shell: """
bash {params.bashscript} {input.bam} {params.outprefix}

# reverse files if method is not dUTP/NSR/NNSR ... ie, R1 in the direction of RNA strand.
strandinfo=`tail -n1 {input.strandinfo} | awk '{{print $NF}}'`
if [ `echo "$strandinfo < 0.25"|bc` -eq 1 ];then
    mv {output.fbw} {output.fbw}.tmp
    mv {output.rbw} {output.fbw}
    mv {output.fbw}.tmp {output.rbw}
fi
"""

shell: """
# Generate RNA QC Dashboard
{params.rwrapper} \
    -m {params.rmarkdown} \
    -r {input.counts} \
    -t {input.tins} \
    -q {input.qc} \
    -o {params.odir} \
    -f RNA_Report.html
"""

shell:"""
set +u
if [ "{params.mode}" == "uge" ]; then 
    module load R
else
    module load R/4.1.3
fi
Rscript {params.rmerger} -tsvs {input} {params.ervs_fam_table} -dir {params.outdir}
"""