Snakefile_wes

import os
import random
import shutil
import glob
import hashlib
import gzip

configfile: "config.yml"

def deterministic_seed(major_seed, index):
	unique_string = f"{major_seed}_{index}"
	hash_object = hashlib.sha256(unique_string.encode())
	seed_int = int(hash_object.hexdigest(), 16)
	# Reduce the seed to a range that is typical for seeds (e.g., 0 to 2^32-1)
	return seed_int % (2**31 - 1)

def generate_runs():
	base_output_dir = config["OUTPUT_DIRECTORY"]
	seed_from_config = config.get("SEED")
	if not seed_from_config:  # This checks for both None and empty string
		major_seed = random.randint(0, 2**31 - 1)
	else:
		try:
			major_seed = int(seed_from_config)  # Ensure major_seed is an integer if provided
		except ValueError:
			major_seed = random.randint(0, 2**31 - 1)  # Generate a random seed if the provided SEED is not a valid integer

	num_tech_replicates = int(config.get("FRAG_REPLICATES", 1))

	os.makedirs(base_output_dir, exist_ok=True)
	seed_log_path = os.path.join(base_output_dir, "seed_log.txt")

	with open(seed_log_path, 'w') as seed_log:
		seed_log.write(f"Major Seed: {major_seed}\n")
		runs = {}
		for i in range(config["NUM_SAMPLES"]):
			sample_dir = os.path.join(base_output_dir, f"sample_{i+1}")
			os.makedirs(sample_dir, exist_ok=True)
			runs[sample_dir] = []
			for j in range(num_tech_replicates):
				tech_rep_dir = os.path.join(sample_dir, f"tech_replicate_{j+1}")
				os.makedirs(tech_rep_dir, exist_ok=True)
				tech_rep_seed = deterministic_seed(major_seed, i * num_tech_replicates + j)
				runs[sample_dir].append((tech_rep_seed, tech_rep_dir))
				seed_log.write(f"Sample {i+1}, Tech Replicate {j+1}: Seed = {tech_rep_seed}\n")
	return runs

# Summoning the runs
runs = generate_runs()


from Bio import SeqIO

def get_genome_info(fasta_file):
	# Determine if the file is gzipped based on its extension
	if fasta_file.endswith('.gz'):
		open_func = gzip.open
	else:
		open_func = open

	with open_func(fasta_file, "rt") as handle:
		# A magical dictionary to hold both IDs and their lengths
		genome_info = {record.id: len(record.seq) for record in SeqIO.parse(handle, "fasta")}
	return genome_info

FASTA_PATH = os.path.join(config["INPUT_DIRECTORY"], config["FASTA_FILE"])
genome_info = get_genome_info(FASTA_PATH)
genome_ids = list(genome_info.keys())
# Function to aggregate PCR files for each GENOME_ID

def aggregate_pcr_files(wildcards):
	pattern = f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/" + \
			  f"{wildcards.GENOME_ID}_batch{{batch_id}}_filtered_pcr.csv"
	return expand(pattern, batch_id=range(config.get("NCORE")))

def aggregate_bedgz_files(wildcards):
	pattern = f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/" + \
			  f"{wildcards.GENOME_ID}_batch{{batch_id}}_filtered.bed.gz"
	return expand(pattern, batch_id=range(config.get("NCORE")))


def should_map_probes():
	return "PROBES_FASTA" in config and config["PROBES_FASTA"] != ""


rule all:
	input:
		expand("{sample_dir}/log_folder/0_Making_SQL_database.log", sample_dir=runs.keys()),
		os.path.join(config["INPUT_DIRECTORY"], "fasta_lengths.csv"),
		expand("{sample_dir}/log_folder/1_Defining_seq_and_mut.log", sample_dir=runs.keys()),
		expand("{sample_dir}/tech_replicate_{tech_rep}/log_folder/2_Fragmentation_prep.log", 
			   sample_dir=runs.keys(), tech_rep=range(1, config["FRAG_REPLICATES"] + 1)),
		[os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes.bed")] if should_map_probes() else [],
		expand("{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_fragmentation_complete.log",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1),
			   GENOME_ID=genome_ids),
		expand("{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_conversion_to_bedgz_complete.log",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1),
			   GENOME_ID=genome_ids),
		expand("{sample_dir}/tech_replicate_{tech_rep}/PCR_filtered.bed.gz", 
			   sample_dir=runs.keys(), tech_rep=range(1, config["FRAG_REPLICATES"] + 1)),
		expand("{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/sequenced_frags.fasta.gz",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, config["FRAG_REPLICATES"] + 1)),
		expand("{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_filtered.log",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1),
			   GENOME_ID=genome_ids),
		expand("{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R1.fastq.gz",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, config["FRAG_REPLICATES"] + 1)),
		expand("{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R2.fastq.gz",
			   sample_dir=runs.keys(),
			   tech_rep=range(1, config["FRAG_REPLICATES"] + 1)),
		os.path.join(config["OUTPUT_DIRECTORY"], "cleanup_1_complete.txt"),
		expand("{sample_dir}/cleanup_2_complete.txt",
		 		sample_dir=runs.keys()),
		os.path.join(config["OUTPUT_DIRECTORY"], "cleanup_3_complete.txt")

		

