无重复样本ChIP-Seq|CUT&Tag的Snakemake自动化流程

   基于snakemake构建了无生物学重复样本的ChIPSeq完整的分析流程。

主要功能是从原始的FASTQ数据开始，通过多个步骤进行质量控制、数据修剪、比对、Peak调用以及差异分析。具体流程如下：

质量控制

：使用FastQC对原始和修剪后的数据进行质量评估，确保数据符合分析标准。

数据修剪

：使用trim_galore修剪FASTQ数据，去除低质量序列。

比对

：通过bowtie2将修剪后的数据比对到参考基因组（如hg38）。

BAM文件处理

：利用samtools对比对后的数据进行排序和索引，生成BAM格式，并生成BigWig文件以便于可视化。

Peak调用

：使用MACS3对比对结果进行Peak调用，提取信号峰值（如narrowPeak、summits等）。

差异分析

：基于不同条件组（处理组与对照组）的Peak文件，执行差异分析，生成差异Peak结果，供后续的生物学解读使用。

### 设置样本和变量 ###SAMPLES = ["KO", "WT"]treat="KO"control="WT"reference="/mnt/8_2/way/datas/refBeds/hg38.bed"index = "/mnt/8_2/way/datas/ref-genome/hg38"

rule all:   input:       # 所有目标文件，包含质量控制报告、比对结果、MACS分析输出、差异分析结果等       expand("qc1/{sample}_1_fastqc.html", sample=SAMPLES),       expand("qc1/{sample}_2_fastqc.html", sample=SAMPLES),       expand("qc2/{sample}_1_val_1_fastqc.html", sample=SAMPLES),       expand("qc2/{sample}_2_val_2_fastqc.html", sample=SAMPLES),       expand("aligned/{sample}_sorted.bam", sample=SAMPLES),       expand("alignment_stats/{sample}_stats.txt", sample=SAMPLES),       expand("aligned/{sample}.bw", sample=SAMPLES),       expand("macs/{sample}_predicted_fragment", sample=SAMPLES),       expand("macs/{sample}_summits.bed", sample=SAMPLES),       expand("macs/{sample}_tags", sample=SAMPLES),       expand("macs/{sample}_control_lambda.bdg", sample=SAMPLES),       expand("macs/{sample}_peaks.narrowPeak", sample=SAMPLES),       expand("macs/{sample}_peaks.xls", sample=SAMPLES),       expand("macs/{sample}_summits.bed", sample=SAMPLES),       expand("macs/{sample}_treat_pileup.bdg", sample=SAMPLES),       expand("macs/{sample}_mean_frag_len", sample=SAMPLES),       "final_output_file.txt",       f"macs/diffPeaks/diff_{treat}vs{control}_c1.0_common.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c1.0_cond1.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c1.0_cond2.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c2.0_common.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c2.0_cond1.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c2.0_cond2.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c3.0_common.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c3.0_cond1.bed",       f"macs/diffPeaks/diff_{treat}vs{control}_c3.0_cond2.bed",       expand("plot/{name}", name=["all.pdf", "all_heatmap.pdf", "all.gz"])

### 质量检查阶段：原始数据 ###rule fastqc:   input:       "{sample}/{sample}_1.fq.gz",  # 输入的第一对FASTQ文件       "{sample}/{sample}_2.fq.gz"   # 输入的第二对FASTQ文件   output:       "qc1/{sample}_1_fastqc.html",  # 生成的FASTQC报告1       "qc1/{sample}_2_fastqc.html"   # 生成的FASTQC报告2   log:         "qc1/{sample}_fastqc"          # 日志文件   threads: 2   params : jobname = "{sample}"      # 设定工作名称   message: "fastqc {input}: {threads}"   shell:       """       # fastqc对.gz文件进行质量检查       fastqc -o qc1 {input[0]} 2> {log}       fastqc -o qc1 {input[1]} 2> {log}       """

