script: "../scripts/bam_to_counts.R"

shell:
     """
     echo "{wildcards.SRR_ID} {wildcards.BASE_FOLDER}/mappings/stringtie/{wildcards.SRR_ID}/{wildcards.SRR_ID}_assembled.gtf" > input.temp

     mv input.temp {wildcards.BASE_FOLDER}/mappings/stringtie/{wildcards.SRR_ID}/{wildcards.SRR_ID}_input.temp

     {input.binary} -i {wildcards.BASE_FOLDER}/mappings/stringtie/{wildcards.SRR_ID}/{wildcards.SRR_ID}_input.temp -g {output.gmat} -t {output.tmat}

     rm -f {wildcards.BASE_FOLDER}/mappings/stringtie/{wildcards.SRR_ID}/{wildcards.SRR_ID}_input.temp
     """

script: "../scripts/txi_to_counts.R"

script: "../scripts/txi_to_counts.R"

run:
    generate_deseq_input(wildcards)

shell:
     """
     {input.fastqc} -t {threads} -o {wildcards.BASE_FOLDER}/quality/{wildcards.SRR_ID} {input.f1} {input.f2}
     """

shell:
     """
     {input.fastqc} -t {threads} -o {wildcards.BASE_FOLDER}/quality/{wildcards.SRR_ID} {input.f1}
     """

shell:
     """
     {input.fastqc} -t {threads} -o {wildcards.BASE_FOLDER}/quality/{wildcards.SRR_ID} {input.f1} {input.f2}
     """

shell:
     """
     {input.fastqc} -t {threads} -o {wildcards.BASE_FOLDER}/quality/{wildcards.SRR_ID} {input.f1}
     """

shell:
     "{params.fasterq} "  # Binary
     "{input} "  # Input SRA File
     "--split-files "  # Split files according to the reads
     "--progress "  # Show progress
     "--skip-technical "  # Skip technical reads
     "--threads {threads} "  # Number of threads
     "--outdir {wildcards.BASE_FOLDER}/fastq"

shell:
     "{params.fasterq} "  # Binary
     "{input} "  # Input SRA File
     "--split-files "  # Split files according to the reads
     "--progress "  # Show progress
     "--skip-technical "  # Skip technical reads
     "--threads {threads} "  # Number of threads
     "--outdir {wildcards.BASE_FOLDER}/fastq"

    shell:
         "{input.salmon} index "  # Salmon index mode
         "-t {input.decoy_index} "  # Decoy transcript
         "-i {wildcards.BASE}/index/salmon "  # Path to index folder
         "--decoys {input.decoy_keys} "  # Decoy keys
         "-k 31 "  # K-mer (31) is recommended by the authors
         "-p {threads}"  # Number of threads


def get_salmon_reads(files):
    if isinstance(files, str):
        return f"-r {files}"

shell:
     """
     rm -rf {params.index_dir}/kvdb

     {params.sortmerna} {params.ref_str} {params.read_str} -workdir {params.index_dir} -fastx -other -a {threads} -num_alignments 1 -paired_in --out2 

     mkdir -p {BASE}/filtered/ 

     mkdir -p {BASE}/logs

     mv {params.index_dir}/out/aligned.log {BASE}/logs/{wildcards.SRR_ID}_filtering.log

     mv {params.index_dir}/out/other_fwd.fastq {BASE}/filtered/{wildcards.SRR_ID}.sra.filtered_1.fastq
     mv {params.index_dir}/out/other_rev.fastq {BASE}/filtered/{wildcards.SRR_ID}.sra.filtered_2.fastq

     rm  -f {params.index_dir}/out/aligned_fwd.fastq
     rm  -f {params.index_dir}/out/aligned_rev.fastq
     """

shell:
     """
     rm -rf {params.index_dir}/kvdb

     {params.sortmerna} {params.ref_str}  -reads {input.reads} -workdir {params.index_dir} -fastx -other -a {threads} -num_alignments 1 -paired_in 

     mkdir -p {BASE}/filtered/ 

     mkdir -p {BASE}/logs

     mv {params.index_dir}/out/aligned.log {BASE}/logs/{wildcards.SRR_ID}_filtering.log

     mv {params.index_dir}/out/other.fastq {BASE}/filtered/{wildcards.SRR_ID}.sra.filtered.fastq 

     rm  -f {params.index_dir}/out/aligned.fastq
     """

