Snakemake处理目录中配对测序reads的方案咨询

1. 自动识别配对文件（核心步骤） #

从目标目录批量抓取所有测序文件，通过正则匹配分组配对reads。Snakemake支持在规则外执行Python代码，直接完成样本与配对文件的映射：

import glob
import re

# 抓取所有兼容格式的测序文件
read_files = glob.glob("path/to/your/reads/*.fastq*")
# 定义配对文件匹配正则（可根据实际命名规则调整，比如替换为_1/_2模式）
pair_pattern = re.compile(r"(.+)_R([12])\.fastq(\.gz)?$")

# 构建样本-配对文件字典
paired_reads = {}
for f in read_files:
    match = pair_pattern.match(f)
    if match:
        sample_id = match.group(1)
        read_type = f"R{match.group(2)}"
        if sample_id not in paired_reads:
            paired_reads[sample_id] = {"R1": None, "R2": None}
        paired_reads[sample_id][read_type] = f

# 过滤缺失配对的无效样本（可选）
valid_samples = [s for s in paired_reads if paired_reads[s]["R1"] and paired_reads[s]["R2"]]

2. 兼容压缩与非压缩文件 #

无需单独区分.fastq和.fastq.gz，在shell命令中用组合命令自动适配：

# 自动识别文件压缩状态并读取
zcat {input.r1} 2>/dev/null || cat {input.r1}

原理：若文件为压缩格式，zcat正常执行；若为非压缩格式，zcat报错后自动切换为cat读取。

3. 动态生成流程规则 #

以fastp质控为例，用样本ID作为wildcard，动态绑定对应配对文件：

流程入口规则 #

rule all:
    input:
        expand("results/clean/{sample}_R1.fastq.gz", sample=valid_samples),
        expand("results/clean/{sample}_R2.fastq.gz", sample=valid_samples)

质控处理规则 #

rule fastp_trim:
    input:
        r1=lambda wildcards: paired_reads[wildcards.sample]["R1"],
        r2=lambda wildcards: paired_reads[wildcards.sample]["R2"]
    output:
        r1_clean="results/clean/{sample}_R1.fastq.gz",
        r2_clean="results/clean/{sample}_R2.fastq.gz",
        report="results/reports/{sample}_fastp.html"
    log: "logs/fastp/{sample}.log"
    shell:
        """
        fastp -i {input.r1} -I {input.r2} \
              -o {output.r1_clean} -O {output.r2_clean} \
              -h {output.report} \
              --thread 4 2>&1 > {log}
        """

4. 额外优化建议 #

• 自定义配对规则：若文件命名不是_R1/_R2，修改pair_pattern即可，比如匹配_1.fastq的正则：r"(.+)_([12])\.fastq(\.gz)?"
• 参数校验：在规则前添加断言，避免无效样本进入流程：

  assert all(paired_reads[s]["R1"] and paired_reads[s]["R2"] for s in valid_samples), "存在单端样本，请检查文件命名"
  

• 临时文件管理：将中间产物标记为temp()，Snakemake会自动清理：

  output:
      r1_clean=temp("results/clean/{sample}_R1.fastq.gz"),
      ...
  

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