### 修剪读取：去除低质量序列 ###rule trim_reads:   input:       R1="{sample}/{sample}_1.fq.gz",  # 原始FASTQ文件1       R2="{sample}/{sample}_2.fq.gz"   # 原始FASTQ文件2   output:       f1="{sample}/{sample}_1_val_1.fq.gz",  # 修剪后的FASTQ文件1       f2="{sample}/{sample}_2_val_2.fq.gz"   # 修剪后的FASTQ文件2   shell:       """       trim_galore -j 4 -q 25 --phred33 --length 35 --stringency 3 --paired --gzip {input.R1} {input.R2} -o {wildcards.sample}       """

### 质量检查阶段：修剪后的数据 ###rule fastqc_2:   input:       "{sample}/{sample}_1_val_1.fq.gz",  # 修剪后的FASTQ文件1       "{sample}/{sample}_2_val_2.fq.gz"   # 修剪后的FASTQ文件2   output:       "qc2/{sample}_1_val_1_fastqc.html",  # 修剪后数据的FASTQC报告1       "qc2/{sample}_2_val_2_fastqc.html"   # 修剪后数据的FASTQC报告2   log:         "qc2/{sample}_fastqc"              # 日志文件   shell:       """       # fastqc对.gz文件进行质量检查       fastqc -o qc2 {input[0]} 2> {log}       fastqc -o qc2 {input[1]} 2> {log}       """

### 比对：使用Bowtie2将修剪后的数据比对到参考基因组 ###rule bowtie:   input:       fastq1="{sample}/{sample}_1_val_1.fq.gz",  # 修剪后的FASTQ文件1       fastq2="{sample}/{sample}_2_val_2.fq.gz"   # 修剪后的FASTQ文件2   output:       "aligned/{sample}.sam"  # 输出的SAM文件   params:       core="20"  # 设置使用的线程数   shell:       "bowtie2 --local --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -x {index} -1 {input.fastq1} -2 {input.fastq2} -p {params.core} > {output}"

### 转换为BAM文件并生成索引，生成大规模数据格式 ###rule samtools:   input:       sam="aligned/{sample}.sam"  # 输入的SAM文件   output:       bam="aligned/{sample}_sorted.bam",  # 输出的排序后的BAM文件       stats="alignment_stats/{sample}_stats.txt",  # BAM文件统计信息       bw="aligned/{sample}.bw"  # 输出的bigwig文件，用于可视化   params:       core="20"  # 设置使用的线程数   shell:       """       samtools view -@ {params.core} -bS {input} | samtools sort -@ {params.core} -o {output.bam}       samtools index -@ {params.core} {output.bam}       samtools flagstat {output.bam} > {output.stats}       bamCoverage -v --binSize 50 --normalizeUsing RPKM --extendReads -p {params.core} -b {output.bam} -o {output.bw}       rm {input.sam}       """

### 使用MACS3预测片段长度 ###rule macs3_pred:   input:       bam="aligned/{sample}_sorted.bam"  # 输入的排序后的BAM文件   output:       fragment="macs/{sample}_predicted_fragment"  # 预测的片段长度   shell:       """       macs3 predictd -i {input.bam} 2>&1 | grep --line-buffered -oP 'predicted fragment length is \K(\d+)' > {output.fragment}       """

### 聚合片段长度数据 ###rule aggregate:   input:       expand("macs/{sample}_predicted_fragment", sample=SAMPLES)  # 输入所有样本的片段长度文件   output:       allok="final_output_file.txt"  # 聚合输出文件   shell:       "echo 1 > {output.allok}"

### 使用MACS3进行Peak调用 ###rule macs3_callpeak:   input:       bam="aligned/{sample}_sorted.bam",  # 输入排序后的BAM文件       fragment="macs/{sample}_predicted_fragment",  # 输入预测的片段长度       allok="final_output_file.txt"  # 聚合的片段长度文件   output:       tags="macs/{sample}_tags",  # 生成的tags文件       con_bdg="macs/{sample}_control_lambda.bdg",  # 控制组的BDG文件       narrowPeak="macs/{sample}_peaks.narrowPeak",  # 生成的narrowPeak文件       peaks="macs/{sample}_peaks.xls",  # 生成的peaks Excel文件       bed="macs/{sample}_summits.bed",  # Summit位置的BED文件       pileup="macs/{sample}_treat_pileup.bdg",  # 处理后的堆积图   shell:       """       macs3 callpeak -t {input.bam} -c {input.allok} --format BAM --gsize hs --nomodel --shift -100 --extsize 200 -n {wildcards.sample} --outdir macs       """

