流程管理工具Snakemake用法总结

介绍

snakemake 原理

image

建议的框架

├── config.yaml
├── environment.yaml
├── scripts
│   ├── script1.py
│   └── script2.R
└── Snakefile

Install

conda create -n py35 python=3.5
source activate py35
conda install -c bioconda -c conda-forge snakemake

用conda在python3环境里安装。

之后,如果python3可以正常用。(which python3)可以被访问,那么snakemake可以直接运行。不需要source到py35环境。

如果没有python3,那么需要source到py35环境才可以用snakemake。

使用

常用参数

  • --snakefile, -s 指定Snakefile,否则是当前目录下的Snakefile
  • --dryrun, -n 不真正执行,一般用来查看Snakefile是否有错
  • --printshellcmds, -p 输出要执行的shell命令
  • --reason, -r 输出每条rule执行的原因,默认FALSE
  • --cores, --jobs, -j 指定运行的核数,若不指定,则使用最大的核数
  • --force, -f 重新运行第一条rule或指定的rule。或指定输出结果(绝对路径)
  • --forceall, -F 重新运行所有的rule,不管是否已经有输出结果。
  • --forcerun, -R 执行指定的rule及下游所有任务.无需target。(但rule all中必须有rule的输出文件)
  • --keep-going, -k Go on with independent jobs if a job fails.(某任务报错,其他不相关的任务可以继续运行)
  • snakemake --dag | dot -Tpdf > dag.pdf 可视化

关于f强制运行

  • output没有wildcards的: -f rulename
  • output有wildcards:-f outputfile

多个output指定一个就可以。要是绝对路径。比如rule的output使用了wildcards,就像下面的

join(outdir0, "{sample}_1.fastq"), join(outdir0, "{sample}_2.fastq")

可以指定(target)即某个样本的输出文件

snakemake -np -f projects/02snakeReq/00.prepare/B_1.fastq

集群投递

  • snakemake --cluster "qsub -V -cwd -q v01" -j 10
  • --cluster 集群运行指令
  • qusb -V -cwd -q 表示输出当前环境变量(-V),在当前目录下运行(-cwd), 投递到指定的队列(-q), 如果不指定则使用任何可用队列
  • --local-cores N: 在每个集群中最多并行N核
  • --cluster-config/-u FILE: 集群配置文件
snakemake -j 99 --debug --immediate-submit --cluster-config cluster.json --cluster 'bsub_cluster.py {dependencies}'

snakemake --cluster "qsub -q res -cwd  -o logs -e logs" --jobs 32

the maximum number of jobs to be queued or executed at the same time

可最大提交32的jobs同时跑。在cpu允许范围内,jobs越大跑的越快

使用例子

常用

# 查看全部流程
snakemake -np
# 指定snakefile
snakemake --snakemakefile myfile.py
# 重跑rule dedup_realn(dry run),注意如果snakemake认为需要跑的step(rull all输出文件推断)还是会跑的
snakemake --snakefile res-snake.py -np -R dedup_realn 

修改代码后,重跑改变的rule

snakemake -n -R `snakemake --list-input-changes`
snakemake -n -R `snakemake --list-params-changes`
snakemake -n -R `snakemake --list-code-changes`

删除snamake 输出的文件

仅删除rule在output中写的文件。最好先和--dry-run跑确认一下。同样的还有--delete-temp-output

 snakemake some_target --delete-all-output
snakemake index --delete-all-output -np --snakefile res-snake.py
Building DAG of jobs...
Would delete /02snakeReq/00.prepare/ref/salmonella.dict
Would delete /02snakeReq/00.prepare/ref/salmonella.fa

workflow太大了,不看output,只看最终结果

snakemake -n --quiet

打开zsh关于snakemake的自动完成:将以下代码放在~/.zshrc中
compdef _gnu_generic snakemake

部分运行

## 所有和这个rule相关的file都会rerun
./pyflow-ChIPseq -R call_peaks_macs2

## 重跑一个样本,只要指定target file
./pyflow-ATACseq -R 04aln/m280.sorted.bam

##只重跑align rule。
./pyflow-ATACseq -f align

-f
强制跑:依赖的前面没有重跑的文件亦可以跳过,只要旧文件存在就可以跑。

