Bioinformatics training course: epigenomics: ChIP-seq数据的Peak Calling

本节课我们主要学习如何使用MACS2来对ChIP-seq数据进行Peak Calling。为了让大家更全面了解该软件的使用方法，以及能够处理历史遗留数据，本课件将以单端数据为例，讲解数据分析的过程。

基本原理

MACS2（Model-based Analysis for ChIP-Seq, version 2）是常用的ChIP-seq数据分析工具，用于识别基因组中蛋白质-DNA相互作用的结合位点（即peaks）。它通过一系列统计模型和算法来区分真正的结合位点与背景噪音。MACS2的基本原理包括以下几个步骤：

模型建立： MACS2基于Poisson分布模型来模拟背景噪音。它假设在没有真正信号的区域，测序reads是随机分布的。因此，结合位点区域相对背景的显著性可以通过偏离Poisson分布的程度来衡量。
数据处理：由于单端ChIP-seq测序只捕获到DNA片段的一小部分，因此MACS2根据实验设定的片段长度（fragment size）对reads进行延伸，以重构完整的DNA片段。
peak calling：在有输入对照（input control）的情况下，MACS2会将ChIP样品与对照进行比较，减少假阳性信号。通过对比ChIP信号与背景噪音的差异，找出显著的结合区域。MACS2使用Poisson模型计算每个潜在peak相对于背景的p值，反映该区域reads富集的显著性。一般情况下，p值越低，peak越有可能是真实的结合位点。
q值校正：由于进行多个假设检验，MACS2会进行FDR（False Discovery Rate）校正，以控制假阳性率，并输出经过校正的q值。用户可以根据q值来过滤不显著的peaks。
输出结果：最终，MACS2输出peak的位置信息、p值、q值和峰值强度等信息，用于后续的生物学分析。

估计片段长度

[参考文献](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0011471)

Figure 1: 参考文献

如上图所示，MACS2通过分析 forward tag和reverse tag的位置关系，估计DNA片段的长度。它的具体做法如下：

统计位移距离：MACS2会分析forward tag和reverse tag的分布情况，计算两者之间的距离。通常情况下，forward tag和reverse tag在蛋白质结合位点的两侧，MACS2通过计算forward tag到reverse tag的距离来估计测序片段的长度。
最大相关性距离：MACS2 会通过将 reverse tags 向 5’ 方向（即向 forward tags 靠近的方向）逐步移动，计算不同位移距离下 forward tags 和 reverse tags 的重叠程度。当两组reads的重叠达到最大时，所对应的位移距离就被认为是片段的长度。

