如何在Snakemake中读取BED文件内容转换为chopchop所需参数格式

你可以通过「提前解析BED文件生成区间列表」+ 适配Snakemake规则的方式完成需求，下面提供两种常用实现方式：

方案1：所有区间统一输出到单个结果文件 #

适合区间数量少、不需要拆分结果的场景：

1. 首先在Snakemake文件开头添加BED解析逻辑：

# 读取BED文件并转换为chr:start-end格式的区间列表
def read_bed_regions(bed_path):
    regions = []
    with open(bed_path, "r") as f:
        for line in f:
            line = line.strip()
            # 跳过空行和注释行
            if not line or line.startswith("#"):
                continue
            cols = line.split()
            regions.append(f"{cols[0]}:{cols[1]}-{cols[2]}")
    return regions

# 替换为你实际的BED文件路径
TARGET_REGIONS = read_bed_regions("/gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/enhancer_dataset/jurkat/test.bed")

2. 修改你的run_chopchop规则，循环调用chopchop：

rule run_chopchop :
    input:
        bed = "/gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/enhancer_dataset/jurkat/test.bed"
    output:
        "/gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/chopchop_output/guide_chopchop.txt"
    params:
        regions = TARGET_REGIONS,
        # 注意你原命令里多了一个多余的`.`前缀，这里修正为正确的绝对路径
        chopchop_path = "/gpfs/home/user/git/chopchop/chopchop.py"
    shell:'''
    set +u; source /gpfs/home/user/Apps/anaconda3/bin/activate chopchop; set -u
    # 提前创建temp目录避免报错
    mkdir -p temp
    # 循环所有区间调用chopchop，结果追加到输出文件
    for region in {params.regions}; do
        python {params.chopchop_path} -G hg38 -o temp -Target $region >> {output}
    done
    '''

方案2：每个区间输出独立结果文件（更推荐） #

符合Snakemake增量运行的特性，中间出错不需要重跑所有区间，方便后续排查：
在前面解析BED的基础上，修改规则为通配符模式：

# 生成每个区间对应的输出文件路径
region_outputs = []
for region in TARGET_REGIONS:
    # 替换区间中的特殊字符为下划线，避免文件名非法
    safe_region_name = region.replace(":", "_").replace("-", "_")
    region_outputs.append(f"/gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/chopchop_output/{safe_region_name}_guide.txt")

# 总入口规则，依赖所有区间的输出结果
rule all:
    input: region_outputs

# 单个区间的处理规则
rule run_chopchop_per_region:
    output:
        "/gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/chopchop_output/{region_name}_guide.txt"
    params:
        # 把通配符中的文件名还原为原始区间格式
        target = lambda wildcards: wildcards.region_name.replace("_", ":", 1).replace("_", "-")
    shell:'''
    set +u; source /gpfs/home/user/Apps/anaconda3/bin/activate chopchop; set -u
    mkdir -p temp
    python /gpfs/home/user/git/chopchop/chopchop.py -G hg38 -o temp -Target {params.target} > {output}
    '''

注意事项 #

• 如果你需要每个区间的结果独立表头，可以在方案1的循环里加过滤逻辑，跳过重复的表头行
• 若chopchop支持单次传入多个-Target参数，你可以直接把所有区间拼在-Target后面，不需要循环，可根据工具实际使用方式调整
• 注意提前创建输出目录chopchop_output，或者在shell开头加mkdir -p /gpfs/home/user/crispr_project/CRISPRi_Enh_TALL/chopchop_output

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