有wildcards的rule不能-f来跑。因为没有跑rule all,无法推断sample。

但可以指定输出文件来跑-f

查看结果

snakemake --summary | sort -k1,1 | less -S

# or detailed summary will give you the commands used to generated the output and what input is used
snakemake --detailed-summary | sort -k1,1 > snakemake_run_summary.txt

output_file date rule version log-file(s) status plan

其他功能

访问其他rule中的参数

rule salmon_quant:
    input:
        r1 = lambda wildcards: FILES[wildcards.sample]['R1'],
        r2 = lambda wildcards: FILES[wildcards.sample]['R2'],
        index = rules.salmon_index.output.index

protected

用来指定某些中间文件是需要保留的,eg.output: protected(“f1.bam”)。

expand :

相当于列表推导式


image

邮件通知
link

config

可在snakemake命令时设置。会覆盖config.yaml中的相同参数。优先认为是number。如果必须识别为string,可以加引号。

snakemake mytarget --config foo=bar
snakemake --config 'version="2018_1"'

priorites

temp

temp: 通过temp方法可以在所有rule运行完后删除指定的中间文件,eg.output: temp(“f1.bam”)。

Further, an output file marked as temp is deleted after all rules that use it as an input are completed:

rule NAME:
    input:
        "path/to/inputfile"
    output:
        temp("path/to/outputfile")
    shell:
        "somecommand {input} {output}"

报错和解决

反义

在rule的shell使用awk注意:大括号要变成双括号;tab要再次反义

grep -v "#"  ${{smp}}_merged.sv.vcf|awk '
            BEGIN{{FS=OFS="\\t"}}
            $10!~/NaN/{{$10=NR}}
            ...

output为文件夹

Outputs of incorrect type (directories when expecting files or vice versa). Output directories must be flagged with directory().

ref_split = directory(join(outdir0, "ref/ref_split"))

当output是一个文件夹?关于snakemake自动建文件夹
好像不行。output只能是文件 才会建文件夹

Unable to set utime on symlink /projects/02snakeReq/00.prepare/ref/salmonella.fa. Your Python build does not support it.

设置output文件来让snakemake了解分析顺序。
output文件必须是文件而不是文件夹。
后rule的input,必须承接前rule的output。目前不知道是否也要写在rule all中
snp_filter的output,和 snp_indel_anno来衔接
annovar_db的output,he snp_indel_anno来衔接

报错

检查rule的output和rule all的input是否一一对应


image

wildcards

shell中{sample}要写成{wildcards.sample}


image

MissingOutputException

Waiting at most 50 seconds for missing files.
MissingOutputException in line 525 of projects/02snakeReq/res-snake.py:
Missing files after 50 seconds:
projects/03svTest/04.sv/D-korps6_svaba/D-korps6.sv.vcf

rule的output文件必须在输出结果文件中。

KeyError

不同的结果文件名不能一模一样。(写在input,output中的文件)

举例:05.sv下有lumpy和manta文件夹,都有A.sv.vcf文件。后面在05.sv文件夹下也有A.sv.vcf文件。而且这三个文件分属不同的rule。sample会识别错误

KeyError: 'A-korps6_pindel/A-korps6'
Wildcards:
sample=A-korps6_pindel/A-korps6

重跑流程
重跑流程,当原来产生的文件要被覆盖时,可能会出现ProtectedOutputException。要用protected()

ProtectedOutputException in line 114 of Resequencing/03.SNP_INDEL_detection_v2/snakefile_v2.py:
Write-protected output files for rule index:
/projects/04liver_cancer/00.prepare/ref/hg19.fa

参考

  1. 生信分析流程构建的几大流派
  2. snakemake框架
  3. snakemake pyflow-ATACseq流程
  4. https://github.com/AAFC-BICoE/snakemake-mothur
  5. 使用GATK和snakemake框架的总结
  6. yaml语言教程
  7. snakemake 官方文档
  8. https://github.com/inodb/snakemake-parallel-bwa/blob/master/Snakefile
  9. mutect流程
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,080评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,422评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,630评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,554评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,662评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,856评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,014评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,752评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,212评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,541评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,687评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,347评论 4 331
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,973评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,777评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,006评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,406评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,576评论 2 349

推荐阅读更多精彩内容