第4节:输入输出文件与通配符归纳总结
说实话,学完前三节,大家应该对 Snakemake 的基本流程有感觉了。但很多人卡在通配符这里——我见过不少同事,明明规则写对了,就是跑不起来,最后发现是输入输出文件里的通配符没配对。这一节咱们集中梳理一下输入输出文件与通配符的关系,把那些容易踩的坑填平。
4.1 通配符到底是什么?
通配符,说白了就是规则里的占位符。它让你的规则能批量处理不同文件。比如你要处理样本 A、B、C,不用写三条规则,写一条带通配符的规则就够了。
它的工作原理是这样的:Snakemake 根据你请求的目标文件(比如 results/A.txt),反向推导通配符的取值(A),然后代入规则里执行。代码里用 {sample} 这种花括号表示。
{sample}、{chrom},别用 {x}、{y}。后续维护时自己都看得懂。
4.2 输入输出文件中的通配符规则
其实就三条核心规则,记住了就行:
- 规则 output 里的通配符定义了这个规则能生成的“变量”。每个 output 文件中的通配符组合就是这个规则的“签名”。
- 规则 input 里可以使用相同的通配符。Snakemake 会自动把通配符的值替换进去。
- output 中没有出现的通配符,也不能在 input 里凭空出现。除非你通过
wildcards对象显式传递(后面会讲)。
举个例子:
rule bwa_align:
input:
"data/{sample}.fastq"
output:
"results/{sample}.bam"
shell:
"bwa mem ref.fa {input} > {output}"
这里的 {sample} 同时出现在 input 和 output 中。Snakemake 根据 output 文件(比如 results/A.bam)提取出 sample=A,然后找 input 文件 data/A.fastq。简单直接。
4.3 通配符匹配的常见陷阱
我在项目中遇到过最坑的一个问题:output 里有两个通配符,但 input 里只用了其中一个。比如:
rule merge_bams:
input:
"data/{sample}.bam",
"data/{sample}_qc.bam"
output:
"results/{sample}_{batch}.merged.bam"
shell:
"samtools merge {output} {input}"
这个规则永远跑不通!因为 output 里有 {sample} 和 {batch} 两个通配符,但 input 里只出现了 {sample}。Snakemake 不知道 {batch} 从哪里来。要修复,要么从 input 里也引入 {batch},要么用 lambda wildcards 在 input 函数里生成路径。
{sample}.txt,input 里写的是 raw/{sample}_R1.fq。跑的时候一直报错说找不到文件。查了半天发现是 input 里的文件名写错了——_R1 应该是 _R1 但我的数据里是 _R1 还是 _R1?唉,通配符匹配不区分大小写?其实不是,是文件名大小写的问题。Snakemake 只是按字符串替换,文件名必须真实存在。
4.4 模式规则与通配符的配合
模式规则(rule with wildcards)加上 CHECKPOINT、expand 函数,就能实现动态批量。下面这张图总结了核心关系:
你想想看,这个流程其实跟“点菜”很像:你告诉服务员你要“宫保鸡丁”(目标文件),服务员根据菜单(规则)知道需要“鸡胸肉、花生、辣椒”(输入文件),然后后厨(shell)开始做。通配符就是菜单上的“空白处”,服务员根据你的要求填上实际食材名。
4.5 expand() 函数的通配符展开
expand() 是 Snakemake 里生成多个文件路径的利器。它会把通配符的所有组合展开成列表。比如:
samples = ["A", "B"]
chroms = ["chr1", "chr2"]
expand("results/{sample}_{chrom}.vcf", sample=samples, chrom=chroms)
这个会生成:
["results/A_chr1.vcf", "results/A_chr2.vcf", "results/B_chr1.vcf", "results/B_chr2.vcf"]
默认是笛卡尔积。但有时候你可能想要一一对应,而不是全组合。那就用 zip 参数:
expand("results/{sample}_{chrom}.vcf", zip, sample=samples, chrom=chroms)
这会生成 ["results/A_chr1.vcf", "results/B_chr2.vcf"]。我在做 paired-end 测序数据处理时常用这个——样本 A 和它的 read1/read2 一一对应。
expand("data/{sample}_R1.fq", sample=["A","B"]) 和 "data/{sample}_R1.fq" 不是一回事。前者生成列表,后者只是一个模板字符串。在 input 里要用 expand 的地方千万别偷懒直接用字符串。
4.6 通配符在输入函数中的高级用法
有时候输入文件路径不是简单的替换能搞定的。比如要根据样本名查一个配置文件,或者对通配符做字符串变换。这时候就要用 输入函数:
def get_input(wildcards):
sample = wildcards.sample
# 假设配置里记录了每个样本的测序平台
platform = config[sample]["platform"]
if platform == "ILLUMINA":
return f"data/{sample}_S1_L001_R1_001.fastq.gz"
else:
return f"data/{sample}.fastq.gz"
rule align:
input: get_input
output: "results/{sample}.bam"
shell:
"bwa mem ref.fa {input} > {output}"
这里 get_input 函数接收 wildcards 对象,可以读取 wildcards.sample,然后根据实际数据命名规则返回不同的输入路径。我们在实际项目中经常用这种方式来处理不同批次数据的命名不一致问题。
print 或日志,先看看实际生成了什么路径再批量跑。
4.7 通配符与 checkpoint:动态批量处理
有心的读者会发现,前面讲的通配符都是在运行前就确定好的。那如果运行过程中才能知道有哪些文件呢?比如先跑一个命令生成文件列表,然后根据这个列表再跑后续规则。这时候就要用 checkpoint。
checkpoint generate_list:
output: "list/{batch}.txt"
shell:
"ls data/{wildcards.batch}/*.fastq > {output}"
rule process_each:
input:
"data/{batch}/{sample}.fastq"
output:
"results/{batch}/{sample}.processed"
shell:
"touch {output}"
rule aggregate:
input:
expand("results/{batch}/{sample}.processed",
batch=checkpoints.generate_list.get(**batch_samples))
output:
"final/{batch}.done"
shell:
"touch {output}"
工作原理:Snakemake 先执行 generate_list,然后 aggregate 的 input 函数通过 checkpoints.generate_list.get() 读取列表文件,动态生成新的通配符组合。这样就能做到真正的“运行时才知道有多少个样本”。
这一块内容比较深,咱们后面会有专门章节讲 checkpoint。这里先记住:通配符不一定要在 Snakefile 里写死,也可以在执行过程中动态生成。
4.8 总结:输入输出文件与通配符的黄金法则
| 场景 | 法则 | 我的建议 |
|---|---|---|
| 静态文件名 | 直接写字符串,不用通配符 | 保持简单 |
| 批量处理样本 | 在 output 中定义通配符 | 通配符名称与变量名一致 |
| 多通配符组合 | 使用 expand() 生成列表 | 注意笛卡尔积 vs zip |
| 复杂输入路径 | 使用输入函数 | 加日志,调试路径 |
| 动态通配符 | 使用 checkpoint | 会单独章节讲解 |
嗯,这一节的内容就到这里。通配符是整个 Snakemake 的重中之重,大部分失败重试问题也跟通配符没配对有关。自己动手写几个规则,把上面的例子跑一遍,比看十遍文档都有用。