Snakemake如何定义规则合并单个通配符对应的所有输入BAM文件

核心解决思路 #

使用Snakemake的动态输入函数，根据当前匹配到的样本通配符值，自动拉取该样本对应的所有测序批次BAM文件，无需硬编码批次信息。

具体修改步骤 #

1. 定义全局公共参数（可选，推荐） #

在流程开头统一配置比对工具、参考基因组参数，避免多处硬编码导致不一致：

# Define samples and runs
RUNS, SAMPLES = glob_wildcards("/labs/jandr/walter/tb/data/Stanford/{run}/{samp}_L001_R1_001.fastq.gz")
# 新增公共参数
MAPPER = "bwa"
REF = "MTB_ancestor_reference"
print("runs are: ", RUNS)
print("samples are: ", SAMPLES)

对应修改rule all中的合并BAM路径生成逻辑：

rule all:
  input:
     trim = expand(['process/trim/{run}_{samp}_trim_1.fq.gz'], zip, run = RUNS, samp = SAMPLES),
     kraken=expand('process/trim/{run}_{samp}_trim_kr_1.fq.gz', zip, run = RUNS, samp = SAMPLES),
     bams=expand('process/bams/{run}_{samp}_bwa_MTB_ancestor_reference_rg_sorted.bam', zip, run = RUNS, samp = SAMPLES),
     # 修改为使用公共参数展开
     combined_bams=expand('process/bams/{samp}_{mapper}_{ref}.merged.rmdup.bam', samp = np.unique(SAMPLES), mapper=MAPPER, ref=REF)

2. 重写combine_bams规则 #

替换原有combine_bams规则，使用动态输入匹配多批次BAM，同时修正原shell命令逻辑错误（sambamba markdup仅支持单个输入BAM，多文件需先合并再去重）：

# 按样本合并所有批次BAM并去重
rule combine_bams:
  input:
    # 动态匹配当前样本对应的所有批次BAM
    bams = lambda wildcards: expand(
        'process/bams/{run}_{samp}_{mapper}_{ref}_rg_sorted.bam',
        run = RUNS,
        samp = wildcards.samp,
        mapper = MAPPER,
        ref = REF
    )
  output: 
    combined_bam = 'process/bams/{samp}_{mapper}_{ref}.merged.rmdup.bam'
  log: 
     'process/bams/{samp}_{mapper}_{ref}_merge_bams.log'
  threads: 8
  shell:
    """
    # 先合并多批次BAM
    sambamba merge -t {threads} /tmp/{wildcards.samp}_merged.tmp.bam {input.bams}
    # 再去重
    sambamba markdup -r -p -t {threads} /tmp/{wildcards.samp}_merged.tmp.bam {output.combined_bam}
    # 清理临时文件
    rm -f /tmp/{wildcards.samp}_merged.tmp.bam
    """

可选适配场景 #

如果存在部分批次没有对应样本测序数据的情况，可以把动态输入改成实时匹配存在的文件，避免生成不存在的输入路径：

bams = lambda wildcards: expand(
    'process/bams/{run}_{samp}_{mapper}_{ref}_rg_sorted.bam',
    # 仅匹配该样本实际存在的run
    run = glob_wildcards(f'process/bams/{{run}}_{wildcards.samp}_{MAPPER}_{REF}_rg_sorted.bam').run,
    samp = wildcards.samp,
    mapper = MAPPER,
    ref = REF
)

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