精华内容
下载资源
问答
  • snakemake使用笔记
    千次阅读
    2018-06-12 11:34:00

    snakemake是一个用来编写任务流程的工具,用python编写的,因此其执行的流程脚本也比较通俗易懂,易于理解。

    一、从一个简单的例子开始

    1、安装snakemake

    安装snakemake的方法有多种,snakemake官方推荐的是conda,安装方法如下:
    conda install -c bioconda snakemake

    2、一个简单的snakemake脚本

    虽然snakemake广泛的应用于生物信息方面的流程编写,但是snakemake的应用并不局限于编写生物信息学的流程,这里以一个简单的合并文件的例子开始介绍snakemake的简单使用。

    #首先我们建立两个文件
    $ echo "Here is hello." > hello.txt
    $ echo "Here is world." > world.txt
    
    #接下来开始编写我们的Snakefile rule concat: # 这里的rule可视为snakemake定义的关键字,concat使我们自定义的这一步任务的名称 input: # input同样是snakemake的关键字,定义了在这个任务中的输入文件 expand("{file}.txt", file=["hello", "world"]) #expand是一个snakemake定义的替换命令 output: # output也是snakemake的关键字,定义输出结果的保存文件 "merged.txt" shell: # 这里表示我们下面的命令将在命令行中执行 "cat {input} > {output}" #最后就可以在Snakefile的路径执行snakemake命令即可 $ snakemake $ cat merge.txt Here is hello. Here is world. 

    在上面的Snakefile脚本中,ruleinputoutputshellexpand均为snakemake中的关键字或者命令。同时Snakefile中的每一个rule其实都可以看作是一个简单的shell脚本,通过Snakefile将多个rule组织在一起并按照我们定义的顺序来执行。另外,在output中的结果文件可以是未存在目录中的文件,这时会自动创建不存在的目录。

    二、snakemake中的一些命令与规则

    1、rule

    rule是Snakefile中最主要的部分。如上面的例子所说,每一个rule定义了一系列pipe中的一步,每一个rule都可以当作一个shell脚本来处理,一般主要包括inputoutputshell3个部分。同时还有许多上面没有列出来的用法:

    1. rule all。不同于其他的rule,在rule all里面一般不会去定义要执行的命令,他一般用来定义最后的输出结果文件。除了rule all中定义的文件外最后输出结果不会保存任何中间文件。例如将上面的脚本改成如下文件则没有输出结果:
    rule all:
        input:
             #"merged.txt"  取消注释后,则能正常输出文件
    rule concat:
        input:
            expand("{file}.txt", file=["hello", "world"]) output: "merge.txt" shell: "cat {input} > {output}" 
    1. wildcards。用来获取通配符匹配到的部分,例如对于通配符"{dataset}/file.{group}.txt"匹配到文件101/file.A.txt,则{wildcards.dataset}就是101,{wildcards.group}就是A。
    2. threads。通过在rule里面指定threads参数来指定分配给程序的线程数,egthreads: 8
    3. resources。可用来指定程序运行的内存,eg. resources: mem_mb=800
    4. message。使用message参数可以指定每运行到一个rule时,在终端中给出提示信息,eg.message: "starting mapping ..."
    5. priority。可用来指定程序运行的优先级,默认为0,eg.priority: 20
    6. log。用来指定生成的日志文件,eg.log: "logs/concat.log"
    7. params。指定程序运行的参数,eg.params: cat="-n",调用方法为{params.cat}
    8. run。在run的缩进区域里面可以输入并执行python代码。
    9. scripts。用来执行指定脚本,eg.scripts: "rm_dup.py"
    10. temp。通过temp方法可以在所有rule运行完后删除指定的中间文件,eg.output: temp("f1.bam")
    11. protected。用来指定某些中间文件是需要保留的,eg.output: protected("f1.bam")
    12. ancient。重复运行执行某个Snakefile时,snakemake会通过比较输入文件的时间戳是否更改(比原来的新)来决定是否重新执行程序生成文件,使用ancient方法可以强制使得结果文件一旦生成就不会再次重新生成覆盖,即便输入文件时间戳已经更新,eg.input: ancient("f1.fastq")
    13. Rule Dependencies。可通过快捷方式指定前一个rule的输出文件为此rule的输入文件:
    rule a:
        input:  "path/to/input"
        output: "path/to/output"
        shell:  ...
    
    rule b:
        input:  rules.a.output   #直接通过rules.a.output 指定rule a的输出
        output: "path/to/output/of/b"
        shell:  ...
    
    1. report。使用snakemake定义的report函数可以方便的将结果嵌入到一个HTML文件中进行查看。

    2. configuration

    每计算一次数据都要重写一次Snakefile有时可能会显得有些繁琐,我们可以将那些改动写入配置文件,使用相同流程计算时,将输入文件的文件名写入配置文件然后通过Snakefile读入即可。
    配置文件有两种书写格式——json和yaml。在Snakefile中读入配置文件使用如下方式:

    configfile: "path/to/config.json" 
    configfile: "path/to/config.yaml" # 也可直接在执行snakemake命令时指定配置 $ snakemake --config yourparam=1.5 

    在shell命令中直接调用config文件中的内容的话,不需要引号,如config[a]而不是config["a"]
    集群计算配置

    三、snakemake的执行

    一般讲所有的参数配置写入Snakefile后直接在Snakefile所在路径执行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 重新执行Snakefile,当更新了rule时候使用此命令
    
    #一些可视化命令
    $ snakemake --dag | dot -Tpdf > dag.pdf
    
    #集群投递
    snakemake --cluster "qsub -V -cwd -q 节点队列" -j 10
    # --cluster /-c CMD: 集群运行指令 # qusb -V -cwd -q, 表示输出当前环境变量(-V),在当前目录下运行(-cwd), 投递到指定的队列(-q), 如果不指定则使用任何可用队列 # --local-cores N: 在每个集群中最多并行N核 # --cluster-config/-u FILE: 集群配置文件 

    参考:snakemake官方文档



    作者:To_2019_1_4
    链接:https://www.jianshu.com/p/14b9eccc0c0e
    來源:简书
    著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
    更多相关内容
  • 使用GATK4进行变体调用管道的Snakemake工作流程以及使用Delly和Tiddit进行结构变体调用的Snakemake工作流程 这是一个使用GATK4进行变体调用管道的Snakemake工作流程,还包括一个对结构变体的故意和轻率的调用(将...
  • snakemake使用小结

    2019-10-06 03:06:15
    自动匹配上了,(注意此时文件夹貌似只能有一个脚本文件),cp了一个好像报错了 接下来,要做排序了,代码最后一起贴 可以使用dag选项和dot命令对“规则的执行和依赖关系”进行可视化, snakemake --dag sorted_...

     首先在linux 里配置conda

    下载

    wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive/Anaconda3-5.3.1-Linux-x86_64.sh

    chmod +x Anaconda3-5.3.1-Linux-x86_64.sh

    bash Anaconda3-5.3.1-Linux-x86_64.sh

    安装完毕,如果忘记选择yes,敲conda命令报错“command not found" 加上source /root/anaconda3/etc/profile.d/conda.sh

    conda env list   得到 /root/anaconda3

    export PATH=~/anaconda3/bin:$PATH

    不然全局无法使用conda命令,(但是重启putty好像就不管用了,还不清楚原因)

    vs code可以不安装

     

    安装tree命令,yum install tree

    tree -af 可以查看树形文件结构 

     

    snakemake是纯python的任务流程工具(基于python3),以前商业环境用过control-M

    https://snakemake.readthedocs.io/en/stable/

     

    首先做个变异检测,就是和标准的序列做对比,有点类似于代码的compare,用过Beyond Compare或svn和git的筒子们应该很熟悉了,

    但是基因序列是个非常大的序列文件,在linux也没有windows那样简便的图形操作见面,而且,这种对比工作是大量重复的,需要脚本化。

    cd snakemake-snakemake-tutorial-623791d7ec6d
    conda env create --name snakemake-tutorial --file environment.yaml

    --------------------------------------------------------------

    export PATH=~/anaconda3/bin:$PATH
    source activate snakemake-tutorial

    --------------------------------------------------------------

     

    # 退出当前环境
    source deactivate

    这里使用到Samtools工具,具体使用方法可以参考https://blog.csdn.net/g863402758/article/details/53081342

    他是一个用于处理sam与bam格式的工具软件,能够实现二进制查看、格式转换、排序及合并等功能,

    结合sam格式中的flag、tag等信息,还可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作命令,

    还可以大大扩展samtools的使用范围与功能。

    conda install snakemake

     conda install samtools

     bowtie2和samtools都是对比工具,bowtie2暂时没安装,安装方法先记录下

    sudo wget https://jaist.dl.sourceforge.net/project/bowtie-bio/bowtie2/2.3.4.1/bowtie2-2.3.4.1-linux-x86_64.zip

    unzip bowtie2-2.3.4.1-linux-x86_64.zip

    vi /etc/environment

    添加 bin 目录的路径,并用 : 隔开

    source /etc/enviroment 使配置生效

    开始写job脚本

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/A.fastq"
        output:
            "mapped_reads/A.bam"
        shell:
            """ bwa mem {input} | samtools view -Sb - > {output} """
    期间一直出一个错误,说Command must be given as string after the shell keyword
    运行snakemake -np mapped_reads/A.bam检查一下是否会出错

    执行这个job,把-n去掉

    可以看到,生成了A.bam文件

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/{sample}.fastq"
        output:
            "mapped_reads/{sample}.bam"
        shell:
            """ bwa mem {input} | samtools view -Sb - > {output} """
    将A改成{sample},在输入命令的时候加上你的参数,自动匹配上了,(注意此时文件夹貌似只能有一个脚本文件),cp了一个好像报错了

    接下来,要做排序了,代码最后一起贴

    可以使用dag选项和dot命令对“规则的执行和依赖关系”进行可视化,

    snakemake --dag sorted_reads/{A,B,C}.bam.bai | dot -Tpdf > dag.pdf  这个命令好像会报错
    snakemake --dag sorted_reads/{A,B,C}.bam.bai | dot -Tsvg > dag.svg

    整合之前的BAM文件,做基因组变异识别
    SAMPLES=["A","B","C"] rule bcftools_call: input: fa="data/genome.fa", bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES), bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES) output: "calls/all.vcf" shell: "samtools mpileup -g -f {input.fa} {input.bam} | " "bcftools call -mv - > {output}"
    其中expand是自动匹配变量求文件路径的语法糖
    检查一下,snakemake -np calls/all.vcf

    最后出report,以上都是在规则里执行shell脚本,snakemake的一个优点就是可以在规则里面写Python脚本,只需要把shell改成run,此外还不需要用到引号。

    测试一下,snakemake -np report.html

    画出流程图

    snakemake --dag report.html | dot -Tsvg > final.svg

     

    执行一下:snakemake -p report.html
    可以看到生成了报告文件

    到此,还有

    rule all:

    log:

    多线程thread:

    -j 指定cpu核心

    params:

    加载configfile: "config.yaml"

    这几个功能没有操作,留个以后有空再处理

    最后,在新建一个snakemake项目时,都先用conda create -n 项目名 python=版本号创建一个全局环境,用于安装一些常用的软件,例如bwa、samtools、seqkit等。然后用如下命令将环境导出成yaml文件

    conda env export -n 项目名 -f environment.yaml

    以后再部署的时候,

    只需要conda env create -f environment.yaml

    这个过程类似于ghost系统,或者打包虚拟机类似

    参考了以下网址,感谢!

    https://www.jianshu.com/p/8e57fd2b81b2

    http://pedagogix-tagc.univ-mrs.fr/courses/ABD/practical/snakemake/snake_intro.html

    转载于:https://www.cnblogs.com/marszhw/p/10572745.html

    展开全文
  • 使用Snakemake搭建分析流程

    万次阅读 多人点赞 2018-12-28 19:28:43
    snakemake使用 threads 来定义当前规则所用的进程数,我们可以对之前的 bwa_map 增加该指令。 rule bwa_map: input: "data/genome.fa", "data/samples/{sample}.fastq" output: "mapped_reads/{sample}.bam" ...

    ## 目前已有的框架

    A review of bioinformatics pipeline framework 的作者对已有的工具进行很好的分类

    工具推荐

    作者的看法:

    1. implicit,也就是Make rule语法更适合用于整合不同执行工具
    2. 基于配置的流程更加稳定,也比较适合用于集群分配任务。

    最后作者建议是:

    • 如果实验室既不是纯粹的生物学试验(不需要workbench这种UI界面),也不需要高性能基于类的流程设计, 不太好选, 主要原则是投入和产出比
    • 如果实验室进行的是重复性的研究,那么就需要对数据和软件进行版本控制, 建议是 configuration-based pipelines
    • 如果实验室做的是探索性的概念证明类工作(exploratory proofs-of-concept),那么需要的是 DSL-based pipeline。
    • 如果实验室用不到高性能计算机(HPC),只能用云服务器,就是server-based frameworks.

    目前已有的流程可以在awesome-pipeline 进行查找。

    就目前来看,pipeline frameworks & library 这部分的框架中 nextflow 是点赞数最多的生物学相关框架。只可惜nextflow在运行时需要创建fifo,而在NTFS文件系统上无法创建,所以我选择 snakemake , 一个基于Python写的DSL流程框架。

    环境准备

    为了能够顺利完成这部分的教程,请准备一个Linux环境,如果使用Windows,则按照biostarhandbook(一)分析环境和数据可重复 部署一个虚拟机,并安装miniconda3。

    如下步骤会下载所需数据,并安装所需要的软件,并且启动工作环境。

    wget https://bitbucket.org/snakemake/snakemake-tutorial/get/v3.11.0.tar.bz2
    tar -xf v3.11.0.tar.bz2 --strip 1
    cd snakemake-snakemake-tutorial-623791d7ec6d
    conda env create --name snakemake-tutorial --file environment.yaml
    source activate snakemake-tutorial
    # 退出当前环境
    source deactivate
    

    当前环境下的所有文件

    ├── data
    │   ├── genome.fa
    │   ├── genome.fa.amb
    │   ├── genome.fa.ann
    │   ├── genome.fa.bwt
    │   ├── genome.fa.fai
    │   ├── genome.fa.pac
    │   ├── genome.fa.sa
    │   └── samples
    │       ├── A.fastq
    │       ├── B.fastq
    │       └── C.fastq
    ├── environment.yaml
    └── README.md
    

    基础:一个案例流程

    如果你编译过软件,那你应该见过和用过make, 但是你估计也没有仔细想过make是干嘛用的。Make是最常用的软件构建工具,诞生于1977年,主要用于C语言的项目,是为了处理编译时存在各种依赖关系,尤其是部分文件更新后,Make能够重新生成需要更新的文件以及其对应的文件。

    Snakemake和Make功能一致,只不过用Python实现,增加了许多Python的特性,并且和Python一样非常容易阅读。下面将使用Snakemake写一个变异检测流程。

    第一步:序列比对

    Snakemake非常简单,就是写各种rule来完成不同的任务。我们第一条rule就是将序列比对到参考基因组上。如果在命令行下就是bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam。 但是按照Snakemake的规则就是下面的写法。

    # 用你擅长的文本编辑器
    vim Snakefile
    # 编辑如下内容
    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/A.fastq"
        output:
            "mapped_reads/A.bam"
        shell:
            """
            bwa mem {input} | samtools view -Sb - > {output}
            """
    

    解释一下:这几行定义了一个规则(rule),在这个规则下,输入(input)有两个,而输出(output)只有一个,在shell中运行命令,只不过里面的文件都用{}形式替代。伪执行一下:snakemake -np mapped_reads/A.bam检查一下是否会出错,真实运行情况如下(不带规则,默认执行第一个规则):

    run snakemake

    第二步:推广序列比对规则

    如果仅仅是上面这样子处理一个文件,还无法体现snakemake的用途,毕竟还不如手动敲代码来的方便。snakemake的一个有点在于它能够使用文件名通配的方式对一类文件进行处理。将上面的A改成{sample},就可以将符合*.fastq的文件处理成*.bam.

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/{sample}.fastq"
        output:
            "mapped_reads/{sample}.bam"
        shell:
            """
            bwa mem {input} | samtools view -Sb - > {output}
            """
    

    那么,用snakemake -np mapped_reads/{A,B,C}.bam,就会发现,他非常机智的就比对了B.fastqC.fastq而不会再比对一遍A.fastq, 也不需要你写一堆的判断语句去手动处理。

    规则统配

    当然,如果你用touch data/samples/A.fastq改变A.fastq的时间戳,他就会认位A.fastq文件发生了改变,那么重复之前的命令就会比对A.fastq。

    第三步:比对后排序

    比对后的文件还需要进一步的排序,才能用于后续分析,那么规则该如何写呢?

    rule samtools_sort:
        input:
            "mapped_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam"
        shell:
            "samtools sort -T sorted_reads/{wildcards.sample}"
            " -O bam {input} > {output}"
    

    以之前的输出作为输出文件名,输出到另一个文件夹中。和之前的规则基本相同,只不过这里用到了wildcards.sample来获取通配名用作-T的临时文件的前缀sample实际名字。

    运行snakemake -np sorted_reads/B.bam,你就会发现他就会非常智能的先比对再排序。这是因为snakemake会自动解决依赖关系,并且按照依赖的前后顺序进行执行。

    第四步: 建立索引和对任务可视化

    这里我们再写一个规则,对之前的排序后的BAM文件建立索引。

    rule samtools_index:
        input:
            "sorted_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam.bai"
        shell:
            "samtools index {input}"
    

    目前已经写了三个规则,那么这些规则的执行和依赖关系如何呢? snakemake提供了--dag选项用于dot命令进行可视化

    snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
    

    运行流程

    第五步:基因组变异识别

    基因组变异识别需要整合之前所有的BAM文件,你可能会打算这样写

    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bamA="sorted_reads/A.bam"
            bamB="sorted_reads/B.bam"
            baiA="sorted_reads/A.bam.bai"
            baiB="sorted_reads/B.bam.bai"
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bamA} {input.bamB} | "
            "bcftools call -mv - > {output}"
    

    这样写的却没有问题,但是以后每多一个样本就需要多写一个输入,太麻烦了。这里就体现出Snakemake和Python所带来的特性了,我们可以用列表推导式的方法搞定。

    ["sorted_reads/{}.bam".format(sample) for sample in ["A","B"]]
    

    进一步,可以在规则外定义SAMPLES=["A","B"],则规则内的输入可以写成bam=["sorted_reads/{}.bam".format(sample) for sample in SAMPLES]. 由于列表推导式比较常用,但是写起来有点麻烦,snakemake定义了expand进行简化, 上面可以继续改写成expand("sorted_reads/{sample}.bam", sample=SAMPLES)

    那么最后的规则就是

    SAMPLES=["A","B"]
    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
            bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv - > {output}"
    

    小练习: 请用snakemake生成当前的DAG图。

    第六步:编写报告

    上面都是在规则里执行shell脚本,snakemake的一个优点就是可以在规则里面写Python脚本,只需要把shell改成run,此外还不需要用到引号。

    rule report:
        input:
            "calls/all.vcf"
        output:
            "report.html"
        run:
            from snakemake.utils import report
            with open(input[0]) as vcf:
                n_calls = sum(1 for l in vcf if not l.startswith("#"))
    
            report("""
            An example variant calling workflow
            ===================================
    
            Reads were mapped to the Yeast
            reference genome and variants were called jointly with
            SAMtools/BCFtools.
    
            This resulted in {n_calls} variants (see Table T1_).
            """, output[0], T1=input[0])
    

    这里还用到了snakemake的一个函数,report,可以对markdown语法进行渲染生成网页。

    第七步:增加目标规则

    之前运行snakemake都是用的snakemake 目标文件名, 除了目标文件名外,snakemake还支持规则名作为目标。通常我们按照习惯定义一个all规则,来生成结果文件。

    rule all:
        input:
            "report.html
    

    基础部分小结:

    总结下学习过程,知识点如下:

    • Snakemake基于规则执行命令,规则一般由input, output,shell三部分组成。
    • Snakemake可以自动确定不同规则的输入输出的依赖关系,根据时间戳来判断文件是否需要重新生成
    • Snakemake以{sample}.fa形式进行文件名通配,用{wildcards.sample}获取sample的实际文件名
    • Snakemake用expand()生成多个文件名,本质是Python的列表推导式
    • Snakemake可以在规则外直接写Python代码,在规则内的run里也可以写Python代码。
    • Snakefile的第一个规则通常是rule all,因为默snakemake默认执行第一条规则

    进阶:对流程进一步修饰

    在基础部分中,我们完成了流程的框架,下一步则是对这个框架进行不断完善,比如说编写配置文件,声明不同rule的消耗资源,记录运行日志等。

    第一步: 声明所需进程数

    对于一些工具,比如说bwa,多进程或者多线程运行能够大大加速计算。snakemake使用threads来定义当前规则所用的进程数,我们可以对之前的bwa_map增加该指令。

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/{sample}.fastq"
        output:
            "mapped_reads/{sample}.bam"
        threads:8
        shell:
            "bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
    

    声明threads后,Snakemake任务调度器就会在程序运行的时候是否并行多个任务。这主要和参数中的--cores相关。比如说

    snakemake --cores 10
    

    由于总体上就分配了10个核心,于是一次就只能运行一个需要消耗8个核心的bwa_map。但是当其中一个bwa_map运行完毕,这个时候snakemaek就会同时运行一个消耗8个核心的bwa_map和没有设置核心数的samtools_sort,来保证效率最大化。因此对于需要多线程或多进程运行的程序而言,将所需的进程单独编码,而不是硬编码到shell命令中,能够更有效的使用资源。

    第二步:配置文件

    之前的SAMPLES写在了snakefile,也就是意味这对于不同的项目,需要对snakefile进行修改,更好的方式是用一个配置文件。配置文件可以用JSON或YAML语法进行写,然后用configfile: "config.yaml"读取成字典,变量名为config。

    config.yaml内容为:

    samples:
        A: data/samples/A.fastq
        B: data/samples/B.fastq
    

    YAML使用缩进表示层级关系,其中缩进必须用空格,但是空格数目不重要,重要的是所今后左侧对齐。上面的YAML被Pytho读取之后,以字典保存,形式为{'samples': {'A': 'data/samples/A.fastq', 'B': 'data/samples/B.fastq'}}

    而snakefile也可以改写成

    configfile: "config.yaml"
    ...
    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bam=expand("sorted_reads/{sample}.bam", sample=config["samples"]),
            bai=expand("sorted_reads/{sample}.bam.bai", sample=config["smaples])
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv - > {output}"
    

    虽然sample是一个字典,但是展开的时候,只会使用他们的key值部分。

    关于YAML格式的教程,见阮一峰的博客:http://www.ruanyifeng.com/blog/2016/07/yaml.html

    第三步:输入函数

    既然已经把文件路径都存入到配置文件中,那么可以进一步的改写之前的bwa_map里的输入部分。也就是从字典里面提取到存放的路径。最开始我就是打算这样写

    rule bwa_map:
        input:
            "data/genome.fa",
            config['samples']["{sample}"]
        output:
            "mapped_reads/{sample}.bam"
        threads:8
        shell:
            "bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
    

    毕竟"{sample}"从理论上应该得到sample的名字。但是snakemake -np显示出现错误

    KeyError in line 11 of /home6/zgxu/snakemake-snakemake-tutorial-623791d7ec6d/Snakefile:
    '{sample}'
    

    这可能是{sample}的形式只能在匹配的时候使用,而在获取值的时候应该用基础第三步的wildcards.sample形式。于是继续改成config["samples"][wildcards.sample]。然而还是出现了错误。

    name 'wildcards' is not defined
    

    为了理解错误的原因,并找到解决方法,我们需要理解Snakemake工作流程执行的一些原理,它执行分为三个阶段

    • 初始化阶段,工作流程会被解析,所有规则都会被实例化
    • DAG阶段,也就是生成有向无环图,确定依赖关系的时候,所有的通配名部分都会被真正的文件名代替。
    • 调度阶段,DAG的任务按照顺序执行

    也就是说在初始化阶段,我们是无法获知通配符所指代的具体文件名,必须要等到第二阶段,才会有wildcards变量出现。也就是说之前的出错的原因都是因为第一个阶段没通过。这个时候就需要输入函数推迟文件名的确定,可以用Python的匿名函数,也可以是普通的函数

    rule bwa_map:
        input:
            "data/genome.fa",
            lambda wildcards: config["samples"][wildcards.sample]
        output:
            "mapped_reads/{sample}.bam"
        threads: 8
        shell:
            "bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
    

    第四步:规则参数

    有些时候,shell命令不仅仅是由input和output中的文件组成,还需要一些静态的参数设置。如果把这些参数放在input里,则会因为找不到文件而出错,所以需要专门的params用来设置这些参数。

    rule bwa_map:
        input:
            "data/genome.fa",
            lambda wildcards: config["samples"][wildcards.sample]
        output:
            "mapped_reads/{sample}.bam"
        threads: 8
        params:
            rg="@RG\tID:{sample}\tSM:{sample}"
        shell:
            "bwa mem -R '{params.rg}' '-t {threads} {input} | samtools view -Sb - > {output}"
    

    写在rule中的params的参数,可以在shell命令中或者是run里面的代码进行调用。

    第五步: 日志文件

    当工作流程特别的大,每一步的输出日志都建议保存下来,而不是输出到屏幕,这样子出错的时候才能找到出错的所在。snakemake非常贴心的定义了log,用于记录日志。好处就在于出错的时候,在log里面定义的文件是不会被snakemake删掉,而output里面的文件则是会被删除。继续修改之前的bwa_map.

    rule bwa_map:
        input:
            "data/genome.fa",
            lambda wildcards: config["samples"][wildcards.sample]
        output:
            "mapped_reads/{sample}.bam"
        params:
            rg="@RG\tID:{sample}\tSM:{sample}"
        log:
            "logs/bwa_mem/{sample}.log"
        threads: 8
        shell:
            "(bwa mem -R '{params.rg}' -t {threads} {input} | "
            "samtools view -Sb - > {output}) 2> {log}"
    

    这里将标准错误重定向到了log中。

    第六步:临时文件和受保护的文件

    由于高通量测序的数据量通常很大,因此很多无用的中间文件会占据大量的磁盘空间。而特异在执行结束后写一个shell命令清除不但写起来麻烦,而且也不好管理。Snakemake使用temp()来将一些文件标记成临时文件,在执行结束后自动删除。

    rule bwa_map:
        input:
            "data/genome.fa",
            lambda wildcards: config["samples"][wildcards.sample]
        output:
            temp("mapped_reads/{sample}.bam")
        params:
            rg="@RG\tID:{sample}\tSM:{sample}"
        log:
            "logs/bwa_mem/{sample}.log"
        threads: 8
        shell:
            "(bwa mem -R '{params.rg}' -t {threads} {input} | "
            "samtools view -Sb - > {output}) 2> {log}"
    

    修改之后的代码,当samtools_sort运行结束后就会把"mapped_reads"下的BAM删掉。同时由于比对和排序都比较耗时,得到的结果要是不小心被误删就会浪费大量计算时间,最后的方法就是用protected()保护起来

    rule samtools_sort:
        input:
            "mapped_reads/{sample}.bam"
        output:
            protected("sorted_reads/{sample}.bam")
        shell:
            "samtools sort -T sorted_reads/{wildcards.sample} "
            "-O bam {input} > {output}"
    

    最后,snakemake就会在文件系统中对该输出文件写保护,也就是最后的权限为-r--r--r--, 在删除的时候会问你rm: remove write-protected regular file ‘A.bam’?.

    进阶部分小结

    • 使用threads:定义不同规则所需线程数,有利于snakemake全局分配任务,最优化任务并行
    • 使用configfile:读取配置文件,将配置和流程分离
    • snakemake在DAG阶段才会知道通配的具体文件名,因此在input和output出现的wildcards就需要推迟到第二步。
    • log里定义的日志文件,不会因任务失败而被删除
    • params定义的参数,可以在shell和run中直接调用
    • temp()中的文件运行结束后会被删除,而protected()中的文件会有写保护,避免意外删除。

    高级:实现流程的自动部署

    上面的分析流程都是基于当前环境下已经安装好要调用的软件,如果你希望在新的环境中也能快速部署你的分析流程,那么你需要用到snakmake更高级的特性,也就是为每个rule定义专门的运行环境。

    全局环境

    我建议你在新建一个snakemake项目时,都先用conda create -n 项目名 python=版本号创建一个全局环境,用于安装一些常用的软件,例如bwa、samtools、seqkit等。然后用如下命令将环境导出成yaml文件

    conda env export -n 项目名 -f environment.yaml
    

    那么当你到了一个新的环境,你就可以用下面这个命令重建出你的运行环境

    conda env create -f environment.yaml
    

    局部环境

    当然仅仅依赖于全局环境或许还不够,对于不同的规则(rule)可能还有Python2和Python3的区别,所以你还得为每个规则创建环境。

    snakemake有一个参数--use-conda,会解析rule中的conda规则,根据其提供的yaml文件安装特定版本的工具,以基础第一步的序列比对为例,

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/A.fastq"
        output:
            "mapped_reads/A.bam"
        conda:
            "envs/map.yaml"
        shell:
            """
            mkdir -p mapped_reads
            bwa mem {input} | samtools view -Sb - > {output}
            """
    

    随后在snakemake执行的目录下创建envs文件夹,增加map.yaml, 内容如下

    name: map
    channels:
      - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
      - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
      - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
      - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
      - defaults
    dependencies:
      - bwa=0.7.17
      - samtools=1.9
    show_channel_urls: true
    

    注意: YAML文件的name行不是必要的,但是建议加上。

    那么当你用snakmake --use-conda执行时,他就会在.snakemake/conda下创建专门的conda环境用于处理当前规则。对于当前项目,该conda环境创建之后就会一直用于该规则,除非yaml文件发生改变。

    如果你希望在实际运行项目之前先创建好环境,那么可以使用--create-envs-only参数。

    由于默认情况下,每个项目运行时只会在当前的.snakemake/conda查找环境或者安装环境,所以在其他目录执行项目时,snakemake又会重新创建conda环境,如果你担心太占地方或者环境太大,安装的时候太废时间,你可以用--conda-prefix指定专门的文件夹。

    代码总结

    最后的代码如下

    configfile: "config.yaml"
    
    
    rule all:
        input:
            "report.html"
    
    
    rule bwa_map:
        input:
            "data/genome.fa",
            lambda wildcards: config["samples"][wildcards.sample]
        output:
            temp("mapped_reads/{sample}.bam")
        params:
            rg="@RG\tID:{sample}\tSM:{sample}"
        log:
            "logs/bwa_mem/{sample}.log"
        threads: 8
        shell:
            "(bwa mem -R '{params.rg}' -t {threads} {input} | "
            "samtools view -Sb - > {output}) 2> {log}"
    
    
    rule samtools_sort:
        input:
            "mapped_reads/{sample}.bam"
        output:
            protected("sorted_reads/{sample}.bam")
        shell:
            "samtools sort -T sorted_reads/{wildcards.sample} "
            "-O bam {input} > {output}"
    
    
    rule samtools_index:
        input:
            "sorted_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam.bai"
        shell:
            "samtools index {input}"
    
    
    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bam=expand("sorted_reads/{sample}.bam", sample=config["samples"]),
            bai=expand("sorted_reads/{sample}.bam.bai", sample=config["samples"])
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv - > {output}"
    
    
    rule report:
        input:
            "calls/all.vcf"
        output:
            "report.html"
        run:
            from snakemake.utils import report
            with open(input[0]) as vcf:
                n_calls = sum(1 for l in vcf if not l.startswith("#"))
    
            report("""
            An example variant calling workflow
            ===================================
    
            Reads were mapped to the Yeast
            reference genome and variants were called jointly with
            SAMtools/BCFtools.
    
            This resulted in {n_calls} variants (see Table T1_).
            """, output[0], T1=input[0])
    

    执行snakemake

    写完Snakefile之后就需要用snakemake执行。snakemake的选项非常多,这里列出一些比较常用的运行方式。

    运行前检查潜在错误:

    snakemake -n
    snakemake -np
    snakemake -nr
    # --dryrun/-n: 不真正执行
    # --printshellcmds/-p: 输出要执行的shell命令
    # --reason/-r: 输出每条rule执行的原因
    

    直接运行:

    snakemake
    snakemake -s Snakefile -j 4
    # -s/--snakefile 指定Snakefile,否则是当前目录下的Snakefile
    # --cores/--jobs/-j N: 指定并行数,如果不指定N,则使用当前最大可用的核心数
    

    强制重新运行:

    snakemake -f
    # --forece/-f: 强制执行选定的目标,或是第一个规则,无论是否已经完成
    snakemake -F
    # --forceall/-F: 也是强制执行,同时该规则所依赖的规则都要重新执行
    snakemake -R some_rule
    # --forecerun/-R TARGET: 重新执行给定的规则或生成文件。当你修改规则的时候,使用该命令
    

    可视化:

    snakemake --dag  | dot -Tsvg > dag.svg
    snakemake --dag  | dit -Tpdf > dag.pdf
    # --dag: 生成依赖的有向图
    snakemake --gui 0.0.0.0:2468
    # --gui: 通过网页查看运行状态
    

    集群执行:

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

    参考资料

    展开全文
  • 使用本包时请引用以下论文:王大鹏。 hppRNA——一种基于 Snakemake 的方便的无参数管道,用于对大量样本进行 RNA-Seq 分析。 生物信息学简报,第 19 卷,第 4 期,2018 年 7 月 20 日,第 622-626 页。
  • Snakemake教程-01基础部分

    千次阅读 2021-07-04 10:42:25
    snakemake简介1.1 初识snakemake2. Snakemake教程2.1安装2.2 Basic: An example workflow2.2.1 Step 1: Mapping readsStep 2: Generalizing the read mapping ruleStep 3: Sorting read alignmentsStep 4: Indexing...

    1. snakemake简介

    目标:流程管理工具,创建reproducible和scalable的数据分析。

    特点:

    • a human readable, Python based language.
    • 可以无缝地拓展到服务器、集群、网格和云环境,不需要修改工作流定义。
    • 可包含所需软件的描述,这些软件可以自动部署到任何执行环境中。

    1.1 初识snakemake

    参考:https://snakemake.github.io/

    特点解释
    Readability(可读性)易于阅读、适应性强、功能强大。每个规则描述了分析中的一个步骤,该步骤定义了如何从输入文件获取输出文件。规则之间的依赖关系是自动确定的。
    Portability(便捷性)通过整合conda包管理工具和container虚拟化技术,所有软件依赖项在执行时都可以自动部署。
    Modularization(模块化)通过整合script和jupyter notebook。很容易创建和使用可重复的工具包装器,并将数据分析分解为良好分割的模块。
    Transparency(透明化)自动的,交互式的,自包含的报告确保了从结果到使用的步骤、参数、代码和软件的完整透明性。
    Scalability(可拓展性)工作流可以无缝地从单核拓展到多核、集群或云,不需要修改工作流定义和自动避免冗余计算。

    2. Snakemake教程

    链接:https://snakemake.readthedocs.io/en/stable/tutorial/tutorial.html#tutorial

    Snakemake提供的符号拓展,维持了流程定义的属性。

    Snakemake整合了包管理工具和容器工具Singularity,使得定义软件堆栈变成了软件的流程的一部分。

    2.1安装

    安装snakemake需要以下软件:

    Python ≥3.5
    Snakemake ≥5.24.1
    BWA 0.7
    SAMtools 1.9
    Pysam 0.15
    BCFtools 1.9
    Graphviz 2.42
    Jinja2 2.11
    NetworkX 2.5
    Matplotlib 3.3
    

    本地安装:

    wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh
    bash Mambaforge-Linux-x86_64.sh  ##注意是否加入PATH时,选yes。
    
    mkdir snakemake-tutorial
    cd snakemake-tutorial
    
    ##开始着手创建流程。
    ##首先下载示例数据。
    $ wget https://github.com/snakemake/snakemake-tutorial-data/archive/v5.24.1.tar.gz
    $ tar --wildcards -xf v5.24.1.tar.gz --strip 1 "*/data" "*/environment.yaml"
    
    ##创建独立的虚拟环境
    $ conda activate base
    $ mamba env create --name snakemake-tutorial --file environment.yaml
    $ conda install -n base -c conda-forge mamba
    $ conda activate snakemake-tutorial
    $ snakemake --help
    $ conda deactivate
    

    报错:

    ERROR   Could not write out repodata file
    '/share/inspurStorage/home1/jialh/mambaforge/pkgs/cache/602d3932.json': No such file or directory
    

    解决办法:https://zhuanlan.zhihu.com/p/380567007

    经多次尝试: rm -rf /Users/Likey/ProgramFiles/mambaforge/pkgs/cache/*

    删除缓存文件,再次安装即可,再次安装时最好开启外网或者配置国内源环境避免其它问题.

    2.2 Basic: An example workflow

    Snakemake流程通过Snakefile中的特定rules来定义。Rules可以将workflow分解为small steps,通过指定如何从input files集合创建output files集合。

    下面的流程来自genome analysis领域。本教程不需要你知道genome analysis。但是我们将在接下来的段落中提供一些背景。

    2.2.1 Step 1: Mapping reads

    Snakemake文件内容如下

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/A.fastq"
        output:
            "mapped_reads/A.bam"
        shell:
            "bwa mem {input} | samtools view -Sb - > {output}"
    

    运行方式为:

    snakemake --cores 4 mapped_reads/A.bam
    

    无论何时执行工作流程,你都可以指定使用的cores数目。本教程,我们只使用单个核。接下来,你将会看到并行化工作。注意,执行完上述命令后,Snakemake将不会再次创建mapped_reads/A.bam。 Snakemake only re-runs jobs if one of the input files is newer than one of the output files or one of the input files will be updated by another job.

    Step 2: Generalizing the read mapping rule

    Snakemake可以generalizing rules by using named wildcards(即通配符). 仅仅使用通配符{sample}替换上述输入文件和输出文件中的A即可。

    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/{sample}.fastq"
        output:
            "mapped_reads/{sample}.bam"
        shell:
            "bwa mem {input} | samtools view -Sb - > {output}"
    

    Step 3: Sorting read alignments

    对于后续步骤,我们需要对读段比对BAM文件排序。这可以通过samtools sort命令来实现。我们可以在bwa_map下面增加如下规则:

    rule samtools_sort:
        input:
            "mapped_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam"
        shell:
            "samtools sort -T sorted_reads/{wildcards.sample} "
            "-O bam {input} > {output}"
    

    该规则会从mapped_reads获取输入文件,将排序后的版本存入sorted_reads目录。注意Snakemake会自动创建缺失的目录。对于排序,samtools需要通过字符T来指定前缀。此处,我们需要通配符sample。Snakemake也可以通过wildcards对象的value属性在shell命令中获取通配符。

    Step 4: Indexing read alignments and visualizing the DAG of jobs

    接下来,我们需要使用samtools来对排序后的读段比对建立索引,这样我们可以快速通过读段比对到基因组上的位置,快速获取读段(reads)。如下所示:

    rule samtools_index:
        input:
            "sorted_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam.bai"
        shell:
            "samtools index {input}"
    

    注意:Snakemake使用Python format mini language来形成shell命令。有时候,你必须使用花括号{}。在那种情形下,你必须避开通过双花括号来避开它们,例如:ls {{A, B}}.txt

    共有三个步骤,可以创建结果的有向无环图(directed acyclic graph,DAG)。 可以执行:

    $ snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
    

    创建visualization of the DAG可以通过Graphviz的dot命令来实现。对于给定的目标文件,Snakemake可以用dot语言指定DAG,将其输入dot命令, 这可以被渲染为SVG格式。渲染的DAG会保存到dag.svg文件,如下图所示:

    image

    Step 5: Calling genomic variants

    接下来,我们将所有mapped reads聚合到一起,联合call genomic variants。对于the variant calling, 我们将结合两个工具samtoolsbcftools。Snakemake提供了helper function for collecting input files, 能够帮助我们来描述这一步的聚合。使用:

    expand("sorted_reads/{sample}.bam", sample=SAMPLES)
    

    我们将获取一个文件列表,其中"sorted_reads/{sample}.bam"模式会通过给定的sample列表SAMPLES进行格式化。例如:

    ["sorted_reads/A.bam", "sorted_reads/B.bam"]
    

    这个功能对于包括多个通配符的模式非常有效,例如:

    expand("sorted_reads/{sample}.{replicate}.bam", sample=SAMPLES, replicate=[0, 1])
    

    将使用所有元素SAMPLES和列表[0, 1]来创建结果,得到:

    ["sorted_reads/A.0.bam", "sorted_reads/A.1.bam", "sorted_reads/B.0.bam", "sorted_reads/B.1.bam"]
    

    此处,我们仅使用expand的简单形式。我们首先让Snakemake知道我们想考虑什么。记住,Snakemake从结果文件回溯,而不是从输入文件。因此,它不推断所有可能的输出形式,例如数据目录中的fastq文件。也要记住,Snakemake主要是增强了声明性语句的Python代码,来定义工作流。因此,我们可以在Snakemake的顶部用普通的Python特别定义的sample列表。

    SAMPLES = ["A", "B"]
    

    接下来,我们可以从config file中学习更多更复杂的形式。但是现在这足以让我们增加下面的规则到Snakefile中:

    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
            bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv - > {output}"
    

    使用多个输入和输出文件,有时在shell命令中分别引用它们是很方便的。这可以通过specifying names for input or output files来实现,例如用fa=...。这些文件能够被shell通过名字引用,例如使用{input.fa}。对于long shell commands, 建议split the string over multiple indented lines。python可以自动将其合并。 进一步,你将注意到input or output file list can contain arbitary Python statements, 只要其放回字符串,或者字符串列表。此处,我们用expand函数来聚合所有样本比对上的reads。

    注意: splitting long shell commands across multiple lines, make sure to include whitespace at the end of each line. 因为每一行的字符串只是简单地连接在一起,否则可能导致错误。例如,不能将当前的例子分割为如下命令:

    "samtools mpileup"
    "-g -f {input.fa} {input.bam}"
    

    这将会连接为:

    "samtools mpileup-g -f {input.fa} {input.bam}"
    

    导致抛出错误:

    [main] unrecognized command 'mpileup-g'
    

    练习:
    生成DAG图,代码为:

    snakemake --dag calls/all.vcf | dot -Tsvg > dag_vcf.svg
    

    结果如下图所示:

    image

    Step 6: Using custom scripts

    通常,流程不仅仅是由各种工具组成,也包含手工代码来计算统计结果或者画图等等。尽管Snakemake允许你直接写python代码,但是通常将这些逻辑移到单独的脚本更合理。出于这个目的,Snakemake提供了script指令。增加如下命令到你的snakefile文件:

    rule plot_quals:
        input:
            "calls/all.vcf"
        output:
            "plots/quals.svg"
        script:
            "scripts/plot-quals.py"
    

    使用这个规则,我们将创建质量打分的柱状图(histogram), 其已经分配到了calls/all.vcf的变异检测中。实际上python代码将通过隐藏的脚本scripts/plot-quals.py来产生柱状图。脚本的路径通常是相对Snakefile的。在脚本中,规则的所有属性像input, output, wildcards等snakemake全局对象属性也可以使用。 使用如下内容,产生scripts/plot-quals.py文件:

    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from pysam import VariantFile
    
    quals = [record.qual for record in VariantFile(snakemake.input[0])]
    plt.hist(quals)
    
    plt.savefig(snakemake.output[0])
    

    除了python代码,可以使用R代码。详情参考:https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#snakefiles-external-scripts

    Step 7: Adding a target rule

    到目前为止,我们总是通过指定命令行中目标文件的名字来执行流程。除了文件名,Snakemake also accepts rule names as targets,如果访问的规则不包括通配符。因此,可以编写收集所需结果的特定子集或所有结果的目标规则。而且,如果命令行没有给目标,Snakemake也可以定义Snakefile的first rule作为目标。 因此最好是在workflow顶部设置all规则,将所有目标输出文件作为输入文件。

    此处,我们可以增加规则:

    rule all:
        input:
            "plots/quals.svg"
    

    到workflow的顶部。 然后执行Snakemake:

    snakemake -n --cores 4
    

    练习:

    • 绘制完整流程的DAG图。
    snakemake --dag plots/quals.svg | dot -Tsvg > dag_quals.svg
    

    0265b54f5d18681deebd5b47acb1bf7

    • 执行完整的流程,看一看结果plots/quals.svg
    snakemake -n --cores 4
    

    quals

    • 查看snakemake --help, 选取部分参数,如下所示:
    -h, --help  ##展示帮助文档。
      --dry-run, --dryrun, -n  ##不执行任何东西,仅仅展示流程会干什么。
      --profile PROFILE        ##用于参数化Snakemake的profile名字。
      --cache [RULE [RULE ...]] ##通过环境变量$SNAKEMAKE_OUTPUT_CACHE, 存储给定规则输出文件的中间cache。
      --snakefile FILE, -s FILE ##snakefile形式的工作流程定义。
      --cores [N], -c [N]   ##最多使用多少个CPU cores/jobs用于并行化。
      --jobs [N], -j [N]    ##最多使用多少个CPU cluster/cloud jobs用于并行化。
      --local-cores N       ##在使用cluster/cloud模式时,主机使用多少个cores来并行化。
      --resources [NAME=INT [NAME=INT ...]], --res [NAME=INT [NAME=INT ...]] ##定义额外的资源。
      --set-threads RULE=THREADS [RULE=THREADS ...] ##重写threads规则。
      --max-threads MAX_THREADS   ##定义任何工作的最大threads数目。
      --set-resources RULE:RESOURCE=VALUE [RULE:RESOURCE=VALUE ...] ##重写rules使用的资源。
      --config [KEY=VALUE [KEY=VALUE ...]], -C [KEY=VALUE [KEY=VALUE ...]] ##设置和重写流程config对象的值。
      --configfile FILE [FILE ...], --configfiles FILE [FILE ...] ##指定或重新流程config file(见文本)。
      --envvars VARNAME [VARNAME ...] ##传递到cloud jobs的环境变量数。
      --directory DIR, -d DIR         ##指定工作路径。
      --touch, -t           ##touch输出文件(更新日期,不需要真正改变它们)。
      --keep-going, -k      ##如果一个jobs失败,与之相关的jobs是否继续执行。
      --force, -f           ##忽略所有输出,从第一个rule或者选定的target来强制执行。
      --forceall, -F        ##忽略所有输出,强制执行选定的(或第一个)rule和所有其依赖的rules。
      --forcerun [TARGET [TARGET ...]], -R [TARGET [TARGET ...]] ##强制重新执行给定的规则或者产生给定的文件。使用这个选项,如果你改变了一个rule, 想要你的流程中所有它的输出都更新。
      --prioritize TARGET [TARGET ...], -P TARGET [TARGET ...] ##告诉调度器为给定目标(及其所有依赖项)的创建分配最高优先级。
      --batch RULE=BATCH/BATCHES  ##为给定RULE的输入文件创建BATCH。
      --until TARGET [TARGET ...], -U TARGET [TARGET ...] ##重新运行此rule, 直到其到达指定的rules或者files。
      --omit-from TARGET [TARGET ...], -O TARGET [TARGET ...] ##防止执行的规则或者创建的文件,是DAG中这些目标的下游。
      --rerun-incomplete, --ri ##重新运行所有输出被认为是不完整的任务。
    

    使用--forcerun重新执行samtools_sort规则,看看会发生什么。

    snakemake --cores 4  --forcerun samtools_sort
    

    流程重新执行了samtools_sort及其后面的所有步骤。

    • 使用--reason,它能够展示每个job执行的理由。尝试使用这个指令以及dry-run和--forcerun,理解snakemake的决策。
    snakemake --cores 1  --forcerun samtools_sort -n --reason
    

    执行结果,会产生reason部分,说明为什么会执行这一步。

    clipboard

    Summary

    总而言之,结果流程如下所示:

    SAMPLES = ["A", "B"]
    
    
    rule all:
        input:
            "plots/quals.svg"
    
    
    rule bwa_map:
        input:
            "data/genome.fa",
            "data/samples/{sample}.fastq"
        output:
            "mapped_reads/{sample}.bam"
        shell:
            "bwa mem {input} | samtools view -Sb - > {output}"
    
    
    rule samtools_sort:
        input:
            "mapped_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam"
        shell:
            "samtools sort -T sorted_reads/{wildcards.sample} "
            "-O bam {input} > {output}"
    
    
    rule samtools_index:
        input:
            "sorted_reads/{sample}.bam"
        output:
            "sorted_reads/{sample}.bam.bai"
        shell:
            "samtools index {input}"
    
    
    rule bcftools_call:
        input:
            fa="data/genome.fa",
            bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
            bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
        output:
            "calls/all.vcf"
        shell:
            "samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv - > {output}"
    
    
    rule plot_quals:
        input:
            "calls/all.vcf"
        output:
            "plots/quals.svg"
        script:
            "scripts/plot-quals.py"
    
    展开全文
  • Snakemake包装器存储库 Snakemake包装器存储库是可重复使用包装器的集合,这些包装器可用于快速使用规则和工作流中流行的命令行工具。 有关更多信息,请访问 。
  • 但这还不是全部,它们确保以有组织的方式捕获所有分析日志,它们明确概述了所使用的软件(和确切的软件版本)、每个步骤的输入和输出文件。 最后,当您的数据不可避免地变成大数据时,您可以轻松地从在笔记本电脑上...
  • shared_snakemake 这是Snakefiles的存储库,其中包含可以插入到各种项目中的常用snakemake规则和工作流。... 使用-o {log.stdout} -e {log.stderr}指定时qsub在呼叫--cluster的参数snakemake 。 添加新规
  • ribosomal_snakemake

    2021-04-22 20:51:02
    ribosomal_snakemake 核糖体蛋白树生成 使用15个核糖体蛋白序列生成核糖体蛋白系统发育树的工作流程。 可以使用蛋白质或DNA序列来构建树,将DNA序列用于更紧密相关的基因组,以及将蛋白质序列用于更趋异的基因组可能...
  • snakePipes是使用构建的灵活而强大的工作流程,可简化NGS数据的分析。 可用的工作流程 DNA映射* 芯片序列* mRNA序列* 非编码RNA-seq * ATAC序列* 核糖核酸序列 嗝 全基因组亚硫酸氢盐Seq / WGBS (*在“等位...
  • Snakemake搭建流程 - 干货级

    千次阅读 2020-08-05 16:20:11
    snakemake 使用说明4.1 定义workflow4.2 配置信息Configuration4.3 运行代码示例参考 1. snakemake 简介 最近碰到一个snakemake搭建的流程,挺好奇,便学习了一波,在此特分享一些体会和心得,仅供想快速入手的骚年...
  • Snakemake管道可处理基于板的scATAC-seq数据 该存储库包含用于处理通过生成的scATAC-seq数据的代码。 与原始Nat相比有什么区别。 通讯出版物? 在此管道中,我们具有: 添加了config.json文件,以使处理更加灵活且...
  • snakemake_align_MAJIQ

    2021-04-12 16:00:58
    Snakemake工作流程:snakemake_align_MAJIQ 这是新的Snakemake工作流程的模板。 用涵盖目的和领域的全面描述代替此文本。 将您的代码插入相应的文件夹,即scripts , rules和envs 。 在Snakefile定义工作流程的入口...
  • snakebids:Snakemake + BIDS

    2021-03-25 23:34:20
    该软件包允许您使用Snakemake构建BIDS应用程序。 它提供: 使用PyBIDS进行灵活的数据获取,只能通过配置文件条目进行配置 帮助程序功能,用于在Snakemake工作流程/规则中创建BIDS路径 命令行调用符合BIDS App的...
  • 文章目录进阶部分:对案例进行进一步修饰Step 1: Specifying the number of used threadsStep 2: Config filesStep 3: Input functionsStep 4: ...参考:https://snakemake.readthedocs.io/en/stable/tutorial/advanced
  • 您也可以从stdin中输入文件,这会将文件打印到屏幕上,或者使用--diff或--check选项。 有关更多详细信息,请参见。 目录 安装 皮皮 pip install snakefmt conda conda install -c bioconda snakefmt 货柜 码头工人 ...
  • Snakemake 常用参数以及进阶用法介绍

    千次阅读 2021-11-12 14:03:50
    上一篇介绍了 Snakemake 入门教程 做了一个简单的示例,具体查看我的上一篇内容 下面会介绍一下 Snakemake的常用参数以及进阶的用法~ 介绍之前大家可以看一个视频了解一下(时长:19min14s, 选择性观看) ...
  • Nanopype是基于snakemake的管道,可提供便捷的牛津纳米Kong技术(ONT)数据处理和存储解决方案。 流水线分为处理部分(包括基本调用和对齐)以及分析部分(涵盖其他下游应用程序)。 部分提供了所有包含的工具的...
  • 这是我们在实验室中使用的工具,用于改进我们的wetlab协议,并提供了一个简单的框架来重现和比较具有不同参数的不同实验。 它使用STAR映射读取。 它可用于使用两个读取的任何单细胞方案,其中第一个读取保存Cell和...
  • 但是都是我一步一步跑过的流程,所以会更有印象,送给想学 Snkaemake 但是一直没有去学的朋友们,这些内容对于有生信基础的人来讲,上手会很快,因为很多的生信软件都使用过,写起来也就没有那么晦涩,下面开始~ ...
  • snakemake使用threads来定义当前规则所用的进程数, samples: A: data/samples/A.fastq B: data/samples/B.fastq configfile: "config.yaml" 虽然sample是一个字典,但是展开的时候,只会使用他们的key值部分。...
  • rules :包含snakemake的规则和所有必需的函数定义。 scripts :生成系统特定脚本的基本脚本。 您可以在此处修改特定的分析或图形渲染。 它们被编写为采用config/config.yaml文件中指定的东西的系统参数,这些参数...
  • snakemake使用
  • snakemake使用 threads 来定义当前规则所用的进程数,我们可以对之前的 bwa_map 增加该指令。 rule bwa_map: input: "data/genome.fa", "data/samples/{sample}.fastq" output: "mapped_reads/{sample}.bam" ...
  • 如果蛇文件开始变得很大,您可以尝试将它们分成较小的子集,并使用include或subworkflows 。 我倾向于将与单个工具关联的规则保留在较小的snakefile中,然后根据需要包含/排除它们。 在这里,您可以有两个文件ruleA...
  • snakemake学习笔记

    2019-06-03 21:07:00
    snakemake 是一个流程搭建的工具,这里主要用来记录一些snakemake使用方法 对于run或者shell部分的需要使用sample变量可以使用wildcards.sample来获取 对于写好的模块可以使用include来载入,然后使用rule all...
  • Snakemake的优点: - 支持多线程 - 支持断点运行 - 支持shell命令,还可以和 Python 库结合使用 行,现在我们知道 Snakemake 比 shell 更好用了。 搓手手,那么如何使用 Snakemake 呢? 来试试栗子 ~ 0. 安装 ...
  • snakemake 学习笔记1

    千次阅读 2019-02-27 14:39:38
    1, snakemake介绍 Snakemake是用Python3写的一个流程化工具, 非常方便....3, 使用cat命令, 将两者合并为hebing.txt echo "hello number1" &a
  • snakemake作为一款主流流程管理软件,对比shell,py,R写的流程它的好处是不言而喻的。基本操作和技巧可以参考官方文档和各种博客上的文章,这里主要介绍如何利用singularity容器配合snakemake。 我们搭建生信流程...

空空如也

空空如也

1 2 3 4 5 ... 11
收藏数 205
精华内容 82
关键字:

snakemake使用