snakemake+Anaconda个人自用入门笔记

最新推荐文章于 2022-12-27 14:41:13 发布

同志我累了

最新推荐文章于 2022-12-27 14:41:13 发布

阅读量1.7k

点赞数

分类专栏： liunx生物信息文章标签： linux python

本文链接：https://blog.csdn.net/weixin_45503019/article/details/108128812

版权

liunx生物信息专栏收录该内容

7 篇文章 0 订阅

订阅专栏

snakemake使用

snakemak各人想法：
是一款能够通过可以按照流程进行计算并将流程过程进行保存的，可以定制重复的流程，比如说有个任务是将一个文件的里的样本进行两两比对，那么这些两两比对是同一个过程的话，那么可以以通过使用snakemake提交后台，并且让这任务同时运行，每个两两比对任务是独特的，他们都会有一个编号ID，这些编号是同时运行的，他们所用的资源也是我们可以分配的，你分配是一个core两个线程，那他们就是跑一个核两个线程，同时任务是互不干扰的，所以即使其中有个两两比对无法正常运行也不会终止其他两两比对的任务，当提交的ID全部完成后（这里我猜测存在着一种判断，当ID全部完成后才可以进行下一个流程也就是rule），关于log是怎么来的，其实是通过捕抓每个shell也就是操i在时候的输出，就是我们操作中间会有一些输出比如说1+1 会出现2而这个2我们通过赋值就给了log中的文本，其实也就是说这个log是我们习惯而来的捕抓，叫pog也是没事的内容不会变。

rule {rulename}:
	input:
	output:
	shell:'sort{input}>{output}'
	    #命令shell,R,python

Anaconda

由于每个人使用的环境不一样，软件版本不一致导致无法在一共地方运行的服务器能在另一个服务器正常运行，所以通过Anaconda来拷贝系统环境，Anaconda相当于虚拟机

安装教程

https://zhuanlan.zhihu.com/p/35711429

查看有哪些虚拟环境

conda env list

创建python3环境

conda create -n py3test python=3*

激活环境

source activate py3test

安装snakemake

conda install snakemake

关闭环境

source deactivate

snakemake

测试Anaconda

 python -v  #查看Python版本
 # /usr/lib64/python2.6/encodings/utf_8.pyc matches /usr/lib64/python2.6/encodings/utf_8.py
import encodings.utf_8 # precompiled from /usr/lib64/python2.6/encodings/utf_8.pyc
Python 2.6.6 (r266:84292, Aug 18 2016, 15:13:37) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-17)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
# 使用snakemake脚本将两个文件夹的内容合并

############################################################
#title:ATAC-seq
#autor---同志我累了
#founction
############################################################
#rule bwt:                           |   +++++++++++++++++++
#       input:"words.txt"            |  ++  =   +++    0   ++    
#       output:"words.sorted.txt"    |   +++++++++++++++++++ 
#       shell:                       |       +++++++++++ 
#             "sort{input}>{output}" |         +++++++
#                                    |           +++
#------------------------------------------------------------



rule concet: #rule snakemake定义语，concet为任务名可以自己更改
        input: #input snakemake 的定义语,定义输入文件
                expand("{file}.txt",file =["hello","world"]) #expand为替换命令，snakemake自带
        output: #output snakemake 的定义语,定义输出结果的保存文件
                "merged.txt"
        shell:  # 这里表示我们下面的命令将在命令行中执行
                "cat {input}>{output}"

snakemake  --snakefile snaketest  --cores

运行结果
在这里插入图片描述

然后更改代码，看看snake有没有权限自己创建文件

在这里插入图片描述

发现确实可以自己创建文件夹，上面有重叠处是因为我是运行了下面代码再上的上一张图。

rule all

在这里插入图片描述

构建简单的snkamake分析ATAC-seq

##########################################################
#2020/8/30
#同志干起来了
#test ATAC-seq pipline with snake
#v 1。00
##########################################################

这是输入的，类似于你使用input一样这是你在运行snkamake前为他准备的参数和一些数据，包才能正常运行，可以看玩后面再来看上面

REP_INDEX = {"rep1","rep2"}
INDEX_BT2 = "" #提前建好的index
PICARD=""#软件包的位置。

这一步是为总的snakemake流程控制输出，也就是运行后保留的数据，没有统计在rule input中的文件将不保存

rule all:
    input:
       expand( "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort.bam",rep=REP_INDEX)
       expand( "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort_rmdup..bam",rep=REP_INDEX)
       expand( "macs2_result/ATAC_seq_{rep}_peaks.narrowPeak",rep=REP_INDEX)

rule cutdadapt:
    input:
         "raw.fastq/ATAC_seq_{rep}_R1.fq.gz",
         "raw.fastq/ATAC_seq_{rep}_R2.fq.gz"
    output:
        "fix.fastq/ATAC_seq_{rep}_R1.fq.gz",
        "fix.fastq/ATAC_seq_{rep}_R2.fq.gz"
    log:"fix.fastq/ATAC_seq_{rep}.cutadapt.log,rep=REP_INDEX"
    shell:
        "cutadapt -j 1 -times 1 -e 0.1 -0 3 --quality-cutoffy 25 \
                -m 50 -a AGATCGGAAGACACGTCTGAACTCCAGTCAG \
                -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTTAGATCTCGGTGGTCGCCGTATCATT \
                -o {utput[0]} -p { output[1]}{input[0]} {input[1]}>{log} 2>&1"

rule bt2_mapping:
    input:
         "fix.fastq/ATAC_seq_{rep}_R1.fq.gz",
         "fix.fastq/ATAC_seq_{rep}_R2.fq.gz"
    output:
        "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38.sam"
    log: "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38.log"
    shell:
        "bowtie2 -x {INDEX_BT2} -p 20 -1 {input[0]} -2 {input[1]} -s {output} > {log} 2>&1"
        #-s表示输出为sam文件，搞定比对了完成比对了

rule bam_file_sort:
        input:
            "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38.sam"
        output:
             "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort.bam"
        log:
            "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort.log"
        shell:
            "samtools sort -0 BAM -o {output} -T {output}.te,p -@ 4 -m 2G {input}"

rule remove_dupliaction:
    input:
         "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort.bam"
    output:
        "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort_rmdup..bam",
        "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort_rmdup..matrix"
    log:
        "bam/fix.fastq/ATAC_seq_{rep}_bt2_hg38_sort_rmdup..log"
    shell:
        "java -Xms5G -Xmx5G -xx:parallelGCThreads=4 \
                -jar {PICARD} MarkDuplicates \
                I={input} 0={output[0]} M={ output[1]} \
                ASO=coordinate REMOVE_DUPLICATES=true 2>{log}"

rule call_peak:
    input:
        "bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.bam"
    output:
        "macs2_result/ATAC_seq_{rep}_peaks.narrowPeak"
    params:
        "ATAC_seq_{rep}",
        "macs2_result"
    log:
    ""
    shell:
        "mac2 callpeak -t {input} -f BAM -g hs --outdir {params[1]} -n {params[0]} -m 2 100 > {log} 2>&1"