zoukankan      html  css  js  c++  java
  • 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	ID:{sample}	PL:illumina	SM:{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

     

    本文版权归作者和博客园共有,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利.
  • 相关阅读:
    English Voice of <<Cups>>
    【线段树】奶牛排队(USACO 2007 January Gold)
    【线段树】买水果
    【线段树】卫星覆盖(NOI97)-矩阵切割
    插入排序 (Insertion Sort)
    选择排序 (Selection Sort)
    springboot整合redis
    redis入门及相关API
    mycat配置文件的详细介绍
    redis常用命令
  • 原文地址:https://www.cnblogs.com/mmtinfo/p/15343697.html
Copyright © 2011-2022 走看看