文章目录
参考： https://snakemake.readthedocs.io/en/stable/tutorial/additional_features.html
额外的特征
接下来，我们将介绍更多上面的案例没有包括的特征。更多细节和更多特征，可以参考 Writing Workflows, Frequently Asked Questions或者使用命令行帮助(snakemake --help)。
1. Benchmarking（基准测试）
使用benchmark指令，Snakemake可以度量一项工作的运行时间。我们可以在bwa_map中激活benchmarking:
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:
repeat("benchmarks/{sample}.bwa.benchmark.txt", 3)
threads: 8
shell:
"(bwa mem -R '{params.rg}' -t {threads} {input} | "
"samtools view -Sb - > {output}) 2> {log}"
可以多次重复基准测试，了解测量量的可变性。 这可以通过另一个注释benchmark文件来实现，例如：repeat("benchmarks/{sample}.bwa.benchmark.txt", 3)。Snakemake会被告知执行三次。重复度量可以作为tab分割的基准文件的后续行。
运行代码为：
##删掉保护的bam文件
rm -rf /share/inspurStorage/home1/jialh/tools/snakemake/snakemake-tutorial/sorted_reads
##重新运行整个流程
snakemake --forceall --cores 4
2. Modularization（模块化）
为了重用构建的模块或者简化结构化的大型流程，有时候需要合理地split a workflow into modules。对此，Snakemake提供了include指令来将另一个Snakefile包含进当前的Snakefile，例如：
include: "path/to/other.snakefile"
另外，Snakemake还允许定义子流程（define sub-workflows， 详情可参考https://snakemake.readthedocs.io/en/stable/snakefiles/modularization.html#snakefiles-sub-workflows）。基本示例如下：
subworkflow otherworkflow:
workdir:
"../path/to/otherworkflow"
snakefile:
"../path/to/otherworkflow/Snakefile"
configfile:
"path/to/custom_configfile.yaml"
rule a:
input:
otherworkflow("test.txt")
output: ...
shell: ...
练习
将读段比对相关的rules作为分开的Snakefile, 然后用include指令，使它们能够在我们的示例流程中工作。
mapping_snakefile文件中写入：
##snakemake --cores 4 -s Snakefile
##snakemake --cores 4 mapped_reads/B.bam
##snakemake --cores 4 mapped_reads/C.bam
def get_bwa_map_input_fastqs(wildcards):
return config["samples"][wildcards.sample]
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:
repeat("benchmarks/{sample}.bwa.benchmark.txt", 3)
threads: 8
shell:
"(bwa mem -R '{params.rg}' -t {threads} {input} | "
"samtools view -Sb - > {output}) 2> {log}"
主snakefile为：
include: "/home1/jialh/tools/snakemake/snakemake-tutorial/mapping_Snakefile" ##注意需要使用引号
#SAMPLES = ["A", "B"]
configfile: "config.yaml"
rule all:
input:
"plots/quals.svg"
##snakemake --cores 4 sorted_reads/B.bam
##snakemake --cores 4 sorted_reads/{A,B,C}.bam
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 --cores 4 sorted_reads/{A,B}.bam.bai
##snakemake --dag sorted_reads/{A,B}.bam | dot -Tsvg > dag_test.svg
rule samtools_index:
input:
"sorted_reads/{sample}.bam"
output:
"sorted_reads/{sample}.bam.bai"
shell:
"samtools index {input}"
##snakemake --cores 4 calls/all.vcf
print(config["pmr"])
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:
pmr=config["pmr"]
log:
"logs/bcftools_call/all_vcf.log"
shell:
"bcftools mpileup -f {input.fa} {input.bam} | "
"bcftools call -mv -P {params.pmr} - > {output}"
##snakemake --dag calls/all.vcf | dot -Tsvg > dag_vcf.svg
##使用自己写的脚本。
##snakemake --cores 4 plots/quals.svg
rule plot_quals:
input:
"calls/all.vcf"
output:
"plots/quals.svg"
script:
"scripts/plot-quals.py"
3.Automatic deployment of software dependencies（自动部署软件依赖项）
为了完全重复数据分析，能够执行每个步骤和指定所有使用的参数是不够的。使用的软件工具和库也需要写清楚。在本教程中，你已经看到了Conda如何为整个流程指定独立的软件环境。使用Snakefile， 你能够更进一步，指定每个rule的Conda环境。这种方式，你可以使用冲突的软件版本（例如，结合Python2和Python3）。
在我们的案例中，除了外部环境，我们可以指定每条rule的环境，例如：
rule samtools_index:
input:
"sorted_reads/{sample}.bam"
output:
"sorted_reads/{sample}.bam.bai"
conda:
"envs/samtools.yaml"
shell:
"samtools index {input}"
其中，envs/samtools.yaml可以被定义为：
channels:
- default
- bioconda
- conda-forge
dependencies:
- samtools =1.9
Snakemake可以被执行为：
snakemake --use-conda --cores 4 --forceall
如果自动创建需要的环境和然后在执行工作前激活它们。最好是指定该环境需要的packages的最低和最高版本。上述运行结果如下：
会尝试下载和激活所需的环境。
4. Tool wrappers（工具包装）
为了简化流行工具的使用，Snakemake提供了包装资源库（ Snakemake wrapper repository）。A wrapper是一段短的脚本（通常包装了命令行应用，使之在Snakemake内部能够直接处理。对此，Snakemake提供了wrapper指令来代替shell, script或者run。例如bwa_map可以被写成如下形式：
rule bwa_mem:
input:
ref="data/genome.fa",
sample=lambda wildcards: config["samples"][wildcards.sample]
output:
temp("mapped_reads/{sample}.bam")
log:
"logs/bwa_mem/{sample}.log"
params:
"-R '@RG\tID:{sample}\tSM:{sample}'"
threads: 8
wrapper:
"0.15.3/bio/bwa/mem"
wrapper指令需要（部分）URL直接指向存储库中的包装器。这些可以在对应的数据库（database）中看到。在调用时，Snakemake将自动下载请求的包装器版本。进一步，结合--use-conda， 所需的软件将自动在执行之前下载。
[注意：按照实际操作看，–use-conda会导致程序执行很慢，而且可能导致奇怪的错误，慎用。]
5. Cluster execution
默认情形下，Snakemake在本地机器上执行工作。然而，其也可以在分布式环境下执行工作，例如：compute cluster or batch system。如果节点共享公共的文件系统，Snakemake支持三种可选的执行模式。
在集群环境中，计算工作通常通过类似qsub的命令提交给shell脚本。Snakemake提供了一种generic mode在clusters上执行。通过调用Snakemake：
snakemake --cluster qsub --jobs 100
此部分暂时用不到， 略过。
6. Using –cluster-status
可用于检测，是否一个cluster job被完全执行成功了。
7.Constraining wildcards
Snakemake使用正则表达式来匹配输入文件和输出文件，进而决定jobs之间的依赖关系。 有时候，限制一个通配符的值非常有用。这可以通过增加描述允许通配符的值的集合的正则表达式来实现。 例如，通配符sample在输出文件"sorted_reads/{sample}.bam"可以限制其只允许字母数字式样的sample name为"sorted_reads/{sample,[A-Za-z0-9]+}.bam"。每个规则都可以定义限制，或者在全局使用wildcard_constrains关键字， 正如Wildcards所展示的。这种机制能够解决两种类型的混淆。
- 其可以避免歧义规则，例如两个或者多个规则被用于产生相同的输出文件。其他处理混淆规则的方法会在 Handling Ambiguous Rules中描述。
- 其有助于基于匹配指导正则表达式，使得通配符被分配到文件名的正确部分。考虑输出文件
{sample}.{group}.txt， 假设目标文件是A.1.normal.txt。不清楚到底是dataset="A.1"和group="normal"， 还是dataset="A"和group="1.normal"是正确的赋值。此处，通过{sample,[A-Z]+}.{group}限制数据集通配符可以解决此问题。
当解决混淆规则时，最好的办法是首先尝试通过合适的文件结构解决混淆，例如，通过将不同步骤的输出文件划分到不同的目录。