Snakefile

import os
import sys
import glob
import re

####################	 FUNCTIONS	############################


# retrieves all forward/reverse fastq files from config file based on the sample name (list of all chunks)
def getFastqs(wildcards):
	for s in config["samples"]:
		if(s["sample"] == wildcards.sampleName):
			chunks=list()
			if(wildcards.read == "R1"):
				for pair in s["files"]:
					chunks.append(pair["forward"]["path"])
			if(wildcards.read == "R2"):
				for pair in s["files"]:
					chunks.append(pair["reverse"]["path"])
			return chunks

def getFastqPair(wildcards):
	for s in config["samples"]:
		if(s["sample"] == wildcards.sampleName):
			if(wildcards.strand=="R1"):
				return s["files"][int(wildcards.pairNumber)]["forward"]["path"]
			elif(wildcards.strand=="R2"):
				return s["files"][int(wildcards.pairNumber)]["reverse"]["path"]


def getCleanedFastq(wildcards):
	fastqs = list()
	for s in config["samples"]:
		if(s["sample"] == wildcards.sampleName):
			nbPairs = len(s["files"])
			if(wildcards.read == "R1"):
				for nb in range(nbPairs):
					fastqs.append(OUTPUTDIR+"/Samples/"+s["sample"]+"/CUTADAPT/"+s["sample"]+"_"+str(nb)+"_R1.fastq.gz")
			if(wildcards.read == "R2"):
				for nb in range(nbPairs):
					fastqs.append(OUTPUTDIR+"/Samples/"+s["sample"]+"/CUTADAPT/"+s["sample"]+"_"+str(nb)+"_R2.fastq.gz")
	fastqs.sort()
	return fastqs

def getConditions():
	sampleList = list()
	for s in config["samplesCondition"]:
		sd = dict()
		sd["samplename"] = s["name"]
		sd["filename"] = s["name"]
		sd["condition"] = s["condition"]
		sampleList.append(sd)
	return sampleList

def getConditionForComp(wildcards):
	s = re.search("(^[\w-]+)__vs__([\w-]+)$", wildcards.comp)
	d = dict()
	d["cond1"] = s.group(1)
	d["cond2"] = s.group(2)
	#return {"cond1":s.group(1),"cond2":s.group(2)}
	return d

## Fonction pour {wildcards.comp}
#def getTargetAll{wildcards.comp}
#	fastqs = list{wildcards.comp}
#	for s in conf{wildcards.comp}
#		nbPairs ={wildcards.comp}
#		for nb in{wildcards.comp}
#			fastq{wildcards.comp}"/Samples/"+s["sample"]+"/CUTADAPT/"+s["sample"]+"_"+str(nb)+"_R1.fastq.gz")
#			fastq{wildcards.comp}"/Samples/"+s["sample"]+"/CUTADAPT/"+s["sample"]+"_"+str(nb)+"_R2.fastq.gz")
#	return fastqs{wildcards.comp}

####################	 CMD to launch pipeline   ############################
## snakemake --config proj="project.json" conf="config.json" -nrp

########## loading configuration variable ##################
configfile: config["proj"]
configfile: config["conf"]

########## define global variables ##########
SNAKEFILEDIR = workflow.basedir
OUTPUTDIR = os.path.abspath(config["outdir"])
SCRIPTPATH = os.path.join(SNAKEFILEDIR,"scripts")

GTF = config["reference"]["gtf"]

FASTA = config["reference"]["fasta"]
FASTABASENAME = os.path.splitext(FASTA)[0]

BIOMART = config["reference"]["biomart"]
REFNAME = config["reference"]["name"]

GENOMEDIR=config["reference"]["STARindexDir"]

FWADAPT=config["cutadapt-forward"]
RVADAPT=config["cutadapt-reverse"]

## Overhang calculation for the generation of the reference files in STAR
OVERHANG = int(config["read-length"]) - 1

# List samples to trigger the target rule to create
SAMPLESALL = list()
for s in config["samples"]:
	SAMPLESALL.append(s["sample"])

# List comparisons to trigger the target rule
COMPARISONS = [c for c in config["comparisons"]]

# List of the direction of the read
READS=["R1","R2"]


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

