【原】snakemake 学习笔记2

一个稍微复杂的案例, 看看snakemake的用法.

过程介绍

• 1, 安装snakemake

• 2, 新建文件

• 3, 新建一个简单的Snakemake参数文件

• 4, 扩展, 去关联输出文件

• 5, 使用全局变量, 关联文件

• 6, 批量运行

1, 安装snakemake

这里需要时python3, 不支持python2

pip3 install --user snakemake pyaml

2, 新建几个FASTQ文件

这里, 我们新建两个配对的RNA-seq数据, 格式是FASTQ的文件, 然后经过下面两步处理:

• 第一步: 数据质量控制

• 第二部: 将基因表达合并为一个文件

创建文件

• 创建genome.fa文件, 使用touch创建空文件即可

• 创建fastq文件夹

• 在fastq文件夹中, 创建Sample1.R1.fastq.gz Sample1.R2.fastq.gz Sample2.R1.fastq.gz Sample2.R2.fastq.gz四个空文件

touch genome.famkdir fastqtouch fastq/Sample1.R1.fastq.gz fastq/Sample1.R2.fastq.gztouch fastq/Sample2.R1.fastq.gz fastq/Sample2.R2.fastq.gz

创建结果, 使用tree查看:

(base) [dengfei@localhost test]$ tree.├── fastq│ ├── Sample1.R1.fastq.gz│ ├── Sample1.R2.fastq.gz│ ├── Sample2.R1.fastq.gz│ └── Sample2.R2.fastq.gz└── genome.fa
1 directory, 5 files

3, 创建snakemake参数文件

将下面代码命名为Snakefile

SAMPLES = ['Sample1', 'Sample2']
rule all: input: expand('{sample}.txt', sample=SAMPLES)
rule quantify_genes: input: genome = 'genome.fa', r1 = 'fastq/{sample}.R1.fastq.gz', r2 = 'fastq/{sample}.R2.fastq.gz' output: '{sample}.txt' shell: 'echo {input.genome} {input.r1} {input.r2} > {output}'

4, 参数解释

我们下面进行代码的讲解:

这里， 定义了一个SAMPLE的数组：

SAMPLES = ['Sample1', 'Sample2']

数组, SAMPLES,里面有两个元素: Sample1和Sample2

定义一个rule， 名称为all， input使用expand函数, 能够将数组的内容解析给{sample}

rule all:    input:     expand('{sample}.txt', sample=SAMPLES)

定义一个rule, 命名为 quantify_genes, 里面有input, output, shell, 其中{sample}是用的rule all里面的name

rule quantify_genes: input: genome = 'genome.fa', r1 = 'fastq/{sample}.R1.fastq.gz', r2 = 'fastq/{sample}.R2.fastq.gz' output: '{sample}.txt' shell: 'echo {input.genome} {input.r1} {input.r2} > {output}'

5, 运行参数

预览命令, 使用命令:

snakemake -np

参数介绍
-n 或者—dryrun, 表示只生成命令, 但是不执行命令, 可以预览一下生成的命令.

-p 或者—printshellcmds, 表示将生成的shell打印出来
注意:

-n 不执行, 只打印命令
-p 执行, 同时打印命令(shell)
两者执行的前提是结果文件还没有生成.

例子:

(snake_test) [dengfei@localhost ex2]$ snakemake -npBuilding DAG of jobs...Job counts: count jobs 1 all 2 quantify_genes 3
[Tue Apr 2 13:49:34 2019]rule quantify_genes: input: genome.fa, fastq/Sample1.R1.fastq.gz, fastq/Sample1.R2.fastq.gz output: Sample1.txt jobid: 1 wildcards: sample=Sample1
echo genome.fa fastq/Sample1.R1.fastq.gz fastq/Sample1.R2.fastq.gz > Sample1.txt
[Tue Apr 2 13:49:34 2019]rule quantify_genes: input: genome.fa, fastq/Sample2.R1.fastq.gz, fastq/Sample2.R2.fastq.gz output: Sample2.txt jobid: 2 wildcards: sample=Sample2
echo genome.fa fastq/Sample2.R1.fastq.gz fastq/Sample2.R2.fastq.gz > Sample2.txt
[Tue Apr 2 13:49:34 2019]localrule all: input: Sample1.txt, Sample2.txt jobid: 0
Job counts: count jobs 1 all 2 quantify_genes 3This was a dry-run (flag -n). The order of jobs does not reflect the order of execution