rule Making_fasta_lengths_info:
	input:
		FASTA_PATH = FASTA_PATH
	output:
		fasta_length_csv = os.path.join(config["INPUT_DIRECTORY"], "fasta_lengths.csv")
	run:
		import os
		# Remove the output file if it already exists
		if os.path.exists(output.fasta_length_csv):
			os.remove(output.fasta_length_csv)
		# Execute the Python script
		shell("python3 /usr/src/app/scripts/wes/fasta_lengths_extraction.py {input.FASTA_PATH} {output.fasta_length_csv}")


rule Making_SQL_database:
	input:
		FASTA_PATH = FASTA_PATH,
		FASTA_LENGTHS=os.path.join(config["INPUT_DIRECTORY"], "fasta_lengths.csv")
	output:
		log = "{sample_dir}/log_folder/0_Making_SQL_database.log"
	params:
		NCORE=config.get("NCORE", "")
	wildcard_constraints:
		sample_dir="|".join(runs.keys())
	run:
		command = "Rscript /usr/src/app/scripts/wes/0making_SQL_db.R"
		if params.NCORE:
			command += f" --cores {params.NCORE}"
		command += f" --fasta_sequence {input.FASTA_PATH}"
		command += f" --fasta_lengths_info {input.FASTA_LENGTHS}"
		command += f" --output {wildcards.sample_dir}"
		command += f" > {output.log}"
		shell(command)

rule defining_seq_and_mut:
	input:
		INPUT_LOG="{sample_dir}/log_folder/0_Making_SQL_database.log",
		FASTA_LENGTHS=os.path.join(config["INPUT_DIRECTORY"], "fasta_lengths.csv"),
		PROBES_TARGETS=os.path.join(config["INPUT_DIRECTORY"], config["PROBES_TARGES"])if config["PROBES_TARGES"] else "",
	output:
		log = "{sample_dir}/log_folder/1_Defining_seq_and_mut.log"
	params:
		seed = lambda wildcards: runs[wildcards.sample_dir][0][0],  
		OUTPUT_DIRECTORY=lambda wildcards: wildcards.sample_dir,
		NCORE=config.get("NCORE", ""),
		MUT_CONTEXT=config.get("MUT_CONTEXT", ""),
		TS_TV_RATIO=config.get("TS_TV_RATIO", ""),
		SUBSTITUTION_PROBABILITY_TABLE=lambda wildcards: os.path.join(config["INPUT_DIRECTORY"], config["SUBSTITUTION_PROBABILITY_TABLE"]) if config["SUBSTITUTION_PROBABILITY_TABLE"] else "",
		NR_COPIES=config.get("NR_COPIES", ""),
		MULTIFASTA_COPIES_TABLE=lambda wildcards: os.path.join(config["INPUT_DIRECTORY"], config["MULTIFASTA_COPIES_TABLE"]) if config["MULTIFASTA_COPIES_TABLE"] else "",
		MUT_GENOME_FRACTION=config.get("MUT_GENOME_FRACTION", ""),
		FRACTION_POSITION=config.get("FRACTION_POSITION", ""),
		MUTATION_RATE=config.get("MUTATION_RATE", ""),
		SBS_SIGNATURES=lambda wildcards: os.path.join(config["INPUT_DIRECTORY"], config["SBS_SIGNATURES"]) if config["SBS_SIGNATURES"] else "",

	wildcard_constraints:
		sample_dir="|".join(runs.keys())
	run:
		command = "Rscript /usr/src/app/scripts/wes/1mut_genomes.R"
		if params.NCORE:
			command += f" --cores {params.NCORE}"
		command += f" --fasta_lengths_info {input.FASTA_LENGTHS}"
		command += f" --target_exome_region {input.PROBES_TARGETS}"
		command += f" --output {params.OUTPUT_DIRECTORY}"
		if params.MUT_CONTEXT:
			command += f" --mutational_context {params.MUT_CONTEXT}"
		if params.TS_TV_RATIO:
			command += f" --Ts_Tv_ratio {params.TS_TV_RATIO}"
		if params.SUBSTITUTION_PROBABILITY_TABLE:
			command += f" --variant_probability_table {params.SUBSTITUTION_PROBABILITY_TABLE}"
		if params.SBS_SIGNATURES:
			command += f" --sbs_signature_table {params.SBS_SIGNATURES}"
		if params.NR_COPIES:
			command += f" --nr_copies {params.NR_COPIES}"
		if params.MULTIFASTA_COPIES_TABLE:
			command += f" --multifasta_copies {params.MULTIFASTA_COPIES_TABLE}"
		if params.MUT_GENOME_FRACTION:
			command += f" --mut_genome_fraction {params.MUT_GENOME_FRACTION}"
		if params.FRACTION_POSITION:
			command += f" --fraction_positions {params.FRACTION_POSITION}"
		if params.MUTATION_RATE:
			command += f" --specific_mut_rate {params.MUTATION_RATE}"
		command += f" --set_seed {params.seed}"
		command += f" > {output.log}"
		shell(command)

