• 欢迎访问daniu,分享一些有趣的事情,推荐使用最新版火狐浏览器和Chrome浏览器访问本网站

snakemake踩坑笔记

生信杂谈 daniu 3个月前 (05-11) 347次浏览 0个评论 扫描二维码

snakemake是什么?

简单来说,snakemake是一个用来编写任务流程的工具,基于python开发,其语法简单,功能强大,易于理解。是居家旅行必备之佳品。

snakemake的安装

snakemake可以通过conda和apt进行安装,下面两种方法二选一即可:

conda install -c bioconda snakemake
sudo apt install snakemake

snakemake快速上手

简单来说snakemake就是由一个个的规则组合而成的,其基本缩进规则与python相同

下面先看一个基本的示例:

SAMPLES = ["A", "B"]


rule all:
    input:
        expand('sorted_reads/{sample}.bam',sample=SAMPLES)


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}"

可以看到,snakemake的代码具有良好的可读性,即使没有任何基础的人也可以很轻易的读懂,其本身还提供了简单却够用的可视化流程帮助我们了解代码的运行逻辑,可以使用下面的方式将流程可视化

snakemake --dag | dot -Tsvg > dag.svg

snakemake执行流程的步骤是从后往前的,我们借助上面的示例代码来了解。

首先我们定义了一个样本列表SAMPLES包含我们所要处理的所有样本,expand是snakemake提供的实用函数,类似于python中的format,它将SAMPLES展开并替换文本中的占位,得到我们最终指定的输出文件:

sorted_reads/A.bam

sorted_reads/B.bam

我们如果希望在后面单独使用sample参数,可以通过{wildcards.sample}进行访问。

当程序执行时,程序会查找当前目录下是否已经存在上述文件,如果不存在则进一步查找能够输出该文件的规则,在这里samtools_sort可以输出目标文件,snakemake会进一步检查samtools_sort的输入文件是否满足,如果不满足则重复此查找过程。

snakemake高级特性

snakemake内置了一些非常使用的高级特性,可以极大的提高我们的工作效率,这些高级特性包括:

临时和写保护文件:

在整个流程中有一些文件可能是我们并不关心的中间产物,有一些文件我们可能不希望他们被轻易改动,snakemake提供了temp( )和protected( )函数实现这一点,比如在上面的例子中,如果我们对bwa_map的结果并不关心的话,则可以将其改为:

SAMPLES = ["A", "B"]


rule all:
    input:
        expand('sorted_reads/{sample}.bam',sample=SAMPLES)


rule bwa_map:
    input:
        "data/genome.fa",
        "data/samples/{sample}.fastq"
    output:
        temp("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}"

这样,bwa_map的输出将作为一个临时文件在流程结束后被删除

多线程:

snakemake对多线程有着良好的支持,我们只需要在响应的规则中使用threads指定线程数(不指定默认为1并在执行的时候使用 –cores int 指定最大线程数量,snakemake即会自动处理。

日志:

在多线程并发处理上,如果将程序的控制台输出都打印在终端上往往不利于查看,我们可以利用snakemake提供的log函数将他们重定向到相应的文件中,注意snakemake仅会自动创建文件,重定向操作仍需要额外在shell中通过重定向操作符来完成。

conda环境:

snakemake对conda环境有着良好的支持,可以通过下面的形式来添加:

conda:
“envs/map.yaml”

这样就为相应的规则指定了一个局部环境。

 

配置文件支持:

snakemake支持将一些配置剥离成独立的yaml或者json文件,并通过在首行指定configfile: “config.yaml”的形式来引入。

Benchmarking

这一个很有趣的功能,它可以将想要的规则执行的时间写到一个文件中。

其他小tips:

shell执行程序的指定:

shell.executable(“/bin/bash”)

一个包含了上面几乎全部内容究极完全体示例:

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"
    benchmark:
        "benchmarks/{sample}.bwa.benchmark.txt"
    threads: 8
    shell:
        "(bwa mem -R '{params.rg}' -t {threads} {input} | "
        "samtools view -Sb - > {output}) 2> {log}"

推荐阅读:

[1]. https://snakemake.readthedocs.io/en/stable/getting_started

[2]. https://blog.csdn.net/u012110870/article/details/85330457

[3]. https://blog.csdn.net/u012110870/article/details/102804222

 

 

 

 


daniu , 版权所有丨如未注明 , 均为原创丨本网站采用BY-NC-SA协议进行授权
转载请注明原文链接:snakemake踩坑笔记
喜欢 (19)
关于作者:
发表我的评论
取消评论
表情 贴图 加粗 删除线 居中 斜体 签到

Hi,您需要填写昵称和邮箱!

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址