shell:
     """
     {input.star} --runMode genomeGenerate --runThreadN {threads} --genomeDir {params.folder} --genomeFastaFiles {input.genome} --sjdbGTFfile {input.gtf}

    if [ -f Log.out ]; then
        mkdir -p {BASE}/logs
        mv Log.out {BASE}/logs/STAR.index.log
    fi 
     """

shell:
     """
     {input.star} --runThreadN {threads} --genomeDir {BASE}/index/star --outFileNamePrefix {params.prefix} --outSAMtype BAM SortedByCoordinate --readFilesIn {input.files}

     rm -f {BASE}/bams/{params.prefix}.Log.progress.out

     """

    shell:
         """
         {input.bowtie2} {input.transcripts} {input.transcripts}
         """

"""
Get statistics from paired-end
Bowtie2:
-p : number of threads
-q : input files are in fastq format
--no-unal: suppress SAM records for unaligned reads
-x: Index files prefix (without X.bt2 extension)

"""

shell:
     """
     {input.bowtie} -p {threads} -q --no-unal -x {input.transcripts} -1 {input.files[0]} -2 {input.files[1]} > {output}
     """

start_time <- Sys.time()

suppressPackageStartupMessages(library("DESeq2"))

# We can relevel our conditions initially so that we do not have to worry
# about them later on
conds <- factor(snakemake@params$conditions)
conds <- relevel(conds, snakemake@params$reference)

colData <- data.frame(sample = snakemake@params$samples,
                      condition = conds)

df <- read.csv(snakemake@input$combined)
rownames(df) <- df$gene_id
df$gene_id <- NULL

dds <- DESeq2::DESeqDataSetFromMatrix(countData = df,
                                      colData = colData,
                                      design = formula
                                      (snakemake@params$design))


# Start the DE analysis
dds <- DESeq2::DESeq(dds)

outputs <- NULL

# Get results from each combination of conditions
for (g in unique(snakemake@params$conditions)) {
  if (g != snakemake@params$reference) {
    contrast <- c(snakemake@params$column, g, snakemake@params$reference)
    res <- results(dds, contrast = contrast)
    fn <- paste(snakemake@wildcards$BASE, "deseq2", "analysis",
                snakemake@wildcards$METHOD, snakemake@wildcards$METHOD,
                sep = "/")
    name <- paste(fn, g, "vs", snakemake@params$reference, sep = "_")
    name <- paste0(name, ".csv")
    df <- as.data.frame(res)
    data.table::setDT(df, keep.rownames = "gene_id")
    write.csv(df, file = name, row.names = FALSE)
    outputs <- c(outputs, name)
  }
}

# Generate file for snakamake by renaming the input file
summary <- paste(snakemake@wildcards$BASE, "deseq2", "analysis",
                 snakemake@wildcards$METHOD, "analysis.log",
                 sep = "/")


fileConn <- file(summary)
fc <- "Analysis \t: DESeq2"
fc <- c(fc, paste0("Start Time \t: ", start_time))
fc <- c(fc, paste0("Finish Time \t: ", Sys.time()))
fc <- c(fc, paste0("Samples \t: ", paste0(snakemake@params$samples, collapse
  = ",")))
fc <- c(fc, paste0("Conditions \t: ", paste0(snakemake@params$conditions,
                                             collapse = ",")))
fc <- c(fc, paste0("Design Formula \t: ", snakemake@params$design))
fc <- c(fc, paste0("Reference \t: ", snakemake@params$reference))
fc <- c(fc, paste0("Outputs \t: ", paste(outputs, sep = "")))
writeLines(fc, fileConn)
close(fileConn)

abort()
{
    echo >&2 '
***************
*** ABORTED ***
***************
'
    echo "An error occurred. Exiting..." >&2
    exit 1
}

trap 'abort' 0

set -e

###############################################
## It assumes awk, bedtools and mashmap is
## available.
## We have tested this script with
## awk 4.1.3, bedtools v2.28.0 and mashmap v2.0
## on an Ubuntu system.
###############################################

threads=1
awk="awk"
bedtools="bedtools"
mashmap="mashmap"