rule Fragmentation_prep:
	input:
		INPUT_LOG="{sample_dir}/log_folder/1_Defining_seq_and_mut.log"
	output:
		log="{sample_dir}/tech_replicate_{tech_rep}/log_folder/2_Fragmentation_prep.log"
	params:
		seed=lambda wildcards: int(runs[wildcards.sample_dir][int(wildcards.tech_rep) - 1][0]),
		OUTPUT_DIRECTORY=lambda wildcards: f"./tech_replicate_{wildcards.tech_rep}",
		INPUT_DIRECTORY=lambda wildcards: wildcards.sample_dir,
		FRAGMENT_FRACTION=config.get("FRAGMENT_FRACTION", "")
	wildcard_constraints:
		sample_dir="|".join(runs.keys()),
		tech_rep="\d+"
	run:
		command = "Rscript /usr/src/app/scripts/wes/2fragmentation_prep_wes.R"
		if params.FRAGMENT_FRACTION:
			command += f" --fragmentation_fraction {params.FRAGMENT_FRACTION}"
		command += f" --set_seed {params.seed}"
		command += f" --input_directory {params.INPUT_DIRECTORY}"
		command += f" --output_directory {params.OUTPUT_DIRECTORY}"
		command += f" > {output.log}"
		shell(command)

rule fragmentation:
	input:
		INPUT_LOG="{sample_dir}/tech_replicate_{tech_rep}/log_folder/2_Fragmentation_prep.log"
	output:
		log="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_fragmentation_complete.log"
	params:
		SELECTED="{sample_dir}/tech_replicate_{tech_rep}/{GENOME_ID}_fragmentation_selected.csv",
		MIN_FRAGMENT_LENGTH=config.get("MIN_FRAGMENT_LENGTH", ""),
		MAX_FRAGMENT_LENGTH=config.get("MAX_FRAGMENT_LENGTH", ""),
		NCORE=config.get("NCORE", ""),
		seed=lambda wildcards: int(runs[wildcards.sample_dir][int(wildcards.tech_rep) - 1][0]),
		OUTPUT_DIRECTORY=lambda wildcards: f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}",
		GENOME_ID=lambda wildcards: wildcards.GENOME_ID,
		GENOME_LENGTH=lambda wildcards: get_genome_info(FASTA_PATH)[wildcards.GENOME_ID]
	wildcard_constraints:
		sample_dir="|".join(runs.keys()),
		tech_rep="\d+",
		GENOME_ID="|".join(get_genome_info(FASTA_PATH).keys())
	run:
		command = f"/usr/src/app/scripts/wes/fragmentation_batches {params.SELECTED} "
		command += f"{params.GENOME_LENGTH} {params.MIN_FRAGMENT_LENGTH} {params.MAX_FRAGMENT_LENGTH} "
		command += f"{params.NCORE} {params.seed} {params.OUTPUT_DIRECTORY}/{params.GENOME_ID}"
		command += f" > {output.log}"
		shell(command)

