Snakemake工作流添加重复样本合并功能：多通配符使用与执行报错排查

针对你遇到的样本读取错误、expand生成重复路径，以及后续要添加的重复合并需求，我整理了分步排查和修复方案：

一、核心错误根源：样本表格读取与expand用法错误 #

你之前的两个关键问题：

1. 设置index_col="sample"后，samples.sample调用的是DataFrame的随机抽样方法（而非sample列），导致rule all生成了奇怪的文件名；
2. expand单独传入三个列的取值，触发了笛卡尔积组合，导致重复生成路径。

二、分步修复现有规则 #

1. 正确读取样本表格 #

推荐使用复合索引（因为sample + bio_unit + tech_unit是唯一标识），既高效又避免列名混淆：

import pandas as pd
import os
import yaml
configfile: "config.yaml"

# 用三个字段作为复合索引，唯一标识每个技术重复
samples = pd.read_table(config["samples"], index_col=["sample", "bio_unit", "tech_unit"])

2. 修复rule all的expand调用 #

要按表格行的对应关系生成路径，需用zip参数避免笛卡尔积：

# 从复合索引中拆分出三个字段的列表
sample_list, bio_list, tech_list = zip(*samples.index)

rule all:
    input:
        expand(
            config["arima_mapping"] + "{sample}_{unit_bio}_{unit_tech}_R1.bam",
            zip,  # 按行打包三个字段，避免笛卡尔积
            sample=sample_list,
            unit_bio=bio_list,
            unit_tech=tech_list
        )

3. 修复r1_mapping的输入读取 #

利用复合索引直接定位行，代码更简洁且不易出错：

rule r1_mapping:
    input:
        # 通过复合索引直接获取对应样本的fq1路径
        read1 = lambda wc: samples.loc[(wc.sample, wc.unit_bio, wc.unit_tech), "fq1"],
        ref = config["PacBio_assembly"],
        linker = config["index"] + "pac_bio_assembly.fna.amb"
    output: config["arima_mapping"] + "unprocessed_bam/({sample},{unit_bio},{unit_tech})_R1.bam"
    params:
    conda: "../envs/arima_mapping.yaml"
    log: config["logs"] + "arima_mapping/mapping/({sample},{unit_bio},{unit_tech})_R1.log"
    threads: 12
    shell: "bwa mem -t {threads} {input.ref} {input.read1} | samtools view --threads {threads} -h -b -o {output} 2> {log}"

三、添加合并重复的规则示例 #

等现有规则正常运行后，你可以添加合并技术/生物学重复的规则，以下是合并同一sample + bio_unit下技术重复的示例：

1. 定义获取待合并文件的函数 #

def get_tech_replicates(wc):
    # 筛选同一sample+bio_unit下的所有tech_unit
    tech_units = samples.loc[(wc.sample, wc.unit_bio), :].index.get_level_values("tech_unit")
    return expand(
        config["arima_mapping"] + "unprocessed_bam/({sample},{unit_bio},{unit_tech})_R1.bam",
        unit_tech=tech_units,
        sample=wc.sample,
        unit_bio=wc.unit_bio
    )

2. 编写合并规则 #

rule merge_tech_replicates:
    input: get_tech_replicates
    output: config["arima_mapping"] + "merged_bam/{sample}_{unit_bio}_merged_R1.bam"
    conda: "../envs/samtools.yaml"
    threads: 8
    shell: "samtools merge -@ {threads} {output} {input}"

3. 更新rule all #

记得把合并后的文件加入rule all的输入，确保Snakemake会触发合并规则：

rule all:
    input:
        expand(
            config["arima_mapping"] + "{sample}_{unit_bio}_{unit_tech}_R1.bam",
            zip,
            sample=sample_list,
            unit_bio=bio_list,
            unit_tech=tech_list
        ),
        # 添加合并后的文件
        expand(
            config["arima_mapping"] + "merged_bam/{sample}_{unit_bio}_merged_R1.bam",
            zip,
            sample=samples.index.get_level_values("sample").unique(),
            unit_bio=samples.index.get_level_values("bio_unit").unique()
        )

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