什么是snakemake

官网地址：https://snakemake.readthedocs.io/en/stable/

使用snakemake的好处

process data with the same pipeline

自动生成拓扑图

ATAC-seq 基本原理

ATAC-seq基本原理

image.png

每个核小体上可以缠绕146-147bp的DNA,大概缠绕两圈左右，但是在核小体上并不是所有的组蛋白都有DNA缠绕；

有些地方转录比较活跃，chromatin accessibility 染色质可接近性比较大，那么这些裸露的DNA的转录活性高；

此时，我们可以用Tn5这种酶，这种酶就可以通过链置换反应，直接把测序建库的接头adapter直接就添加到裸露的DNA上，简单来说就是Tn5把裸露的DNA链打断，并通过链置换反应加上测序建库需要的adapter。

然后再对整个基因组测序，就能够捕获这样的信号，通过mapping到genome上，就知道哪一段DNA序列是裸露的了。

因为越裸露的区域，最终捕获的reads数目就越多。

如果某一段DNA紧紧地缠绕在核小体周围，那它基本上就是没有信号的。

那最后我们通过鉴定atac-seq的峰在哪里，我们就可以知道染色质哪里是比较开放的。

image.png

那么一个基因如果处于开放染色质，那么这个基因就是高表达的；

那么如果是一个promoter在开放空间的话，那么它就很可能导致下游的基因高表达；

那么如果是一个enhancer在atac-seq信号比较弱的区域，那么它可能不太可能发生近端的相互作用；

ATAC-seq的数据分析pipeline

image.png

1. 去除接头，cutadapt；

2. 对于切除adapter的fastq文件，我们需要使用类似于bowtie2的软件去mapping到参考基因组上；

3. 把比对完的结果要进行排序，并把它记录成BAM文件；

4. 对BAM文件进行remove PCR duplication;

5. 找peak，peak calling

snakemake基础教学与上机实践

snakenake

1. snake的基本单元是rule，一个rule接一个rule；

2. 每一个rule包括input，output和你要在shell里面执行的命令；

3. 然后是下一步的输入，下一步的输出；

4. 最后的输入和最后的输出；

5. 同时还可以写python运行；

Snakemake Rule

image.png

rule sort: #关键词+rule的名字；

4个空格 #换起一行前面一定要有4个空格；

input #然后就是input、output和shell

output #然后就是input、output和shell

shell #然后就是input、output和shell

花括号里是可被替换的内容

rule bwt

使用anaconda 来管理软件和运行环境

image.png

可以用 “which python” 来查看，是在虚拟环境下的python还是外部的python

image.png

image.png

查看是否下载完好，用 “snakemake -h”，如果出现说明文档就证明已经安装好了；

image.png

然后退出环境即可

snakemake分析atac-seq的流程

每个raw fastq文件分别有两个atac-seq的数据，分别是rep1和rep2；因为是双端测序，每个rep有reads1和reads2的数据

reads分别是一一对应的

用sublime text写pipeline

######################################################

# 2020/05/15

# Yue Qu

# test ATAC-seq pipeline with snakemake

######################################################

## rep 和 index_BT2 一直没有填,可以最后填

REP_INDEX = {"rep1","rep2"}

INDEX_BT2 = "~/menghw_HD/reference/bowtie2_hg38/hg38_only_chromosome" # 每个人的index的引用路径不一样，自行修改

## bowtie2 软件包也要写明软件的路径

#1. $ cd ~/menghw_HD/software_package/

#2. $ ls

#3. $ bbmap Bismark_v0.20.0 ... ...

#4. $ cd bin/

#5. $ ls

#6. $ picard.jar VarScan.v2.3.9.jar ... ...

#7. $ pwd

#8. $ /home/menghaowei/menghw_HD/software_package/bin # picard.jar的路径

## picard 这样引用不容易报错

PICARD = "/home/menghaowei/menghw_HD/software_package/bin/picard.jar"

## rule all 应该写要保留的中间结果和要保留的最后结果

rule all:

input:

expand("bam/ATAC_seq_{rep}_bt2_hg38_sort.bam",rep=REP_INDEX) # 保留step 3的output结果-bam 文件的保存；需要用expand把rep填充进去；

expand("bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.bam",rep=REP_INDEX) # 保留step 4的output结果-bam 文件的保存；

expand("macs2_result/ATAC_seq_{rep}_peaks.narrowPeak",rep=REP_INDEX) # 保存call_peaks的结果

## step 1

rule cutadapt:

input:

# "ATAC_seq_{rep}_R{index}.fq.gz"

"raw.fastq/ATAC_seq_{rep}_R1.fq.gz" #因为ATAC_seq再raw fastq里面，所以前面要加路径

"raw.fastq/ATAC_seq_{rep}_R2.fq.gz"

output:

"fix.fastq/ATAC_seq_{rep}_R1.fq.gz" # 如果fix.fastq文件存在，则把ATAC_seq 文件放到文件夹中，如果不存在则创建该文件夹再放入；

"fix.fastq/ATAC_seq_{rep}_R2.fq.gz"

log:

"fix.fastq/ATAC_seq_{rep}_cutadapt.log" # 要保存的log文件

