Snakemake+RMarkdown定制你的分析流程和报告

生信技能树

发布于 2022-06-27 20:51:40

2.6K1

发布于 2022-06-27 20:51:40

文章被收录于专栏：生信技能树生信技能树

下面是温州医科大硕士“何物昂”的笔记

前言

之前在健明老师的安排下，做了几次兼职项目，体验了一把“数字游民”。

数字游民第三波有你吗 https://mp.weixin.qq.com/s/q864LQvsOOmd9nUyxk939w

数字游民从学徒作业开始 https://mp.weixin.qq.com/s/b3rR--dUwAZSvibF07-WQQ

常见的组学分析的基本流程比较固定，两三次手动的分析后，开始尝试使用snakemake搭建分析流程，以及发现配合RMarkdown可以自动化分析数据然后生成对应的分析报告。虽然使用shell脚本或者其他编程语言也是能实现一个分析流程的。不过这样的话，需要考虑的细节问题就有许多，比如：

路径问题，结果或日志文件的输出，需要提前创建好对应的父目录
需要自行编写特定命令实现并行运算
总线程数控制，内存资源控制
调用其他语言的脚本运行任务，还得考虑如何进行参数传递
断点运行，要是程序中断，得考虑从程序从哪里中断的，然后从哪里重新开始运行
......

不过更主要的是，我想要一个直接分析完然后直接生成结果报告的流程。因为一开始提供给用户分析结果时，我都是手动将部分内容复制到Typora里，然后生成pdf/html的，这很麻烦，而且容易出错。snakemake里是提供了report 功能。不过日常分析中，我们常用R语言，不少文档也用Rmarkdown写出来，可能用Rmarkdown起来更熟悉和方便一些。

这里使用snakemake 来实现一个ATAC-Seq的分析流程，同时采用Rmarkdown 来生成一个简单的分析报告。大致包含以下内容：

fastq质控
fastq比对
bam过滤
callpeak
peak注释
peak邻近基因功能富集
差异peak寻找

ATAC-Seq

ATAC-Seq 介绍和教程参考：

入门ATAC-seq，你需要知道什么？https://mp.weixin.qq.com/s/txb22VjztZhCiEyKMI6GAw
ATAC-seq数据分析实战 https://www.bilibili.com/video/BV1C7411C7ez

Snakemake流程

Snakemake简介

Snakemake是一个工作流引擎系统，提供了基于Python的可读性流程定义语言,可重现,可扩展的数据分析的工具和强大的执行环境，无需流程更改就可从单核环境迁移到集群，云服务环境上运行。

snakemake workflow 由一系列的rules 组成，每个rule为一个分析步骤，用于执行特定的功能。snakemake 流程是以输出为导向的。输出为导向是相对于输出导向的流程，我们平常在linux 写的shell 脚本是以输入为导向的。

比如，我们像把data/的文件的fastaq.gz 改为fq.gz, 如果是平常shell 脚本的写法则为：

