snakmake 小练习

最近在学习snakemake 用于生信流程管理，现在用一个snakemake 来完成小任务：将在某一文件夹下的多个bam文件截取一部分，然后建立索引，在提取出fastq序列，最后比对回基因组。

需要两个文件，一个配置文件config.yaml和snakemake文件。

 

config.yaml 文件内容用空格分隔，主要是要用到的软件，也可以包含一些流程中要用的参数

TMP_DIR: /home/share/work/TMP 

HG19: /home/share/BioSoft/hg19_bwa/ucsc.hg19.fasta

BWA: /home/share/BioSoft/bin/bwa

 

snakemake文件：我一般都以 .py 做后缀名，主流程写在这里

configfile: "./config.yaml"

hg19=config["HG19"]
tmp=config["TMP_DIR"]
Bwa=config["BWA"]

SAMPLES,=glob_wildcards("bam/{sample}.05.mapped.bam")

rule all:
    input:
        expand("out/{sample}.bam.bai",sample=SAMPLES),
        expand("out/bwa/{sample}.sam",sample=SAMPLES)

rule split:
    input:
        "bam/{sample}.05.mapped.bam"
    output:
        "out/{sample}.bam"
    log:
        "log/{sample}.log"
    shell:
        "samtools view {input} chr7:3240000-5260000  -b -o {output} 2> {log}"

rule index:
    input:
        "out/{sample}.bam"
    output:
        "out/{sample}.bam.bai"
    shell:
        "samtools index {input} {output}"

rule fq:
    input:
        "out/{sample}.bam"
    output:
        r1="out/{sample}_r1.fq",
        r2="out/{sample}_r2.fq",
    params:
        null="/dev/null"
    shell:
        "samtools fastq -N {input} -1 {output.r1} -2 {output.r2} -0 {params.null} -s {params.null}"

rule bwa:
    input:
        r1="out/{sample}_r1.fq",
        r2="out/{sample}_r2.fq"
    output:
        "out/bwa/{sample}.sam"
    params:
        RG=r"'@RG\tID:{sample}\tPL:illumina\tSM:{sample}'",
        ref=hg19,
        bwa=Bwa
    threads: 3
    log:
        "log/bwa/{sample}.log"
    shell:
        "{params.bwa} mem -t {threads} -M -R {params.RG} {params.ref} {input.r1} {input.r2} -o {output}"

 

写完后伪运行下，检查各个步骤代码

snakemake -np  -s snakemake.py

 

检查没有问题之后投递计算节点，任务结束后还可以做个流程图：

snakemake --dag -s snakemake.py |dot -Tpdf > pipeline.pdf