如何在Snakemake中用通配符合并不同通配符的测序文件

问题背景 #

使用Snakemake处理数百个基因组测序样本，部分样本因测序质量问题进行了二次测序，每个样本对应4或8个fastq文件（分属不同run/lane）。已完成trimming和mapping步骤，但合并不同run/lane的BAM文件时，原通配符无法匹配多组文件，此前用shell条件判断的实现方式不够优雅且适配性差，需优化通配符使用逻辑并设计通用合并规则。

文件结构示例 #

multiple_subdirectories/BM1P149_A2/111353_BY/
├── 111353_BY_run757_AATGGTAG_S238_L001_R1_001.fastq.gz
├── 111353_BY_run757_AATGGTAG_S238_L001_R2_001.fastq.gz
├── 111353_BY_run757_AATGGTAG_S238_L002_R1_001.fastq.gz
├── 111353_BY_run757_AATGGTAG_S238_L002_R2_001.fastq.gz
├── 111353_BY_run758_AATGGTAG_S85_L001_R1_001.fastq.gz
├── 111353_BY_run758_AATGGTAG_S85_L001_R2_001.fastq.gz
├── 111353_BY_run758_AATGGTAG_S85_L002_R1_001.fastq.gz
└── 111353_BY_run758_AATGGTAG_S85_L002_R2_001.fastq.gz

现有Trimming阶段代码 #

(PLAQUES,SAMPLES,IDS,PAIRED,) = glob_wildcards(config["dataDirInput"]+"{plaque}/{sample}/{id}_R{paire}_001.fastq.gz")

rule all:
  input: 
    expand(config["dataDirOutput"]+"{plaque}/{sample}/{id}_R{paire}_001_unpaired.fastq.gz", zip, plaque = PLAQUES, sample = SAMPLES, id = IDS, paire = PAIRED)

rule trimmomatic:
  input: 
    read1=config["dataDirInput"]+"{plaque}/{sample}/{id}_R1_001.fastq.gz",
    read2=config["dataDirInput"]+"{plaque}/{sample}/{id}_R2_001.fastq.gz"
  output: 
    read1_paired=config["dataDirOutput"]+"{plaque}/{sample}/{id}_R1_001_paired.fastq.gz",
    read1_unpaired=config["dataDirOutput"]+"{plaque}/{sample}/{id}_R1_001_unpaired.fastq.gz",
    read2_paired=config["dataDirOutput"]+"{plaque}/{sample}/{id}_R2_001_paired.fastq.gz",
    read2_unpaired=config["dataDirOutput"]+"{plaque}/{sample}/{id}_R2_001_unpaired.fastq.gz"
  threads: 20  
  resources:
    mem_mb = 20,
    time_min = 300 
  params: 
    env=config["trimmomatic_env"] # 使用conda环境调用trimmomatic
  shell: "{params.env} PE {input.read1} {input.read2} {output.read1_paired} {output.read1_unpaired} {output.read2_paired} {output.read2_unpaired} ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:True LEADING:3 TRAILING:3 MINLEN:36"

配置文件片段 #

trimmomatic_env: "/directory_of_the_conda_environment"
dataDirInput: "/multiple_subdirectories/Sample_GData"  
dataDirOutput: "/multiple_subdirectories/Trimmed"    

优化方案 #

1. 重新解析通配符，拆分样本核心标识与Run/Lane信息 #

原{id}包含run、lane等冗余信息，导致无法按样本分组。需拆分文件名字段，提取样本核心ID（如111353_BY）、run号、lane号等独立通配符：

# 拆分文件名中的可变字段，匹配plaque、sample、样本核心ID、run、lane、配对端信息
(PLAQUES, SAMPLES, SAMPLE_CORES, RUNS, LANES, PAIRED) = glob_wildcards(
    config["dataDirInput"]+"{plaque}/{sample}/{sample_core}_run{run}_*_L{lane}_R{paired}_001.fastq.gz"
)

# 去重得到唯一的样本组合（plaque+sample+样本核心ID），用于合并阶段
UNIQUE_SAMPLES = sorted(list(set(zip(PLAQUES, SAMPLES, SAMPLE_CORES))))

2. 调整Mapping规则（按Run/Lane生成独立BAM） #

确保每个run/lane的测序数据生成独立BAM文件，方便后续合并：

rule mapping:
    input:
        read1=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_*_L{lane}_R1_001_paired.fastq.gz",
        read2=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_*_L{lane}_R2_001_paired.fastq.gz"
    output:
        bam=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_L{lane}.bam",
        bai=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_L{lane}.bam.bai"
    threads: 20
    params:
        env=config["mapping_env"] # 替换为你的mapping软件conda环境
    shell:
        """
        {params.env} bwa mem -t {threads} reference.fasta {input.read1} {input.read2} | samtools sort -@ {threads} -o {output.bam}
        samtools index {output.bam}
        """

3. 编写通用BAM合并规则 #

利用lambda函数动态收集同一样本下的所有run/lane BAM文件，无需手动指定数量：

rule merge_bams:
    input:
        # 根据通配符匹配当前样本的所有run/lane BAM
        lambda wildcards: expand(
            config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_L{lane}.bam",
            run=RUNS, lane=LANES,
            plaque=wildcards.plaque, sample=wildcards.sample, sample_core=wildcards.sample_core
        )
    output:
        merged_bam=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_merged.bam",
        merged_bai=config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_merged.bam.bai"
    threads: 20
    shell:
        """
        samtools merge -@ {threads} {output.merged_bam} {input}
        samtools index {output.merged_bam}
        """

4. 更新rule all，纳入合并后的BAM #

确保Snakemake识别最终需要生成的合并文件：

rule all:
    input:
        # 保留原trimmed文件的输出要求
        expand(config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_run{run}_*_L{lane}_R{paired}_001_unpaired.fastq.gz",
               zip, plaque=PLAQUES, sample=SAMPLES, sample_core=SAMPLE_CORES, run=RUNS, lane=LANES, paired=PAIRED),
        # 添加合并后的BAM文件
        expand(config["dataDirOutput"]+"{plaque}/{sample}/{sample_core}_merged.bam",
               zip, plaque=[p for p,s,sc in UNIQUE_SAMPLES],
               sample=[s for p,s,sc in UNIQUE_SAMPLES],
               sample_core=[sc for p,s,sc in UNIQUE_SAMPLES])

关键注意事项 #

• 通配符拆分要精准：确保文件名中每个可变字段（run、lane）都被拆分为独立通配符，这是按样本分组的核心前提。
• 动态输入生成：merge规则中用lambda+expand自动匹配所有符合条件的BAM，适配不同样本的run/lane数量差异。
• 去重唯一样本：通过set(zip(...))获取不重复的样本组合，避免重复生成合并任务。

内容的提问来源于stack exchange，提问作者Benoît Madec