Snakemake配置样本组合问题：指定特定样本配对合并

问题场景 #

配置文件中样本按分组和SRA编号组织：

samples:
    group1:
        sra1: 
            sample: "SRR14724462"
            cell_line: "NA24385"
            exome_bedfile: "/bedfiles/truseq.sorted.bed"
        sra2: 
            sample: "SRR14724472"
            cell_line: "NA24385"
            exome_bedfile: "/bedfiles/idt.sorted.bed"
    group2:
        sra1: 
            sample: "SRR14724463"
            cell_line: "NA12878"
            exome_bedfile: "/bedfiles/truseq.sorted.bed"
        sra2:
            sample: "SRR14724473"
            cell_line: "NA12878"
            exome_bedfile: "/bedfiles/idt.sorted.bed"

需求是按SRA编号配对合并：

• group1的sra1 ↔ group2的sra1（SRR14724462 和 SRR14724463）
• group1的sra2 ↔ group2的sra2（SRR14724472 和 SRR14724473）

但当前rule all使用expand时，默认生成样本列表的笛卡尔积，导致出现不需要的组合。

现有代码问题 #

当前规则生成了所有样本的笛卡尔积：

rule combine:       
    output:
        r1 = TRIMMED_DIR + "/{sample1}_{sample2}_R1.fastq",
        r2 = TRIMMED_DIR + "/{sample1}_{sample2}_R2.fastq"
    params:
        trimmed_dir = TRIMMED_DIR,
        a = "{sample1}", 
        b = "{sample2}"
    shell: 
        cd {params.trimmed_dir}
        /combine.sh {params.a}_R1_trimmed.fastq {params.a}_R2_trimmed.fastq {params.b}_R1_trimmed.fastq {params.b}_R2_trimmed.fastq

rule all:
     expand(TRIMMED_DIR + "/{sample1}_{sample2}_R1.fastq", sample1=list_a, sample2=list_b),
     expand(TRIMMED_DIR + "/{sample1}_{sample2}_R2.fastq", sample1=list_a, sample2=list_b)

解决方案 #

方法1：直接指定配对组合 #

如果配对关系固定，可直接构造配对元组列表，再生成目标文件路径：

# 定义需要的配对
PAIRS = [("SRR14724462", "SRR14724463"), ("SRR14724472", "SRR14724473")]

rule all:
    input:
        [TRIMMED_DIR + f"/{s1}_{s2}_R1.fastq" for s1, s2 in PAIRS],
        [TRIMMED_DIR + f"/{s1}_{s2}_R2.fastq" for s1, s2 in PAIRS]

或者利用expand的zip参数，按列表顺序配对：

rule all:
    input:
        expand(TRIMMED_DIR + "/{sample1}_{sample2}_R1.fastq", zip, sample1=list_a, sample2=list_b),
        expand(TRIMMED_DIR + "/{sample1}_{sample2}_R2.fastq", zip, sample1=list_a, sample2=list_b)

注：zip参数要求sample1和sample2对应的列表长度一致，会按索引一一配对。

方法2：从配置文件动态提取配对（推荐） #

如果配置文件结构固定，可通过遍历SRA键动态生成配对，避免硬编码样本ID：

# 从config中自动提取配对
PAIRS = []
for sra_key in config["samples"]["group1"].keys():
    sample1 = config["samples"]["group1"][sra_key]["sample"]
    sample2 = config["samples"]["group2"][sra_key]["sample"]
    PAIRS.append( (sample1, sample2) )

rule all:
    input:
        [TRIMMED_DIR + f"/{s1}_{s2}_R1.fastq" for s1, s2 in PAIRS],
        [TRIMMED_DIR + f"/{s1}_{s2}_R2.fastq" for s1, s2 in PAIRS]

这种方式的优势是：后续新增SRA样本（如sra3）时，无需修改代码，只要配置文件保持对应结构，就能自动生成正确配对。

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