2_2_simulation.sm

localrules: all_simulate_data, all_apply_ggm_simulation, all_simulation, all_simulation2, all_simulation3,
            validate_all, validate_minimal, validate_all_subsets

#------------------------------------------------------------------------------
# Simulate ground truth and data for simulation study
#------------------------------------------------------------------------------
RUNS = 100
ITERATIONS = range(1,RUNS+1)
ITERATIONS_TEST = range(1,51)
SUBSETS = [str(n) for n in range(50,650,50)]
#SUBSETS = [str(n) for n in range(10,210,10)]
#SUBSETS.append("all")

# for the rs60626639 sentinel, simulation seems to take too long at the moment, so we should recalculate
# on a specialiced run (e.g. more runtime or threads)
# This is mostly an issue for the "all" subset, as we do 11 inference cycles (for each noisey matrix).
MEQTL_SIM = MEQTL.sentinel
MEQTL_SIM.remove("rs60626639")

rule simulate_data:
        input:
                data=DCOHORT_DATA + "lolipop/{sentinel}_meqtl.rds",
                ranges=DRANGES + "{sentinel}_meqtl.rds",
                priors=DPRIORS + "{sentinel}_meqtl.rds"
        output:
                DSIM_DATA + "{sentinel}.RData"
        threads: 10
        resources:
                mem_mb=15000
        params:
                sentinel="{sentinel}",
                runs=RUNS,
                time="1-00:00:00"
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/simulate_data/{sentinel}.log"
        benchmark:
                "benchmarks/simulation/simulate_data/{sentinel}.bmk"
        script:
                "../scripts/simulation/simulate_data.R"

rule all_simulate_data:
	input:
		expand(DSIM_DATA + "{sentinel}.RData", sentinel=MEQTL.sentinel)

#------------------------------------------------------------------------------
# Apply ggm on simulated data
#------------------------------------------------------------------------------
rule apply_ggm_simulation:
        input:
                data=DSIM_DATA + "{sentinel}.RData",
                ppi_db=PPI_DB,
		cpg_context = "results/current/cpg_context.rds"
        output:
                DSIM_FITS + "{sentinel}-iter{iteration}-subset{subset}.RData"
        params:
                iteration="{iteration}",
                time="2-00:00:00",
		minimal_subset_size = 70
        threads: 16
        wildcard_constraints:
                subset = "all|minimal"
        resources:
                mem_mb=40000
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/apply_ggm/{sentinel}-iter{iteration}-subset{subset}.log"
        benchmark:
                "benchmarks/simulation/apply_ggm/{sentinel}-iter{iteration}-subset{subset}.bmk"
        script:
                "../scripts/simulation/run_ggm.R"

rule apply_ggm_simulation_subsets:
        input:
                data=DSIM_DATA + "{sentinel}.RData",
                ppi_db=PPI_DB,
		cpg_context = "results/current/cpg_context.rds"
        output:
                DSIM_FITS + "subset{subset}/{sentinel}-iter{iteration}.RData"
        params:
                iteration="{iteration}",
                time="2-00:00:00",
		minimal_subset_size = 70
        threads: 16
        wildcard_constraints:
                subset = "[0-9]+0*"
        resources:
                mem_mb=40000
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/apply_ggm/subset{subset}/{sentinel}-iter{iteration}.log"
        benchmark:
                "benchmarks/simulation/apply_ggm/subset{subset}/{sentinel}-iter{iteration}.bmk"
        script:
                "../scripts/simulation/run_ggm.R"

rule apply_ggm_simulation_prior_completeness:
        input:
                data=DSIM_DATA + "{sentinel}.RData",
                ppi_db=PPI_DB,
		cpg_context = "results/current/cpg_context.rds"
        output:
                DSIM_FITS + "prior_completeness/{sentinel}-iter{iteration}.RData"
        params:
                iteration="{iteration}",
                time="2-00:00:00"
        threads: 16
        resources:
                mem_mb=40000
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/apply_ggm/prior_completeness/{sentinel}-iter{iteration}.log"
        benchmark:
                "benchmarks/simulation/apply_ggm/prior_completeness/{sentinel}-iter{iteration}.bmk"
        script:
                "../scripts/simulation/run_ggm_prior_completeness.R"


