Snakemake流水线中根据配置参数选择Shell命令的优化方法

原问题与现有实现 #

需求是根据配置文件中是否存在excluded_regions参数，执行不同的bamCoverage命令。现有实现如下：

rule deeptools_bamCoverage_pe:
    input:
        bam="DATA/BAM/{sample_sp}_pe.bam",
        bai="DATA/BAM/{sample_sp}_pe.bam.bai"
    output:
        "DATA/BIGWIG/{sample_sp}_pe_RPGC.bw"
    log:
        "snakemake_logs/deeptools_bamCoverage/{sample_sp}_pe.log"
    run:
        if config["excluded_regions"]:
            shell("bamCoverage --normalizeUsing RPGC -bl "+config["excluded_regions"]+
                    " --effectiveGenomeSize $((2913022398-"+config["excluded_regions"].split(".")[0].split("_")[-1]+")) -e -b {input.bam} -o {output} 2>{log}")
        else:
            shell("bamCoverage --normalizeUsing RPGC --effectiveGenomeSize 2913022398 -e -b {input.bam} -o {output} 2>{log}")

配置文件中excluded_regions参数为类似DATA/GENOMES/DAC_excl_HG38_71570285.bed的路径，其中71570285是排除区域总长度，用于计算有效基因组大小。

更简洁的实现方式 #

方法1：利用Snakemake参数格式化（推荐） #

通过params模块提前定义动态参数，结合Snakemake的条件格式化语法，无需拆分shell命令：

rule deeptools_bamCoverage_pe:
    input:
        bam="DATA/BAM/{sample_sp}_pe.bam",
        bai="DATA/BAM/{sample_sp}_pe.bam.bai"
    output:
        "DATA/BIGWIG/{sample_sp}_pe_RPGC.bw"
    log:
        "snakemake_logs/deeptools_bamCoverage/{sample_sp}_pe.log"
    params:
        # 若配置存在excluded_regions则生成-bl参数，否则为空
        blacklist = config.get("excluded_regions", ""),
        # 计算有效基因组大小：存在配置时用总大小减去排除区域长度，否则用默认值
        effective_genome_size = (
            f"$((2913022398-{config['excluded_regions'].split('.')[0].split('_')[-1]}))"
            if "excluded_regions" in config
            else "2913022398"
        )
    shell:
        """
        bamCoverage --normalizeUsing RPGC {params.blacklist: -bl {}} \
            --effectiveGenomeSize {params.effective_genome_size} -e -b {input.bam} -o {output} 2>{log}
        """

这里{params.blacklist: -bl {}}是Snakemake的条件格式化语法：当params.blacklist非空时，输出-bl 参数值；为空时则不输出该部分，自动适配两种场景。

方法2：简化run块的字符串拼接 #

如果偏好使用run块，可提前拼接动态参数，只写一次核心shell命令：

rule deeptools_bamCoverage_pe:
    input:
        bam="DATA/BAM/{sample_sp}_pe.bam",
        bai="DATA/BAM/{sample_sp}_pe.bam.bai"
    output:
        "DATA/BIGWIG/{sample_sp}_pe_RPGC.bw"
    log:
        "snakemake_logs/deeptools_bamCoverage/{sample_sp}_pe.log"
    run:
        # 生成黑名单参数
        blacklist_arg = f"-bl {config['excluded_regions']}" if "excluded_regions" in config else ""
        # 计算有效基因组大小
        if "excluded_regions" in config:
            excluded_size = config["excluded_regions"].split(".")[0].split("_")[-1]
            genome_size = f"$((2913022398-{excluded_size}))"
        else:
            genome_size = "2913022398"
        # 统一执行shell命令
        shell(
            f"bamCoverage --normalizeUsing RPGC {blacklist_arg} "
            f"--effectiveGenomeSize {genome_size} -e -b {input.bam} -o {output} 2>{log}"
        )

额外优化建议 #

从文件名中提取排除区域长度的做法易出错（比如文件名格式变动），建议直接将长度写入配置文件：

# config.yaml
excluded_regions: DATA/GENOMES/DAC_excl_HG38_71570285.bed
excluded_regions_size: 71570285

之后在rule中直接引用config['excluded_regions_size']，替换字符串分割的逻辑，提升可靠性。

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