使用Snakemake构建按染色体拆分BAM流程遇Wildcards错误求助

Let's break down what's going wrong here and fix it step by step:

1. Core Issues Identified #

• {OutputDir} is misused as a wildcard: If OutputDir is just a fixed output path (not a variable sample/pipeline parameter), wrapping it in curly braces tells Snakemake to treat it as a wildcard. Since your output files don't have a variable OutputDir component, Snakemake can't resolve what value it should take.
• Incorrect expand() usage in split_bam: Snakemake rules are designed to handle single, atomic tasks (like splitting one chromosome at a time). Using expand() here creates a full list of all split BAMs as outputs for a single rule run, which breaks the wildcard matching logic—Snakemake can't map one input BAM to dozens of outputs in a single rule instance.
• Missing chr wildcard definition: Your shell command uses {chr} but the rule doesn't define this wildcard in its inputs/outputs, so Snakemake has no idea where to pull that value from.

2. Step-by-Step Fixes #

a. Define OutputDir as a global variable #

First, set OutputDir (and other static variables like name, nt) at the top of your Snakefile as regular variables, not wildcards. For example:

# At the top of your Snakefile, after loading config
OutputDir = config.get("output_dir", "results")  # Use config value or default to "results"
name = config.get("sample_name", "my_sample")
nt = config.get("threads", 8)
OntSoftware = config.get("singularity_image", "ont_software.sif")
ref_mmi = config.get("ref_mmi", "reference.mmi")
ref_fasta = config.get("ref_fasta", "reference.fasta")
TargetFastq = config.get("fastq", "input.fastq.gz")

b. Rewrite split_bam to handle single chromosomes #

Instead of expand(), use a {chr} wildcard in the output to represent one chromosome per rule run. This lets Snakemake spawn a separate job for each chromosome:

rule split_bam:
    input:
        bam = rules.minimap2.output.Sortbam
    output:
        # Double curly braces escape the wildcard for Python's f-string
        splitBam = f"{OutputDir}/split/{name}.{{chr}}.bam",
        splitBamBai = f"{OutputDir}/split/{name}.{{chr}}.bam.bai"
    # Ensure only valid chromosomes are processed
    wildcard_constraints:
        chr = "|".join(ref_chr)
    resources:
        mem_mb = 30000
    threads: nt
    singularity:
       OntSoftware
    shell:
        """
        # Use wildcards.chr to access the current chromosome for this job
        samtools view -@ {threads} -b {input.bam} {wildcards.chr} > {output.splitBam}
        samtools index -@ {threads} {output.splitBam}
        """

• The double {{chr}} in the f-string ensures Python outputs a single {chr} that Snakemake recognizes as a wildcard.
• wildcard_constraints prevents invalid chromosome values from triggering jobs.

c. Keep rule all as-is #

Your rule all is correct—it uses expand() to list all split BAMs as final targets, which tells Snakemake to spawn the necessary split_bam jobs for every chromosome in ref_chr.

3. Full Corrected Snakefile #

# Load config file (adjust path as needed)
configfile: "config.yaml"

# Global static variables
OutputDir = config.get("output_dir", "results")
name = config.get("sample_name", "my_sample")
nt = config.get("threads", 8)
OntSoftware = config.get("singularity_image", "ont_software.sif")
ref_mmi = config.get("ref_mmi", "reference.mmi")
ref_fasta = config.get("ref_fasta", "reference.fasta")
TargetFastq = config.get("fastq", "input.fastq.gz")

# Define reference chromosomes based on config
if config['ref'] == 'hg38':
    ref_chr = []
    for i in range(1,23):
        ref_chr.append('chr'+str(i))
    ref_chr.extend(['chrX','chrY'])
    
elif config['ref'] == 'b37':
    ref_chr = []
    for i in range(1,23):
        ref_chr.append(str(i))
    ref_chr.extend(['X','Y'])

rule all:
    input:
        expand(f"{OutputDir}/split/{name}.{{chr}}.bam", chr=ref_chr)

rule minimap2:
    input:
        TargetFastq
    output:
        Sortbam = f"{OutputDir}/{name}.sorted.bam",
        Sortbai = f"{OutputDir}/{name}.sorted.bam.bai"
    resources:
        mem_mb = 40000
    threads: nt
    singularity:
       OntSoftware
    shell:
        """
        minimap2 -ax map-ont -d {ref_mmi} --MD -t {threads} {ref_fasta} {input} | samtools sort -O BAM -o {output.Sortbam}
        samtools index {output.Sortbam}
        """

rule split_bam:
    input:
        bam = rules.minimap2.output.Sortbam
    output:
        splitBam = f"{OutputDir}/split/{name}.{{chr}}.bam",
        splitBamBai = f"{OutputDir}/split/{name}.{{chr}}.bam.bai"
    wildcard_constraints:
        chr = "|".join(ref_chr)
    resources:
        mem_mb = 30000
    threads: nt
    singularity:
       OntSoftware
    shell:
        """
        samtools view -@ {threads} -b {input.bam} {wildcards.chr} > {output.splitBam}
        samtools index -@ {threads} {output.splitBam}
        """

4. Why This Works #

• OutputDir is now a fixed path, so Snakemake no longer treats it as an unresolved wildcard.
• split_bam now runs one job per chromosome, which aligns with Snakemake's task-oriented design. The expand() in rule all triggers all necessary jobs automatically.
• wildcards.chr correctly passes the chromosome value for each individual job to the samtools commands.

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