Snakemake工作流试运行Python索引错误求助

环境信息 #

• Snakemake版本：7.32.4
• Python版本：3.10.14
• 运行系统：Linux Ubuntu

问题描述 #

运行Snakemake工作流时触发KeyError：('ShS-DIF4-rep1', 'FC2')，该样本实际仅与units表中的'FC1'单元对应，但试运行时Snakemake尝试匹配不存在的组合ShS-DIF4-rep1 + FC2，完整错误栈如下：

InputFunctionException in rule cutAdapt_pe in file /home/x/snakemake_longRNAseq/workflow/rules/trim.smk, line 1:
Error:
  KeyError: ('ShS-DIF4-rep1', 'FC2')
Wildcards:
  sample=ShS-DIF4-rep1
  unit=FC2
Traceback:
  File "/home/x/snakemake_longRNAseq/workflow/rules/common.smk", line 27, in cutAdapt_rawFastq
  File "/home/x/snakemake_longRNAseq/workflow/rules/common.smk", line 23, in single_end
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1184, in __getitem__
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1368, in _getitem_tuple
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1041, in _getitem_lowerdim
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1153, in _getitem_nested_tuple
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1431, in _getitem_axis
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexing.py", line 1381, in _get_label
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/generic.py", line 4293, in xs
  File "/home/x/mambaforge/envs/longRNAseq/lib/python3.11/site-packages/pandas/core/indexes/multi.py", line 3222, in _get_loc_level (rule cutAdapt_pe, line 1, /home/x/snakemake_longRNAseq/workflow/rules/trim.smk)

问题根因 #

1. rule all中的expand命令错误生成了所有样本与所有单元的笛卡尔积，而非units表中实际存在的(sample, unit)配对组合，导致Snakemake尝试访问不存在的索引对。
2. cutAdapt_se规则的shell命令存在变量错误：输出引用了不存在的output.forwardPaired，实际输出变量名为single_end。

解决方案与代码修正 #

1. 修正rule all的expand逻辑 #

直接基于units表的索引生成有效配对，避免无效的笛卡尔积：

# Snakefile中的rule all修改后
rule all:
    input:
        # 生成配对端结果路径：仅保留units表中fq2非空的有效组合
        expand(
            "results/cutAdapt/{sample}-{unit}.{read}.fastq.gz",
            sample=units.index.get_level_values("sample"),
            unit=units.index.get_level_values("unit"),
            read=["1", "2"],
            condition=lambda sample, unit: not pd.isnull(units.loc[(sample, unit), "fq2"])
        ),
        # 生成单端结果路径：仅保留units表中fq2为空的有效组合
        expand(
            "results/cutAdapt/{sample}-{unit}.fastq.gz",
            sample=units.index.get_level_values("sample"),
            unit=units.index.get_level_values("unit"),
            condition=lambda sample, unit: pd.isnull(units.loc[(sample, unit), "fq2"])
        )

2. 修正cutAdapt_se规则的shell命令 #

将输出变量改为正确的output.single_end：

# trim.smk中的cutAdapt_se修改后
rule cutAdapt_se:
    input:
        cutAdapt_rawFastq
    output:
        single_end="results/cutAdapt/{sample}-{unit}.fastq.gz"
    threads:
        4
    shell:
        """
        cutadapt -q 20 -m 30 --quality-base=33 -j {threads} -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o {output.single_end} {input}
        """

3. 可选优化：提前拆分单端/配对端单元 #

在common.smk中提前拆分单端和配对端单元集合，进一步简化rule all逻辑：

# common.smk中新增代码
# 拆分单端和配对端的(sample, unit)索引
se_units = units[pd.isnull(units["fq2"])].index
pe_units = units[~pd.isnull(units["fq2"])].index

# 简化后的rule all
rule all:
    input:
        expand(
            "results/cutAdapt/{sample}-{unit}.{read}.fastq.gz",
            sample=pe_units.get_level_values("sample"),
            unit=pe_units.get_level_values("unit"),
            read=["1", "2"]
        ),
        expand(
            "results/cutAdapt/{sample}-{unit}.fastq.gz",
            sample=se_units.get_level_values("sample"),
            unit=se_units.get_level_values("unit")
        )

验证方法 #

运行snakemake -n进行试运行，确认不再出现KeyError，且生成的目标文件路径均为units表中存在的有效配对。

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