[参考文献](https://genome.cshlp.org/content/22/9/1813.full)

Figure 2: 参考文献

估计背景参数

假设一个窗口中期望的Reads是\(\lambda\)，那么观测到不同数目的Reads的可能性：

[参考文献](https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html)

Figure 3: 参考文献

在MACS2中，背景\(\lambda\)参数用于估计ChIP-seq数据中的背景噪音水平，以帮助区分真正的信号峰值（peaks）和随机噪音。背景lambda的估计过程旨在提供对基因组区域中期望的背景reads分布的合理估计，并基于此来计算每个潜在结合区域的富集程度。当ChIP-seq实验提供了输入对照（input control），通过在不同尺度上（1 kb、5kb，10 kb，基因组）计算背景reads水平，并选取最大值最作为\(\lambda\)的估计。

Figure 4: 参考文献

测试数据

我们从SRA下载了GSM3112103（PU.1 ChIP-seq）和GSM3112107（Input）数据，物种为人，单端测序。先进行数据的预处理：去除接头序列、序列比对。

ln -s /sibcb1/bioinformatics/biocourse/biocourse01/epigenomics/FASTQ .
ln -s /sibcb1/bioinformatics/biocourse/biocourse01/epigenomics/trimmed .
ln -s /sibcb1/bioinformatics/biocourse/biocourse01/epigenomics/output .
ln -s /sibcb1/bioinformatics/biocourse/biocourse01/dnaseqAlignment/reference/GRCh38.fasta .

ls -la FASTQ | grep fastq.gz

-rw-r--r-- 1 biocourse01 ods   529375283 Nov 13  2023 GSM3112103_PU1.fastq.gz
-rw-r--r-- 1 biocourse01 ods   538523573 Nov 13  2023 GSM3112107_input.fastq.gz

去除接头序列

for i in $(ls FASTQ/*.fastq.gz | awk -F'/' '{print $NF}' | awk -F'.' '{print $1}'); do
  echo $i
  fastqc  -o FASTQ FASTQ/$i.fastq.gz
  /sibcb1/bioinformatics/training_shared/software/fastp \
    -i FASTQ/$i.fastq.gz \
    -o trimmed/$i.fastq.gz \
    -l 25 -w 12 -t 4 \
    --adapter_fasta FASTQ/illumina_multiplex.fa \
    --json trimmed/$i.json \
    --html trimmed/$i.html \
    2>trimmed/$i.fastp.log
    fastqc -o trimmed trimmed/$i.fastq.gz
done

比对

sample='GSM3112103_PU1 GSM3112107_input'

for sam in $sample; do
        bwa mem -t 16 \
        /sibcb1/bioinformatics/biocourse/biocourse01/dnaseqAlignment/reference/GRCh38 \
        trimmed/$sam.fastq.gz \
        2>output/$sam.log | samtools view -bS - -o output/$sam\_bwa.bam
        samtools sort -@ 4 -o output/$sam.bam output/$sam\_bwa.bam
        samtools index output/$sam.bam
        rm output/$sam\_bwa.bam
done

Peak calling

命令行

使用一下的命令，鉴定q < 0.05的peaks：

macs2 callpeak -B --SPMR -t output/GSM3112103_PU1.bam -c output/GSM3112107_input.bam -g hs -n PU1 -f BAM --keep-dup 1 --outdir macs2

可以检查fragment size的估计：

source('macs2/PU1_model.r')

-g: 有效基因组大小。可以是1.0e+9或1000000000，或者使用快捷方式：‘hs’ 表示人类 (2.7e9)，‘mm’ 表示小鼠 (1.87e9)，‘ce’ 表示线虫 (9e7)，‘dm’ 表示果蝇 (1.2e8)，默认值为：hs（人类）。
--keep-dup: 如何处理重复片段。可以是只保留一个（1）,或者所有（all），或者自动判定（auto，对高覆盖度数据推荐使用）。注意：MACS2能够自动判定那些Reads是重复序列，不需要预先使用samtools/picard定义。
--SPMR: 配合-B使用，计算’fragment pileup per million reads’。

结果

在本例子中，PU1_peaks.xls中可以看到产生的peak的区间和显著性信息，包括：-log10(pvalue)、Fold enrichment、-log10(qvalue)等。在该文件中，我们还可以看到每个样本的library size：tags after filtering in treatment: xxxx。当我们使用macs2 bdgdiff的时候，需要用到library size。阅读此参考资料获取更多信息。

注释到基因组

我们可以使用ChIPseeker方便地将富集区间注释到基因组上：

suppressPackageStartupMessages(library(ChIPseeker))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene))
suppressPackageStartupMessages(library(org.Hs.eg.db))

peak <- readPeakFile('macs2/PU1_peaks.narrowPeak')
anno <- annotatePeak(peak, TxDb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb='org.Hs.eg.db', verbose=F)
out <- as.data.frame(anno)
head(anno@anno)

GRanges object with 6 ranges and 19 metadata columns:
      seqnames      ranges strand |          V4        V5          V6
         <Rle>   <IRanges>  <Rle> | <character> <integer> <character>
  [1]     chr1 21210-21387      * |  PU1_peak_1        42           .
  [2]     chr1 26044-26228      * |  PU1_peak_2        47           .
  [3]     chr1 28966-29386      * |  PU1_peak_3        42           .
  [4]     chr1 33070-33346      * |  PU1_peak_4        90           .
  [5]     chr1 75648-75899      * |  PU1_peak_5        93           .
  [6]     chr1 96395-96936      * |  PU1_peak_6       510           .
             V7        V8        V9       V10             annotation
      <numeric> <numeric> <numeric> <integer>            <character>
  [1]   4.87062   6.38964   4.29801        67 Intron (ENST00000488..
  [2]   4.81400   6.91255   4.78891        59 Intron (ENST00000488..
  [3]   4.87062   6.38964   4.29801       342       Promoter (<=1kb)
  [4]   7.37812  11.36856   9.03715       151       Promoter (2-3kb)
  [5]   7.19848  11.66750   9.31923       109      Distal Intergenic
  [6]  22.94265  54.12161  51.00166       310 Intron (ENST00000466..
        geneChr geneStart   geneEnd geneLength geneStrand      geneId
      <integer> <integer> <integer>  <integer>  <integer> <character>
  [1]         1     17369     17436         68          2   102466751
  [2]         1     29554     31097       1544          1   107985730
  [3]         1     29554     31097       1544          1   107985730
  [4]         1     30366     30503        138          1   100302278
  [5]         1     65419     71585       6167          1       79501
  [6]         1     65419     71585       6167          1       79501
           transcriptId distanceToTSS         ENSEMBL      SYMBOL
            <character>     <numeric>     <character> <character>
  [1] ENST00000619216.1         -3774 ENSG00000278267   MIR6859-1
  [2] ENST00000473358.1         -3326            <NA> MIR1302-2HG
  [3] ENST00000473358.1          -168            <NA> MIR1302-2HG
  [4] ENST00000607096.1          2704 ENSG00000284332   MIR1302-2
  [5] ENST00000641515.2         10229 ENSG00000186092       OR4F5
  [6] ENST00000641515.2         30976 ENSG00000186092       OR4F5
                    GENENAME
                 <character>
  [1]        microRNA 6859-1
  [2]    MIR1302-2 host gene
  [3]    MIR1302-2 host gene
  [4]        microRNA 1302-2
  [5] olfactory receptor f..
  [6] olfactory receptor f..
  -------
  seqinfo: 87 sequences from hg38 genome

plotAnnoPie(anno)

Motif富集分析

# 添加路径
export PATH="/sibcb/program/src/homer/bin:$PATH"

# 为了演示的方法，我们这里只提取前1000个区间进行分析，在实际操作中，请分析整个文件
head -1000  macs2/PU1_peaks.narrowPeak > macs2/PU1_peaks_1000.narrowP
findMotifsGenome.pl macs2/PU1_peaks_1000.narrowPeak GRCh38.fasta homer2 -mask -p 4

findMotifsGenome.pl同时进行de novo motif和known motif的分析。注意：该命令需要读取参考序列，并在同一目录下建立文件夹，所以确保你有写的权限。

从已知motif的富集结果看，PU.1的motif排在第一。排在前五的motif都是该motif的冗余形式。

从de novo motif的结果来看，排在第一motif也是正确的结果。在实际应用中，我们可以报告两个分析结果。

数据展示

信号文件

macs2 callpeak产生了两个bedGraph文件，在设置--SPMR参数后，这两个信号文件是可比的。还可以进一步，得到Fold-Enrichment的track：

# 染色体长度
samtools faidx GRCh38.fasta
cut -f1,2 GRCh38.fasta.fai > GRCh38.length

# ChIP信号
cat macs2/PU1_treat_pileup.bdg | sort -k1,1 -k2,2n > macs2/PU1_treat_pileup_sorted.bdg
bedGraphToBigWig macs2/PU1_treat_pileup_sorted.bdg GRCh38.length macs2/PU1_treat.bw

# control信号
cat macs2/PU1_control_lambda.bdg | sort -k1,1 -k2,2n > macs2/PU1_control_lambda_sorted.bdg
bedGraphToBigWig macs2/PU1_control_lambda_sorted.bdg GRCh38.length macs2/PU1_control_lambda.bw

# Fold-enrichment track
macs2 bdgcmp -t macs2/PU1_treat_pileup.bdg -c macs2/PU1_control_lambda.bdg -o macs2/PU1_FE.bdg -m FE
cat macs2/PU1_FE.bdg | sort -k1,1 -k2,2n > macs2/PU1_FE_sorted.bdg
bedGraphToBigWig macs2/PU1_FE_sorted.bdg GRCh38.length macs2/PU1_FE.bw

信号分布

# 加载软件环境
source /sibcb1/bioinformatics/training_shared/software/python3.6_profile

computeMatrix scale-regions -p 5 --binSize 50 -b 2000 -a 2000 -m 6000 -S macs2/PU1_FE.bw \
  -R /sibcb1/bioinformatics/biocourse/biocourse01/bulkRNAseq/idx/gencode.v34.primary_assembly.annotation.gtf \
  -o macs2/PU1_global.mat.gz

plotProfile -m macs2/PU1_global.mat.gz --samplesLabel PU.1 -o macs2/PU1_global.mat.pdf

IGV

我们可以把ChIP-seq peaks、信号、motif同时在IGV上展示。前两种文件前面已经获得，这里准备motif instance文件：

cat homer2/knownResults/known1.motif

>AGAGGAAGTG PU.1(ETS)/ThioMac-PU.1-ChIP-Seq(GSE21512)/Homer 7.613173    -651.698381 0   T:432.0(43.20%),B:2278.5(4.89%),P:1e-283
0.643   0.001   0.149   0.207
0.122   0.171   0.705   0.002
0.830   0.012   0.157   0.001
0.001   0.001   0.997   0.001
0.001   0.001   0.997   0.001
0.997   0.001   0.001   0.001
0.990   0.001   0.001   0.008
0.024   0.074   0.901   0.001
0.001   0.005   0.001   0.993
0.079   0.097   0.773   0.051

# 扫描Peak中的TFBS
findMotifsGenome.pl macs2/PU1_peaks_1000.narrowPeak GRCh38.fasta homer2 -mask -p 4 \
  -find homer2/knownResults/known1.motif > homer2/Instance_motif1.txt

# 扫描全基因组的TFBS
scanMotifGenomeWide.pl homer2/knownResults/known1.motif GRCh38.fasta -p 16 > homer2/Instance_motif1_GW.txt
awk '{print $2, $9-1, $10}' OFS='\t' homer2/Instance_motif1_GW.txt > homer2/Instance_motif1_GW.bed

使用IGV展示：

结束语

本课程主要学些了单端ChIP-seq数据的分析方法，重点需要掌握MACS2 peak calling的方法。在ATAC-seq数据分析部分，我们将学习MACS2如何处理双端数据。