https://www.bilibili.com/video/BV1jb411i76T?p=2
什么是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
网友评论