使用Snakemake规则重命名原始reads时遇索引越界错误求助

问题核心是Snakemake默认会生成{population}和{individual}的所有笛卡尔积组合，但样本表中很多组合并不存在，导致get_sequences查询不到对应数据触发索引错误。以下是保留原函数逻辑的修复方案：

1. 提取有效种群-个体配对 #

从样本表中提取所有实际存在的种群-个体组合，避免Snakemake生成无效配对：

import pandas as pd

samples = pd.read_csv(parameters["reference"]["samples_csv"], dtype="str")
# 获取所有有效的(population, individual)配对，转为二维列表
valid_pop_ind_pairs = samples[["pop", "ind"]].drop_duplicates().values.tolist()

2. 优化get_sequences函数 #

添加空值判断，抛出更明确的错误信息，同时保留原查询逻辑：

def get_sequences(wildcards):
    # 筛选当前wildcard对应的样本记录
    match_row = samples[(samples["pop"] == wildcards.population) & (samples["ind"] == wildcards.individual)]
    # 若未找到匹配记录，抛出明确错误
    if match_row.empty:
        raise ValueError(f"无匹配样本：种群={wildcards.population}，个体={wildcards.individual}")
    # 提取R1和R2路径
    R1, R2 = match_row[["R1", "R2"]].values[0]
    return R1, R2

3. 指定规则的有效处理目标 #

添加rule all，通过expand配合zip生成所有需要处理的输出文件，确保Snakemake只处理有效组合：

rule all:
    input:
        # 生成所有R1和R2的目标路径
        expand("data/raw/{pop}.{ind}_1.fq.gz", zip, pop=[pair[0] for pair in valid_pop_ind_pairs], ind=[pair[1] for pair in valid_pop_ind_pairs]),
        expand("data/raw/{pop}.{ind}_2.fq.gz", zip, pop=[pair[0] for pair in valid_pop_ind_pairs], ind=[pair[1] for pair in valid_pop_ind_pairs])

rule symlinks:
    input:
        get_sequences
    output:
        R1 = "data/raw/{population}.{individual}_1.fq.gz",
        R2 = "data/raw/{population}.{individual}_2.fq.gz",
    log:
        R1 = "logs/make_symlinks/{population}.{individual}_1.log",
        R2 = "logs/make_symlinks/{population}.{individual}_2.log",
    shell:
        """
        ln --symbolic $(readlink --canonicalize {input[0]}) {output.R1} > {log.R1} 2>&1
        ln --symbolic $(readlink --canonicalize {input[1]}) {output.R2} > {log.R2} 2>&1
        """

修复原理 #

• 不再单独提取去重的种群和个体列表，而是直接获取两者的有效配对，避免生成无效组合。
• 通过rule all明确告诉Snakemake需要处理的所有目标文件，让工作流只针对存在的样本执行symlink操作。
• 优化后的函数增加了错误判断，便于快速定位异常情况。

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