rule all_apply_ggm_simulation_subsetall:
	input:
		expand(DSIM_FITS + "{sentinel}-iter{iteration}-subsetall.RData", sentinel=MEQTL_SIM, iteration=ITERATIONS)

rule all_apply_ggm_simulation_subsets:
	input:
		expand(DSIM_FITS + "subset{subset}/{sentinel}-iter{iter}.RData", sentinel=MEQTL_SIM, iter=ITERATIONS_TEST, subset=SUBSETS)

rule all_apply_ggm_simulation_subsetminimal:
	input:
		expand(DSIM_FITS + "{sentinel}-iter{iter}-subsetminimal.RData", sentinel=MEQTL_SIM, iter=ITERATIONS,)

rule all_apply_ggm_simulation_prior_completeness:
	input:
		expand(DSIM_FITS + "prior_completeness/{sentinel}-iter{iter}.RData", sentinel=MEQTL_SIM, iter=ITERATIONS_TEST)

rule apply_ggm_simulation_iterations:
	input:
		expand(DSIM_FITS + "{sentinel}-iter{iter}-subsetall.RData", iter=ITERATIONS, sentinel=MEQTL_SIM),
		expand(DSIM_FITS + "{sentinel}-iter{iter}-subsetminimal.RData", iter=ITERATIONS, sentinel=MEQTL_SIM)

#------------------------------------------------------------------------------
# Validate a simulation run for the 'all' subset
#------------------------------------------------------------------------------
rule validate_ggm_simulation_subsets:
        input:
                fits=expand(DSIM_FITS + "subset{{subset}}/{{sentinel}}-iter{iteration}.RData", iteration=ITERATIONS_TEST)
        output:
                DSIM_VALIDATION + "subset{subset}/{sentinel}.txt"
        threads: 1
        wildcard_constraints:
                subset = "[0-9]+0*"
        resources:
                mem_mb=1000
        params:
                time="00:10:00"
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/validate_ggm/{sentinel}_subset{subset}.log"
        script:
                "../scripts/simulation/validate.R"

rule validate_ggm_simulation:
        input:
                fits=expand(DSIM_FITS + "{{sentinel}}-iter{iteration}-subset{{subset}}.RData", iteration=ITERATIONS)
        output:
                DSIM_VALIDATION + "{sentinel}-subset{subset}.txt"
        threads: 1
        wildcard_constraints:
                subset = "all|minimal"
        resources:
                mem_mb=1000
        params:
                time="00:10:00"
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/validate_ggm/{sentinel}_subset{subset}.log"
        script:
                "../scripts/simulation/validate.R"

rule validate_ggm_simulation_prior_completeness:
        input:
                fits=expand(DSIM_FITS + "prior_completeness/{{sentinel}}-iter{iteration}.RData", iteration=ITERATIONS_TEST)
        output:
                DSIM_VALIDATION + "prior_completeness/{sentinel}.txt"
        threads: 1
        resources:
                mem_mb=1000
        params:
                time="00:10:00"
        conda:
                "../envs/bioR.yaml"
        log:
                "logs/simulation/validate_ggm/prior_completeness/{sentinel}.log"
        script:
                "../scripts/simulation/validate.R"

#------------------------------------------------------------------------------
# 3 Target rule to validate all simulation runs: all, minimal and subsets
#------------------------------------------------------------------------------
rule validate_all:
	input: expand(DSIM_VALIDATION + "{sentinel}-subsetall.txt", sentinel=MEQTL_SIM)
	output:
		DSIM + "validation-subsetall.txt"		
	resources:
		mem_mb = 1000
	threads: 16
	log: "logs/validate_subsetall.log"
	shell:
		"""
		sort --parallel={threads} -r {input} | uniq > {output}
		"""

rule validate_minimal:
	input: 
		results = expand(DSIM_VALIDATION + "{sentinel}-subsetminimal.txt", sentinel=MEQTL_SIM)
	output:
		combined = DSIM + "validation-subsetminimal.txt"
	resources:
		mem_mb = 1000
	threads: 1
	log: "logs/validate_subsetminimal.log"
	script:
		"../scripts/simulation/reorder_columns.R"

