Snakemake 6.3.0：如何按样本自动深度遍历DAG并完成聚合？

需求说明 #

使用Snakemake 6.3.0，要求对生成相近数量文件的规则采用深度优先遍历DAG，直至执行到按通配符wildcard聚合多文件的规则。

已通过设置步骤优先级，让Snakemake在sort_and_index_binning规则前以深度模式运行，但排序后的文件体积仍较大，因此需尽快完成每个样本的深度聚合步骤（summarize_contig_depth），以清理前期步骤的临时文件。最终目标是自动按通配符src1逐个样本执行流程。

期望执行流程（以2个样本sample_1和sample_2为例） #

• 第一阶段（src1=sample_2）：
  • binning_mapping：执行sample_1_to_sample_2、sample_2_to_sample_2（src为sample_1和sample_2）
  • 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_depth：生成depth_sample_2
• 第二阶段（src1=sample_1）：
  • binning_mapping：执行sample_1_to_sample_1、sample_2_to_sample_1（src为sample_1和sample_2）
  • 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_depth：生成depth_sample_1
• 第三阶段：调用依赖depth_sample_1和depth_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:
    '''
        将源reads文件比对到组装的contigs文件，评估contigs的丰度。
    '''
    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:
    """
    根据比对质量和一致性过滤reads。
    输出为临时文件，因为后续会被排序。
    """
    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:
    '''
        计算reads覆盖深度，用于后续分箱。
    '''
    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}"

内容的提问来源于stack exchange，提问作者Hugo Lefeuvre