library(DESeq2)
library(readr)

parallel <- FALSE
if (snakemake@threads > 1) {
    library("BiocParallel")
    # setup parallelization
    register(MulticoreParam(snakemake@threads))
    parallel <- TRUE
}

# Loading the parameters
project=snakemake@params[["project"]]
samples=snakemake@params[["samples"]]
ref_level=snakemake@params[["ref_level"]]
normalized_counts_file=snakemake@output[["normalized_counts_file"]]

# Rename column name of the count matrix as coldata
# colData and countData must have the same sample order
cts <- as.matrix(read.table(snakemake@input[["cts"]], header=T, row.names = 1))
coldata_read <- read.delim(snakemake@input[["coldata"]], header=TRUE, comment.char="#", quote="")
colnames(cts) <- coldata_read[,1]

coldata <- coldata_read[,-1]
rownames(coldata) <- coldata_read[,1]
coldata$condition <- factor(coldata_read$condition)
coldata$type <- factor(coldata_read$type)

rmproj_list = as.list(strsplit(snakemake@params[["rmproj_list"]], ",")[[1]])

if(length(rmproj_list)!=0){
  for (i in 1:length(rmproj_list)) {
      name <- rmproj_list[[i]]
      coldata <- coldata[-match((name), table = rownames(coldata)), ]
  }
}

# Check that sample names match in both files
if (all(colnames(cts) %in% rownames(coldata)) & all(colnames(cts) == rownames(coldata))){
  # Create the DESeq2 object
  dds <- DESeqDataSetFromMatrix(countData = cts,
                                colData = coldata,
                                design = ~ condition)
} else {
  print("sample names doesn't match in both files")
}

# Remove uninformative columns (to do when filter not already done with the CPM threshold)
#dds <- dds[ rowSums(counts(dds)) > 10, ]

# Specifying the reference level
dds$condition <- relevel(dds$condition, ref = ref_level)

# DESeq : Normalization and preprocessing (counts divided by sample-specific size factors
# determined by median ratio of gene counts relative to geometric mean per gene)
dds <- DESeq(dds, parallel=parallel)
# To save the object in a file for later use
saveRDS(dds, file=snakemake@output[["rds"]])

# Already done in the DESeq function
dds <- estimateSizeFactors( dds)
print(sizeFactors(dds))
# Save the normalized data matrix
normalized_counts <- counts(dds, normalized=TRUE)
write.table(normalized_counts, file=normalized_counts_file, sep="\t", quote=F, col.names=NA)

shell:
  "fastqc --outdir 05_Output/01_fastqc/ {input}"

shell:
  """
  sample=({input.sample})
  sample_trimmed=({output.sample_trimmed})
  sample_untrimmed=({output.sample_untrimmed})
  len=${{#sample[@]}}
  for (( i=0; i<$len; i=i+2 ))
  do trimmomatic PE -threads 4 ${{sample[$i]}} ${{sample[$i+1]}} ${{sample_trimmed[$i]}} ${{sample_untrimmed[$i]}} ${{sample_trimmed[$i+1]}} ${{sample_untrimmed[$i+1]}} LEADING:20 TRAILING:15 SLIDINGWINDOW:4:15 MINLEN:36
  done
  """

shell: 
  """
  multiqc -n {output.trim_multi_html} {input.trim_qc} --force #run multiqc
  rm -rf {params.multiqc_output_trim} #clean-up
  """

shell:
  """
  input={input.ref}
  output={params.index}
  if [ ${{input}} == '*.gz' ];then
    gunzip ${{input}}
    input="${{input%.*}}"
  fi
  hisat2-build ${{input}} ${{output}}
  """

shell:
  """
  index={params.index}
  sample_trimmed=({input.sample_trimmed})
  len=${{#sample_trimmed[@]}}
  reads=({params.reads})
  sam=({params.sam})
  bam=({output.bam})
  flag=0
  if [ ${{reads}} == 'paired' ];then
    for (( i=0; i<$len; i=i+2 ))
      do hisat2 -p 12 -x ${{index}} -1 ${{sample_trimmed[$i]}} -2 ${{sample_trimmed[$i+1]}} -S ${{sam[$flag]}}
        samtools sort ${{sam[$flag]}} > ${{bam[$flag]}}
        rm ${{sam[$flag]}}
        flag=$((${{flag}}+1))
    done
  elif [ ${{reads}} == 'unpaired' ];then
    for (( i=0; i<$len; i++ ))
      do hisat2 -p 12 -x ${{index}} -U ${{sample_trimmed[$i]}} -S ${{sam[$i]}}
        samtools sort ${{sam[$i]}} > ${{bam[$i]}}
        rm ${{sam[$i]}}
    done
  else
    echo "Your fastq files have to be paired or unpaired. Please check the config.yaml"
  fi
  """

shell:
  """
  reads=({params.reads})
  geneid=({params.geneid})
  annotation={input.annotation}
  bam=({input.bam})
  countmatrices=({output.countmatrices})
  len=${{#bam[@]}}
  if [ ${{reads}} == 'paired' ];then
    for (( i=0; i<$len; i++ ))
      do featureCounts -T 12 -p -t exon -g ${{geneid}} -a ${{annotation}} -o ${{countmatrices[$i]}} ${{bam[$i]}}
    done
  elif [ ${{reads}} == 'unpaired' ];then
    for (( i=0; i<$len; i++ ))
      do featureCounts -T 12 -t exon -g ${{geneid}} -a ${{annotation}} -o ${{countmatrices[$i]}} ${{bam[$i]}}
    done
  fi
  """