rule all:
	input:
	 	qc=OUTPUTDIR+"/Report/data/general/fastqc/multiqc_report.html",
	 	#fastqc=expand(OUTPUTDIR+"/Samples/{sampleName}/FASTQC_cleaned/{sampleName}_{read}.concat_fastqc.html",sampleName=SAMPLESALL,read=READS)
	 	index=expand(OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam.bai",sampleName=SAMPLESALL),
		#htseq=expand(OUTPUTDIR+"/DESEQ2/counts/{sampleName}",sampleName=SAMPLESALL)
	 	deg=expand(OUTPUTDIR+"/DESEQ2/results/{suffixe}",suffixe=["NormalizedCountMatrix.txt","NormalizedCountMatrixFiltered.txt","PCAplot.png","sampletosampledistance.jpeg"]),
		rep=OUTPUTDIR+"/Report/report.html"

rule fastqc:
	input: 	OUTPUTDIR+"/Samples/{sampleName}/FASTQC/{sampleName}_{read}.fastq.gz"
	output: OUTPUTDIR+"/Samples/{sampleName}/FASTQC/{sampleName}_{read}_fastqc.html"
	params:	outFolder = OUTPUTDIR+"/Samples/{sampleName}/FASTQC"
	shell: 	"""
			fastqc -o {params.outFolder} {input}
			"""

rule fastqc_cleaned:
	input:	OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{read}.concat.fastq.gz"
	output: OUTPUTDIR+"/Samples/{sampleName}/FASTQC_cleaned/{sampleName}_{read}.concat_fastqc.html"
	params:	outFolder = OUTPUTDIR+"/Samples/{sampleName}/FASTQC_cleaned"
	shell: 	"""
			fastqc -o {params.outFolder} {input}
			"""

rule concat_fastq_chunks:
	input: 	getFastqs
	output: temp(OUTPUTDIR+"/Samples/{sampleName}/FASTQC/{sampleName}_{read}.fastq.gz")
	run:
		chunks = ' '.join(input)
		shell("cat {chunks} > {output}")

rule multiQC:
	input: 	f1=expand(OUTPUTDIR+"/Samples/{sampleName}/FASTQC/{sampleName}_{read}_fastqc.html", sampleName=SAMPLESALL, read=READS),
			f2=expand(OUTPUTDIR+"/Samples/{sampleName}/FASTQC_cleaned/{sampleName}_{read}.concat_fastqc.html", sampleName=SAMPLESALL, read=READS),
			bam=expand(OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam", sampleName=SAMPLESALL)
	output: OUTPUTDIR+"/Report/data/general/fastqc/multiqc_report.html"
	shell: 	"""
			multiqc -f -e general_stats -e tophat -e bowtie2 {OUTPUTDIR}/Samples -o {OUTPUTDIR}/Report/data/general/fastqc
			"""

rule dezipFastq:
	input:	getFastqPair
	output:	temp(OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_{strand}.temp.fastq")
	shell:	"""
			gzip -dc {input} > {output}
			"""

rule prinseq:
	input: 	expand(OUTPUTDIR+"/Samples/{{sampleName}}/PRINSEQ/{{sampleName}}_{{pairNumber}}_{strand}.temp.fastq", strand=READS)
	output:	temp(expand(OUTPUTDIR+"/Samples/{{sampleName}}/PRINSEQ/{{sampleName}}_{{pairNumber}}_{suffixe}.fastq", suffixe=["bad_1","bad_2","good_1","good_2","good_1_singletons","good_2_singletons"])),
			log = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}.log"
	params:	preGood = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_good",
			preBad = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_bad"
	shell:	"""
			perl `which prinseq-lite.pl` -min_qual_mean 30 -no_qual_header -fastq {input[0]} -fastq2 {input[1]} -out_good {params.preGood} -out_bad {params.preBad} -log {output.log}
			if(! test -f {params.preBad}_1.fastq);then touch {params.preBad}_1.fastq {params.preGood}_1_singletons.fastq;fi
			if(! test -f {params.preBad}_2.fastq);then touch {params.preBad}_2.fastq {params.preGood}_2_singletons.fastq;fi
			"""

rule concatPrintseqOutput:
	input:	bad = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_bad_{strand}.fastq",
			good = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_good_{strand}.fastq",
			singleton = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_good_{strand}_singletons.fastq"
	output:	temp(OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_R{strand}.fastq.gz")
	params: tmpdir = OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/"
	shell:	"""
			awk '{{if(NR%2==1){{print}} else{{lineLength=length();if(NR%4==0){{curCar="#"}}else{{curCar="N"}};for(c=0;c<lineLength;c++)printf(curCar); printf "\\n"}}}}' {input.bad} | cat {input.good} {input.singleton} - | paste - - - - | LC_ALL=C sort -k1,1 -t " " -T {params.tmpdir} | tr "\\t" "\\n" | gzip --best - > {output}
			"""

rule cutadaptR1:
	input:	OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_R1.fastq.gz"
	output:	fastq = temp(OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{pairNumber}_R1.fastq.gz"),
			report = OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{pairNumber}_R1.fastq.report"
	shell:	"""
			cutadapt -a {FWADAPT} 2> {output.report} {input} | java -jar {SCRIPTPATH}/pademptyfastq-fat.jar 2> /dev/null | gzip --best -c > {output.fastq}
			"""

rule cutadaptR2:
	input:	OUTPUTDIR+"/Samples/{sampleName}/PRINSEQ/{sampleName}_{pairNumber}_R2.fastq.gz"
	output:	fastq = temp(OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{pairNumber}_R2.fastq.gz"),
			report = OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{pairNumber}_R2.fastq.report"
	shell:	"""
			cutadapt -a {RVADAPT} 2> {output.report} {input} | java -jar {SCRIPTPATH}/pademptyfastq-fat.jar 2> /dev/null | gzip --best -c > {output.fastq}
			"""

rule concat_cleanedfastq_chunks:
	input: 	getCleanedFastq
	output: OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_{read}.concat.fastq.gz"
	run:
		chunks = ' '.join(input)
		shell("cat {chunks} > {output}")


rule genomeFiles:
	input:	GTF,
    		FASTA
	output: GENOMEDIR+"/Genome"
	params: cpu=config["align-cpu"]
	shell: 	"""
			STAR --runThreadN {params.cpu} --runMode genomeGenerate --genomeDir {GENOMEDIR} --genomeFastaFiles {FASTA} --sjdbGTFfile {GTF} --sjdbOverhang {OVERHANG} --limitGenomeGenerateRAM 168632718037
			"""

rule star_PE:
	input:
        	R1=OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_R1.concat.fastq.gz",
	        R2=OUTPUTDIR+"/Samples/{sampleName}/CUTADAPT/{sampleName}_R2.concat.fastq.gz",
	        genome=GENOMEDIR+"/Genome"
	output: OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam"
	params:	cpu=config["align-cpu"]
	resources: parallel_star=1
	shell:	"""
			STAR --runThreadN {params.cpu} --genomeDir {GENOMEDIR} --readFilesCommand zcat --outFileNamePrefix {OUTPUTDIR}/Samples/{wildcards.sampleName}/STAR/{wildcards.sampleName} --readFilesIn {input.R1} {input.R2} --outSAMtype BAM SortedByCoordinate --outBAMsortingThreadN 6 --outSAMattrRGline ID:"{wildcards.sampleName}" SM:"{wildcards.sampleName}" LB:"{wildcards.sampleName}" CN:"PFBirdNantes"
			"""


rule createBamIndex:
	input: 	OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam"
	output: OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam.bai"
	shell: 	"""
			samtools index {input} {output}
			"""


rule htseq:
	input: 	OUTPUTDIR+"/Samples/{sampleName}/STAR/{sampleName}Aligned.sortedByCoord.out.bam"
	output: OUTPUTDIR+"/DESEQ2/counts/{sampleName}"
	params:	library_type=config["library-type"]
	shell:  """
			if [ {params.library_type} = "reverse" ]
			then
				htseq-count -r pos -s reverse -f bam {input} {GTF} > {output}
			elif [ {params.library_type} = "yes" ]
			then
				htseq-count -r pos -s yes -f bam {input} {GTF} > {output}
			else
				htseq-count -r pos -s no -f bam {input} {GTF} > {output}
			fi
			"""

rule deseq2_conditions:
	output: tab = OUTPUTDIR+"/DESEQ2/DESEQ2_CONDITIONS.tab"
	run:
		condArray = getConditions()
		with open(output.tab, 'w') as condfile:
			condfile.write("samplename"+"\t"+"filename"+"\t"+"condition"+"\n")
			for s in condArray:
				condfile.write(s["samplename"]+"\t"+s["filename"]+"\t"+s["condition"]+"\n")

rule deseq2:
	input:	conditions=OUTPUTDIR+"/DESEQ2/DESEQ2_CONDITIONS.tab",
			counts=expand(OUTPUTDIR+"/DESEQ2/counts/{sampleName}",sampleName=SAMPLESALL)
	output: expand(OUTPUTDIR+"/DESEQ2/results/{suffixe}",suffixe=["NormalizedCountMatrix.txt","NormalizedCountMatrixFiltered.txt","PCAplot.png","sampletosampledistance.jpeg"]),
			OUTPUTDIR+"/DESEQ2/rData/dds.RData"
	params:	outdir = OUTPUTDIR+"/DESEQ2/results",
			countsdir = OUTPUTDIR+"/DESEQ2/counts"
	shell:	"""
			cat {SCRIPTPATH}/run_deseq2.R | R --slave --args {input.conditions} {params.countsdir} {params.outdir} {BIOMART} 
			"""

rule deg:
	input:	rdata=OUTPUTDIR+"/DESEQ2/rData/dds.RData"
	output:	expand(OUTPUTDIR+"/DESEQ2/results/{{comp}}/{{comp}}_{file}",file=["DEseqRes.tsv","All_DEG.tsv","UP_DEG.tsv","DOWN_DEG.tsv","Volcano-plot.png","MA-plot.png","clustDEgene.png"])
	params:	cond1=lambda wildcards : getConditionForComp(wildcards)["cond1"],
			cond2=lambda wildcards : getConditionForComp(wildcards)["cond2"]
	shell:	"""
			cat {SCRIPTPATH}/deg.R | R --slave --args {input.rdata} {params.cond1} {params.cond2}
			"""
			
rule annot:
	input:	expand(OUTPUTDIR+"/DESEQ2/results/{{comp}}/{{comp}}_{file}", file=["All_DEG.tsv","DEseqRes.tsv"]),
			corrAnnotations=(SCRIPTPATH+"/corresIDorg.txt"),
			rdata=OUTPUTDIR+"/DESEQ2/rData/dds.RData"
	output:	expand(OUTPUTDIR+"/DESEQ2/results/{{comp}}/{{comp}}_{file}",file=["gseGo.txt","gseKegg.txt"])
	params:	cond1=lambda wildcards : getConditionForComp(wildcards)["cond1"],
			cond2=lambda wildcards : getConditionForComp(wildcards)["cond2"]
	shell:	"""
			cat {SCRIPTPATH}/annot.R | R --slave --args {input.rdata} {params.cond1} {params.cond2} {REFNAME} {input.corrAnnotations}
			"""

rule copyDataForReport:
	input:
			expand(OUTPUTDIR+"/DESEQ2/results/{suffixe}",suffixe=["NormalizedCountMatrix.txt","NormalizedCountMatrixFiltered.txt","PCAplot.png","sampletosampledistance.jpeg"]),
			expand(OUTPUTDIR+"/DESEQ2/results/{comp}/{comp}_{file}",file=["gseGo.txt","gseKegg.txt","DEseqRes.tsv","All_DEG.tsv","UP_DEG.tsv","DOWN_DEG.tsv","Volcano-plot.png","MA-plot.png","clustDEgene.png"], comp=COMPARISONS)
	output:
			temp(touch(OUTPUTDIR+"/Report/TEMPLATE/templateCopied.txt"))
	params:
			templateFolder = OUTPUTDIR+"/Report/TEMPLATE",
			deseqFolder = OUTPUTDIR+"/Report/DESEQ2"
	shell:
			"""
			if test -d {params.templateFolder}; then rm -rf {params.templateFolder};fi
			mkdir {params.templateFolder}
			cp -r {SNAKEFILEDIR}/TEMPLATE/assets {SNAKEFILEDIR}/TEMPLATE/images {params.templateFolder}
			if test -d {params.deseqFolder}; then rm -rf {params.deseqFolder};fi
			mkdir {params.deseqFolder}
			cp -r {OUTPUTDIR}/DESEQ2/results {params.deseqFolder}
			"""

rule buildReport:
	input:	OUTPUTDIR+"/Report/TEMPLATE/templateCopied.txt"
	output:	report = OUTPUTDIR+"/Report/report.html"
	shell:
			"""
			python {SNAKEFILEDIR}/scripts/make_html.py -c {config[conf]} -s {config[proj]} -t {SNAKEFILEDIR}/TEMPLATE/index.html -o {output.report}
			"""