连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中一个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！
下面是温州医科大硕士“何物昂”的笔记
前言
之前在健明老师的安排下，做了几次兼职项目，体验了一把“数字游民”。
常见的组学分析的基本流程比较固定，两三次手动的分析后，开始尝试使用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://www.bilibili.com/video/BV1C7411C7ez
Snakemake流程
Snakemake简介
Snakemake是一个工作流引擎系统，提供了基于Python的可读性流程定义语言,可重现,可扩展的数据分析的工具和强大的执行环境，无需流程更改就可从单核环境迁移到集群，云服务环境上运行。
snakemake workflow 由一系列的rules 组成，每个rule为一个分析步骤，用于执行特定的功能。snakemake 流程是以输出为导向的。输出为导向 是相对于输出导向的流程，我们平常在linux 写的shell 脚本是以输入为导向的。
data/$ 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_fqrule 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.gzdata/ENCFF035OMK.fastq.gz data/ENCFF288CVJ.fastq.gz data/ENCFF820PVO.fastq.gz data/ENCFF883SEZ.fastq.gzdata/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└── scripts7 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/
data/$ tree datadata├── ENCFF035OMK.fastq.gz├── ENCFF279LMU.fastq.gz├── ENCFF288CVJ.fastq.gz├── ENCFF518FYP.fastq.gz├── ENCFF820PVO.fastq.gz├── ENCFF823XXU.fastq.gz├── ENCFF883SEZ.fastq.gz└── ENCFF888ZZV.fastq.gz0 directories, 8 files
重命名
我们手上分析的数据可能来源于不同地方，其命名方式也可能多种多样的。流程分析中先要规范下文件命名。所以本文分析流程的第一步是文件的重命名， 重命名，我们不采用提前手动更改命名的方式，而是直接集成至到分析流程中。
# 创建一个配置文件$ touch config/config.yaml
我们将文件的样本信息写到 config/config.yaml里
workdir: ./resultPE: truesample:liver_rep1:r1: /disk/public/test/data/ENCFF288CVJ.fastq.gzr2: /disk/public/test/data/ENCFF888ZZV.fastq.gzliver_rep2:r1: /disk/public/test/data/ENCFF883SEZ.fastq.gzr2: /disk/public/test/data/ENCFF035OMK.fastq.gzheart_rep1:r1: /disk/public/test/data/ENCFF279LMU.fastq.gzr2: /disk/public/test/data/ENCFF820PVO.fastq.gzheart_rep2:r1: /disk/public/test/data/ENCFF823XXU.fastq.gzr2: /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值
Snakefiletouch 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_fqoutput:"01seq/raw/{sample}_{end}.fq.gz"shell:"ln -s {input} {output}"
config.yamlrename rule 使用 shell命令为"ln -s {input} {output}"， 采用了ln 来生成对应文件的软链接来实现文件的重命名。
fastqc质控
workflow/Snakefile$ touch workflow/rules/fastqc.smk使用fastqc进行质控， 在fastqc.smk写入以下内容
## workflow/rules/fastqc.smkrule raw_fq:input:raw = rules.rename.output,output:"02fqc/raw/{sample}_{end}_fastqc.zip",conda:"../envs/test.yaml"threads: 8log:"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.smkif 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: 8log:"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: 8log:"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.smkif config["PE"]:rule btw2_map:input:rules.trim.outputoutput:"03align/{sample}.bam"conda:"../envs/test.yaml"threads: 8log:"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.outputoutput:"03align/{sample}.bam"conda:"../envs/test.yaml"threads: 8log:"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.yamlbwt2_idx: /home/public/genome/mm/release_M25/btw2/mm10
在rule btw2_map 使用params 来指示idx 参数
bam过滤
创建bamfilter.smk, 写入以下内容，使用samtools进行过滤操作
## workflow/rules/bamfilter.smkrule rmdup:input:rules.btw2_map.outputoutput:"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.outputoutput:"03align/{sample}.filter.bam"conda:"../envs/test.yaml"threads: 8log:"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.smkrule macs2_callpeak:input:rules.filter_bam.outputoutput:"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.yamlgenome: mm
peak注释
peak注释，我们借助R里的包进行注释，创建文件
$ touch workflow/scripts/peak_annotate.R$ touch workflow/rules/peak_anno.smk
在peak_anno.smk 写入以下内容
## workflow/rules/peak_anno.smkrule 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$OrgDblibrary(OrgDb, character.only = T)txdb <- snakemake@config$txdblibrary(txdb, character.only = T)txdb <- eval(parse(text = txdb))sn <- snakemake@wildcards$sample[1]outdir <- snakemake@params$outdirpeak_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.dbtxdb: TxDb.Mmusculus.UCSC.mm10.knownGene
差异peak寻找
差异peak 需要先确定实验组，对照组，同样，我们将这些信息写入在config.yaml 中
## config/config.yaml......diffGroup:lh:control:- liver_rep1- liver_rep2case:- 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 bamsdef get_diff_case_bam(wildcards):sn = get_diff_case_sn(wildcards)bams = [f"03align/{i}.filter.bam" for i in sn]return bamsdef get_diff_control_peak(wildcards):sn = get_diff_control_sn(wildcards)peaks = [f"06callpeak/{i}_peaks.narrowPeak" for i in sn]return peaksdef get_diff_case_peak(wildcards):sn = get_diff_case_sn(wildcards)peaks = [f"06callpeak/{i}_peaks.narrowPeak" for i in sn]return peaksrule 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_peakoutput:"10diffpeak/{dgroup}_result.csv",params:control_sn = get_diff_control_sn,case_sn = get_diff_case_snconda:"../envs/test.yaml"script:"../scripts/peak_diff.R"
在peak_diff.R 写入以下内容, 使用DiffBind进行差异分析
library(DiffBind)output <- snakemake@output[[1]]control_sn <- snakemake@params$control_sncase_sn <- snakemake@params$case_snsns <- 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_testchannels:- defaults- conda-forge- biocondadefault_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/msys2custom_channels:conda-forge: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloudbioconda: https://mirrors.tuna.tsinghua.edu.cn/anaconda/clouddependencies:- 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 -dresult├── 01seq│ ├── clean│ └── raw├── 02fqc│ └── raw├── 03align├── 06callpeak├── 08peakanno│ └── GO├── 10diffpeak└── logs├── bamfilter├── bowtie2├── callpeak├── fastqc│ └── raw└── trim17 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.smkrule make_report:input:diff_peak_resultoutput:"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: truetoc_float: true---```{r setup, include=FALSE}library(ggplot2)library(ChIPseeker)library(enrichplot)samples <- snakemake@params$samplesdiffgroups <- snakemake@params$diffgroupspeak_dir <- snakemake@params$peak_dirdiffpeak_dir <- snakemake@params$diffpeak_dirsp1.Rdata <- paste(peak_dir, "/", samples[1], "_anno.Rdata", sep="")load(sp1.Rdata)```### Peak CallingPeak 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) <- samplespeak_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) <- samplespeak_count_df <- t(as.data.frame(peak_counts))colnames(peak_count_df) <- samplesrownames(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(内含子)或DistalIntergenic(远端基因间区)。```{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://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
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
