通过 bowtie2 + macs3 分析 ChIP-Seq 数据

前言

以下是相关的官网链接和参考资料，如果想要得到除本文内容以外更全面深入的了解，建议直接跳转相关链接进行查阅：

MACS3 github:

https://github.com/macs3-project/MACS

MACS3 documentation:

https://macs3-project.github.io/MACS/

本文比对及去重等流程参考自以下资料：

Alignment and filtering from Harvard Chan Bioinformatics Core (HBC)

Histone ChIP-seq pipeline from ENCODE project (Version 2.0 with bowtie2)

本文流程根据以下资料进行过优化：

sambamba -F "not duplicate" processed bam still have duplicated (sambamba issue#477) (Update 2024/4/24)

Encode blacklist (Update 2024/4/24)

以下内容在本文中不会涉及到：

数据的质控
除 macs3 以外其他软件的安装

你可以考虑参考以下 github repository 安装用于 ChIP-seq 数据处理的 conda 环境，可补齐下文中需要的上游分析软件：

https://github.com/JuseTiZ/Blog-Snakemake-pipeline/

正文

使用 pip 安装 macs3：

1	$ pip install macs3

如果安装过程中无法安装依赖 cykhash，可以尝试使用 mamba 进行安装：

1	$ mamba install -c conda-forge cykhash

安装后可以输入 macs3 --version 查看是否成功安装。

比对

比对前，使用 bowtie2 构建基因组的索引：

1	$ bowtie2-build genome.fa indexname

将以上代码中的 genome.fa 替换为自己的基因组文件，indexname 替换为想要设定的索引名称，后续比对将用到 indexname。

构建后，对数据进行比对得到 sam 文件：

1
2
3

$ bowtie2 --threads 4 \
       -q -x indexname -U chip-seq.fastq.gz \
       -S output.sam > bowtie2.log 2>&1

--threads 使用的线程数量
-X 指定索引
-U 质控后的测序数据文件
-S 输出的 sam 文件

将 indexname 替换为之前建立的索引名称，将 chip-seq.fastq.gz 替换为 chip-seq 的数据，如果是双端则需进行一定修改，考虑到 chip-seq 数据一般都是单端数据此处不做拓展。该过程的屏幕输出将记录在 bowtie2.log 中（例如比对率等）。

以上命令在 slurm 调度系统下的一个并行示例：

#!/bin/bash
#SBATCH -J alignment
#SBATCH -o alignment.out
#SBATCH -e alignment.err
#SBATCH -n 24

input_dir="./fastq"
output_dir="./bowtie2_alignment"

mkdir -p "${output_dir}"

bowtie2 --version > bowtie2.version 2>&1

align_fastq() {
    local fastq_file="$1"
    local sam_file="${output_dir}/$(basename "${fastq_file}" .fq.gz).sam"
    local log_file="${output_dir}/$(basename "${fastq_file}" .fq.gz).log"

    bowtie2 --threads 4 -q -x indexname -U "$fastq_file" -S "$sam_file" > "$log_file" 2>&1
}

export -f align_fastq
export input_dir
export output_dir

parallel --jobs 6 align_fastq ::: "${input_dir}"/*.fq.gz

上述脚本将从当前目录下的 fastq 文件夹里读取所有 .fq.gz 结尾的文件，并通过 bowtie2 进行比对，索引为 indexname。该过程将并行 6 个比对任务，每个比对任务使用 4 个线程，通过 slurm 使用 24 个 CPU。

请根据个人情况（如空闲内存、集群信息等）进行修改以保证正常运行。在比对结束后，请对相关 log 文件进行检查，确保比对率等信息未出现问题。

排序及去重

由于 PCR 过程中可能部分片段有 bias（被过度扩增），将这些片段去重将有助于降低 call peak 时的假阳性。

单个运行时的示例：

$ samtools view -h -S -b -o [output_unsorted_bam] [input_sam_file] # 转为 bam
$ sambamba sort -t 4 -o [output_sorted_bam] [input_unsorted_bam] # 对 bam 排序
$ sambamba markdup -t 4 [output_markdup_bam] [input_sorted_bam] # 标记 duplicate
$ sambamba view -h -t 4 -f bam -F "[XS] == null and not unmapped and not duplicate" [input_markdup_bam] > [output_filter_bam] # 过滤及去重

slurm 调度系统下的批量运行示例：

#!/bin/bash
#SBATCH -J sam2bam
#SBATCH -o sam2bam.out
#SBATCH -e sam2bam.err
#SBATCH -n 24

# Define directories
input_dir="./bowtie2_alignment"
output_dir="./bamfiles"

# Create output directory
mkdir -p "${output_dir}"

# Record versions
samtools --version > samtools.version 2>&1
sambamba --version > sambamba.version 2>&1

# Function to process each SAM file into BAM format
process_sam_to_bam() {
    local sam_file="$1"
    local base_name="${output_dir}/$(basename "${sam_file}" .sam)"
    local unsorted_bam="${base_name}.unsorted.bam"
    local sorted_bam="${base_name}.sorted.bam"
    local marked_bam="${base_name}.marked.bam"
    local final_bam="${base_name}.bam"

    # Convert SAM to BAM
    samtools view -h -S -b -o "${unsorted_bam}" "${sam_file}"

    # Sort BAM file
    sambamba sort -t 4 -o "${sorted_bam}" "${unsorted_bam}"
    rm "${unsorted_bam}"
    
    # Mark duplicate
    sambamba markdup -t 4 "${sorted_bam}" "${marked_bam}"
    rm "${sorted_bam}"

    # Filter BAM file
    sambamba view -h -t 4 -f bam -F "[XS] == null and not unmapped and not duplicate" "${marked_bam}" > "${final_bam}"
    rm "${marked_bam}"
}

# Export function and variables for parallel execution
export -f process_sam_to_bam
export input_dir
export output_dir

# Execute processing in parallel
parallel --jobs 6 process_sam_to_bam ::: "${input_dir}"/*.sam

以上注释均通过 ChatGPT 添加，请确保可以通过环境变量直接调用 samtools 和 sambamba。该脚本主要做了以下工作：

samtools 将 sam 转换为 bam 格式。
sambamba 标记重复 reads。
sambamba 对文件进行排序和索引，此后对其进行去重。

以下是过滤规则：

[XS] == null：
- 过滤掉具有 XS 标签的比对，选出比对位置唯一或最佳的比对（primary alignments），即没有其他次优的比对位置。
not unmapped：
- 过滤未成功比对到参考基因组的序列。
not duplicate：
- 过滤掉标记为重复的比对。在使用 PCR 方法扩增样本时，会产生重复的 reads。

通过上述命令，就能得到所有后续分析中 macs3 所需要使用到的 bam 文件。如果想要将不同 replicate 的数据放在一起获得最终 bedGraph 文件，可在该步后对同一处理下所有 bam 文件进行合并。

Filtering out Blacklisted Regions

详细可见：The ENCODE Blacklist: Identification of Problematic Regions of the Genome

基因组中部分区域的序列高度重复，此外现有基因组版本也可能存在基因组组装错误所导致的 artifact。这些因素都可能导致 signal 出现异常，从而增加假阳性。解决这一问题的方法之一是过滤掉这些区域中的 alignment。

一些模式物种的 blacklist 文件可以在以下 github repositories 获得：

https://github.com/Boyle-Lab/Blacklist/

https://github.com/dozmorovlab/excluderanges

请注意选择正确的基因组版本，如果目标物种的 blacklist 文件并不存在于以上 repositories 中，可以自行通过 Encode Blacklist 软件识别。

以下是使用 deeptools alignmentSieve 过滤 blacklist 的方法：

1	$ alignmentSieve --blackListFileName blacklist.bed.gz --filterMetrics filterMetricBlacklist.txt -p 8 -b input.bam -o output.bam

--blackListFileName 指定 blacklist region file。
-p 指定使用线程数，-b 和 -o 指定输入和输出 bam 文件。

如果过滤的为 ATAC-seq 数据，则还需添加 --ATACshift 参数。在 filterMetricBlacklist.txt 中可以看到总 reads 数和过滤后留存的 reads 数。

由于该步骤并非必需，因此个人看需求选择是否加入到 排序及去重 部分的 Pipeline 中。

macs3 call peak

一般而言，CHIP-seq（染色质免疫沉淀测序）实验中，通常包括 treat 组和 control 组，后者用于提供实验的背景信号，帮助区分特异性的信号和非特异性的背景噪音。通常有两种主要的对照方式：

Input control：在这种对照中，不进行抗体沉淀，直接使用同一样本的一部分进行 DNA 提取和测序。
Mock IP control：这种对照使用不特异的抗体（如与目标蛋白无关的抗体或前清液）来进行免疫沉淀。

因此，假设目前已经有了 treat 和 control 的 bam files，那么可以通过以下命令来进行 callpeak：

1	$ macs3 callpeak -t [treat bam] -c [control bam] -f BAM -n [outname] --broad -g [species] -B --broad-cutoff 0.1

以上命令适用于 broad peak calling on Histone Mark ChIP-seq，如果是其他 ChIP-seq 数据，请根据 github 上的 usage 示例进行相应修改。以 TF ChIP-seq 为例，示例命令如下：

1	$ macs3 callpeak -t [treat bam] -c [control bam] -f BAM -g [species] -n [outname] -B -q 0.01

由于 macs3 的 bdgdiff 只能处理一个重复下的情况，因此如果单个样本有多个重复，在 callpeak 步骤中可以一起输入（即 -t 和 -c 可以接多个 bam 文件）。或者也可以选择每个重复单独进行 callpeak，然后使用其他工具例如 DESeq2 或 edgeR 进行后续的比较分析。

这里将对应参数后的输入依据自己的实际情况和需求进行修改，着重提几个参数：

-n 输出文件的前缀（可包括输出路径）。
-g 有效基因组大小，目前 macs3 内置有：

hs for human (2,913,022,398)

mm for mouse (2,652,783,500)

ce for C. elegans (100,286,401)

dm for fruitfly (142,573,017)

Default:hs

-B 生成 bedgraph 文件，该文件可用于后续计算得到 signal p-value 或 fold change over control 的 bigWig 文件。
-f 指定输入的文件类型，该选项可留空以让 macs3 自行检查，不过如果输入的文件为双端配对数据则需自行指定（BAMPE 或 BEDPE）。

如果分析的物种有效基因组大小在 macs3 中没有内置，且在 deeptools 上也没给出，可以通过以下方法进行计算：

下载 faCount：wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faCount
统计物种非 N 碱基数量：faCount genome.fa -summary > genome.faCount.report
计算：awk 'NR==2 {print $3+$4+$5+$6}' genome.faCount.report

将 faCount 后的 genome.fa 替换为自己的基因组文件后，就能得到相应的信息，例如 rheMac10 基因组的非 N 碱基数量为 2936892725。此外也可以 uniquely mappable regions 确定。关于这两种方法的区别，可见 Effective Genome Size。

由于 macs3 的默认值均根据非 N 碱基数量确定，因此这里只讨论 faCount。

确定相关参数后，macs3 callpeak 将产生相关的 peak 信息文件，若仅需要 peak 信息则进行到该步即完成任务。不添加 -B 参数时会生成的一些文件：

*_peaks.narrowPeak：包含了每个 peak 的详细信息，包括峰的起始和终止位置、峰的名称、峰的得分、以及峰顶相对于峰起始位置的偏移值。
*_summits.bed：包含了每个 peak 的峰顶位置，即信号最强的点。每一行通常表示一个峰顶，包含染色体、峰顶的起始位置、终止位置（起始位置+1），以及峰的强度值。
*_peaks.xls：详细结果文件，包含产生该文件所使用的命令，各参数的情况，以及峰的各个信息等。

Signal P-value & Fold change over control

如果在 call peak 步骤中添加了 -B 参数，那么此时应该就已经生成了 treat 组和 control 组的 bed 文件，通过比较这两个 bed 文件，可以得到基因组上不同位置的 signal p-value 或者 `fold change over control，该结果适用于分析和比较整个基因组的信号（而不仅限于 peak 区域）。

以下是一个运行示例，该例子将产生 Signal P-value 文件，其中第四列为 -log10(pvalue)，因此第四列越大，该处的 signal 越显著：

1	$ macs3 bdgcmp -t treatment.bedGraph -c control.bedGraph -m ppois -p 1.0 -S 1.0 -o output.bedGraph

各参数含义详细可见 documentation，这里提几个重要的：

-m 计算第四列 score 时所用的方法，ppois 即泊松 p 值，选择 FE 则能得到 fold change over control 信息。
-p 为每个计数添加 x 个伪计数，尤其在对数转换时可以避免计数为 0 造成的计算问题，此外也有助于减少因样本稀疏导致的统计波动。
-S 设置处理组（treatment）和对照组（control）数据的缩放因子，如果在 call peak 时未指定 --SPMR 参数，则该项设定为 1 即可。

对于输出后的 bedgraph 文件，可以通过 bedGraphToBigWig 转为 bigWig 文件：

1 2	$ wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/bedGraphToBigWig # 下载工具 $ bedGraphToBigWig output.bedGraph genome.chrom.sizes output.bw

每个物种的 genome.chrom.sizes 文件可在相应的 UCSC Genome Browser goldenpath 中得到，如果所分析物种尚未在基因组浏览器上记载，可自行通过 samtools 或 Biopython 得到。例：

1 2	samtools faidx reference.fasta cut -f1,2 reference.fasta.fai > chrom.sizes

如果此处存在一个样本下有多个重复的情况，并且已经对每个重复都单独计算得到了 bed 文件，你也可以通过 macs3 的 cmbreps 对不同的重复进行合并。例如以 signal p-value 为例：

1	$ macs3 cmbreps -i replicate1.bedGraph replicate2.bedGraph replicate3.bedGraph -o combined.bedGraph --method fisher

其通过费舍尔联合概率检验合并 p 值并计算得到新 p 值。

对于 fold change over control，可将 -m 后的 fisher 换为 mean 进行。

后记

以上流程已被打包为一个完整的 Snakemake Pipeline，可在以下链接中获取：

https://github.com/JuseTiZ/Blog-Snakemake-pipeline/tree/main/ChIP-Pipeline

该 Pipeline 的理解和使用可能需要一定的 Snakemake 语法基础。

考虑到 ChIP-seq 下游分析能弄的东西五花八门，所以这篇文章就只介绍到如何产生结果，对于相关结果的解释和分析就各看需求了。