rule validate_all_subsets:
	input: 
		results = expand(DSIM_VALIDATION + "subset{subset}/{sentinel}.txt", sentinel=MEQTL_SIM, subset=SUBSETS)
	output:
		combined = DSIM + "validation-subsets.txt"		
	resources:
		mem_mb = 1000
	threads: 1
	log: "logs/validate_subsets.log"
	script:
		"../scripts/simulation/reorder_columns.R"

rule validate_all_prior_completeness:
	input: expand(DSIM_VALIDATION + "prior_completeness/{sentinel}.txt", sentinel=MEQTL_SIM)
	output:
		DSIM + "validation-prior-completeness.txt"
	resources:
		mem_mb = 2000
	threads: 16
	log: "logs/validate_prior_completeness.log"
	shell:
		"""
		sort --parallel={threads} -r {input} | uniq > {output}
		"""

#------------------------------------------------------------------------------
# 2 rules to check recovery of mixed edges (SNP-gene)
#------------------------------------------------------------------------------
rule summarize_snp_gene_recovery:
	output:
		summary = DSIM_VALIDATION + "snp_gene_recovery/{sentinel}.tsv"
	params:
		dresults = DSIM_FITS,
		time = "01:00:00"
	log:
		"logs/simulation/snp_gene_recovery/{sentinel}.log"
	benchmark:
		"benchmarks/simulation/snp_gene_recovery/{sentinel}.bmk"
	threads: 6
	resources:
		mem_mb = 6000
	script:
		"../scripts/simulation/summarize_snp_gene_recovery.R"

rule all_summarize_snp_gene_recovery:
	input:
		expand(DSIM_VALIDATION + "snp_gene_recovery/{sentinel}.tsv", sentinel = MEQTL.sentinel)
	output:
		DSIM_VALIDATION + "snp_gene_recovery/all.tsv"
	shell:
		"""
		cat {input} | sort -r | uniq > {output}
		"""

# -----------------------------------------------------------------------------
# Overall target rule to run everything regarding loci based simulation
# -----------------------------------------------------------------------------
rule all_simulation:
	input:
		DSIM_VALIDATION + "snp_gene_recovery/all.tsv",
		DSIM + "validation-subsets.txt",
		DSIM + "validation-subsetall.txt",
		DSIM + "validation-subsetminimal.txt"

# -----------------------------------------------------------------------------
# Runtime benchmarking using simulated data
# -----------------------------------------------------------------------------
BENCHMARK_ITERATIONS = 50
BENCHMARK_SAMPLE_SIZES = [200, 400, 600]
BENCHMARK_NUMBER_OF_NODES = [30, 60, 90, 120]

rule benchmark_model:
	output:
		result_table = "results/current/benchmark/result_table_{model}_N{sample_size}_P{number_nodes}.rds"
	params:
		benchmark_number_iterations = BENCHMARK_ITERATIONS,
		time = "3-00:00:00"
	log:
		"logs/benchmark/{model}_N{sample_size}_P{number_nodes}.log"
	benchmark:
		"benchmarks/benchmark/{model}_N{sample_size}_P{number_nodes}.bmk"
	resources:
		mem_mb = 3000
	threads: 1
	script:
		"../scripts/benchmark.R"

MODELS = ["glasso", "genenet", "bdgraph", "genie3", "irafnet"]
rule all_benchmark:
	input: 
		expand("results/current/benchmark/result_table_{model}_N{sample_size}_P{number_nodes}.rds", model = MODELS, sample_size = BENCHMARK_SAMPLE_SIZES, number_nodes = BENCHMARK_NUMBER_OF_NODES)
	output:
		overview_plot = "results/current/benchmark/summary.pdf"
	params:
		time = "00:10:00",
		benchmark_number_iterations = BENCHMARK_ITERATIONS,
	resources:
		mem_mb = 1000
	threads: 1
	log:
		"logs/benchmark/all_benchmark.log"
	script:
		"../scripts/benchmark_gather_results.R"