shell:
	"unset TMPDIR; module load star/2.7.5c; set -euo pipefail;  " +
	"STAR --readFilesCommand zcat {params.extra} --runThreadN {threads} --genomeDir {params.index} " +
	"--readFilesIn {input.fastq1} {input.fastq2} --outSAMtype BAM SortedByCoordinate " +
	"--outFileNamePrefix {params.star_prefix} --outStd Log  > {log} 2>&1"

shell:
	"""
	module load R;
	Rscript scripts/deseq2-init.R {input.counts} {output} {params.samples} {log} {threads} &> {log}
	"""

shell:
	"unset TMPDIR; module load samtools; samtools index {input}"

shell:
	"unset TMPDIR; module load cutadapt; cutadapt {params} -o {output.fastq1} -p {output.fastq2} {input} > {output.qc}"

import pysam
import pandas as pd
import numpy as np


def get_sample_counts(bam_file, reference):

	bamfile_obj = pysam.AlignmentFile(bam_file,'rb')

	ref = open(reference + "/chrName.txt",'r')

	curr_file_counts = []
	ref_list = []
	for seq in ref:
		seq = seq.strip()
		curr_seq_reads = bamfile_obj.fetch(seq)
		read_counter = 0
		for curr_read in curr_seq_reads:
			read_counter = read_counter + 1
		curr_file_counts.extend([read_counter])
		ref_list.extend([seq])

	counts = pd.DataFrame(curr_file_counts)
	counts = counts.transpose()
	counts.columns=ref_list
	counts = counts.transpose()
	return counts

############################################################


bam_files = snakemake.input.bams
counts = [get_sample_counts(f, snakemake.params.ref)
			for f in bam_files]

samples = snakemake.params.samples
print("checking sample order in the python script")
print(samples)
for t, sample in zip(counts, samples):
	t.columns = [sample]
print("checking sample order in the python script after the zip")
for t in counts:
	print(t.columns)
matrix = pd.concat(counts, axis=1)
matrix.index.name = "gene"
print("Checking columns in matrix before groupby")
print(matrix.columns)
matrix = matrix.groupby(matrix.columns, axis=1).sum()
print("Checking columns in matrix")
print(matrix.columns)
matrix.to_csv(snakemake.output[0], sep="\t")

args=commandArgs(trailingOnly=TRUE)
counts_file=args[1]
output_file=args[2]
sample_file=args[3]
log_file=args[4]
threads=args[5]



log <- file(log_file, open="wt")
sink(log)
sink(log, type="message")

library("DESeq2")
library(dplyr)

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

# colData and countData must have the same sample order, but this is ensured
# by the way we create the count matrix
cts <- read.table(counts_file, header=TRUE, row.names="gene", check.names=FALSE)
coldata <- read.table(sample_file, header=TRUE, row.names="sample", check.names=FALSE)

sorted_cts <- cts[order(names(cts))]


dds <- DESeqDataSetFromMatrix(countData= cts[order(names(cts))],
                              colData=coldata[order(row.names(coldata)),],
                              design=~ condition)

# remove uninformative columns
dds <- dds[ rowSums(counts(dds)) > 1, ]
# normalization and preprocessing
dds <- DESeq(dds, parallel=parallel)

saveRDS(dds, file=output_file)

args=commandArgs(trailingOnly=TRUE)
rds=args[1]
output_table=args[2]
output_ma_plot=args[3]
params.config_yaml=args[4]
log=args[5]
threads=args[6]

log <- file(log, open="wt")
sink(log)
sink(log, type="message")

library("DESeq2")

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

dds <- readRDS(rds)

library("yaml")
config_yaml <- read_yaml(params.config_yaml)
contrast <- c("condition", as.vector(config_yaml$diffexp$contrasts$"MJC_cambium-vs-MJ_cambium"))
res <- results(dds, contrast=contrast, parallel=parallel)
# shrink fold changes for lowly expressed genes
res <- lfcShrink(dds, contrast=contrast, res=res,type="normal")
# sort by p-value
res <- res[order(res$padj),]
# TODO explore IHW usage


# store results
svg(output_ma_plot)
plotMA(res, ylim=c(-8,8))
dev.off()

write.table(as.data.frame(res), file=output_table)

import pysam
import pandas as pd
import numpy as np
from scipy.stats import variation
import re

def get_number_of_reads_in_cat(star_log_file, target_string):
	star_log = open(star_log_file,'r')
	for line in star_log:
		if(target_string in line):
			regex = r"\s\|\s(\d+)"
			match_list = re.findall(regex, line, flags=0)
			return match_list[0]

def get_sample_counts(star_log_file):

	#print("in get sample counts")
	#print("this is star log file")
	#print(star_log_file)

	curr_file_counts = []
	target_strings = ["Number of input reads",
	"Uniquely mapped reads number",
	"Number of reads mapped to multiple loci",
	"Number of reads mapped to too many loci",
	"Number of reads unmapped: too many mismatches",
	"Number of reads unmapped: too short",
	"Number of reads unmapped: other"]

	for string in target_strings:
		curr_file_counts.append(
			get_number_of_reads_in_cat(star_log_file,string))

	counts = pd.DataFrame(curr_file_counts)
	counts = counts.transpose()
	counts.columns=target_strings
	counts = counts.transpose()

	return counts

############################################################
star_logs = np.unique(snakemake.input)
#print("Passing this to get_sample_counts")
#print(star_logs)
#print(len(star_logs))
counts = [get_sample_counts(f)
			for f in star_logs]

samples = snakemake.params.samples
#print("checking sample order in the python script")
#print(samples)
for t, sample in zip(counts, samples):
	t.columns = [sample]
#print("checking sample order in the python script after the zip")
#for t in counts:
#	print(t.columns)
matrix = pd.concat(counts, axis=1)
matrix.index.name = "gene"
#print(matrix.columns)
matrix = matrix.groupby(matrix.columns, axis=1).sum()
matrix = matrix.transpose()
#print("Checking columns in matrix")
#print(matrix.columns)
matrix["proportion_of_library_reads_out_of_all_reads"] = \
	matrix["Number of input reads"] / matrix["Number of input reads"].sum()
matrix["percent_of_library_reads_mapped"] = \
	(matrix["Uniquely mapped reads number"] + matrix["Number of reads mapped to multiple loci"]) / matrix["Number of input reads"]

matrix.to_csv(snakemake.output[0], sep="\t")