Версия Snakemake: 6.3.0
Я хочу применить DAG в глубину к моему змейке для правил, создающих одинаковое количество файлов, вплоть до правила, объединяющего несколько файлов в соответствии с подстановочный знак.
На данный момент я применил совет, данный в этой теме: Snakemake: Tranverse DAG в глубину?
Я расставил приоритеты шаги, чтобы мой змейка находился в режиме глубины до правила sort_and_index_binning. Однако в конце этапа сортировки файлы по-прежнему довольно велики, и мне хотелось бы достичь этапа, на котором суммируется глубина (summarize_contig_depth) как можно быстрее для каждого образца, чтобы у меня не осталось временных файлов из предыдущих шагов. В конце концов, идея заключалась бы в том, чтобы запускать образец с подстановочным знаком src1, но я не понимаю, как это сделать автоматически.
Например, с двумя образцами sample_1 и sample_2 Я хотел бы сделать:
binning_mapping: sample_1_to_sample_2, sample_2_to_sample_2 (здесь sample_1 и sample_2 — это подстановочные знаки src, а sample_2 — это подстановочные знаки src1)filter_bam : sample_1_to_sample_2, sample_2_to_sample_2
sort_and_index_binning : sample_1_to_sample_2, sample_2_to_sample_2
summarize_contig_length : deep_sample_2
THEN
binning_mapping : sample_1_to_sample_1, sample_2_to_sample_1 (здесь sample_1 и sample_2 — это подстановочные знаки src, а sample_1 — это подстановочные знаки src1)
filter_bam : sample_1_to_sample_1, sample_2_to_sample_1
sort_and_index_binning : sample_1_to_sample_1, sample_2_to_sample_1
summarize_contig_length : deep_sample_1
THEN< /p>
правило с использованием глубины_sample_1 и глубины_sample_2
def input_cmd(wildcards):
if wildcards.assembly == "single_assembly":
list_reads = []
for run in reads2use[wildcards.src]:
list_reads.extend(reads2use[wildcards.src][run])
return list_reads
elif wildcards.assembly == "co_assembly":
if simka_type is "None":
return os.path.join(tmpdir, 'samples.txt')
return os.path.join(intermediate_results_dir, "assembly/co_assembly/clusters/{src}.txt")
else:
raise ValueError
rule binning_mapping:
'''
Align the source reads files against the assembled contigs file to assess contigs' abundance.
'''
output:
temp(os.path.join(tmpdir, "{assembly}/{src}_to_{src1}_" + f"{index}_bin_filtering.sam"))
input:
assembly = os.path.join(intermediate_results_dir, "assembly/{assembly}", assembler, "{src1}", f"contigs/{assembly}"),
index1 = expand(os.path.join(intermediate_results_dir, "assembly/{{assembly}}", assembler, "{{src1}}/index/{{src1}}_" + index + "_filtering.{id}.bt2l"), id=range(1, 4)),
index2 = expand(os.path.join(intermediate_results_dir, "assembly/{{assembly}}", assembler, "{{src1}}/index/{{src1}}_" + index + "_filtering.rev.{id}.bt2l"), id=range(1,2)),
reads = input_cmd,
finished_assembly = os.path.join(tmp, "assembly.checkpoint")
params:
prefix = os.path.join("intermediate_results/assembly/{assembly}", assembler, "{src1}/index/{src1}_" + f"{index}_filtering"),
input_reads = lambda wildcards, input : cmdparser.cmd(wildcards.src, input.reads, reads2use, "bowtie2").cmd,
cmd = lambda wildcards : conf.mapping_cmd(config, wildcards.assembly),
threads: 5
priority: 1
conda:
os.path.join(CONDAENV, "bowtie2.yaml")
shell:
"bowtie2 "
"-p {threads} "
"--no-unal "
"-x {params.prefix} "
"{params.input_reads} "
"{params.cmd} "
"-S {output} "
rule filter_bam:
"""
Filter reads based on mapping quality and identity.
Output is temporary because it will be sorted.
"""
output:
temp(os.path.join(intermediate_results_dir, "assembly/{assembly}", assembler, "{src1}/mapped_reads/{src}.filtered.sam")),
input:
os.path.join(tmpdir, "{assembly}/{src}_to_{src1}_" + f"{index}_bin_filtering.sam")
conda:
os.path.join(CONDAENV, "bamutils.yaml")
priority: 2
params:
min_mapq = config["bam_filtering_before_binning"]["min_quality"],
min_idt = config["bam_filtering_before_binning"]["min_identity"],
min_len = config["bam_filtering_before_binning"]["min_len"],
pp = config["bam_filtering_before_binning"]["properly_paired"],
script:
"../scripts/bamprocess.py"
rule sort_and_index_binning:
output:
temp(os.path.join(intermediate_results_dir, "assembly/{assembly}", assembler, "{src1}/mapped_reads/{src}_to_{src1}.sorted.bam"))
input:
os.path.join(intermediate_results_dir, "assembly/{assembly}", assembler, "{src1}/mapped_reads/{src}.filtered.sam"),
threads: 1
priority: 3
conda:
os.path.join(CONDAENV, "samtools.yaml")
shell:
"samtools view -u {input} | "
"samtools sort "
"-@ {threads} "
"-o {output[0]} "
def aggregate_bam_input(wildcards):
if "CASB" in strategies or "CACB" in strategies:
checkpoint_output_simka = checkpoints.cluster_simka.get(**wildcards).output[0]
assembly_dict["co_assembly"] = glob_wildcards(os.path.join(checkpoint_output_simka, "{clusterid}.txt")).clusterid
assembly_request = "co_assembly"
if "SASB" in strategies or "SACB" in strategies:
assembly_dict["single_assembly"] = list(samples.keys())
assembly_request = "single_assembly"
inputs = expand(os.path.join(intermediate_results_dir,
"assembly",
assembly_request,
assembler,
wildcards.src,
"mapped_reads",
"{src1}_to_" + wildcards.src + ".sorted.bam"
), src1=assembly_dict.get(assembly_request))
rule summarize_contig_depth:
'''
Compute reads coverage depth to perform binning.
'''
output:
os.path.join(intermediate_results_dir, "binning/{binning_strategy}/{src1}/depth.txt"),
input:
aggregate_bam_input,
params:
bams = aggregate_bam_input,
conda:
os.path.join(CONDAENV, "metabat.yaml")
threads: 5
priority: 4
shell:
"jgi_summarize_bam_contig_depths --outputDepth {output} {params.bams}"
Подробнее здесь: https://stackoverflow.com/questions/789 ... at-the-end