rule conversion_to_bedgz:
	input:
		LOG="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_fragmentation_complete.log",
		FASTA_LENGTHS=os.path.join(config["INPUT_DIRECTORY"], "fasta_lengths.csv")
	output:
		LOG="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_conversion_to_bedgz_complete.log"
	params:
		NCORE=config.get("NCORE", ""),
		OUTPUT_DIRECTORY=lambda wildcards: f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}"
	wildcard_constraints:
		sample_dir="|".join(runs.keys()),
		tech_rep="\d+",
		GENOME_ID="|".join(get_genome_info(FASTA_PATH).keys())
	run:
		command = f"/usr/src/app/scripts/wes/csv_to_bed_gz {params.OUTPUT_DIRECTORY} {input.FASTA_LENGTHS} "
		command += f"{params.OUTPUT_DIRECTORY} {params.NCORE}"
		command += f" > {output.LOG}"
		shell(command)




# Conditional inclusion of the map_probes rule
if should_map_probes():
	rule map_probes:
		input:
			probes_fasta=os.path.join(config["INPUT_DIRECTORY"], config["PROBES_FASTA"]),
			ref_genome=os.path.join(config["INPUT_DIRECTORY"], config["FASTA_FILE"])
		output:
			bed=os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes.bed")
		params:
			ref_index=os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "ref_index"),
			probes_sam=os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes.sam"),
			probes_bam=os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes.bam"),
			probes_sorted_bam=os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes_sorted.bam"),
		shell:
			"""
			# Build Bowtie2 index
			bowtie2-build {input.ref_genome} {params.ref_index}

			# Map probe sequences to GRCh38_latest_genomic_modified
			bowtie2 -f -a -x {params.ref_index} -U {input.probes_fasta} -S {params.probes_sam}

			# Convert SAM to BAM
			samtools view -bS {params.probes_sam} > {params.probes_bam}

			# Sort BAM file
			samtools sort {params.probes_bam} -o {params.probes_sorted_bam}

			# Convert sorted BAM to BED
			bedtools bamtobed -split -i {params.probes_sorted_bam} > {output.bed}
			"""


rule filter_bed_files:
	input:
		LOG="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_conversion_to_bedgz_complete.log"
	output:
		log="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_filtered.log"
	params:
		exome_bed = (os.path.join(config["INPUT_DIRECTORY"], "mapping_reference", "probes.bed") 
					 if config.get("PROBES_FASTA", "") 
					 else os.path.join(config["INPUT_DIRECTORY"], config["PROBES_BED"])),
		NCORE=config.get("NCORE", ""),
		MATCHING_LENGTH=config.get("MATCHING_LENGTH", "")
	wildcard_constraints:
		sample_dir="|".join(runs.keys()),
		tech_rep="\d+",
		GENOME_ID="|".join(get_genome_info(FASTA_PATH).keys())
	shell:
		"""
		list_file="{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/bed_file_list_{wildcards.GENOME_ID}.txt"
		find {wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep} -name "{wildcards.GENOME_ID}_batch*.bed.gz" > $list_file
		bash /usr/src/app/scripts/wes/bedtools_parallels.sh $list_file {params.exome_bed} {params.NCORE} {params.MATCHING_LENGTH} > {output.log}
		rm -f $list_file
		"""

rule prep_pcr_file:
	input:
		input_log = "{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_filtered.log"
	output:
		log="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_converted_to_csv.log"
	params:
		NCORE=config.get("NCORE", "")
	wildcard_constraints:
		sample_dir="|".join(runs.keys()),
		tech_rep="\d+",
		GENOME_ID="|".join(get_genome_info(FASTA_PATH).keys())
	shell:
		"""
		list_file="{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/bed_file_list_{wildcards.GENOME_ID}.txt"
		find {wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep} -name "{wildcards.GENOME_ID}_batch*_filtered.bed.gz" > $list_file
		bash /usr/src/app/scripts/wes/csv_prep_pcr.sh $list_file {params.NCORE} > {output.log}
		rm -f $list_file
		"""

rule merge_pcr_files:
	input:
		log="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_converted_to_csv.log"
	output:
		merged_pcr = "{sample_dir}/tech_replicate_{tech_rep}/{GENOME_ID}_filtered_pcr.csv"
	params:
		pcr_files = aggregate_pcr_files
	shell:
		"cat {params.pcr_files} > {output.merged_pcr} && rm {params.pcr_files}"

