6.3.3 实验分辨率
-
Peak caller 通常将邻近的
Peak合并到较大的区域中，这可能导致分辨率损失。因此，将两种类型的区域分离开来是很重要的：Broad与sharp区域。一方面，来自组蛋白 ChIP-seq 的数据显示信号的广泛富集，并且需要大区域的识别（ 示例代码见章节 6.5）。因此，组蛋白修饰类型的Peak finder应将重点放在确定区域边界上，而不是确定Peak summits。 这是一项非平凡的任务，仍然是一项挑战，目前的Peak caller仍然难以在样本之间找到一致的组蛋白结合区域边界。这也可能是由于这种类型的信号的性质。 - 另一方面，来自于转录因子类型的
Sharp信号类型富集，需要鉴定narrow区域的分析方法（示例代码见章节 6.4）。由于调节区由多个TF结合位点组成（ref8、ref9）， ChIP-seq 数据的分辨率对于区分它们是至关重要的。在最近的方法改进中，如 ChIP-exo （ref10）和 ChIP-nexus （ref11），实验分辨率增加到几个碱基对。因此，应用于TFChIP-seq 数据的Peak caller应适用于这些没有control的新方法，并通过聚焦于Peak summits而不是整个Peak regions来利用实验分辨率。要确定哪种方法对TF数据最有效，最好的指标是Peak summits之间的最小距离，而不是Peak regions的大小（ref12）。
6.3.4 新一代的 Peak caller
- 最近开发的方法确实考虑了上面讨论的一些细节。MACS 的升级版
MACS2，既能鉴定Broad Peaks有你鉴定Sharp Peaks（ 示例代码见章节 6.5）。Peakzilla（ref12）是专门开发用于在高分辨率下从转录因子 ChIP-seq 数据中鉴定Peak（ 示例代码见章节 6.4）。HOMER （ref13）, 最初设计为在Peak区域重新识别motif的工具（ 示例代码见章节 9）。也有其他类型的工具，如 findPeaks，用于鉴定不同类型的 ChIP-seq 的Peak。JAMM（ref14）使用重复的样本来提高Peak宽度的分辨率和精度。GEM（ref15）通过包括关于TF motif的信息作为附加输入来识别高分辨率的Peak，然而，我们更喜欢使用Motif信息来验证Peak，而不是在识别步骤。已经专门开发了其他方法来比较不同的样本，在 章节 8 将对此进行讨论。
6.3.5 Post-processing
- 可以
Post-processing处理步骤来去除在Peak calling过程中不能过滤的人为造成的Peak。这样的Peak出现在基因组blacklisted区域中，这些区域在任何ChIP-seq数据集（ ChIP 和 Input ）中显示非常高的 reads 丰度，通常位于着丝粒和端粒，并已被ENCODE定义为几个物种（ 见章节 2.2.2）。我们还选择去除位于线粒体染色体（ chrM ）上的Peak。
6.4 Peakzilla：转录因子类型数据
-
Peakzilla 专为
Sharp的 TF ChIP-seq 数据而设计，以便在高分辨率下call peak，即解析由 TF homotypic 结合产生的紧密间隔的Peak summits。Peakzilla 可以使用于没有对照control的ChIP-exo数据。简而言之，它使用来自高度富集区域的正向和反向reads 的双重分布来估计片段大小。然后，它使用滑动窗口直接扫描沿基因组的 reads 的双重分布来对区域进行评分。为此，它首先计算 ChIP 样本中的 reads 数减去 control 样本中的 reads 数（ 通过文库中比对的 reads 总数进行归一化 ）。此方法允许比变化倍数更改更好的排序。然后，它用p-value对这个原始分数进行加权，p-value值检查数据与预期的正向和反向 reads 的双高斯分布的匹配程度。这允许在不去除重复 reads 的情况下用 PCR 人工产物过滤掉位置，以及更好地识别准确的Peak summits位置。最后，通过交换 ChIP 和 control 样本计算经验 FDR，计算每个区域的富集倍数，并通过多重检验进行校正。默认情况下，它返回最小得分为 1 的峰值和 2 的富集倍数。唯一需要的输入是两个 ChIP 和 control 文件。# Create directory for peaks mkdir peaks sample=NRF1_CHIP_WT_1 control=NRF1_INPUT_WT ## Break down of the commands (many intermediate files) # Run Peakzilla bamToBed -i reads/${sample}.bam > reads/${sample}.bed bamToBed -i reads/${control}.bam > reads/${control}.bed peakzilla.py reads/${sample}.bed reads/${control}.bed -l peaks/${sample}_peakzilla_report.txt > peaks/${sample}_peakzilla.tsv # Removal of peaks on chrM (density is usually abnormally high along the complete chrM) # Reorganisation of output into BED-like format: chromosome / start (0-based) / end / summit / score / fold_enrichment awk -v OFS="\t" ’(NR>1&&$1!="chrM"){print $1,$2-1,$3,$5,$6,$9}’ peaks/${sample}_peakzilla.tsv > peaks/${sample}_peaks_peakzilla_tmp.bed # Removal of peaks in blacklisted regions intersectBed -v -a peaks/${sample}_peaks_peakzilla_tmp.bed -b mm10_blacklist.bed.gz > peaks/${sample}_peakzilla.bed # Remove intermediate files rm reads/${sample}.bed reads/${control}.bed peaks/${sample}_peakzilla.tsv peaks/${sample}_peaks_peakzilla_tmp.bed ## 一行命令 (˜6 min) peakzilla.py < (bamToBed -i reads/${sample}.bam) < (bamToBed -i reads/${control}.bam) -l peaks/${sample}_peakzilla_report.txt | awk -v OFS="\t" ’(NR>1&&$1!="chrM"){print $1,$2-1,$3,$5,$6,$9}’ | intersectBed -v -a stdin -b mm10_blacklist.bed.gz > peaks/${sample}_peaks_peakzilla.bed
6.5 MACS2：组蛋白修饰类型数据
-
MACS是最早用来鉴定 ChIP-seq 数据的Peaks软件之一，并且仍然是最受欢迎的一种。最初，它被开发用来识别sharp的Peak，但定义的Peak相对较宽。最新版本MACS2现在两者都可以鉴定。简而言之，它首先删除所有重复的 reads 。然后，使用来自高度富集区域的正向和反向 reads 的双重分布来估计片段大小。它将所有reads 扩展到估计的片段大小，并将 ChIP 和 Input 样本扩展到相同的测序深度，即标准化为文库中比对上的 reads 的总数。然后，它使用滑动窗口扫描沿基因组的片段分布，对区域进行评分。为此，它将 ChIP 中的富集程度与对照样本进行比较，并使用局部泊松分布计算显著性得分。最后，它使用 Benjamini-Hochberg 过程对多重检验进行校正。默认情况下，它返回最小q-value得分为 0.05 的Peak。输入参数包括 ChIP 和 Input 文件的路径、文件格式、基因组大小、输出目录和样本名称。--broad是用来鉴定broad peaks。# Run MACS2 (˜7 min) # q-value threshold can be adapted with the option --broad-cutoff (default 0.1) macs2 callpeak -t reads/${sample}.bam -c reads/${control}.bam -f BAM -g mm --outdir peaks -n ${sample} --broad 2> peaks/${sample}_macs2_report.txt # Reorganisation of output into bed-like format: chromosome / start (0-based) / end / summit / score / fold_enrichment cat peaks/${sample}_peaks.xls | grep -v "#" | awk -v OFS="\t" '(NR>2&&$1!="chrM"){print $1,$2,$3,$4,$6,$7}' | intersectBed -v -a stdin -b mm10_blacklist.bed.gz > peaks/${sample}_peaks_macs.bed # Removal of intermediate files rm peaks/${sample}_peaks.xls peaks/${sample}_model.r peaks/${sample}_peaks.broadPeak peaks/${sample}_peaks.gappedPeak
6.6 饱和度分析
-
为了检查样品的测序深度是否足以识别大多数结合区域，建议进行饱和分析。它涉及对 reads 的数量进行随机二次取样，并计算使用这些子集识别的峰值数量。然后根据使用的 reads 数绘制峰值数量。如果峰的数量显示饱和并达到一个平台，那么样品的测序足够深。NRF1 ChIP-seq 样本在
2000万 reads时显示饱和（ 图 6.2 ）。图 6.2 NRF1 ChIP-seq 饱和度分析
sample=NRF1_CHIP_WT_1
control=NRF1_INPUT_WT
## Number of peaks found
wc -l peaks/${sample}_peaks_peakzilla.bed
## Saturation analysis
# Does the number of peaks saturate when using all available reads
NB=$(samtools idxstats reads/${sample}.bam | awk '{t=t+$3}END{print t}')
L=$sample
for subset in ‘seq 1000000 1000000 $NB‘;
do
L=$L"\t"$(peakzilla.py < (bamToBed -i reads/${sample}.bam | shuf -n $subset
--random-source=reads/${sample}.bam) < (bamToBed -i
reads/${control}.bam) -l /dev/null | wc -l)
done
echo -e $L > peaks/${sample}_peakzilla_saturation_table.txt
# Open R
R
# Load saturation table
d = read.table(paste("peaks/",sample,
"_peakzilla_saturation_table.txt",sep=""),row.names=1)
# Generate the saturation plot
pdf(paste("peaks/",sample,"_peakzilla_saturation_plot.pdf", sep=""))
par(bg="white")
plot(seq(1,ncol(d),1),
d[1,],
type="l",
xlab="Number of reads in million",
ylab="Number of peakscalled",
main=c("Saturation analysis",
rownames(d)[1],
paste("Total",d[1,ncol(d)],"peaks")))
dev.off()
# Close R
q()
ref1：Practical guidelines for the comprehensive analysis of ChIP-seq data.
ref2：How does multiple testing correction work?
ref3：Evaluation of algorithm performance in ChIP-seq peak detection
ref4：Model-based analysis of ChIP- Seq (MACS).
ref5：Design and analysis of ChIP-seq experiments for DNA-binding proteins.
ref6：Systematic evaluation of factors influencing ChIP-seq fidelity
ref7：Genome-wide mapping of in vivo protein-DNA interactions.
ref10：Comprehensive genomewide protein-DNA interactions detected at single-nucleotide resolution.
ref11：ChIP-nexus enables improved detection of in vivo transcription factor binding footprints
ref12：Identification of transcription factor binding sites from ChIP-seq data at high resolution
ref14：JAMM: a peak finder for joint analysis of NGS replicates
网友评论