documentaries的中文
Snakemake⼊门教程(创建⼀个简单的⼯作流)
写在前⾯
既然写了教程就需要具有普适性,能适合⼤多数⼈的胃⼝,我这部分的内容以及⽰例主要还是参考了官⽅教程,但是都是我⼀步⼀步跑过的流程,所以会更有印象,送给想学 Snkaemake 但是⼀直没有去学的朋友们,这些内容对于有⽣信基础的⼈来讲,上⼿会很快,因为很多的⽣信软件都使⽤过,写起来也就没有那么晦涩,下⾯开始~
Snakemake 定义
Snakemake ⼯作流管理系统是⼀种⽤于创建可重复和可扩展的数据分析的⼯具。
⼯作流是通过⼀种⼈类可读的、基于 Python 的语⾔来描述的。它们可以⽆缝扩展到服务器、集、grid和云环境,⽆需修改⼯作流。
最后,Snakemake ⼯作流可能需要对所需软件的准备,这些软件将⾃动部署到任何执⾏环境。
安装
Snakemake 可在 PyPi 上以及通过 Bioconda 和源代码获得。可以使⽤任意⽅法安装 Snakemake,我们这⾥仅介绍使⽤ Conda 安装
通过 Conda/Mamba 安装
这是安装 Snakemake的推荐⽅式,因为Conda安装⽐较简单。
⾸先,必须已经安装了⼀个基于 Conda 的 Python3 发⾏版。推荐的选择是,它不仅提供所需的 Python 和 Conda 命令,⽽且包括 ,它是强烈推荐的Conda包管理器的极其快速和强⼤的替代品。默认的 conda 求解器有点慢,有时在选择最新的软件包版本时会出现问题。因此,建议在任何情况下都使⽤Mamba。
如果不安装 Mambaforge ,也可以直接安装 Mamba
conda install -n base -c conda-forge mamba
evalua单独安装到⼀个⼩环境
$ 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 /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映射到给定的参考基因组。为此,我们将使⽤⼯具,特别
是 subcommand 。在⼯作⽬录中,创建⼀个使⽤您选择的编辑器调⽤的新⽂件。官⽅建议使⽤,因为它为 Snakemake 提供了开箱即⽤的语法突出显⽰。在 Snakefile 中,定义以下命令:bwa mem Snakefile
第⼀个 Snakemake 命令将给定样本的reads映射到给定的参考基因组。 为此,我们将使⽤⼯具 ,特别是⼦命令 bwa mem。 在⼯作⽬录中,使⽤编辑器创建⼀个名为 Snakefile 的新⽂件。 官⽅建议使⽤ ,因为它为 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)和许多指令,这⾥是 input 、 output 和 shell 。
在 input 和 output 指令之后是那些预计将在命令中使⽤或创建的⽂件列表。
在最简单的情况下,这些只是 Python 字符串。 shell 指令后跟⼀个包含要执⾏的 shell 命令的 Python 字符串。
在 shell 命令字符串中,我们可以通过**⼤括号表⽰法(类似于 Python 格式函数)**引⽤命令的元素。
在这⾥,我们通过指定 {output} 来引⽤输出⽂件,通过指定 {input} 来引⽤输⼊⽂件。由于命令有多个输⼊⽂件,Snakemake 将连接它们,⽤空格分隔。换句话说,Snakemake 会在执⾏命令之前将 {input} 替换为 data/genome.fa data/samples/A.fastq。
shell 命令使⽤参考基因组和reads调⽤ bwa mem,并将输出通过管道传输到 samtools,后者创建包含⽐对的压缩 ⽂件。 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:python入门教程 下载
"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
⼀些语法特别⽅便。例如,可以选择在⼀次组合多个⽬标
$ 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 ⽂件中的读取对齐进⾏排序。这可以通过 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}"xss平台开发
Snakemake 使⽤Python 格式来格式化 shell 命令。在 shell 命令中使⽤⼤括号 ({}) 来表⽰其他内容。在这种情况下,必须加倍对我们上⾯提到的是bash括号扩展依靠时逃避它们,
例如:
ls {{A,B}}.txt
已经完成了三个步骤,现在可以查看⽣成的有向⽆环图 (DAG)
$ snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
Snakemake 使⽤ 提供的 dot 命令创建 DAG 的可视化。 对于给定的⽬标⽂件,Snakemake 以 dot 语⾔指定 DAG 并将其通过管道传输到 dot 命令,该命令将定义呈现为 。 渲染的 DAG 通过管道传输到⽂件 dag.svg 中,如下所⽰:
5. Calling genomic variants
我们⼯作流程的下⼀步将聚合所有样本的映射reads,并共同调⽤它们的基因组变异。 对于变体调⽤,我们将结合两个实⽤程序 和 。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。