rule merge_bed_files:
	input:
		log_convert_csv="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_converted_to_csv.log",
		log_filter_bed="{sample_dir}/tech_replicate_{tech_rep}/log_folder/{GENOME_ID}_all_filtered.log",
		merged_pcr_done = "{sample_dir}/tech_replicate_{tech_rep}/{GENOME_ID}_filtered_pcr.csv"
	params:
		bed_files = aggregate_bedgz_files,
	output:
		merged_bed = "{sample_dir}/tech_replicate_{tech_rep}/{GENOME_ID}_filtered.bed.gz"
	shell:
		"gunzip -c {params.bed_files} | gzip > {output.merged_bed} && rm {params.bed_files}"

rule merge_all_files:
	input:
		pcr_csv_files = lambda wildcards: [f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/{genome_id}_filtered_pcr.csv" 
										   for genome_id in genome_ids],
		bed_files = lambda wildcards: [f"{wildcards.sample_dir}/tech_replicate_{wildcards.tech_rep}/{genome_id}_filtered.bed.gz"
									   for genome_id in genome_ids]
	output:
		merged_pcr_csv = "{sample_dir}/tech_replicate_{tech_rep}/PCR_filtered.csv",
		merged_bed = "{sample_dir}/tech_replicate_{tech_rep}/PCR_filtered.bed.gz"
	shell:
		"""
		cat {input.pcr_csv_files} > {output.merged_pcr_csv} && rm {input.pcr_csv_files}
		gunzip -c {input.bed_files} | gzip > {output.merged_bed} && rm {input.bed_files}
		"""


rule pcr_sim:
	input:
		frag_coordinates_lengths = "{sample_dir}/tech_replicate_{tech_rep}/PCR_filtered.csv"
	output:
		csv_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/csv_list.txt",
		fasta_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/fasta_list.txt",
	params:
		seed = lambda wildcards: int(runs[wildcards.sample_dir][int(wildcards.tech_rep) - 1][0]),
		NCORE=config.get("NCORE", ""),
		NUM_CYCLES=config.get("NUM_CYCLES", ""),
		MIDPOINT_CYCLE=config.get("MIDPOINT_CYCLE", ""),
		K_PARAMETER_pcr=config.get("K_PARAMETER_pcr", ""),
		NUM_READS=config.get("NUM_READS", ""),
		OPT_FRAG_MODE=config.get("OPT_FRAG_MODE", ""),
		K_PARAMETER_seq=config.get("K_PARAMETER_seq", ""),
	log:
		log = "{sample_dir}/tech_replicate_{tech_rep}/log_folder/PCR_reaction/PCR_reaction.log"
	run:
		command = "/usr/src/app/scripts/wes/PCR_wes.sh"
		if params.seed:
			command += f" --set_seed {params.seed}"
		if params.NCORE:
			command += f" --cores {params.NCORE}"
		if params.NUM_READS:
			command += f" --num_reads {params.NUM_READS}"
		if params.OPT_FRAG_MODE:
			command += f' --optimal_length_mode "{params.OPT_FRAG_MODE}"'
		if params.K_PARAMETER_seq:
			command += f" --k_parameter_seq {params.K_PARAMETER_seq}"	
		if params.NUM_CYCLES:
			command += f" --num_cycles {params.NUM_CYCLES}"
		if params.MIDPOINT_CYCLE:
			command += f" --midpoint_cycle {params.MIDPOINT_CYCLE}"
		if params.K_PARAMETER_pcr:
			command += f" --k_parameter_pcr {params.K_PARAMETER_pcr}"
		command += f" --input_file {input.frag_coordinates_lengths}"
		command += f" > {log.log}"
		shell(command)



rule simulate_reads:
	input:
		csv_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/csv_list.txt",
		fasta_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/fasta_list.txt"
	output:
		genereted_reads1 = "{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R1.fastq.gz",
		genereted_reads2 = "{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R2.fastq.gz"
	params:
		ERROR_MODEL=config.get("ERROR_MODEL", ""),
		MODE=config.get("MODE", ""),
		NCORE=config.get("NCORE", ""),
		FASTA_GZ_OUTPUT=config.get("FASTA_GZ_OUTPUT", ""),
		GC_BIAS=config.get("GC_BIAS", ""),
		output_prefix = "{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR",
		seed = lambda wildcards: int(runs[wildcards.sample_dir][int(wildcards.tech_rep) - 1][0]),
	run:
		command = "iss generate"
		command += f" --fasta_list {input.fasta_list}"
		command += f" --csv_reads_list {input.csv_list}"
		if params.MODE:
			command += f" --mode {params.MODE}"
		if params.ERROR_MODEL:
			command += f" --model {params.ERROR_MODEL}"
		if params.FASTA_GZ_OUTPUT:
			command += f" {params.FASTA_GZ_OUTPUT}"
		if params.NCORE:
			command += f" --cpus {params.NCORE}"
		if params.GC_BIAS:
			command += f" {params.GC_BIAS}"
		command += f" --seed {params.seed}" 
		command += f" --output {params.output_prefix}"
		shell(command)

