前言

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

MACS3 github:

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

MACS3 documentation:

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

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

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

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

  • 数据的质控
  • 除 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 调度系统下的一个并行示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#!/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 时的假阳性。

单个运行时的示例:

1
2
3
4
$ 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 调度系统下的批量运行示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#!/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 添加,请确保可以通过环境变量直接调用 samtoolssambamba。该脚本主要做了以下工作:

  • samtoolssam 转换为 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,然后使用其他工具例如 DESeq2edgeR 进行后续的比较分析。

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

  • -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 文件。

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

  • 下载 faCountwget 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 文件,可以得到基因组上不同位置中 peak 的 signal p-value 或者 fold change over control

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

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 中得到,如果所分析物种尚未在基因组浏览器上记载,可自行通过 samtoolsBiopython 得到。例:

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 下游分析能弄的东西五花八门,所以这篇文章就只介绍到如何产生结果,对于相关结果的解释和分析就各看需求了。