用外门邪道的方式使用snakemake实现流程控制

Snakemake 确实是非常好用的流程开发及管理工具. 但是在特定的场景下…它也会带来一些问题, 然后在机缘巧合下, 我发现了一种非常蛋疼的使用方式: 只使用snakemake的依赖处理和人物管理, 执行脚本另外生成.

其实起因是因为让我写流程的时候要在特定位置放一份实际的执行脚本, 这样使用流程的人可以读这些脚本之后可以对特定的步骤进行修改后再次投递执行.

初衷是非常好的: 让使用流程, 但不了解snakemake的人也能轻松的对特定样本的分析参数作修改. 可是问题是…snakemake并不支持将实际执行的脚本单独导出形成文件, 只能使用-p参数将执行的内容直接输出在日志里, 并且之前在找解决方案的时候我还不知道可以在cluster模式的投递命令中使用wildcard…所以几经周折, 我找到了如下这个非常歪门邪道的解决方案.

1. 将整个流程划分为若干部分, 然后自写程序生成所有需要执行的脚本
2. 在投递到SGE时, 使用qsub的-sync y参数, 使qsub投递后并不直接结束程序, 而是的等待到投递的完成/失败, 如果人物失败, 会返回非0数值
3. 使用snakemake构建流程依赖, input / output照常指定, 在input中加上执行脚本, 然后shell:后只写qsub的投递命令

这样一来, 可以曲线解决不能直接生成脚本的问题 –> 我自己生成, 给你投递就好
同时也可以利用snakemake的依赖解决和任务管理特性 –> 可以进度监测, 可以断点继续
但同时也挖了一个巨坑: 脚本执行内容和各文件名未在snakemake内指定, 如果不一致流程会出错, 且这种问题不会有什么比较好的方法能定位, 全靠写的人自己小心…

比如: 你写了这样一个脚本test.sh:

1
2

#$ -sync y
cat test.txt | cut -f 2 > tset.col2.txt

然后写的snakemakefile是这样的:

1
2
3
4
5
6
7
8
9
10
11
12
13

rule all:
    input:
        test.col2.txt

rule test:
    input:
        {sample}.txt
    output:
        {sample}.col2.txt
    shell:
        '''
        qsub test.sh
        '''

那么在使用snakemake运行的时候, 是永远不会运行成功的…因为脚本生成的结果名字和snakemake的目的文件名不一致…
然后在默认情况下, 因为判定rule test执行失败, 其他被用于判定的结果文件(如果存在的话)会被删除…

也就是说, 用这种方式的话, 得自己保证目标文件的一致性. 排除这个缺点之外, 这种方式能充分利用snakemake在任务管理上的优势, 又可以避免初期上手snakemake时在理解wildcard机制时所必须要花的时间.