rule merging_fastas_and_read_csvs:
	input:
		csv_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/csv_list.txt",
		fasta_list="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/fasta_list.txt",
		genereted_reads1 = "{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R1.fastq.gz",
		genereted_reads2 = "{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R2.fastq.gz"
	output:
		merged_fasta="{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/sequenced_frags.fasta.gz",
	shell:
		"""
		/usr/src/app/scripts/wes/merging_sequenced_fasta_and_csv.sh {input.csv_list} {input.fasta_list} {output.merged_fasta}
		rm {input.csv_list} {input.fasta_list}
		"""


rule cleanup_1:
	input:
		merged_fastas = expand("{sample_dir}/tech_replicate_{tech_rep}/PCR_reaction/sequenced_frags.fasta.gz",
									   sample_dir=runs.keys(),
									   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1)),
		generated_reads = expand("{sample_dir}/generated_reads/tech_replicate_{tech_rep}_PCR_R1.fastq.gz",
									   sample_dir=runs.keys(),
									   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1)),
		pcr_filtered_files = expand("{sample_dir}/tech_replicate_{tech_rep}/PCR_filtered.csv",
									sample_dir=runs.keys(),
									tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1)),
		frag_selected_dirs = expand("{sample_dir}/tech_replicate_{tech_rep}",
							sample_dir=runs.keys(), 
							tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1))
	params:
		NCORE=config.get("NCORE", "")
	output:
		touch(os.path.join(config["OUTPUT_DIRECTORY"], "cleanup_1_complete.txt"))
	shell:
		"""
		rm -f {input.pcr_filtered_files}
		for dir in {input.frag_selected_dirs}; do
			find $dir -name '*_fragmentation_selected.csv' -exec rm {{}} +
			if [ -f $dir/PCR_reaction/read_pairs_per_fragment.csv ]; then
				rm -f $dir/PCR_reaction/read_pairs_per_fragment.csv
			fi
		done
		touch {output}
		"""
rule cleanup_2:
	input:
		clean_up_1_log = os.path.join(config["OUTPUT_DIRECTORY"], "cleanup_1_complete.txt"),
		sample_tables = "{sample_dir}/"
	output:
		clean_up_2_log = touch("{sample_dir}/cleanup_2_complete.txt")
	wildcard_constraints:
		sample_dir="|".join(runs.keys())
	params:
		NCORE=config.get("NCORE", 1)
	shell:
		"""
		python3.11 /usr/src/app/scripts/wes/aggregate_mut_parallel_v3.py {input.sample_tables} --cores {params.NCORE}
		python3.11 /usr/src/app/scripts/wes/merging_sample_mut_v2.py {input.sample_tables}
		touch {output.clean_up_2_log}
		"""

rule cleanup_3:
	input:
		clean_up_2_log = expand("{sample_dir}/cleanup_2_complete.txt",
								sample_dir=runs.keys()),
		pcr_files = expand("{sample_dir}/tech_replicate_{tech_rep}",
						   sample_dir=runs.keys(), 
						   tech_rep=range(1, max([len(replicates) for replicates in runs.values()]) + 1))
	output:
		clean_up_3_log = touch(os.path.join(config["OUTPUT_DIRECTORY"], "cleanup_3_complete.txt"))
	wildcard_constraints:
		sample_dir="|".join(runs.keys())
	shell:
		"""
		for dir in {input.pcr_files}; do
			if [ -f $dir/PCR_reaction/seq_mutation_counts.csv ]; then
				python3.11 /usr/src/app/scripts/wes/seq_mut_filter.py $dir/PCR_reaction/seq_mutation_counts.csv
			fi
		done

		for dir in {input.pcr_files}; do
			if [ -f $dir/PCR_reaction/PCR_filtered_amp_fragments.csv ]; then
				rm -f $dir/PCR_reaction/PCR_filtered_amp_fragments.csv
			fi
		done

		touch {output.clean_up_3_log}
		"""