精华内容
下载资源
问答
  • hppRNA包专用于从头到尾同时对大量样本进行RNA-Seq分析,在Snakemake管道管理系统中制定。 它从fastq文件开始,结合最先进的软件,将生成基因/异构体表达矩阵、差异表达基因、样本簇以及SNP和融合基因的检测。 第一...
  • 使用Snakemake和R进行经济学研究的可重现工作流 介绍 我们为使用Snakemake和R编程语言的可重复研究项目提供了模板。 我们使用Snakemake构造一组规则,这些规则形成一个DAG,该DAG实施了整个研究流程,首先进行一些...
  • Snakemake工作流管理系统是创建可重现和可扩展的数据分析的工具。 Snakemake非常受欢迎,平均每周有6次以上新引用,下载量超过20万次。 通过基于Python的人类可读语言来描述工作流。 它们可以无缝扩展到服务器,群集...
  • 使用Snakemake管理工作流有以下好处: 配置文件将所有参数包含在一个文件中,因此您可以查看工作流程在做什么,并为以后的运行进行更改。 每次电气石运行的目录结构都是相同的,因此您始终知道输出的位置。 按需...
  • ta_dna_snakemake_pbs-源码

    2021-03-16 20:46:22
    ta_dna_snakemake_pbs ================================================== ==== Snakemake工作流程铺道芹DNA用PBS-扭矩批作业提交系统的MACH2 HPC群集上分析。 该存储库最初基于和。 如果您要运行自己的自定义...
  • MAG Snakemake工作流程 内容 概述 该管道可用于从短读,与宿主相关的宏基因组学数据集中恢复和评估原核MAG。 使用两个文件指定了分析的数据,该文件详细描述了要分析的联合样品和单次运行。 文件run.txt具有单次运行...
  • 蛇汉堡 这是模块化Snakemake的初学者教程,其灵感来自教程。 warm_bun pickle roasted_onion tomato cheese patty
  • 在Conda环境中安装R,Python和Snakemake 通过Conda使用R,Python和Snakemake进行可重现的分析 该存储库包含安装脚本,用于自动安装R,Python和Snakemake以及在conda环境中用于生物信息学和数据科学项目的其他软件。...
  • shared_snakemake 这是Snakefiles的存储库,其中包含可以插入到各种项目中的常用snakemake规则和工作流。 目录 假设条件 此仓库中的规则是根据以下假设编写的。 如果您的项目工作流程违反了这些假设,则此处的某些...
  • Snakemake中的基本散装RNA-seq管线 目录 描述 该存储库包含两种基本形式的基本批量RNA-seq管道的演示,即Snakemake (在workflow/目录中)和bash脚本(在bash_workflow/目录中)。 这两个工作流程执行相同的分析,...
  • 用于Jupyter Notebook的类似于Snakemake的管道,生成交互式管道报告,如下所示: 安装和一般备注 这些仍是该软件的早期产品,因此请记住,该软件尚未准备好投入生产。 注意:为简单起见,我假设您使用的是安装了git...
  • Snakemake包装器存储库 Snakemake包装器存储库是可重复使用包装器的集合,这些包装器可用于快速使用规则和工作流中流行的命令行工具。 有关更多信息,请访问 。
  • snakemake_workflows 使用snakemake的NGS工作流程
  • snakemake_multiqc 这是一个简单,可重现的示例,其运行fastqc,然后使用测试fastq文件在Snakemake工作流程中运行multiqc。 要运行该示例,请从工作目录中运行以下命令,并确保在本地系统或蛇形特定的conda环境中...
  • UAHPC_Snakemake-源码

    2021-02-17 09:50:55
    UAHPC_Snakemake
  • 上一篇介绍了 Snakemake 入门教程 做了一个简单的示例,具体查看我的上一篇内容 下面会介绍一下 Snakemake的常用参数以及进阶的用法~ 介绍之前大家可以看一个视频了解一下(时长:19min14s, 选择性观看) ...

    写在前面

    上一篇介绍了 Snakemake 入门教程 做了一个简单的示例,具体查看我的上一篇内容
    公众号:生信技术

    下面会介绍一下 Snakemake的常用参数以及进阶的用法

    介绍之前大家可以看一个视频了解一下(时长:19min14s, 选择性观看)

    Snakemake的简单介绍

    参数介绍

    命令行参数

    1. 内核数调用
    $ snakemake --cores 1
    
    # 指定多个可用内核
    $ snakemake --cores 4
    

    Snakemake 执行在同一目录中 名为 Snakefile 的文件中指定的工作流(Snakefile 可以通过参数 -s 给出)。

    1. 试运行
    $ snakemake -n
    
    # 试运行并打印试行内容
    $ snakemake -n -r
    

    可以进行试运行,测试工作流是否正确定义以及估计所需计算量。

    Snakemake 工作流通常定义某些规则的使用线程数。 有时,覆盖工作流定义中给出的默认值可以通过使用--set-threads参数来完成,例如,

    $ snakemake --cores 4 --set-threads myrule=2
    

    将覆盖为rule myrule 定义的任何线程数,并使用2代替。 类似地,可以覆盖规则中的其他资源定义,通过

    $ snakemake --cores 4 --set-resources myrule:partition="foo"
    

    当与集群执行结合使用时,这两种机制都特别方便。

    处理非常大的工作流程

    Snakemake 允许批量运行大型工作流。这样,一次需要评估的文件更少,因此可以更快地推断出作业 DAG

    # 指示仅计算规则myrule 的三批输入中的第一批
    $ snakemake --cores 4 --batch myrule=1/3
    
    # 生成第二批
    $ snakemake --cores 4 --batch myrule=2/3
    
    # 生成第三批
    $ snakemake --cores 4 --batch myrule=3/3
    
    一些参数介绍具体内容
    target要建立的目标,可能是rule或文件
    –dry-run, --dryrun, -n不执行任何操作,并显示将要执行的操作。如果是一个非常大的工作流程,使用 –dry-run –quiet 仅打印作业 DAG 的摘要
    –profileSnakemake 的配置文件的名称
    –snakefile, -s指定Snakefile,否则是当前目录下的Snakefile
    –cores, -c , -j并行使用 N 个 CPU cores/jobs
    –resources, --res指定程序运行的内存
    –config设置或覆盖工作流配置对象中的值
    –directory, -d指定工作目录
    –touch, -t将某文件标记为最新,不真正更改它们
    –force强制执行某一条
    –forceall强制执行某条Rule及它的依赖
    –report创建包含结果和统计信息的 HTML 报告
    –list, -l显示给定 Snakefile 中的可用rule
    –dag生成整个流程的有向无环图,不进行分析
    –quiet, -q不输出任何进度或规则信息
    –all-temp将所有输出文件标记为临时文件
    –bash-completion文件名、规则名和参数的 bash 补全

    进阶用法

    1. 指定使用的线程数

    对于某些工具,建议使用多个线程以加快计算速度。Snakemake可以通过threads指令了解规则所需的“threads”。在示例工作流中,为规则bwa_map使用多个线程:

    $ snakemake --cores 4 --set-threads myrule=2
    
    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})将线程数传到 shell 命令。如果没有threads指令,则默认使用 1 个线程。

    执行工作流时,Snakemake 调度程序会考虑作业所需的线程数。 调度程序会确保同时运行的所有作业的线程总和不超过给定的可用 CPU 内核数。

    这个数字是通过 --cores 命令行参数给出的,这对于实际运行工作流的 snakemake 调用是必需的。 例如

    $ snakemake --cores 10
    

    将使用 10 个内核执行工作流。 由于rule bwa_map 需要 8 个线程,因此一次只能运行该rule的一项作业,

    Snakemake 调度程序将尝试用其他作业(例如 samtools_sort )来使剩余的内核饱和。

    当提供的内核数少于线程数时,规则使用的线程数将减少到给定的内核数。

    如果--cores 没有给出数字,则使用所有可用的核心。

    使用标志--forceall,可以强制完全重新执行工作流。 将此标志与 --cores 的不同值结合起来,并检查调度程序如何选择要并行运行的作业。

    除了非常常见的线程资源外,Snakemake 还提供了一个resources指令,可用于指定任意资源,例如内存使用或辅助计算设备(如 GPU)。与线程类似,当使用命令行参数提供该资源的可用量时,调度程序可以考虑 --resources( 参考资料

    2. 配置文件

    Snakemake 提供了一个配置文件机制。配置文件可以用JSONYAML编写,并与configfile指令一起使用。在我们的示例工作流程中,我们添加了以下行

    configfile: "config.yaml"
    

    写到 Snakefile 的顶部

    Snakemake 将加载配置文件并将其内容存储到名为 config 的全局可用字典中。 在我们的例子中,将 config.yaml 中的示例指定为

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

    现在,我们可以从 Snakefile 中删除定义 SAMPLES 的语句,并将规则 bcftools_call 更改为

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

    3. 输入函数

    由于我们已经在配置文件中存储了 FASTQ 文件的路径,我们也可以在 rule bwa_map 来使用这些路径。这种情况与我们上面修改的规则“bcftools_call”不同。

    要理解这一点,重要的是要知道 Snakemake 工作流分三个阶段执行。

    1. 初始化阶段,定义工作流的文件被解析并实例化所有规则。

    2. DAG 阶段,通过填充通配符并将输入文件与输出文件匹配来构建所有作业的有向无环依赖图。

    3. 调度阶段,执行作业的 DAG,根据可用资源启动作业。

    rule bcftools_call 的输入文件列表中的扩展函数在初始化阶段执行。

    在这个阶段,我们不知道作业、通配符值和规则依赖关系。因此,在此阶段,我们无法从配置文件中确定rule bwa_map 的 FASTQ 路径,因为我们甚至不知道将从该rule生成哪些作业。相反,我们需要将输入文件的确定推迟到 DAG 阶段。这可以通过指定输入函数而不是输入指令内部的字符串来实现。对于规则bwa_map,其工作方式如下:

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

    任何正常的功能都可以正常工作。 输入函数将 wildcards 对象作为单个参数,该对象允许通过属性(此处为 wildcards.sample )访问通配符值。 它们必须返回一个字符串或一个字符串列表,这些字符串被解释为输入文件的路径(这里,我们返回在配置文件中为 示例存储的路径 )。 一旦确定了作业的通配符值,就会评估输入函数。

    当添加新的输入文件时,Snakemake 不会自动重新运行作业,如下所示。 但是,可以使用 snakemake --list-input-changes 获取受此类更改影响的输出文件列表。 要触发重新运行,这一点 bash 有帮助:

    示例

    data/samples 文件夹中,还有一个文件 C.fastq

    将该文件添加到配置文件中,并查看当使用 snakemake -n --reason --forcerun bcftools_call 调用时,Snakemake 如何重新计算属于新示例的工作流部分。

    4. rule 参数

    有时,shell 命令不仅由inputoutput 文件以及一些静态标志组成。

    特别是,可能需要根据作业的通配符值设置其他参数。 为此,Snakemake 允许使用 params 指令为 rule 定义任意参数。

    有时,shell 命令不仅由输入和输出文件以及一些静态标志组成。 特别是,可能需要根据作业的通配符值设置其他参数。 为此,Snakemake 允许使用 params 指令为规则定义任意参数。

    在我们的工作流程中,使用所谓的read groups注释对齐的reads是合理的,其中包含像sample名称这样的元数据。 我们相应地修改规则bwa_map

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

    inputoutput 类似,可以从 shell 命令、基于 Python 的 run 块或脚本指令访问 params

    params 指令也可以使用第 3 步中的函数将初始化推迟到 DAG 阶段。 与 input 函数相比,这些函数可以选择采用额外的参数 inputoutputthreadsresources

    5. log文件记录

    在执行大型工作流时,通常希望将每个作业的日志输出存储到单独的文件中,而不是仅仅将所有日志输出打印到终端(当多个作业并行运行时,这会导致输出混乱)

    为此,Snakemake 允许为rule指定日志文件。 日志文件是通过 log 指令定义的,处理方式与输出文件类似,但它们不受rule匹配的影响,并且在作业失败时不会被清除。 修改规则 bwa_map 如下:

    rule bwa_map:
        input:
            "data/genome.fa",
            get_bwa_map_input_fastqs
        output:
            "mapped_reads/{sample}.bam"
        params:
            rg=r"@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}"
    

    以上内容通过修改 shell 命令以收集 bwasamtools 的 STDERR 输出,并将其通过管道传输到 {log} 引用的文件中。 日志文件必须包含与输出文件完全相同的通配符,以避免相同规则的不同作业之间的文件名冲突

    最佳做法是将所有日志文件存储在子目录logs/中,并以 rule 或工具名称为前缀。

    6. 临时文件和受保护文件

    在工作流程中,为每个样本创建两个 BAM 文件,即rule bwa_mapsamtools_sort 的输出。 在不处理示例时,底层数据通常是巨大的。 因此,生成的 BAM 文件需要大量磁盘空间,并且它们的创建需要一些时间。

    为了节省磁盘空间,您可以将输出文件标记为临时文件。 一旦所有消耗作业(需要它作为输入)都已执行,Snakemake 将为您删除标记的文件。 我们将这种机制用于规则 bwa_map 的输出文件:

    rule bwa_map:
        input:
            "data/genome.fa",
            get_bwa_map_input_fastqs
        output:
            temp("mapped_reads/{sample}.bam")
        params:
            rg=r"@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 作业,将导致BAM 文件被删除 。 由于通过读取映射和排序创建 BAM 文件的计算成本很高,因此保护最终的 BAM 文件免遭意外删除或修改是合理的。

    修改规则 samtools_sort 以将其输出文件标记为 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 会对文件系统中的输出文件进行写保护,以免被意外覆盖或删除。
    

    重新执行整个工作流程并观察 Snakemake 如何处理临时文件和受保护文件。
    使用目标运行 Snakemake mapped_reads/A.bam。虽然该文件被标记为临时文件,但您会看到 Snakemake 并没有删除它,因为它被指定为目标文件。
    尝试使用试运行选项再次重新执行整个工作流程。您将看到它失败,因为 Snakemake 无法覆盖受保护的输出文件。

    总结

    对于本教程,我们学习创建了一个config.yaml配置文件:

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

    有了这个,我们工作流程的最终版本Snakefile看起来像这样:

    configfile: "config.yaml"
    
    rule all:
        input:
            "plots/quals.svg"
    
    def get_bwa_map_input_fastqs(wildcards):
        return config["samples"][wildcards.sample]
    
    rule bwa_map:
        input:
            "data/genome.fa",
            get_bwa_map_input_fastqs
        output:
            temp("mapped_reads/{sample}.bam")
        params:
            rg=r"@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"
        params:
            rate=config["prior_mutation_rate"]
        log:
            "logs/bcftools_call/all.log"
        shell:
            "(samtools mpileup -g -f {input.fa} {input.bam} | "
            "bcftools call -mv -P {params.rate} - > {output}) 2> {log}"
    
    rule plot_quals:
        input:
            "calls/all.vcf"
        output:
            "plots/quals.svg"
        script:
            "scripts/plot-quals.py"
    
    展开全文
  • fastq_to_bam_paired_snakemake 遵循GATK最佳做法,将成对的fastq文件转换为可分析的bams的snakemake。 该存储库包含运行snakemake所需的所有文件夹。 由于GitHub不允许空文件夹,所以results /和logs / cluster /...
  • 写在前面 既然写了教程就需要具有普适...Snakemake 定义 Snakemake 工作流管理系统是一种用于创建可重复和可扩展的数据分析的工具。 工作流是通过一种人类可读的、基于 Python 的语言来描述的。它们可以无缝扩展到服务

    写在前面

    既然写了教程就需要具有普适性,能适合大多数人的胃口,我这部分的内容以及示例主要还是参考了官方教程,但是都是我一步一步跑过的流程,所以会更有印象,送给想学 Snkaemake 但是一直没有去学的朋友们,这些内容对于有生信基础的人来讲,上手会很快,因为很多的生信软件都使用过,写起来也就没有那么晦涩,下面开始~

    Snakemake 定义

    Snakemake 工作流管理系统是一种用于创建可重复和可扩展的数据分析的工具。

    工作流是通过一种人类可读的、基于 Python 的语言来描述的。它们可以无缝扩展到服务器、集群、grid和云环境,无需修改工作流。

    最后,Snakemake 工作流可能需要对所需软件的准备,这些软件将自动部署到任何执行环境。

    安装

    Snakemake 可在 PyPi 上以及通过 Bioconda 和源代码获得。可以使用任意方法安装 Snakemake,我们这里仅介绍使用 Conda 安装

    通过 Conda/Mamba 安装

    这是安装 Snakemake的推荐方式,因为Conda安装比较简单。

    首先,必须已经安装了一个基于 Conda 的 Python3 发行版。推荐的选择是Mambaforge,它不仅提供所需的 Python 和 Conda 命令,而且包括 Mamba,它是强烈推荐的Conda包管理器的极其快速和强大的替代品。默认的 conda 求解器有点慢,有时在选择最新的软件包版本时会出现问题。因此,建议在任何情况下都使用Mamba。

    如果不安装 Mambaforge ,也可以直接安装 Mamba

    conda install -n base -c conda-forge mamba
    

    单独安装到一个小环境

    $ conda activate base
    $ mamba create -c conda-forge -c bioconda -n snakemake snakemake
    

    安装到单独环境中比较好,为了避免与其他软件包冲突,并通过如下方式进行激活

    $ conda activate snakemake
    $ snakemake --help
    

    仅安装必备软件的 snakemake 版本

    可以安装仅依赖基本必需软件的 minimal 版本 Snakemake

    $ conda activate base
    $ mamba create -c bioconda -c conda-forge -n snakemake snakemake-minimal
    

    基础示例

    首先先激活 Snakemake ,看各自下载的环境,我是单独创建了一个小环境 所以我进行如下操作:

    $ conda activate snakemake
    

    Snakemake 工作流是通过在 Snakefile 中指定命令来定义的命令通过指定如何从输入文件集创建输出文件集将工作流分解为小步骤(例如,单个工具的应用)。Snakemake通过匹配文件名自动确定命令之间的依赖关系

    下面通过创建示例工作流来介绍 Snakemake 语法。 工作流程来自基因组分析领域。 它将测序reads映射到参考基因组,并在映射reads上调用变体。 本教程不需要您知道这是关于什么的。 尽管如此,我们在以下段落中提供了一些背景知识。

    测试数据下载

    git clone https://bitbucket.org/snakemake/snakemake-tutorial.git
    
    cd snakemake-tutorial
    
    ├── 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
    
    # 开启 Snakemake 学习之旅
    

    1. Mapping reads

    第一个 Snakemake 命令将给定样本的reads映射到给定的参考基因组。为此,我们将使用工具bwa,特别是 subcommand 。在工作目录中,创建一个使用您选择的编辑器调用的新文件。官方建议使用Atom编辑器,因为它为 Snakemake 提供了开箱即用的语法突出显示。在 Snakefile 中,定义以下命令:bwa mem Snakefile

    第一个 Snakemake 命令将给定样本的reads映射到给定的参考基因组。 为此,我们将使用工具 bwa,特别是子命令 bwa mem。 在工作目录中,使用编辑器创建一个名为 Snakefile 的新文件。 官方建议使用 Atom 编辑器,因为它为 Snakemake 提供了开箱即用的语法突出显示。 在 Snakefile 中,定义以下命令:

    首先创建一个名称为 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}"
    
    # 一个常见的错误是忘记输入或输出项之间的逗号。 由于 Python 连接后续字符串,这可能会导致抱错
    

    Snakemake rule有一个名称(这里是 bwa_map)和许多指令,这里是 inputoutputshell

    inputoutput 指令之后是那些预计将在命令中使用或创建的文件列表。

    在最简单的情况下,这些只是 Python 字符串。 shell 指令后跟一个包含要执行的 shell 命令的 Python 字符串。

    在 shell 命令字符串中,我们可以通过**大括号表示法(类似于 Python 格式函数)**引用命令的元素。

    在这里,我们通过指定 {output} 来引用输出文件,通过指定 {input} 来引用输入文件。由于命令有多个输入文件,Snakemake 将连接它们,用空格分隔。换句话说,Snakemake 会在执行命令之前将 {input} 替换为 data/genome.fa data/samples/A.fastq

    shell 命令使用参考基因​​组和reads调用 bwa mem,并将输出通过管道传输到 samtools,后者创建包含比对的压缩 BAM 文件。 samtools 的输出被重定向到命令定义的输出文件中,并带有 >

    执行工作流时,Snakemake 会尝试生成给定的目标文件。可以通过命令行指定目标文件。通过执行

    # 试运行
    $ snakemake -np mapped_reads/A.bam
    

    在包含 Snakefile 的工作目录中, Snakemake 生成目标文件 mapping_reads/A.bam。 由于我们使用了 -n(或 --dry-run)标志,Snakemake 将只显示执行计划而不是实际执行步骤-p 标志指示 Snakemake 打印出生成的 shell 命令以供说明。

    为了生成目标文件,Snakemake 以自上而下的方式应用 Snakefile 中给出的命令。 应用命令来生成一组输出文件称为作业。 对于作业的每个输入文件,Snakemake 再次(即递归地)确定可应用于生成它的命令。 这产生了作业的有向无环图 (DAG),其中边代表依赖关系。 到目前为止,我们只有一个命令,作业的 DAG 由单个节点组成。 尽管如此,我们可以执行我们的工作流程

    $ snakemake --cores 1 mapped_reads/A.bam
    

    无论何时执行工作流,都需要指定要使用的核心数。 对于本教程,现在将使用单个内核。 稍后介绍并行化是如何工作的。 完成上述命令后,Snakemake 将不会再次尝试创建mapped_reads/A.bam,因为它已经存在于文件系统中。 Snakemake 仅在输入文件之一比输出文件之一新 或 输入文件之一将被另一个作业更新时重新运行作业

    2. Generalizing the read mapping rule

    前面的rule仅适用于在文件 data/samples/A.fastq 中读取。 但是,Snakemake 允许使用命名通配符。 只需用通配符 {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}"
    

    当 Snakemake 通过用适当的值替换输出文件中的通配符 {sample} 确定可以应用此命令来生成目标文件时,它将该值传播到输入文件中所有出现的 {sample},从而确定 结果工作的必要输入。

    注意,您的文件路径中可以有多个通配符,但是,为了避免与同一命令的其他作业发生冲突,命令的所有输出文件必须包含完全相同的通配符

    $ snakemake -np mapped_reads/B.bam
    
    # 运行之后输出
    rule bwa_map:
        input: data/genome.fa, data/samples/B.fastq
        output: mapped_reads/B.bam
        jobid: 0
        wildcards: sample=B
        resources: tmpdir=/tmp
    # 可以看到内容随着输入命令变化匹配到了B.bam
    

    Snakemake 将通过将通配符 {sample} 替换为值 B 来确定可以应用命令 bwa_map 来生成目标文件。在试运行的输出中,可以看到通配符值如何传播到输入文件和 shell 命令中的所有文件名。 还可以指定多个目标,例如:

    $ snakemake -np mapped_reads/A.bam mapped_reads/B.bam
    

    一些Bash语法特别方便。例如,可以选择在一次组合多个目标

    $ snakemake -np mapped_reads/{A,B}.bam
    # Bash 只是将其大括号扩展应用于集合 {A,B},为每个元素创建给定的路径并用空格分隔结果路径。
    
    # snakemake -np mapped_reads/{1..10}.bam
    # 会匹配1.bam; 2.bam; ... ; 10.bam 
    

    在这两种情况下, Snakemake 只创建输出文件 mapping_reads/B.bam

    这是因为之前已经执行过 mapping_reads/A.bam 并且没有比输出文件更新的输入文件。

    可以更新输入文件data/samples/A.fastq的文件修改日期

    $ touch data/samples/A.fastq
    

    并运行 Snakemake 重新运行作业来创建文件 mapping_reads/A.bam

    $ snakemake -np mapped_reads/A.bam mapped_reads/B.bam
    

    3. Sorting read alignments

    对于后面的步骤,我们需要对 BAM 文件中的读取对齐进行排序。这可以通过 samtools sort 命令实现。我们在 bwa_map rule下添加以下rule:

    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}"
    
    
    # 在上面的 shell 命令中,我们将字符串分成两行,但是 Python 会自动将它们连接成一行。
    # 分行写的话可以避免 shell 命令行过长。使用它时,需要在每行但最后一行中都有一个尾随空格,以避免参数无法正确分隔。
    

    此命令将从mapped_reads 文件夹中获取输入文件,并将排序后的版本存储在sorted_reads 目录中。

    注意,Snakemake 会在作业执行前自动创建丢失的目录。 对于排序,samtools 需要使用标志 -T 指定的前缀。

    在这里,我们需要通配符sample的值。 Snakemake 允许通过 wildcards 对象访问 shell 命令中的通配符,该对象具有带有每个通配符值的属性。

    wildcards指通配符,学过类 LINUX 系统的,应该都知道什么是通配符。
    * 代表任意多个字符
    ? 代表任意单个字符
    [ ] 代表“[”和“]”之间的某一个字符,比如[0-9]可以代表0-9之间的任意一个数字,[a-zA-Z]可以代表a-z和A-Z之间的任意一个字母,字母区分大小写
    代表一个字符
    ~ 用户的根目录

    $ snakemake -np sorted_reads/B.bam
    

    看到 Snakemake 首先运行bwa_map,然后运行samtools_sort来创建所需的目标文件:如前所述,依赖项通过匹配文件名自动解析。

    4. Indexing read alignments and visualizing the DAG of jobs

    接下来,我们需要再次使用 samtools 来索引已排序的读取比对,以便我们可以通过它们映射到的基因组位置快速访问读取。这可以通过以下命令来完成:

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

    Snakemake 使用Python 格式mini language来格式化 shell 命令。在 shell 命令中使用大括号 ({}) 来表示其他内容。在这种情况下,必须加倍对我们上面提到的是bash括号扩展依靠时逃避它们,
    例如:
    ls {{A,B}}.txt

    已经完成了三个步骤,现在可以查看生成的有向无环图 (DAG)

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

    Snakemake 使用 Graphviz 提供的 dot 命令创建 DAG 的可视化。 对于给定的目标文件,Snakemake 以 dot 语言指定 DAG 并将其通过管道传输到 dot 命令,该命令将定义呈现为 SVG 格式。 渲染的 DAG 通过管道传输到文件 dag.svg 中,如下所示:

    在这里插入图片描述

    5. Calling genomic variants

    我们工作流程的下一步将聚合所有样本的映射reads,并共同调用它们的基因组变异。 对于变体调用,我们将结合两个实用程序 samtoolsbcftools。 Snakemake 提供了一个辅助函数来收集输入文件,帮助我们描述这一步中的聚合。

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

    获取文件列表,其中给定模式sorted_reads/{sample}.bam被格式化为给定样本列表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 文件。

    Snakefiles 原则上是 Python 代码,通过一些声明性语句来定义工作流。 因此,我们可以在 Snakefile 顶部的纯 Python 中临时定义示例列表:

    SAMPLES = ["A", "B"]
    

    可以将以上命令添加到 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 命令中分别引用它们会很方便。 这可以通过指定输入或输出文件的名称来完成,例如使用fa=....然后可以在shell命令中通过名称引用这些文件,例如使用{input.fa}

    对于像这样的长 shell 命令,建议将字符串拆分为多个缩进行。 Python 会自动将其合二为一。 此外,您会注意到输入或输出文件列表可以包含任意 Python 语句,只要它返回一个字符串或字符串列表。 在这里,我们调用我们的 expand 函数来聚合所有样本的对齐reads。

    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"
    

    使用此规则,我们最终将生成已分配给文件 calls/all.vcf 中的variant calls的质量分数的直方图。 生成绘图的实际 Python 代码在脚本 scripts/plot-quals.py 中。 在脚本中,命令的所有属性,如 inputoutputwildcards 等,都可以作为全局 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脚本;有关详细信息和示例,可以阅读官方教程中的外部脚本部分

    7. Adding a target rule

    到目前为止,我们总是通过在命令行指定目标文件来执行工作流。除了文件名,如果请求的规则没有通配符,Snakemake还接受规则名作为目标。

    因此,可以编写目标规则来收集所需结果或所有结果的特定子集。此外,如果命令行中没有给出目标,Snakemake会将 Snakefile 的第一条规则定义为目标。因此,最好的做法是在工作流的顶部有一个规则all,该rule将所有通常需要的目标文件作为输入文件

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

    把这个rule添加到我们工作流程的顶部。执行 Snakemake 时

    $ snakemake -n
    
    # 可以在 Snakefile 的顶部添加多个目标rules。虽然 Snakemake 将默认执行第一个,但可以通过命令行(例如,snakemake -n mytarget)定位其中的任何一个。
    

    执行此命令将显示用于创建文件的执行计划plots/quals.svg,其中包含并总结了我们所有的结果。

    除了 Snakemake 将工作流的第一条规则视为默认目标之外,Snakefile 中的规则顺序是任意的,不会影响作业的 DAG。

    总结

    生成的工作流程如下所示:

    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"
    
    
    展开全文
  • Moeller Lab宏基因组学处理流程 Snakemake管道用于基本处理实验室中的宏基因组数据。
  • 使用Snakemake搭建分析流程

    万次阅读 多人点赞 2018-12-28 19:28:43
    snakemake-snakemake-tutorial-623791d7ec6d conda env create --name snakemake-tutorial --file environment.yaml source activate snakemake-tutorial # 退出当前环境 source deactivate 当前环境下...

    ## 目前已有的框架

    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: 集群配置文件
    

    参考资料

    展开全文
  • Snakemake管道可处理基于板的scATAC-seq数据 该存储库包含用于处理通过生成的scATAC-seq数据的代码。 与原始Nat相比有什么区别。 通讯出版物? 在此管道中,我们具有: 添加了config.json文件,以使处理更加灵活且...
  • snakemake 学习笔记3

    千次阅读 2019-04-02 21:25:40
    Complete log: /home/dengfei/test/snakemake/ex4/.snakemake/log/2019-04-02T211109.153566.snakemake.log 查看*add_a.txt文件: (snake_test) [dengfei@localhost ex4]$ ls *add_a.txt 1_add_a.txt 2_add_a....

    目标

    这次, 我要实现这个路程图.

    在这里插入图片描述

    目标介绍

    • 第一: 生成1.txt , 2.txt, 3.txt
    • 第二: 向每个文件中加入"add a"字符, 命名为:1_add_a.txt, 2_add_a.txt, 3_add_a.txt
    • 第三: 向文件中增加"add b", 命名为:1_add_a_add_b.txt, 2_add_a_add_b.txt, 3_add_a_add_b.txt
    • 第四: 向文件中增加"add c", 命名为: 1_add_a_add_b_add_c.txt, 2_add_a_add_b_add_c.txt, 3_add_a_add_b_add_c.txt
    • 第五: 将1_add_a_add_b.txt, 2_add_a_add_b.txt, 1_add_a_add_b_add_c.txt, 2_add_a_add_b_add_c.txt, 3_add_a_add_b_add_c.txt 合并为hebing.txt文件

    1. 生成三个文件

    (snake_test) [dengfei@localhost ex4]$ ls *txt
    1.txt  2.txt  3.txt
    (snake_test) [dengfei@localhost ex4]$ cat *txt
    this is 1.txt
    this is 2.txt
    this is 3.txt
    
    

    2. 在每个文件中增加"add a"

    对应的Snakefile内容如下:

    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    
    

    预览一下命令:snakemake -np {1,2,3}_add_a.txt

    注意: 这里要把生成的文件{1,2,3}_add_a.txt写出来, 命令才可以运行.

    (snake_test) [dengfei@localhost ex4]$ snakemake -np {1,2,3}_add_a.txt
    Building DAG of jobs...
    Job counts:
    	count	jobs
    	3	adda
    	3
    
    [Tue Apr  2 21:09:19 2019]
    rule adda:
        input: 3.txt
        output: 3_add_a.txt
        jobid: 2
        wildcards: file=3
    
    cat 3.txt |xargs echo add a >3_add_a.txt
    
    [Tue Apr  2 21:09:19 2019]
    rule adda:
        input: 2.txt
        output: 2_add_a.txt
        jobid: 0
        wildcards: file=2
    
    cat 2.txt |xargs echo add a >2_add_a.txt
    
    [Tue Apr  2 21:09:19 2019]
    rule adda:
        input: 1.txt
        output: 1_add_a.txt
        jobid: 1
        wildcards: file=1
    
    cat 1.txt |xargs echo add a >1_add_a.txt
    Job counts:
    	count	jobs
    	3	adda
    	3
    This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.
    
    

    执行命令:

    snakemake  {1,2,3}_add_a.txt
    
    
    Building DAG of jobs...
    Using shell: /usr/bin/bash
    Provided cores: 1
    Rules claiming more threads will be scaled down.
    Job counts:
    	count	jobs
    	3	adda
    	3
    
    [Tue Apr  2 21:11:09 2019]
    rule adda:
        input: 3.txt
        output: 3_add_a.txt
        jobid: 0
        wildcards: file=3
    
    [Tue Apr  2 21:11:09 2019]
    Finished job 0.
    1 of 3 steps (33%) done
    
    [Tue Apr  2 21:11:09 2019]
    rule adda:
        input: 1.txt
        output: 1_add_a.txt
        jobid: 1
        wildcards: file=1
    
    [Tue Apr  2 21:11:09 2019]
    Finished job 1.
    2 of 3 steps (67%) done
    
    [Tue Apr  2 21:11:09 2019]
    rule adda:
        input: 2.txt
        output: 2_add_a.txt
        jobid: 2
        wildcards: file=2
    
    [Tue Apr  2 21:11:09 2019]
    Finished job 2.
    3 of 3 steps (100%) done
    Complete log: /home/dengfei/test/snakemake/ex4/.snakemake/log/2019-04-02T211109.153566.snakemake.log
    
    

    查看*add_a.txt文件:

    (snake_test) [dengfei@localhost ex4]$ ls *add_a.txt
    1_add_a.txt  2_add_a.txt  3_add_a.txt
    (snake_test) [dengfei@localhost ex4]$ cat *add_a.txt
    add a this is 1.txt
    add a this is 2.txt
    add a this is 3.txt
    
    

    搞定.

    3. 在每个文件中增加"add b"

    对应的Snakefile内容如下:

    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    rule addb:
        input:
            "{file}_add_a.txt"
        output:
            "{file}_add_a_add_b.txt"
        shell:
            "cat {input} | xargs echo add b >{output}"
    
    
    

    预览一下命令:snakemake -np {1,2,3}_add_a_add_b.txt

    (snake_test) [dengfei@localhost ex4]$ snakemake  {1,2,3}_add_a_add_b.txt
    Building DAG of jobs...
    Using shell: /usr/bin/bash
    Provided cores: 1
    Rules claiming more threads will be scaled down.
    Job counts:
    	count	jobs
    	3	addb
    	3
    
    [Tue Apr  2 21:13:57 2019]
    rule addb:
        input: 2_add_a.txt
        output: 2_add_a_add_b.txt
        jobid: 0
        wildcards: file=2
    
    [Tue Apr  2 21:13:57 2019]
    Finished job 0.
    1 of 3 steps (33%) done
    
    [Tue Apr  2 21:13:57 2019]
    rule addb:
        input: 1_add_a.txt
        output: 1_add_a_add_b.txt
        jobid: 1
        wildcards: file=1
    
    [Tue Apr  2 21:13:57 2019]
    Finished job 1.
    2 of 3 steps (67%) done
    
    [Tue Apr  2 21:13:57 2019]
    rule addb:
        input: 3_add_a.txt
        output: 3_add_a_add_b.txt
        jobid: 2
        wildcards: file=3
    
    [Tue Apr  2 21:13:57 2019]
    Finished job 2.
    3 of 3 steps (100%) done
    Complete log: /home/dengfei/test/snakemake/ex4/.snakemake/log/2019-04-02T211357.666661.snakemake.log
    
    
    

    执行命令:

    snakemake  {1,2,3}_add_a_add_b.txt
    
    

    查看流程图

    命令:

    snakemake --dag {1,2,3}_add_a_add_b.txt |dot -Tpdf >a.pdf
    
    

    这里生成的a.pdf如下:

    4. 在每个文件中增加"add c"

    Snakemake命令:

    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    rule addb:
        input:
            "{file}_add_a.txt"
        output:
            "{file}_add_a_add_b.txt"
        shell:
            "cat {input} | xargs echo add b >{output}"
    
    rule addc:
        input:
            "{file}_add_a_add_b.txt"
        output:
            "{file}_add_a_add_b_add_c.txt"
        shell:
            "cat {input} | xargs echo add c >{output}"
    
    

    流程图:

    命令:

    snakemake --dag {1,2,3}_add_a_add_b_add_c.txt |dot -Tpdf >a1.pdf
    
    

    在这里插入图片描述

    5. 将文件合并

    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    rule addb:
        input:
            "{file}_add_a.txt"
        output:
            "{file}_add_a_add_b.txt"
        shell:
            "cat {input} | xargs echo add b >{output}"
    
    rule addc:
        input:
            "{file}_add_a_add_b.txt"
        output:
            "{file}_add_a_add_b_add_c.txt"
        shell:
            "cat {input} | xargs echo add c >{output}"
    
    rule hebing:
        input:
           a=expand("{file}_add_a_add_b_add_c.txt",file=["1","2","3"]),
           b=expand("{file}_add_a_add_b.txt",file=["1","2"])
        output:"hebing.txt"
        shell:"cat {input.a} {input.b} >{output}"
    
    

    执行命令:

    snakemake hebing.txt
    

    执行结果:

    Building DAG of jobs...
    Using shell: /usr/bin/bash
    Provided cores: 1
    Rules claiming more threads will be scaled down.
    Job counts:
    	count	jobs
    	3	addc
    	1	hebing
    	4
    
    [Tue Apr  2 21:21:04 2019]
    rule addc:
        input: 1_add_a_add_b.txt
        output: 1_add_a_add_b_add_c.txt
        jobid: 1
        wildcards: file=1
    
    [Tue Apr  2 21:21:04 2019]
    Finished job 1.
    1 of 4 steps (25%) done
    
    [Tue Apr  2 21:21:04 2019]
    rule addc:
        input: 3_add_a_add_b.txt
        output: 3_add_a_add_b_add_c.txt
        jobid: 3
        wildcards: file=3
    
    [Tue Apr  2 21:21:04 2019]
    Finished job 3.
    2 of 4 steps (50%) done
    
    [Tue Apr  2 21:21:04 2019]
    rule addc:
        input: 2_add_a_add_b.txt
        output: 2_add_a_add_b_add_c.txt
        jobid: 2
        wildcards: file=2
    
    [Tue Apr  2 21:21:04 2019]
    Finished job 2.
    3 of 4 steps (75%) done
    
    [Tue Apr  2 21:21:04 2019]
    rule hebing:
        input: 1_add_a_add_b_add_c.txt, 2_add_a_add_b_add_c.txt, 3_add_a_add_b_add_c.txt, 1_add_a_add_b.txt, 2_add_a_add_b.txt
        output: hebing.txt
        jobid: 0
    
    [Tue Apr  2 21:21:04 2019]
    Finished job 0.
    4 of 4 steps (100%) done
    Complete log: /home/dengfei/test/snakemake/ex4/.snakemake/log/2019-04-02T212104.719887.snakemake.log
    
    

    流程图:
    在这里插入图片描述

    搞定

    欢迎关注我的公众号: R-breeding
    在这里插入图片描述

    相关阅读

    snakemake 学习笔记1 - CSDN博客
    snakemake-学习笔记2 - CSDN博客

    后记1

    今天测试了一下rule all的功能, 它是定义输出文件的, 如果没有定义, 需要在命令行中书写.

    因为最后的输出文件是hebing.txt, 所以我们这里在Snakefile中定义一下输出文件.

    rule all:
        input:"hebing.txt"
    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    rule addb:
        input:
            "{file}_add_a.txt"
        output:
            "{file}_add_a_add_b.txt"
        shell:
            "cat {input} | xargs echo add b >{output}"
    
    rule addc:
        input:
            "{file}_add_a_add_b.txt"
        output:
            "{file}_add_a_add_b_add_c.txt"
        shell:
            "cat {input} | xargs echo add c >{output}"
    
    rule hebing:
        input:
           a=expand("{file}_add_a_add_b_add_c.txt",file=["1","2","3"]),
           b=expand("{file}_add_a_add_b.txt",file=["1","2"])
        output:"hebing.txt"
        shell:"cat {input.a} {input.b} >{output}"
    
    

    执行命令:

    snakemake
    
    

    结果如下:

    (base) [dengfei@localhost ex4]$ snakemake
    Provided cores: 1
    Rules claiming more threads will be scaled down.
    Job counts:
    	count	jobs
    	3	adda
    	3	addb
    	3	addc
    	1	all
    	1	hebing
    	11
    
    rule adda:
        input: 1.txt
        output: 1_add_a.txt
        jobid: 7
        wildcards: file=1
    
    Finished job 7.
    1 of 11 steps (9%) done
    
    rule adda:
        input: 2.txt
        output: 2_add_a.txt
        jobid: 9
        wildcards: file=2
    
    Finished job 9.
    2 of 11 steps (18%) done
    
    rule adda:
        input: 3.txt
        output: 3_add_a.txt
        jobid: 10
        wildcards: file=3
    
    Finished job 10.
    3 of 11 steps (27%) done
    
    rule addb:
        input: 3_add_a.txt
        output: 3_add_a_add_b.txt
        jobid: 8
        wildcards: file=3
    
    Finished job 8.
    4 of 11 steps (36%) done
    
    rule addb:
        input: 1_add_a.txt
        output: 1_add_a_add_b.txt
        jobid: 3
        wildcards: file=1
    
    Finished job 3.
    5 of 11 steps (45%) done
    
    rule addb:
        input: 2_add_a.txt
        output: 2_add_a_add_b.txt
        jobid: 6
        wildcards: file=2
    
    Finished job 6.
    6 of 11 steps (55%) done
    
    rule addc:
        input: 3_add_a_add_b.txt
        output: 3_add_a_add_b_add_c.txt
        jobid: 5
        wildcards: file=3
    
    Finished job 5.
    7 of 11 steps (64%) done
    
    rule addc:
        input: 2_add_a_add_b.txt
        output: 2_add_a_add_b_add_c.txt
        jobid: 2
        wildcards: file=2
    
    Finished job 2.
    8 of 11 steps (73%) done
    
    rule addc:
        input: 1_add_a_add_b.txt
        output: 1_add_a_add_b_add_c.txt
        jobid: 4
        wildcards: file=1
    
    Finished job 4.
    9 of 11 steps (82%) done
    
    rule hebing:
        input: 1_add_a_add_b_add_c.txt, 2_add_a_add_b_add_c.txt, 3_add_a_add_b_add_c.txt, 1_add_a_add_b.txt, 2_add_a_add_b.txt
        output: hebing.txt
        jobid: 1
    
    Finished job 1.
    10 of 11 steps (91%) done
    
    localrule all:
        input: hebing.txt
        jobid: 0
    
    Finished job 0.
    11 of 11 steps (100%) done
    
    

    查看结果:

    (base) [dengfei@localhost ex4]$ cat hebing.txt 
    add c add b add a this is 1.txt
    add c add b add a this is 2.txt
    add c add b add a this is 3.txt
    add b add a this is 1.txt
    add b add a this is 2.txt
    
    

    后记2

    snakemake如果是默认的名称, 为Snakefile, 但是这样写没有高亮, 可以写为a.py, 然后用snakemake -s a.py运行即可.

    
    rule all:
        input:"hebing.txt"
    rule adda:
        input: "{file}.txt"
        output: "{file}_add_a.txt"
        shell: "cat {input} |xargs echo add a >{output}"
    rule addb:
        input:
            "{file}_add_a.txt"
        output:
            "{file}_add_a_add_b.txt"
        shell:
            "cat {input} | xargs echo add b >{output}"
    
    rule addc:
        input:
            "{file}_add_a_add_b.txt"
        output:
            "{file}_add_a_add_b_add_c.txt"
        shell:
            "cat {input} | xargs echo add c >{output}"
    
    rule hebing:
        input:
           a=expand("{file}_add_a_add_b_add_c.txt",file=["1","2","3"]),
           b=expand("{file}_add_a_add_b.txt",file=["1","2"])
        output:"hebing.txt"
        shell:"cat {input.a} {input.b} >{output}"
    
    

    执行结果:

    (base) [dengfei@localhost ex4]$ snakemake -s a.py 
    Provided cores: 1
    Rules claiming more threads will be scaled down.
    Job counts:
    	count	jobs
    	3	adda
    	3	addb
    	3	addc
    	1	all
    	1	hebing
    	11
    
    rule adda:
        input: 1.txt
        output: 1_add_a.txt
        jobid: 8
        wildcards: file=1
    
    Finished job 8.
    1 of 11 steps (9%) done
    
    rule adda:
        input: 3.txt
        output: 3_add_a.txt
        jobid: 10
        wildcards: file=3
    
    Finished job 10.
    2 of 11 steps (18%) done
    
    rule adda:
        input: 2.txt
        output: 2_add_a.txt
        jobid: 9
        wildcards: file=2
    
    Finished job 9.
    3 of 11 steps (27%) done
    
    rule addb:
        input: 3_add_a.txt
        output: 3_add_a_add_b.txt
        jobid: 7
        wildcards: file=3
    
    Finished job 7.
    4 of 11 steps (36%) done
    
    rule addb:
        input: 2_add_a.txt
        output: 2_add_a_add_b.txt
        jobid: 4
        wildcards: file=2
    
    Finished job 4.
    5 of 11 steps (45%) done
    
    rule addb:
        input: 1_add_a.txt
        output: 1_add_a_add_b.txt
        jobid: 3
        wildcards: file=1
    
    Finished job 3.
    6 of 11 steps (55%) done
    
    rule addc:
        input: 3_add_a_add_b.txt
        output: 3_add_a_add_b_add_c.txt
        jobid: 2
        wildcards: file=3
    
    Finished job 2.
    7 of 11 steps (64%) done
    
    rule addc:
        input: 2_add_a_add_b.txt
        output: 2_add_a_add_b_add_c.txt
        jobid: 5
        wildcards: file=2
    
    Finished job 5.
    8 of 11 steps (73%) done
    
    rule addc:
        input: 1_add_a_add_b.txt
        output: 1_add_a_add_b_add_c.txt
        jobid: 6
        wildcards: file=1
    
    Finished job 6.
    9 of 11 steps (82%) done
    
    rule hebing:
        input: 1_add_a_add_b_add_c.txt, 2_add_a_add_b_add_c.txt, 3_add_a_add_b_add_c.txt, 1_add_a_add_b.txt, 2_add_a_add_b.txt
        output: hebing.txt
        jobid: 1
    
    Finished job 1.
    10 of 11 steps (91%) done
    
    localrule all:
        input: hebing.txt
        jobid: 0
    
    Finished job 0.
    11 of 11 steps (100%) done
    
    
    展开全文
  • 使用snakemake的好处process data with the same pipeline自动生成拓扑图ATAC-seq 基本原理ATAC-seq基本原理image.png每个核小体上可以缠绕146-147bp的DNA,大概缠绕两圈左右,但是在核小体上并不是所有的组蛋白都...
  • snakemake使用
  • Snakemake搭建流程 - 干货级

    千次阅读 2020-08-05 16:20:11
    snakemake 简介2. snakemake 安装3. snakemake 参数简介4. snakemake 使用说明4.1 定义workflow4.2 配置信息Configuration4.3 运行代码示例参考 1. snakemake 简介 最近碰到一个snakemake搭建的流程,挺好奇,便...
  • snakemake_align_MAJIQ-源码

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

    2021-03-18 08:33:55
    这是snakemake内置的工作流程,用于处理来自10X scRNA-seq实验的功能条形码库。从原始测序读取开始,管道将映射到功能条形码的参考,量化细胞条形码和UMI,并生成分析报告。 作者 罗宾·迈耶斯(@robinmeyers) 用法...
  • Nanopype是基于snakemake的管道,可提供便捷的牛津纳米Kong技术(ONT)数据处理和存储解决方案。 流水线分为处理部分(包括基本调用和对齐)以及分析部分(涵盖其他下游应用程序)。 部分提供了所有包含的工具的...
  • 蛇形 该存储库提供文件的格式。 它遵循的设计和规格。...docker run -it snakemake/snakefmt snakefmt --help 您可以在找到所有可用的标签。 奇点 URI= " docker://snakemake/snakefmt " singularity exec " $
  • 从零开始用snakemake搭建完整的甲基化生信分析流程 阅读文章后的收获:专业的snakemake编程技能;甲基化分析pipeline 作者比利时根特大学生物信息学硕士,5年生信分析师从业经验,快速响应读者问题,提供技术支持,...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 382
精华内容 152
关键字:

snakemake