### 差异Peak分析 ###rule macs3_diffpeaks_c1:   input:       treatment="macs/{treat}_peaks.xls",  # 处理组的Peak文件       control="macs/{control}_peaks.xls",  # 对照组的Peak文件       allok="final_output_file.txt"  # 聚合的片段长度文件   output:       "macs/diffPeaks/diff_{treat}vs{control}_c1.0_common.bed"  # 生成的差异Peak文件   shell:       """       macs3 callpeak -t {input.treatment} -c {input.control} --format xls -n {wildcards.treat}vs{wildcards.control} -B --SPMR --outdir macs/diffPeaks       """

如何运行：

准备环境：

安装Snakemake（可以通过conda或pip安装）。

安装所需的工具，如FastQC,Trim Galore,Bowtie2,Samtools,MACS3和MACS2。

准备数据：

将你的测序数据放入指定的文件夹结构中，确保数据文件格式正确（_1.fq.gz和_2.fq.gz）。

执行工作流：

使用命令行进入包含Snakefile的目录。

运行以下命令来启动分析：

snakemake --cores <number_of_cores>

4. 结果分析：

输出的结果将包括：

QC 报告（HTML格式）

比对后的BAM文件及索引文件

峰值调用结果

差异峰分析文件（例如：diff_1.0_common.bed等）

最后再总结下前面流程定义的所有的规则：

1. all规则

功能

：定义所有输出文件的依赖关系，最终的目标文件包括原始数据的质量控制报告、比对结果、MACS分析输出、差异峰值分析结果等。

2. fastqc规则

功能

：对原始FASTQ数据进行质量检查，生成质量控制报告（qc1/{sample}_1_fastqc.html和qc1/{sample}_2_fastqc.html）。

输入

：原始FASTQ文件。

输出

：FastQC生成的HTML报告。

3. trim_reads规则

功能

：对原始FASTQ数据进行修剪，去除低质量的序列并保留高质量部分。

输入

：原始的FASTQ文件。

输出

：修剪后的FASTQ文件。

4. fastqc_2规则

功能

：对修剪后的FASTQ数据进行再次质量检查。

输入

：修剪后的FASTQ文件。

输出

：第二次FastQC质量报告。

5. bowtie规则

功能

：使用bowtie2将处理后的FASTQ数据比对到参考基因组，生成SAM文件。

输入

：修剪后的FASTQ文件。

输出

：比对结果的SAM文件。

6. samtools规则

功能

：使用samtools将SAM文件转换为BAM文件并进行排序、索引，计算BAM文件的统计信息，并生成bigwig文件用于可视化。

输入

：SAM文件。

输出

：BAM文件、BAM文件的统计信息、bigwig文件。

7. macs3_pred规则

功能

：使用macs3预测片段长度，为后续的Peak调用准备数据。

输入

：BAM文件。

输出

：预测的片段长度文件。

8. aggregate规则

功能

：聚合所有样本的片段长度数据。

输入

：多个样本的片段长度预测文件。

输出

：汇总的输出文件。

9. macs3_callpeak规则

功能

：使用macs3进行Peak调用，识别ChIP-Seq实验中的富集区域，生成多个相关文件，包括Peak、summits、tags等。

输入

：BAM文件和预测的片段长度。

输出

：Peak调用的多个输出文件（如narrowPeak、xls等）。

10. macs3_diffpeaks_c1、macs3_diffpeaks_c2、macs3_diffpeaks_c3规则

功能

：进行差异Peak分析，比较处理组和对照组的Peak差异，生成不同条件下的差异Peak文件（common、cond1、cond2等）。

输入

：MACS分析生成的tags、control和treatment的文件，片段长度文件。

输出

：差异Peak文件（根据不同的条件和标准进行过滤）。

这些规则定义了ChIP-Seq数据分析中的各个关键步骤，从数据的质量检查、修剪、比对、Peak调用，到差异分析，最终为研究者提供了差异化的Peak区域及其他分析结果。