安装
由于服务器没有连外网只能访问学校自己的镜像站，在镜像中发现snakemake的python最多只有3.6版本于是利用conda建立了一个新环境：
conda -n snakemake python=3.6
进入该环境后安装snakemake
conda activate snakemake
conda install snakemake
此时产生报错CondaValueError: Malformed version string ‘~’ : invalid character(s)
搜了半天有说要把conda更新的，但是我更新的时候仍然显示该错误。
于是将~/.condarc中镜像源的 …/clode/conda-forge删除
顺利更新
之后成功安装
示例程序
步骤1：映射读取
下载连接如下：
wget https://bitbucket.org/snakemake/snakemake-tutorial/get/v5.2.3.tar.bz2
snakemake使用手册地址：snakemake手册
按照手册中安装environment中的软件:
- snakemake-minimal =5.2.4
- python =3.6
- jinja2 =2.10
- networkx =2.1
- matplotlib =2.2.3
- graphviz =2.38.0
- bcftools =1.9
- samtools =1.9
- bwa =0.7.17
- pysam =0.15.0
默认的文件是Snakefile
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}"
#测试pipeline
snakemake -np mapped_reads/A.bam
#执行pipeline
snakemake --cores 1 mapped_reaas/A.bam
为什么要加 mapped_reads/A.bam?
samtools报错：缺少libtinfow.so 安装相关软件包解决
conda install -c conda-forge ncurses
2:通用化读取映射规则
rule bwa_map:
input:
"data/genome.fa",
"data/samples/{sample}.fastq"
output:
"mapped_reads/{sample}.bam"
shell:
"bwa mem {input} | samtools view -Sb - > {output}"
使用-np测试
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
bwa mem data/genome.fa data/samples/B.fastq | samtools view -Sb - > mapped_reads/B.bam
Job counts:
count jobs
1 bwa_map
1
所以 这个参数是target 用来匹配通配符？
(snakemake) a@test02:~/study/snakemake/tutorial$ snakemake -np mapped_reads/A.bam
Nothing to be done.
文件夹中如果已经有A.bam则不会产生信息（没有输入文件比输出文件更新）
还可以使用bash命令中的{ }进行通配符的匹配
snakemake mapped_reads/{A,B,C}.bam
3.对读取对齐方式进行排序
Snakemake 在执行作业之前自动创建缺少的目录
较长的shell命令也可以分为两行
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访问通配符（input和output已经定义好了，其他通配符没有声明）
#如果删除wildcards
NameError: The name 'sample' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}}
4.索引读取对齐并可视化作业的 DAG
重新运行第2步的pipieline获得bam文件
我把命令放在了 /tutorial/S1 中
--snakefile FILE, -s FILE
The workflow definition in a snakefile.
snakemake -s S1 -np mapped_reads/{A,B,C}.bam
rule bwa_map:
input: data/genome.fa, data/samples/A.fastq
output: mapped_reads/A.bam
jobid: 0
wildcards: sample=A
bwa mem data/genome.fa data/samples/A.fastq | samtools view -Sb - > mapped_reads/A.bam
rule bwa_map:
input: data/genome.fa, data/samples/B.fastq
output: mapped_reads/B.bam
jobid: 1
wildcards: sample=B
bwa mem data/genome.fa data/samples/B.fastq | samtools view -Sb - > mapped_reads/B.bam
rule bwa_map:
input: data/genome.fa, data/samples/C.fastq
output: mapped_reads/C.bam
jobid: 2
wildcards: sample=C
bwa mem data/genome.fa data/samples/C.fastq | samtools view -Sb - > mapped_reads/C.bam
Job counts:
count jobs
3 bwa_map
3
snakemake -s S1 mapped_reads/{A,B,C}.bam
...
3 of 3 step(100%) done.
同理调用第三步中的pipeline /tutorial/S2
在绘制DAG图时发现只有一个==>所有流程都在一个文件中
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tpng > dag.png
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tpdf > dag.pdf
snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg
DAG 包含每个作业的节点，连接它们的边缘表示依赖项。 不需要运行的作业框架（因为它们的输出是最新的）是虚线的。 对于具有通配符的规则，特定作业的通配符值将显示在作业节点中。
删除掉mapped_reads和sort_reads中的文件，重新生成DAG图
分别测试各个部分的输出文件：
snakemake -np mapped_reads/{A,B,C}.bam
Job counts:
count jobs
3 bwa_map
3
snakemake -np sorted_reads/{A,B,C}.bam
Job counts:
count jobs
3 bwa_map
3 samtools_sort
6
snakemake -np sorted_reads/{A,B,C}.bam.bai
Job counts:
count jobs
3 bwa_map
3 samtools_index
3 samtools_sort
9
snakemake sorted_reads/{A,B,C}.bam.bai
#全部生成
5.调用基因组变异
先安装使用到的软件
#在我的snakemake环境中
conda install -c bioconda bcftools
Snakemake 提供了一个用于收集输入文件的辅助函数：
expand:该函数类似Python的format
SAMPLES = ["A", "B"]
expand("sorted_reads/{sample}.bam", sample=SAMPLES)
#用SAMPLES中的值分别替换掉simple
这里定义了三个输入文件的变量，在命令中可以指定对于变量名称进行访问，input.fa input.bam，整体函数如下：
SAMPLES = ["A", "B"]
rule bwa_map:
..
rule samtools_sort:
..
rule samtools_index:
..
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:
"bcftools mpileup -f {input.fa} {input.bam} | "
"bcftools call -mv - > {output}"
生成DAG图
snakemake --dag calls/all.vcf | dot -Tpng > dag.png
6.使用自定义脚本
Snakemake提供script来指定需要运行的脚本，脚本路径始终相对于引用的Snakemake文件
这些脚本在运行时都会传入一个 snakemake 对象，包含属性：input、output、params、 wildcards、 log、 threads、 resources、 config。
mkdir scripts && vim 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])
在这里使用了snakemake中的input对象中的第一个元素snakemake.input[0]
除了使用python外还可以使用R脚本—名为snakemake的S4对象，想访问input对象中的第一个元素可以使用snakemake@input[[1]]或snakemake@input[["myfile"]]
rule plot_quals:
input:
"calls/all.vcf"
output:
"plots/quals.svg"
script:
"scripts/plot-quals.py"
查看当前环境下python中的包–没有matplotlib和pysam
conda install matplotlib pysam
安装失败，conda中搜索pysam，发现只有py2的环境
# Name Version Build Channel
pysam 0.6 py26_0 pkgs/free
pysam 0.6 py27_0 pkgs/free
在在官方手册中发现bioconda中有最新的。
conda install -c bioconda pysam
成功安装
snakemake plots/quals.svg
得到直方图
7.添加目标规则
因此，可以编写目标规则来收集所需结果或所有结果的特定子集。 此外，如果在命令行中没有给出目标，Snakemake会将Snakefile的第一条规则定义为目标。因此，最佳做法是在工作流的顶部设置一个规则，该规则将所有通常需要的目标文件作为输入文件。
在顶部添加 relu all后就可以不指定输出文件来规定任务，而是按照rule all中的input进行自动规划，Snakefile 中的规则顺序是任意的，不会影响作业的 DAG。
snakemake -n
Job counts:
count jobs
1 all
1 bcftools_call
2 bwa_map
1 plot_quals
2 samtools_index
2 samtools_sort
9
564
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
