rule bam_to_sam:
	input: "04_bam_files/{sample}.sorted.dedup.bam"
	output: "03_sam_files/{sample}.sorted.dedup.sam"
	log: "00_logs/{sample}_bam_to_sam.log"
	shell: "samtools view -@ THREADS {input} > {output} 2> {log}"
		
rule preprocess_data:
	input: "03_sam_files/{sample}.sorted.dedup.sam"
	output: directory("06_peakseq_files/{sample}_preprocess")
	log: "00_logs/{sample}_data_preprocess.log"
	run:
		shell("mkdir -p {output}")
		shell("PeakSeq -preprocess SAM {input} {output}")

def sample_group(wildcards):
	group = Samples[wildcards.sample]
	return expand("04_bam_files/{sample}.sorted.dedup.bam", sample = group)

def control_group(wildcards):
	group = Controls[wildcards.control]
	return expand("04_bam_files/{sample}.sorted.dedup.bam", sample = group)

rule make_configfile:
	input:
		s = sample_group,
		c = control_group
	output: "06_peakseq_files/{sample}_{control}_config.txt"
	run:
		shell('''
			echo \'Experiment_id {wildcards.sample}_{wildcards.control}\' >> {output}
			echo \'Enrichment_mapped_fragment_length 200\' >> {output}
			echo \'target_FDR 0.05\' >> {output}
			echo \'N_Simulations 50\' >> {output}
			echo \'Minimum_interpeak_distance 200\' >> {output}
			echo \'Mappability_map_file MAPPABILITY_MAP_FILE\' >> {output}
			echo \'ChIP_Seq_reads_data_dirs\' >> {output}
		''')
		for sample in input.s: shell('echo \'\t{sample}\' >> {output}')
		shell('echo \'Input_reads_data_dirs\' >> {output}')
		for control in input.c: shell('echo \'\t{control}\' >> {output}')
		shell('''
			echo \'narrowPeak_output_file_path 06_peakseq_peaks/{wildcards.sample}_{wildcards.control}.narrowPeak\' >> {output}
			echo \'Background_model Simulated\' >> {output}
			echo \'max_Qvalue 0.05\' >> {output}
		''')
	
rule call_peak:
	input: "06_peakseq_files/{sample}_{control}_config.txt"
	output: "06_peakseq_peaks/{sample}_{control}.narrowPeak"
	log: "00_logs/{sample}_{control}_peakseq_peaks.log"
	run:
		shell("PeakSeq -peak_select {input} 2> {log}")
		shell("rm beacon.log {wildcards.sample}_{wildcards.control}.txt seed.txt")
