正文
snakmake 小练习
提示：扫一扫查出行【扫一扫了解最新限行尾号】
复制提示
最近在学习snakemake 用于生信流程管理，现在用一个snakemake 来完成小任务：将在某一文件夹下的多个bam文件截取一部分，然后建立索引，在提取出fastq序列，最后比对回基因组。
需要两个文件，一个配置文件config.yaml和snakemake文件。
config.yaml 文件内容用空格分隔，主要是要用到的软件，也可以包含一些流程中要用的参数
1 TMP_DIR: /home/share/work/TMP
2
3 HG19: /home/share/BioSoft/hg19_bwa/ucsc.hg19.fasta
4
5 BWA: /home/share/BioSoft/bin/bwa
snakemake文件：我一般都以 .py 做后缀名，主流程写在这里
1 configfile: "./config.yaml"
2
3 hg19=config["HG19"]
4 tmp=config["TMP_DIR"]
5 Bwa=config["BWA"]
6
7 SAMPLES,=glob_wildcards("bam/{sample}.05.mapped.bam")
8
9 rule all:
10 input:
11 expand("out/{sample}.bam.bai",sample=SAMPLES),
12 expand("out/bwa/{sample}.sam",sample=SAMPLES)
13
14 rule split:
15 input:
16 "bam/{sample}.05.mapped.bam"
17 output:
18 "out/{sample}.bam"
19 log:
20 "log/{sample}.log"
21 shell:
22 "samtools view {input} chr7:3240000-5260000 -b -o {output} 2> {log}"
23
24 rule index:
25 input:
26 "out/{sample}.bam"
27 output:
28 "out/{sample}.bam.bai"
29 shell:
30 "samtools index {input} {output}"
31
32 rule fq:
33 input:
34 "out/{sample}.bam"
35 output:
36 r1="out/{sample}_r1.fq",
37 r2="out/{sample}_r2.fq",
38 params:
39 null="/dev/null"
40 shell:
41 "samtools fastq -N {input} -1 {output.r1} -2 {output.r2} -0 {params.null} -s {params.null}"
42
43 rule bwa:
44 input:
45 r1="out/{sample}_r1.fq",
46 r2="out/{sample}_r2.fq"
47 output:
48 "out/bwa/{sample}.sam"
49 params:
50 RG=r"'@RG\tID:{sample}\tPL:illumina\tSM:{sample}'",
51 ref=hg19,
52 bwa=Bwa
53 threads: 3
54 log:
55 "log/bwa/{sample}.log"
56 shell:
57 "{params.bwa} mem -t {threads} -M -R {params.RG} {params.ref} {input.r1} {input.r2} -o {output}"
写完后伪运行下，检查各个步骤代码
1 snakemake -np -s snakemake.py
检查没有问题之后投递计算节点，任务结束后还可以做个流程图：
1 snakemake --dag -s snakemake.py |dot -Tpdf > pipeline.pdf