# Argument Parsing
print_usage_and_exit () {
    echo "Usage: $0 [-j <N> =1 default] [-b <bedtools binary path> =bedtools default] [-m <mashmap binary path> =mashmap default] -a <gtf file> -g <genome fasta> -t <txome fasta> -o <output path>"
    exit 1
}

echo "****************"
echo "*** getDecoy ***"
echo "****************"
while getopts ":a:b:o:j:h:g:t:m:" opt; do
    case $opt in
        b)
            bedtools=`realpath $OPTARG`
            echo "-b <bedtools binary> = $bedtools"
            ;;
        m)
            mashmap=`realpath $OPTARG`
            echo "-m <mashmap binary> = $mashmap"
            ;;
        a)
            gtffile=`realpath $OPTARG`
            echo "-a <Annotation GTF file> = $gtffile"
            ;;
        o)
            outfolder="$OPTARG"
            echo "-o <Output files Path> = $outfolder"
            ;;
        j)
            threads="$OPTARG"
            echo "-j <Concurrency level> = $threads"
            ;;
        g)
            genomefile=`realpath $OPTARG`
            echo "-g <Genome fasta> = $genomefile"
            ;;
        t)
            txpfile=`realpath $OPTARG`
            echo "-t <Transcriptome fasta> = $txpfile"
            ;;
        h)
            print_usage_and_exit
            ;;
        \?)
            echo "Invalid option: -$OPTARG"
            print_usage_and_exit
            ;;
        :)
            echo "Option -$OPTARG requires an argument."
            print_usage_and_exit
            ;;
    esac
done

# Required arguments
if [ -z "$gtffile" -o -z "$outfolder" -o -z "$genomefile" -o -z "$txpfile" -o -z "$mashmap" -o -z "$awk" -o -z "$bedtools" -o -z "$threads" ]
then
    echo "Error: missing required argument(s)"
    print_usage_and_exit
fi

mkdir -p $outfolder
cd $outfolder

# extracting all the exonic features to mask
echo "[1/10] Extracting exonic features from the gtf"
$awk -v OFS='\t' '{if ($3=="exon") {print $1,$4,$5}}' $gtffile > exons.bed

# masking the exonic regions from the genome
echo "[2/10] Masking the genome fasta"
$bedtools maskfasta -fi $genomefile -bed exons.bed -fo reference.masked.genome.fa

# aligning the transcriptome to the masked genome
echo "[3/10] Aligning transcriptome to genome"
$mashmap -r reference.masked.genome.fa -q $txpfile -t $threads --pi 80 -s 500

# extracting the bed files from the reported alignment
echo "[4/10] Extracting intervals from mashmap alignments"
$awk -v OFS='\t' '{print $6,$8,$9}' mashmap.out | sort -k1,1 -k2,2n - > genome_found.sorted.bed

# merging the reported intervals
echo "[5/10] Merging the intervals"
$bedtools merge -i genome_found.sorted.bed > genome_found_merged.bed

# extracting relevant sequence from the genome
echo "[6/10] Extracting sequences from the genome"
$bedtools getfasta -fi reference.masked.genome.fa -bed genome_found_merged.bed -fo genome_found.fa

# concatenating the sequence at per chromsome level to extract decoy sequences
echo "[7/10] Concatenating to get decoy sequences"
$awk '{a=$0; getline;split(a, b, ":");  r[b[1]] = r[b[1]]""$0} END { for (k in r) { print k"\n"r[k] } }' genome_found.fa > decoy.fa

# concatenating decoys to transcriptome
echo "[8/10] Making gentrome"
cat $txpfile decoy.fa > gentrome.fa

# extracting the names of the decoys
echo "[9/10] Extracting decoy sequence ids"
grep ">" decoy.fa | $awk '{print substr($1,2); }' > decoys.txt

# removing extra files
echo "[10/10] Removing temporary files"
rm exons.bed reference.masked.genome.fa mashmap.out genome_found.sorted.bed genome_found_merged.bed genome_found.fa decoy.fa reference.masked.genome.fa.fai

trap : 0
echo >&2 '
**********************************************
*** DONE Processing ...
*** You can use files `$outfolder/gentrome.fa`
*** and $outfolder/decoys.txt` with
*** `salmon index`
**********************************************
'