$ cp data/ENCFF035OMK.fastq.gz data/ENCFF035OMK.fq.gz
# 对于多个文件可以写for 循环
$ for i in $(ls data/*.fastq.gz);do cp $i data/$(basename $i fastq.gz)fq.gz; done

输入导向的运行方式，需要先确定输入文件.

如果是在输出导向的snakemake 中，则需要先确定输出文件。

new_fq = [
        "data/ENCFF035OMK.fq.gz",
        "data/ENCFF279LMU.fq.gz",
        "data/ENCFF288CVJ.fq.gz",
        "data/ENCFF518FYP.fq.gz",
        "data/ENCFF820PVO.fq.gz",
        "data/ENCFF823XXU.fq.gz",
        "data/ENCFF883SEZ.fq.gz",
        "data/ENCFF888ZZV.fq.gz"
]

rule all:
    input:
        new_fq

rule change_suffix:
    input:
        "data/{sample}.fastq.gz"
    output:
        "data/{sample}.fq.gz"
    shell:
        "cp {input} {output}"

在shell 命令中的cp 命令，在snakemake中，写成一个rule change_suffix，rule中的input, output，则由wildcards "sample"表示组成的字符表达式。rule all 用来确定流程最后的输出文件是哪些。

snakemake wildcards ，类似于linux 的通配符,用来匹配对应的字符，这里用来匹配样本名

$ ls data/*.fastq.gz
data/ENCFF035OMK.fastq.gz  data/ENCFF288CVJ.fastq.gz  data/ENCFF820PVO.fastq.gz  data/ENCFF883SEZ.fastq.gz
data/ENCFF279LMU.fastq.gz  data/ENCFF518FYP.fastq.gz  data/ENCFF823XXU.fastq.gz  data/ENCFF888ZZV.fastq.gz

在snakemake 中，先通过rule all input 确定了输出文件new_fq，继而在其他rule output中寻找可以匹配的字符表达式。对应上面的小例子中，在rule change_suffix 找到了对应匹配的output。即new_fq 可以匹配 "data/{sample}.fq.gz"，确定了{sample}实际值，进而确定input

额，不要嫌原来shell 命令只要一行就能解决的问题，改成了snakemake要多了许多行代码。对于简单的日常任务是shell要方便许多。但是对于一个稍显复杂分析流程而言，使用snakemake 会更合适。

软件安装

$ conda install -c conda-forge -c bioconda snakemake -y

项目结构

$ mkdir {config,data,result,workflow}
$ mkdir workflow/{rules,envs,scripts}
$ tree -d
.
├── config
├── data
├── result
└── workflow
    ├── envs
    ├── rules
    └── scripts

7 directories

其中：

config 存放配置文件
data 存放原始数据
result 存放分析的结果
workflow 分析流程相关代码相关文件
- envs 存放环境配置文件
- rules 存放smk文件
- scripts 存放代码脚本

数据

这里用ENCODE里小鼠的 heart tissue embryo (11.5 days)， liver tissue embryo (11.5 days) 两个数据作为演示例子：

https://www.encodeproject.org/experiments/ENCSR785NEL/
https://www.encodeproject.org/experiments/ENCSR820ACB/

下载好fastq数据放在data/ 文件夹内,

$ tree data
data
├── ENCFF035OMK.fastq.gz
├── ENCFF279LMU.fastq.gz
├── ENCFF288CVJ.fastq.gz
├── ENCFF518FYP.fastq.gz
├── ENCFF820PVO.fastq.gz
├── ENCFF823XXU.fastq.gz
├── ENCFF883SEZ.fastq.gz
└── ENCFF888ZZV.fastq.gz

0 directories, 8 files

重命名

我们手上分析的数据可能来源于不同地方，其命名方式也可能多种多样的。流程分析中先要规范下文件命名。所以本文分析流程的第一步是文件的重命名，重命名，我们不采用提前手动更改命名的方式，而是直接集成至到分析流程中。

# 创建一个配置文件
$ touch config/config.yaml

我们将文件的样本信息写到 config/config.yaml里

workdir: ./result
PE: true
sample:
  liver_rep1:
    r1:  /disk/public/test/data/ENCFF288CVJ.fastq.gz
    r2: /disk/public/test/data/ENCFF888ZZV.fastq.gz
  liver_rep2:
    r1:  /disk/public/test/data/ENCFF883SEZ.fastq.gz
    r2:  /disk/public/test/data/ENCFF035OMK.fastq.gz
  heart_rep1:
    r1:  /disk/public/test/data/ENCFF279LMU.fastq.gz
    r2: /disk/public/test/data/ENCFF820PVO.fastq.gz
  heart_rep2:
    r1:  /disk/public/test/data/ENCFF823XXU.fastq.gz
    r2:  /disk/public/test/data/ENCFF518FYP.fastq.gz

配置文件采用yaml 格式, 这是一种简单的格式。

YAML 语言教程: http://ruanyifeng.com/blog/2016/07/yaml.html

目前配置文件中，目前定义了3个对象:

workdir: 设置工作目录
PE: 用来确定是否为paired-end 测序数据
sample 样本信息，其下一级为样本名：
- liver_rep1 样本名自定义，再下一级为read1.read2样本数据
- r1: read1的文件
- r2: read2的文件
- se,如果是单端的，我们使用se 作为key值

然后编写代码进行文件的更名，创建Snakefile 文件，snakemake默认运行该文件的内容

touch workflow/Snakefile

## workflow/Snakefile

# 设置配置文件
configfile: "config/config.yaml"
# 设置工作目录
workdir: config["workdir"]

## 获取配置文件中的样本名
SAMPLES = config["sample"].keys()

## 单端双端的一些配置
if config["PE"]:
    ENDS = ["r1", "r2"]
else:
    ENDS = ["se"]


def get_fq(wildcards):
    return config["sample"][wildcards.sample][wildcards.end]

rule rename:
    input:
        get_fq
    output:
        "01seq/raw/{sample}_{end}.fq.gz"
    shell:
        "ln -s {input} {output}"

rename 是自定义的文件重命名rule, rule 中的output 为"01seq/raw/{sample}_{end}.fq.gz"，而input 是一个函数，因为原文件是类似ENCFF279LMU.fastq.gz，这样的命名，没法直接推导出input 文件，所以这里借用一个函数，来获取匹配到的{sample}和{end}，通过{sample}和{end}实际值，来获取config.yaml 中定义的样本文件。

rename rule 使用 shell命令为"ln -s {input} {output}"，采用了ln 来生成对应文件的软链接来实现文件的重命名。

fastqc质控

流程下一步进行fastq质控, 创建一个对应smk文件来执行质控功能，所有rules都可以直接写在workflow/Snakefile里, 但这里我们将不同功能分别写在不同文件里，进行模块分离，将其写在单独的文件里。

$ touch workflow/rules/fastqc.smk

使用fastqc进行质控，在fastqc.smk写入以下内容

## workflow/rules/fastqc.smk

rule raw_fq:
    input: 
        raw = rules.rename.output,
    output:
        "02fqc/raw/{sample}_{end}_fastqc.zip",
    conda:
        "../envs/test.yaml"           
    threads: 8        
    log:
        "logs/fastqc/raw/{sample}_{end}.log",
    shell:
        "fastqc -o 02fqc/raw -f fastq -t {threads} --noextract {input} 2> {log}"

rule raw_fq 用于raw fastq的qc, 其中input 可以是"01seq/raw/{sample}_{end}.fq.gz"，也可以直接用rules.rename.output 来引用rule output的输出。用conda 来指定特定conda环境，用threads 来限定线程数， log 来指定输出日志。

去除adapter

创建trimming.smk

$ touch workflow/rules/trimming.smk

采用trim_galore进行接头去除，在trimming.smk 写入以下内容

## workflow/rules/trimming.smk 

if config["PE"]:
    rule trim:
        input:
            "01seq/raw/{sample}_r1.fq.gz",
            "01seq/raw/{sample}_r2.fq.gz"
        output:
            "01seq/clean/{sample}_r1_val_1.fq.gz",
            "01seq/clean/{sample}_r2_val_2.fq.gz"
        conda:
            "../envs/test.yaml"   
        threads: 8 
        log:
            "logs/trim/{sample}.log"          
        shell:
            "trim_galore -j {threads} -q 25 --phred33 --length 36 --stringency 3 -o 01seq/clean --paired {input} 2> {log}"
else:
    rule trim:
        input:
            "01seq/raw/{sample}_se.fq.gz"
        output:
            "01seq/clean/{sample}_se_trimmed.fq.gz"
        conda:
            "../envs/test.yaml"              
        threads: 8 
        log:
            "logs/trim/{sample}.log"          
        shell:
            "trim_galore -j {threads} -q 25 --phred33 --length 36 --stringency 3 -o 01seq/clean {input} 2> {log}"

这里通过if else 来根据来返回不同的trim rule，以适应不同的场景。snakemake 是基于Python扩展的，Python原来的语法照样可以在snakmake里使用。

比对

创建bowtie2.smk,

$ touch workflow/rules/bowtie2.smk

写入以下内容，采用bowtie2进行比对

## workflow/rules/bowtie2.smk 


if config["PE"]:
    rule btw2_map:
        input:
            rules.trim.output
        output:
            "03align/{sample}.bam" 
        conda:
            "../envs/test.yaml"
        threads: 8 
        log:
            "logs/bowtie2/{sample}.log"
        params:
            idx = config['bwt2_idx']
        shell:
            "bowtie2  -X2000 --mm -x {params} -p {threads} -1 {input[0]} -2 {input[1]} 2> {log} | "
            "samtools sort -O bam -o {output}"
else:
    rule btw2_map:
        input:
            rules.trim.output
        output:
            "03align/{sample}.bam" 
        conda:
            "../envs/test.yaml"
        threads: 8 
        log:
            "logs/bowtie2/{sample}.log"
        params:
            idx = config['bwt2_idx']
        shell:
            "bowtie2 -X2000 --mm -x {params} -p {threads} -U {input} 2> {log} | "
            "samtools sort -O bam -o {output} && samtools index {output} -@ {threads}"

bowtie2 比对需要index, 在config.yaml 设置bowtie2的index

## config/config.yaml

bwt2_idx: /home/public/genome/mm/release_M25/btw2/mm10

在rule btw2_map 使用params 来指示idx 参数

bam过滤

创建bamfilter.smk, 写入以下内容，使用samtools进行过滤操作

## workflow/rules/bamfilter.smk 


rule rmdup:
    input:
        rules.btw2_map.output
    output:
        "03align/{sample}.rmdup.bam"
    conda:
        "../envs/test.yaml"
    log:
        "logs/bamfilter/{sample}.rmdup.log"
    params:
        mode = "" if config["PE"] else "-s"
    shell:
        "samtools rmdup {params.mode} --output-fmt BAM {input} {output} 2> {log}"

rule filter_bam:
    input:
        rules.rmdup.output
    output:
        "03align/{sample}.filter.bam"
    conda:
        "../envs/test.yaml"
    threads: 8
  log: 
        "logs/bamfilter/{sample}.filter.log"
    params:
        mode = "" if config["PE"] else "-s"
    shell:
        "samtools view -@ {threads} -q 30  -h {input} | grep -v -P '\tchrM\t' |samtools view -Sb -o {output} "

callpeak

创建callpeak.smk, 写入以下内容，采用macs2，进行callpeak

## workflow/rules/callpeak.smk 

rule macs2_callpeak:
    input:
        rules.filter_bam.output
    output:
        "06callpeak/{sample}_peaks.narrowPeak"
    conda:
        "../envs/test.yaml"
    params:
        np = "06callpeak/{sample}", 
        fp = "BAMPE" if config["PE"] else "",
        gp = config["genome"]
    log:
        "logs/callpeak/{sample}.log"
    shell:
        """
        macs2 callpeak -t {input[0]} -g {params.gp} -f {params.fp} -n {params.np} -B -q 0.05 2> {log}
        """

这里，添加了一个新的参数config["genome"]，在config.yaml中添加该参数

## config/config.yaml

genome: mm

peak注释

peak注释，我们借助R里的包进行注释，创建文件

$ touch workflow/scripts/peak_annotate.R
$ touch workflow/rules/peak_anno.smk

在peak_anno.smk 写入以下内容

## workflow/rules/peak_anno.smk 

rule peak_annoate:
    input:
         "06callpeak/{sample}_peaks.narrowPeak"
    output:
        "08peakanno/{sample}_peak_anno.csv",
    conda:
        "../envs/test.yaml"
    params:
        outdir = "08peakanno"
    script:
        "../scripts/peak_annotate.R"

其中，"../scripts/peak_annotate.R" 是用来进行注释的R脚本, 里面采用ChIPseeker进行peak注释，clusterProfiler进行邻近基因GO注释

在snakemake rule中除了利用 shell 运行shell命令外，还可以通过script来直接调用脚本，当前支持Python, R, R Markdown, Julia 和Rust。

在peak_annotate.R 写入以下内容：

library(ChIPseeker)
library(clusterProfiler)
library(ggplot2)

OrgDb <- snakemake@config$OrgDb
library(OrgDb, character.only = T)


txdb <- snakemake@config$txdb
library(txdb, character.only = T)
txdb <- eval(parse(text = txdb))

sn <- snakemake@wildcards$sample[1]
outdir <- snakemake@params$outdir


peak_path <- file.path(getwd(), snakemake@input[1])
peak_anno <- annotatePeak(peak_path, tssRegion=c(-3000, 3000), TxDb=txdb)
peak_anno_df <- as.data.frame(peak_anno)

peak_gene_df <- AnnotationDbi::select(eval(parse(text = OrgDb)),
                                      keytype = "ENTREZID",
                                      keys = peak_anno_df$geneId,
                                      columns = c("ENTREZID", "SYMBOL", "GENENAME"))


coln <- c("seqnames", "start", "end", "width", "V4", "V7", "annotation", "geneStart", "geneEnd", "geneLength",
          "geneStrand", "geneId", "transcriptId", "distanceToTSS", "SYMBOL")

new_peak_anno_df <- cbind(peak_anno_df, peak_gene_df)[, coln]
colnames(new_peak_anno_df)[c(5, 6)] <- c("peak_Id", "peak_score")



### 保持peak注释文件为CSV文件
peak.anno.fn <- sprintf("%s_peak_anno.csv", sn)
write.csv(new_peak_anno_df, file = file.path(outdir,  peak.anno.fn))    


### Peak 邻近基因 GO 富集分析
all_ego <- enrichGO(
  gene          = new_peak_anno_df$geneId,
  keyType       = "ENTREZID",
  OrgDb         = OrgDb,
  ont           = "ALL",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.1,
  readable      = TRUE,
  pool          = FALSE)
dir.create(file.path(outdir, "GO"), showWarnings = F)            
save.image(sprintf("%s/GO/%s_GO_BP.RData", outdir, sn))            
if (dim(as.data.frame(all_ego))[1] > 0){
  p <- dotplot(all_ego, split="ONTOLOGY", showCategory = 10)  + facet_grid(ONTOLOGY~., scale="free")
  ggsave(sprintf("%s/GO/%s_GO_all.pdf", outdir, sn), width=10, height = 10, plot=p)
  write.csv(as.data.frame(all_ego), file = sprintf("%s/GO/%s_GO_all.csv", outdir, sn))
  p
}

BP_ego <- enrichGO(
  gene          = new_peak_anno_df$geneId,
  keyType       = "ENTREZID",
  OrgDb         = OrgDb,
  ont           = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.1,
  readable      = TRUE,
  pool          = FALSE)


if (dim(as.data.frame(BP_ego))[1]>0){
  BP_p <- dotplot(BP_ego, showCategory = 12) 
  ggsave(sprintf("%s/GO/%s_GO_BP.pdf", outdir, sn), width=10, height = 10, plot=BP_p)
  write.csv(as.data.frame(BP_ego), file = sprintf("%s/GO/%s_GO_BP.csv", outdir, sn))
}

CC_ego <- enrichGO(
  gene          = new_peak_anno_df$geneId,
  keyType       = "ENTREZID",
  OrgDb         = OrgDb,
  ont           = "CC",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.1,
  readable      = TRUE,
  pool          = FALSE)
if (dim(as.data.frame(CC_ego))[1] > 0){
  CC_p <- dotplot(CC_ego, showCategory = 12) 
  ggsave(sprintf("%s/GO/%s_GO_CC.pdf", outdir, sn), width=10, height = 10, plot=CC_p)
  write.csv(as.data.frame(CC_ego), file = sprintf("%s/GO/%s_GO_CC.csv", outdir, sn))
}

MF_ego <- enrichGO(
  gene          = new_peak_anno_df$geneId,
  keyType       = "ENTREZID",
  OrgDb         = OrgDb,
  ont           = "MF",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.1,
  qvalueCutoff  = 0.1,
  readable      = TRUE,
  pool          = FALSE)
if (dim(as.data.frame(MF_ego))[1]>0){
  MF_p <- dotplot(MF_ego, showCategory = 12) 
  ggsave(sprintf("%s/GO/%s_GO_MF.pdf", outdir, sn), width=10, height = 10, plot=MF_p)
  write.csv(as.data.frame(MF_ego), file = sprintf("%s/GO/%s_GO_MF.csv", outdir, sn))
}

Rdata.out <- sprintf("%s/%s_anno.Rdata", outdir, sn)

save.image(Rdata.out)

代码有点长，先把它折叠起来了。

在config.yaml中添加一些注释使用的参数

OrgDb: org.Mm.eg.db
txdb: TxDb.Mmusculus.UCSC.mm10.knownGene

差异peak寻找

差异peak 需要先确定实验组，对照组，同样，我们将这些信息写入在config.yaml 中

## config/config.yaml

......

diffGroup:
    lh:
      control:
        - liver_rep1
        - liver_rep2
      case:
        - heart_rep1
        - heart_rep2

将差异组别信息保存在diffGroup下，其中:

lh 是差异分析组名
control 差异分析组的对照组
case 差异分析组的实验组

差异peak寻找，也借助R里的包进行分析，创建文件

$ touch workflow/scripts/peak_diff.R
$ touch workflow/rules/diffpeak.smk

在diffpeak.smk 写入以下内容

def get_diff_control_sn(wildcards):
    return config["diffGroup"][wildcards.dgroup]["control"]

def get_diff_case_sn(wildcards):
    return config["diffGroup"][wildcards.dgroup]["case"]


def get_diff_control_bam(wildcards):
    sn = get_diff_control_sn(wildcards)
    bams = [f"03align/{i}.filter.bam" for i in  sn]
    return bams

def get_diff_case_bam(wildcards):
    sn = get_diff_case_sn(wildcards)
    bams = [f"03align/{i}.filter.bam" for i in  sn]
    return bams

def get_diff_control_peak(wildcards):
    sn = get_diff_control_sn(wildcards)
    peaks = [f"06callpeak/{i}_peaks.narrowPeak" for i in  sn]
    return peaks

def get_diff_case_peak(wildcards):
    sn = get_diff_case_sn(wildcards)
    peaks = [f"06callpeak/{i}_peaks.narrowPeak" for i in  sn]
    return peaks


rule peak_diff:
    input:
        control_bam = get_diff_control_bam,
        case_bam = get_diff_case_bam,
        control_peak = get_diff_control_peak,
        case_peak = get_diff_case_peak
    output:
        "10diffpeak/{dgroup}_result.csv",
    params:
        control_sn = get_diff_control_sn,
        case_sn = get_diff_case_sn
    conda:
        "../envs/test.yaml"        
    script:
        "../scripts/peak_diff.R"

在peak_diff.R 写入以下内容, 使用DiffBind进行差异分析

library(DiffBind)
output <- snakemake@output[[1]]
control_sn <- snakemake@params$control_sn
case_sn <- snakemake@params$case_sn
sns <- c(control_sn, case_sn)
bams <- c(snakemake@input$control_bam, snakemake@input$case_bam)
peaks <- c(snakemake@input$control_peak, snakemake@input$case_peak)
samples <- data.frame("SampleID"  = sns,
                      "Tissue"    = "all",
                      "Factor"    = c(rep("A", length(control_sn)), rep("B", length(case_sn))),
                      "Replicate" = c(seq(1, length(control_sn)), seq(1, length(case_sn))),
                      "bamReads"  = bams,
                      "Peaks"     = peaks,
                      "PeakCaller" = "narrow")
dbObj <- dba(sampleSheet=samples)
dbObj <- dba.count(dbObj , bUseSummarizeOverlaps=TRUE)
dbObj <- dba.normalize(dbObj)
dbObj <- dba.contrast(dbObj, categories=DBA_FACTOR, minMembers = 2)
dbObj <- dba.analyze(dbObj, method = DBA_DESEQ2)
dbObj.report <- dba.report(dbObj)
write.csv(dbObj.report, file=output)

确定最终输出文件

snakemake 使用all rule 来收集所有最终输出文件。

raw_fq_qc_zips = expand("02fqc/raw/{sample}_{end}_fastqc.zip", sample=SAMPLES, end=ENDS)
peak_anno = expand("08peakanno/{sample}_peak_anno.csv", sample=SAMPLES)
diff_peak_result = expand("10diffpeak/{dgroup}_result.csv", dgroup=DIFFGROUPS)

rule all:
    input:
        raw_fq_qc_zips,
        peak_anno,
        diff_peak_result

我们只需在all input 里，写入程序测最终输出。没有后续程序依赖的输出，而中间步骤的输出，会有snakemake自动运行生成。

raw_fq_qc_zips 由于是fastqc.zip文件，没有后续程序依赖，索要生成它，需要指定为最终输出
peak_anno 也是，peak_anno.csv 没有后续程序依赖，索要生成它，需要指定为最终输出
diff_peak_result 为主要的最终输出，它之前上面的peak, bam 文件不要指定，因为diff_peak_result 的生成依赖于它们提前运行生成结果

conda 环境

上面中通过conda 设置conda环境为"../envs/test.yaml"，然后rule中运行的程序会自动激活conda环境，使用环境中的程序来运行。该分析流程中，所需的软件都能通过conda 安装，包括R包。

name: smk_test
channels:
  - defaults
  - conda-forge
  - bioconda  
default_channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/r
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/msys2
custom_channels:
  conda-forge: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud
  bioconda: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud

dependencies:
 - fastqc
 - trim-galore
 - samtools
 - bowtie2
 - macs2
 - meme=5.4.1
 - r-base=4.1.1
 - r-ggplot2=3.3.6
 - bioconductor-do.db=2.9.0=r41hdfd78af_12
 - bioconductor-txdb.mmusculus.ucsc.mm10.knowngene=3.10.0
 - bioconductor-clusterprofiler=4.2.0
 - bioconductor-chipseeker=1.30.0
 - bioconductor-org.mm.eg.db=3.14.0
 - bioconductor-diffbind
 - r-rmarkdown
 - pandoc

不同的rule可以设置不同conda环境，同时还可以利用本地的其他conda环境，只需把yaml文件地址改成conda 环境名就行了。

snakemake运行

snakemake流程运行

$ snakemake -c 24 -p --use-conda

-c 指定运行cpu核数
-p 打印出运行shell命令
-- use-conda 使用conda 环境

运行完后，查看下生成的目录结构

$ tree result -d
result
├── 01seq
│   ├── clean
│   └── raw
├── 02fqc
│   └── raw
├── 03align
├── 06callpeak
├── 08peakanno
│   └── GO
├── 10diffpeak
└── logs
    ├── bamfilter
    ├── bowtie2
    ├── callpeak
    ├── fastqc
    │   └── raw
    └── trim

17 directories

RMarkdown生成报告

如前面提到，snakemake可以直接运行RMarkdown 脚本。

那我们编写一个rule 来调用生成报的的Rmd程序。

创建文件

$ touch workflow/rules/make_report.smk
$ touch workflow/scripts/make_report.Rmd

在make_report.smk 写入以下内容

### workflow/rules/make_report.smk
rule make_report:
    input:
        diff_peak_result
    output:
        "ATACseq_report.html",
    conda:
        "../envs/test.yaml"
    params:
        samples = SAMPLES,
        diffgroups = DIFFGROUPS,
        peak_dir = "08peakanno",
        diffpeak_dir = "10diffpeak"
    script:
        "../scripts/make_report.Rmd"

input 来确定rule make_report能在其他rule运行完后再运行,
output 来指定生成的报告文件名，
conda指定运行环境
params 确定一些参数，让make_report.Rmd里的程序寻找生成报告所需要的文件
script Rmd脚本路径

再workflow/scripts/make_report.Rmd，写入以下内容

---
title: "ATAC-seq report"
output:
    html_document:
        toc: true
        toc_float: true
    
---

```{r setup, include=FALSE}

library(ggplot2)
library(ChIPseeker)
library(enrichplot)


samples <- snakemake@params$samples
diffgroups <- snakemake@params$diffgroups
peak_dir <- snakemake@params$peak_dir
diffpeak_dir <- snakemake@params$diffpeak_dir


sp1.Rdata <- paste(peak_dir, "/", samples[1], "_anno.Rdata", sep="")
load(sp1.Rdata)
```

### Peak Calling
Peak Calling 是指使用统计学方法计算出参考基因组上比对上的 Reads 显著富集的区域（称为 Peak），这些区域在不同的实验中表
示不同的研究重点，在组蛋白/转录因子的 ChIPseq 和 CUT&Tag 实验中为组蛋白的修饰位点/转录因子结合位点，ATACseq 实验中为开
放染色质区。获得 Peak 后，可对 Peak 进行后续的 Peak 序列 Motif 分析、 Peak 注释等分析，进一步挖掘感兴趣的方向。
本次分析使用主流的 Peak Calling 软件 MACS2。MACS2 在 Call Peak 时，可以选择使用 narrow 或 broad 参数进行分析，这两个参
数寻找 Peak 的方法略有不同，找到的 Peak 的峰形也不同，narrow 峰形较窄，broad 峰形较宽，转录因子和一些组蛋白如 H3K27ac 的
Peak 的峰形是窄的，一些组蛋白如 H3K36me3、H3K9me3 等的 Peak 是宽的。我们默认使用 narrow 参数进行分析。

所有样本 Peak 数目统计如下：
```{r , echo=FALSE, message=F}
peak_anno.csv <- paste(peak_dir, "/", samples, "_peak_anno.csv", sep="")

names(peak_anno.csv) <- samples

peak_anno.ls <- lapply(peak_anno.csv, function(file){
    read.csv(file)
})

peak_counts <- sapply(peak_anno.ls, function(df){
    dim(df)[1]
})

names(peak_counts) <- samples

peak_count_df <- t(as.data.frame(peak_counts))

colnames(peak_count_df) <- samples
rownames(peak_count_df) <- "Peaksum"
knitr::kable(peak_count_df, caption = "Sample peak number")
```

### 功能性区域 Peak 分布
对 Peak 注释分析结果中 Peak 在所覆盖的基因的功能性区域上的分布进行统计。一般来说，大部分调控因子结合的位置在启动子区域。
而基因间区上一般也有微弱的信号，因为基因间区在基因组上占比极大，所以检测到的 Peak 相对其他区域来说可能比较多，
但这种 Peak 一般不是真的调控因子结合位点。该分析使用 ChIPseeker软件实现。

下图为柱状图示例

```{r, echo=FALSE, message=F}
plotAnnoBar(peak_anno)
```

下图为饼图示例

```{r, echo=FALSE, message=F}
plotAnnoPie(peak_anno)
```


Peak 相对于TSS 数量分布比例                                                                                                                                                                                                                                 

```{r, echo=FALSE, message=F}
plotDistToTSS(peak_anno,
                title="Distribution of transcription factor-binding loci\nrelative to TSS")
```

### Peak 注释

转录因子一般结合在基因的转录起始位点附近，对其上下游基因的转录进行调控。而结合在基因外显子和内含子区域的调控因子
则可能影响该基因的可变剪切行为。Peak 覆盖的基因及邻近位置的基因均有可能是转录因子或其他调控因子调控的靶基因，获取这些
基因并进行注释，有助于研究目标蛋白调控的基因功能或代谢通路。我们使用 R包ChIPseeker(Yu G, Wang L, He Q 2015)的 annotatePeak
函数注释Peak所在的基因组特征。

本分析注释的功能性区域包括Promoter(启动子)、Exon(外显子)、UTR(非翻译区，有5’UTR和3’UTR)、Intron(内含子)或Distal
Intergenic(远端基因间区)。


```{r, echo=FALSE, message=F}
### 显示peak注释表格
knitr::kable(new_peak_anno_df[1:10, -c(5, 6,10,11,12,13)], caption = "peak annotation")
```

### Peak 邻近基因 GO 富集分析

 Peak 邻近基因 GO 富集分析的目的是得出  Peak 邻近基因 显著富集的 GO 条目， 从而发现不同发育时期、不同细胞状态
下可能被转录因子调控的基因功能 。对每个  Peak 邻近基因 进行 GO 注释，得到其所有的 GO 功能条目，将所有  Peak 邻近
基因 分别从 MF（Molecular Function）、BP（Biological Process）、 CC（Cellular Component）这三个层面分配到不同的 GO 功能分类
中，采用 R 包 clusterProfiler 计算每个 GO 分类中基因富集的显著性水平，从而筛选出  Peak 邻近基因 显著富集的 GO 分类。(pval<0.1, qval<0.1)

```{r, echo=FALSE, message=F}


if (dim(as.data.frame(all_ego))[1] > 0){
    p <- dotplot(all_ego, split="ONTOLOGY", showCategory = 10)  + facet_grid(ONTOLOGY~., scale="free")
    p
}
```

###  组间差异 Peak 分析

ATAC-seq 实验一般都会设计多组样本，以研究不同发育时期、不同细胞状态下的特异性染色质开放区域。通过找到样本间的差异
Peak，可以得到某种状态下特异性染色质开放区域，进而找到特异性结合的 Motif、靶基因等，推测表观遗传是否在发育过程及疾病中
起到一定作用，从而进行后续机制研究。分析方法为，首先将每个样本的 Peak 文件合并，然后使用 bedtools 工具对合并之后的 Peak
文件进行处理，如果两个 Peak 有重叠区域，则合并成一个新的 Peak。计算每个样本在每个合并的新 Peak 区域上的 Read 数目，最后
使用 DESeq2 进行差异分析，得到样本间的差异 Peak 即差异染色质开放区域。

差异 Peak 计算结果:

```{r, echo=FALSE, message=F}


diff_peak_file <- paste(diffpeak_dir, "/", diffgroups[1], "_result.csv", sep ="")
diff_peak_df <- read.csv(diff_peak_file)


knitr::kable(diff_peak_df[1:10, ], caption = "diff peak ")

```

内容有点长，先折叠了。

将 "ATACseq_report.html" 添加到rule all input里

rule all:
    input:
        raw_fq_qc_zips,
        peak_anno,
        diff_peak_result,
        "ATACseq_report.html"

运行

$ snakemake -c 24 -p --use-conda

运行后生成ATACseq_report.html，其内容如下

其他

上述流程，我是成功运行了一遍的。理论上对读者来说是非常友好的，前提是你具备基础的计算机知识，

我把它粗略的分成基于R语言的统计可视化，以及基于Linux的NGS数据处理：

因为大家使用时，可能遇到一些问题：

使用conda 环境时，可能遇到类似这样的错误

subprocess.CalledProcessError: Command 'conda export -n 'xx'' returned non-zero exit status

该问题参考该解决方案

https://stackoverflow.com/questions/69001097/conda-4-10-3-and-snakemake-5-conda-exe-problem

使用yaml配置安装conda环境时，自动安装的依赖包可能用不了，可以更换环境或者手动重新安装
一些snakemake 错误提示，具体问题具体分析了
也不排除上文代码，我从本地复制粘贴到这里时，出现问题。

文末

如大家所见，上文中不管背景又或者代码都没有详细的介绍，同时本文流程及分析报告也稍显简陋，该文抛砖引玉。不管ATAC-Seq或者snakemake，还是Rmarkdown网上都有许多优秀的教程，相信大家能创建出更好的流程报告来~

参考

《R数据科学》

https://snakemake.readthedocs.io/en/stable/index.html

https://mp.weixin.qq.com/s/txb22VjztZhCiEyKMI6GAw

https://www.bilibili.com/video/BV1C7411C7ez

https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf

https://stackoverflow.com/questions/69001097/conda-4-10-3-and-snakemake-5-conda-exe-problem

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2022-05-31，如有侵权请联系 cloudcommunity@tencent.com 删除

https

网络安全

数据分析

shell

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！