shell:

# 分行写后面加斜线“\”，后面不能有任何字符

"cutadapt -j 4 --times 1 -e 0.1 -0 3 --quality-cutoff 25 \

-m 50 -a AGATCGCGATTGCGTAGTCACGTACGTTGCATGAGCTTA \

-A AGATCGCGATTGCGTAGTCACGTACGTTGCATGAGCTTA \

-o {output[0]} -p {output[1]} {inpuit[0]} {input[1]} > {log} 2>&1"

# input在命令行中的引用，input[0]代表第一个文件，input[1]代表第二个文件

# 输出的话，看cutadapt说明，用 “-o” “-p” 进行输出

# log 文件的输出按照linux的标准格式输出即可

## step 2

rule bt2_mapping:

input:

"fix.fastq/ATAC_seq_{rep}_R1.fq.gz" # 这一步的input就是上一步的output

"fix.fastq/ATAC_seq_{rep}_R2.fq.gz"

output:

# output得换一个新的文件夹

"bam/ATAC_seq_{rep}_bt2_hg38.sam" # hg38，命名看个人习惯，这里作者写上了参考基因组的版本

log:

"bam/ATAC_seq_{rep}_bt2_hg38.log"

shell:

"bowtie2 -x {INDEX_BT2} -p 4 -1 {input[0]} -2 {input[1]} -S {output} > {log} 2>&1"

# -x 是index；-p 后面是需要的核心的数目；-1和-2分别是两个input；-S 是指输出成sam文件

## step 3

rule bam_file_sort:

input:

"bam/ATAC_seq_{rep}_bt2_hg38.sam"

output:

"bam/ATAC_seq_{rep}_bt2_hg38_sort.bam"

log:

"bam/ATAC_seq_{rep}_bt2_hg38_sort.log"

shell:

"samtools sort -O BAM -o {output} -T {output}.temp -@ 4 -m 2G {input}" # 输出成BAM格式；输出文件是 “-o”；“-T” 是临时文件；“-@ 4”代表核心数4个；“-m 2G” 代表内存使用2G

## step 4

rule remove_duplication:

input:

"bam/ATAC_seq_{rep}_bt2_hg38_sort.bam"

output:

"bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.bam" # rmdup这步骤通常用picard软件来处理，这个软件处理以后会输出两个文件，一个是“bam”，一个是“matrix”

"bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.matrix"

log:

"bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.log"

shell:

# Xms5g 表示用多少g内存，最小5g，最大5g；-XX以后是设置它的核心数，如果不设置的话会把服务器跑满，这样可能反而会变慢

"java -Xms5g -Xmx5g -XX:ParallelGCThreads=4 \

-jar {PICARD} MarkDuplicates \

I={input} O={output[0]} M={output[1]} \

ASO=coordinate REMOVE_DUPLICATES=true 2>{log}"

# output[0]和output[1]分别指第一个和第二个文件

## step 5

rule call_peak:

input:

"bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.bam"

output:

# 输出一个call_peak的结果

# 我们一般用macs2进行call_peak

"macs2_result/ATAC_seq_{rep}_peaks.narrowPeak" # rep1的call peak结果

params:

# 此外我们需要生成一个参数，首先要写清楚前缀名，把它当作一个变量进行保存

“ATAC_seq_{rep}”,

# 第二个的话，要把上面输出的路径当作一个参数传进去

"macs2_result"

# 这个是macs2本身shell的问题，不是说所有地方都是这么操作的

log:

"bam/ATAC_seq_{rep}_bt2_hg38_sort_rmdup.log"，

shell：

# -c = control组；-t = treatment组/input；用的是BAM cell；用的是hs = homo sapensis的基因组;

# --outdir {params[1]} 就是把params里面的1(params中的第二个文件)输出到这个目录下；

# 输出的项目的名字呢就叫做params[0],也就是params中的第一个文件“ATAC_seq_{seq}”;

# -m 2 100 = 我们只输出大于2倍和小于100倍的数集；

# 最后我们进行call summits --- 这个加不加都可以；

# 最后的最后 > {log} 2>&1

"macs2 callpeak -t {input} -f BAM -g hs --outdir {params[1]} -n {params[0]} -m 2 100 --call-summits > {log} 2>&1"

然后再把pipeline上传到服务器上面

上传到服务器上

image.png

然后用snakemake运行一遍

image.png

“-n” 代表不是真的运行它，而是检查它的逻辑是不是有错误;

"-p" 代表是把每一行的命令行展现出来

分为那几个部分呢

第一部分：job counts

第一部分

告诉你一共有多少个任务

第二部分

第二部分

告诉你每一步的命令行是什么；输入输出是什么；

告诉你job ID，是从小到大进行排序的；

image.png

第三部分--拓扑图的生成

image.png

生成的拓扑图过程，

image.png

输出为pdf，“dot” 是linux画图的命令

image.png

输出内容的效果，

image.png

运行

snakemake -s xxx.py -p -j 2 & # -p: 输出的命令；-j: 同时运行的任务的数目

image.png

查看是否已经在运行了

1

2

一步一步运行

中断进程

1

2

3

批量kill

常用的命令

image.png