Snakemake按染色体拆分文件后无法识别输入/输出文件的问题

嘿，这个问题我太熟了——在Snakemake里处理动态生成的文件（比如按染色体拆分BAM这种没法提前知道输出数量的场景）确实容易踩坑，不过咱们用checkpoint功能就能完美解决！下面给你一步步的具体方案：

核心思路 #

Snakemake的checkpoint专门用来处理规则生成未知数量输出的情况。我们先通过checkpoint执行BAM拆分，再用自定义函数收集拆分后的文件，让后续的变异检测和VCF合并规则能自动识别这些动态生成的输入输出。

具体实现步骤 #

1. 用Checkpoint拆分BAM文件 #

先写一个checkpoint规则，把拆分后的BAM统一放到单独目录里，方便后续收集：

checkpoint split_bam:
    input:
        "sample.bam"
    output:
        directory("split_bams")  # 创建目录存放拆分后的BAM
    shell:
        """
        mkdir -p split_bams
        # bamtools split会自动在-out指定的前缀后加上染色体名和.bam后缀
        bamtools split -in {input} --reference -out split_bams/sample.REF_
        """

2. 编写函数收集拆分后的BAM #

接下来写一个辅助函数，从checkpoint的输出目录里收集所有拆分后的BAM文件，并提取对应的染色体名称：

import os
import glob

def get_split_bams(wildcards):
    # 获取checkpoint的输出目录
    checkpoint_dir = checkpoints.split_bam.get_output(wildcards)[0]
    # 遍历目录下所有BAM文件，提取染色体名
    bam_files = glob.glob(os.path.join(checkpoint_dir, "*.bam"))
    chroms = [os.path.basename(f).replace("sample.REF_", "").replace(".bam", "") for f in bam_files]
    # 返回所有拆分BAM的路径列表
    return expand("split_bams/sample.REF_{chrom}.bam", chrom=chroms)

3. 变异检测规则 #

现在可以针对每个拆分后的BAM文件写变异检测规则，用wildcard匹配染色体：

rule call_variants:
    input:
        lambda wildcards: f"split_bams/sample.REF_{wildcards.chrom}.bam"
    output:
        "split_vcfs/sample.REF_{chrom}.vcf"
    shell:
        """
        mkdir -p split_vcfs
        # 替换成你实际使用的变异检测命令，比如FreeBayes、GATK等
        freebayes -f reference.fasta {input} > {output}
        """

记得把freebayes -f reference.fasta {input} > {output}换成你真实的变异检测命令哦。

4. 合并VCF文件 #

最后写合并规则，收集所有拆分后的VCF文件，用vcf-concat合并：

def get_split_vcfs(wildcards):
    # 收集所有拆分后的VCF文件
    vcf_files = glob.glob("split_vcfs/*.vcf")
    chroms = [os.path.basename(f).replace("sample.REF_", "").replace(".vcf", "") for f in vcf_files]
    return expand("split_vcfs/sample.REF_{chrom}.vcf", chrom=chroms)

rule merge_vcfs:
    input:
        get_split_vcfs
    output:
        "sample.vcf"
    shell:
        """
        vcf-concat {input} > {output}
        """

可选方案：提前提取染色体列表 #

如果你不想用checkpoint，也可以先写一个规则提取BAM里的染色体列表，再基于这个列表生成所有输入输出：

rule get_chromosomes:
    input:
        "sample.bam"
    output:
        "chroms.txt"
    shell:
        """
        # 用samtools提取染色体名，过滤掉未比对的序列（*）
        samtools idxstats {input} | cut -f1 | grep -v '*' > {output}
        """

# 读取染色体列表
CHROMS = open("chroms.txt").read().splitlines()

rule split_bam:
    input:
        "sample.bam"
    output:
        expand("sample.REF_{chrom}.bam", chrom=CHROMS)
    shell:
        """
        bamtools split -in {input} --reference
        """

# 后续的call_variants和merge_vcfs规则可以直接用CHROMS变量来生成输入输出

不过这个方案的缺点是，如果BAM文件更新了，需要手动重新运行get_chromosomes规则，而checkpoint的方式更灵活，能自动适配BAM的变化。

这样处理后，Snakemake就能自动识别所有动态生成的文件，再也不会报unknown output/input files的错误啦！

内容的提问来源于stack exchange，提问作者Wouter De Coster