ChIP-Seq 分析流程-下游 (1)

写在前面的话：

参考使用的文件资料是由哈佛陈生物信息学核心 (HBC) 教学团队成员开发的。另外也看了多个公众号文章和书籍。

Website: https://hbctraining.github.io/Intro-to-ChIPseq/

Github: https://github.com/hbctraining/Intro-to-ChIPseq

接上游分析结果，这次我们继续探索下游流程。

上游分析：https://mp.weixin.qq.com/s/0KOk_ij29xrqdQzNaeVVvA

$ tree
.
|-- Nanog-rep1-macs2.log
|-- Nanog-rep1_model.r
|-- Nanog-rep1_peaks.narrowPeak
|-- Nanog-rep1_peaks.xls
|-- Nanog-rep1_summits.bed
|-- Nanog-rep2-macs2.log
|-- Nanog-rep2_model.r
|-- Nanog-rep2_peaks.narrowPeak
|-- Nanog-rep2_peaks.xls
|-- Nanog-rep2_summits.bed
|-- Pou5f1-rep1-macs2.log
|-- Pou5f1-rep1_model.r
|-- Pou5f1-rep1_peaks.narrowPeak
|-- Pou5f1-rep1_peaks.xls
|-- Pou5f1-rep1_summits.bed
|-- Pou5f1-rep2-macs2.log
|-- Pou5f1-rep2_model.r
|-- Pou5f1-rep2_peaks.narrowPeak
|-- Pou5f1-rep2_peaks.xls
`-- Pou5f1-rep2_summits.bed

上游分析呢，我们共得到了四组文件， Nanog-rep1 , Nanog-rep2 , Pou5f1-rep1 , Pou5f1-rep2 ，每组数据又有五个不同类型的文件。

• -macs2.log : 日志文件

• _peaks.narrowPeak ：BED6+4 格式文件，包含峰值位置以及峰值峰值、pvalue 和 qvalue

• _peaks.xls ：包含有关峰值信息的表格文件。其他信息包括堆积和折叠富集

• _summits.bed ：每个峰的峰顶位置。要找到结合位点的基序，建议使用此文件

• _model.R ：一个 R 脚本，你可以使用它根据数据和互相关图生成有关模型的 PDF 图像

创建目录结构

新的开始，先创建目录结构，将需要的数据和未来归档的数据放在该放的位置一定是一个好习惯。这一步在Rstudio中完成，Rstudio中如何创建文件和移动文件就不多说了，网上说的很详细了，这里直接展示成品号，哈哈哈哈。

.
|-- chipseq.Rproj
|-- data
|   |-- bams
|   `-- peakcalls
|-- figures
|-- meta
`-- results

创建完目录就是将上游分析结果移动过来了，另外我们还需要对我们的样本根据需求编写一个 meatdata 文件，原始信息链接在这里：https://raw.githubusercontent.com/hbctraining/Intro-to-ChIPseq/master/samplesheet_chr12.csv 我们下载过来就ok。防止链接失效就也直接粘贴过来吧。请严格按照这个格式创建信息！

SampleID,Factor,Replicate,bamReads,ControlID,bamControl,Peaks,PeakCaller,Tissue,Condition
Nanog.Rep1,Nanog,1,data/bams/H1hesc_Nanog_Rep1_aln.bam,Nanog-Input1,data/bams/H1hesc_Input_Rep1_aln.bam,data/peakcalls/Nanog-rep1_peaks.narrowPeak,narrow,NA,NA
Nanog.Rep2,Nanog,2,data/bams/H1hesc_Nanog_Rep2_aln.bam,Nanog-Input2,data/bams/H1hesc_Input_Rep2_aln.bam,data/peakcalls/Nanog-rep2_peaks.narrowPeak,narrow,NA,NA
Pou5f1.Rep1,Pou5f1,1,data/bams/H1hesc_Pou5f1_Rep1_aln.bam,Pou5f1-Input1,data/bams/H1hesc_Input_Rep1_aln.bam,data/peakcalls/Pou5f1-rep1_peaks.narrowPeak,narrow,NA,NA
Pou5f1.Rep2,Pou5f1,2,data/bams/H1hesc_Pou5f1_Rep2_aln.bam,Pou5f1-Input2,data/bams/H1hesc_Input_Rep2_aln.bam,data/peakcalls/Pou5f1-rep2_peaks.narrowPeak,narrow,NA,NA

表格每一列代表信息如下

• SampleID ：样本的标识符字符串
• Factor、Tissue,Condition ：最多三个不同因素的标识符字符串（需要列出所有列。如果您没有信息，则将值设置为 NA）
• Replicate本的重复次数
• bamReads ：包含 ChIP 样本对齐reads的 BAM 文件的文件路径
• ControlID ：对照样本的标识符字符串
• bamControl ：包含对照样本对齐reads的 bam 文件的文件路径
• Peaks ：包含样本峰值的文件路径
• PeakCaller ：所用峰值调用器的标识符字符串。可能的值包括“raw”、“bed”、“narrow”、“macs”

最后我们的文件和文件结构就是这样的

|-- chipseq.Rproj
|-- data
|   |-- bams
|   |   |-- H1hesc_Input_Rep1_aln.bam
|   |   |-- H1hesc_Input_Rep1_aln.bam.bai
|   |   |-- H1hesc_Input_Rep2_aln.bam
|   |   |-- H1hesc_Input_Rep2_aln.bam.bai
|   |   |-- H1hesc_Nanog_Rep1_aln.bam
|   |   |-- H1hesc_Nanog_Rep1_aln.bam.bai
|   |   |-- H1hesc_Nanog_Rep2_aln.bam
|   |   |-- H1hesc_Nanog_Rep2_aln.bam.bai
|   |   |-- H1hesc_Pou5f1_Rep1_aln.bam
|   |   |-- H1hesc_Pou5f1_Rep1_aln.bam.bai
|   |   |-- H1hesc_Pou5f1_Rep2_aln.bam
|   |   `-- H1hesc_Pou5f1_Rep2_aln.bam.bai
|   `-- peakcalls
|-- figures
|-- meta
|   `-- samplesheet.csv
`-- results

peak 质量评估

## Load libraries
library(ChIPQC)

## Load sample data
samples 'meta/samplesheet.csv')
View(samples)

## Create ChIPQC object
register(SerialParam())
chipObj "hg19") 

## Create ChIPQC report
ChIPQCreport(chipObj, reportName="ChIP QC report: Nanog and Pou5f1", reportFolder="ChIPQCreport")

Summary

生成的报告首先是一个表格

• ID - 样本的唯一标识符。
• Tissue/Factor/Condition - 与样本相关的元数据。
• Replicate - 样本组内的重复编号。
• Reads - 在分析的染色体内的样本reads数。
• Dup% - 通过 MapQ 过滤的 reads 中标记为重复的百分比。
• FragLen - 通过交叉覆盖方法估计的片段长度。
• SSD - SSD 分数（htSeqTools）。
• FragLenCC - 在片段长度处的交叉覆盖分数。
• RelativeCC - 在片段长度处的交叉覆盖分数与reads长度处的交叉覆盖分数的比值。
• RIP% - 峰内reads的百分比。
• RIBL% - 黑名单区域内reads的百分比。

RiP（峰内reads）

这个指标表示与“已识别峰”重叠的reads的百分比。它可以被视为一个“信噪比”度量，用于衡量文库中由结合位点的片段组成的比例与背景reads的比例。

RiP（也称为 FRiP）值会根据感兴趣的蛋白质而有所不同：

• 一个典型的高质量转录因子（具有尖锐/窄峰）在成功富集的情况下会表现出大约 5% 或更高的 RiP 值。
• 一个高质量的 RNA 聚合酶 II（具有尖锐/窄峰和分散/宽峰的混合）会表现出 30% 或更高的 RiP 值。
• 也有一些已知的高质量数据集的 RiP 值低于 1%（例如 RNA 聚合酶 III 或结合位点较少的蛋白质）。

良好富集的 ChIP 样本将有更高比例的 reads 与已识别的峰重叠。尽管 Nanog 的 RiP 较高，但 Nanog 样本的箱线图显示出与 Pou5f1 相比，不同重复之间有相当大的差异，这可能可以通过 reads 长度和测序深度的差异来解释。

SSD

对于 IP 样本，我们预期在某些区域会有reads的富集或高覆盖率，而其他区域则覆盖率较低。对于对照样本，我们预期基因组覆盖率差异较小。一个“好”的或富集的样本通常在某些区域有显著的reads堆积（覆盖率差异较大），因此更高的 SSD 更能表明更好的富集。

SSD 分数对人工高信号区域和真正的 ChIP 富集区域都很敏感。因此，高 SSD 可能是 ChIP 富集的结果，也可能是黑名单区域中某些人工高信号的结果。

在我们的数据集中，Pou5f1 重复样本的分数比 Nanog 重复样本的分数更高。这可能表明 Pou5f1 样本中有更大的富集，但我们需要仔细查看 ChIPQC 的其余输出，以确保 Pou5f1 中的高 SSD 不是由于某些未知的伪影造成的。

可以使用报告中的“覆盖直方图”来可视化 SSD 探索的覆盖均匀性。x 轴表示在碱基对位置的reads堆积高度，y 轴表示有多少位置具有这种堆积高度。该图是对数刻度。

1. 富集效果好的样本：

• 合理的“尾部” ：在y轴上有一个合理的“尾部”意味着在图表中，表示测序深度的y轴上有一些较高的值。这些较高的y轴值表示在这些位置上有更多的reads，也就是说这些位置上有更高的堆积。
• 更多的位置具有较高的测序深度 ：这意味着在基因组的多个位置上，测序深度较高。这是富集效果好的一个标志，因为它表明在这些位置上有更多的目标蛋白质结合。

• 主要由背景reads组成 ：这类样本中的reads主要是背景噪音，而不是目标蛋白质结合的reads。
• 大多数位置堆积较低 ：在这种情况下，基因组中的大多数位置（y轴上的高值）将会有较低的堆积（x轴上的低值）。这意味着在这些位置上，reads的数量较少，测序深度较低。

在我们的数据集中，Nanog 样本相比 Pou5f1 样本有相当重的尾部，特别是重复样本 2。“重尾”指的是曲线比指数曲线更重，曲线下面积更大。这表明 Nanog 样本在基因组中有更多位置具有更高的测序深度。然而，SSD 分数对 Pou5f1 来说更高。当 SSD （平方差之和） 很高但覆盖率看起来较低时，这可能是由于存在大面积高深度区域，且可能是基因组区域被列入黑名单的标志。

RiBL（reads 黑名单区域中的重叠部分）

1. RiBL 分数 ：

• 定义 ：RiBL（Read in Blacklisted regions）分数表示reads重叠已知具有高信号的区域的百分比。
• 重要性 ：较低的 RiBL 分数表示样本中较少的reads落在这些已知的高信号区域，这通常是更好的，因为这些区域的信号可能是异常的或非特异性的。

• 特性 ：这些区域通常显示为唯一可映射区域，这意味着传统的可映射性过滤器无法有效去除它们。
• 常见位置 ：这些区域通常位于特定类型的重复序列中，如着丝粒、端粒和卫星重复序列。

• RiBL 分数的作用 ：RiBL 分数提供了样本中背景信号水平的一个指标。尽管这些黑名单区域仅占基因组的 0.5%，但它们可能贡献了超过 10% 的总信号。
• 干扰分析 ：来自这些黑名单区域的信号已被证明会干扰峰值调用和片段长度估计。因此，跟踪并过滤映射到这些区域的reads是必要的，以确保数据分析的准确性。

Strand cross-correlation （链互相关）

交叉相关图通常会产生两个峰：一个对应主要片段长度的富集峰（最高相关值），另一个对应读取长度的峰（“幻影”峰）。其实这里不理解也没关系，主要涉及测序相关知识，后续可补充说明吧。

1. 要片段长度峰 ：这个峰值对应于 ChIP-seq 实验中实际富集的 DNA 片段的长度，反映了真实的生物学信号。
2. 幻影峰 ：这个峰值对应于读取长度，并不是由于真实的 DNA 片段富集引起的，而是由于读取方向的对称性导致的假象信号。这个峰值通常出现在读取长度的位置上。

以下是我们的结果出图：

这里也给了三种示例：

强信号

高质量的 ChIP-seq 数据集通常会有一个比读取长度峰更大的片段长度峰。下面的例子展示了在人类细胞中使用 CTCF（锌指转录因子）数据的强信号。使用好的抗体，转录因子通常会产生 45,000 到 60,000 个峰。红色垂直线显示了在真实峰位移处的主要峰，而蓝色垂直线处的小突起代表读取长度。

弱信号

以下是一个使用 Pol2 数据的弱信号示例。在这里，这种特定的抗体效率不高，因此这些峰是宽泛且分散的。我们在交叉相关图中观察到两个峰：一个在真实峰位移处（约 185-200 bp），另一个在读取长度处。对于弱信号的数据集，读取长度峰将开始占主导地位。

无信号

一次失败的实验会类似于仅使用输入数据的交叉相关图，在这种情况下，我们几乎观察不到片段长度的峰值或根本没有峰值。请注意，在下面的示例中，最强的峰值是蓝线（读取长度），并且在图谱中基本没有其他显著的峰值。峰值的缺失是预料之中的，因为在特定目标位点周围不应该有显著的片段聚集（除了可能在开放染色质区域中存在的弱偏差，这取决于所使用的实验方案）。

报告中还有一些其他信息，可自行探索！！！

Diffbind 差异 peak 鉴定

DiffBind 是一个 R 语言的 Bioconductor 软件包，用于识别在两个或多个样本组之间差异富集的位点。它主要处理峰值调用集（‘peaksets’），这些是代表每个样本的候选蛋白质结合位点的基因组区间集合。DiffBind 包含支持处理峰值集的函数，包括在整个数据集中重叠和合并峰值集、在峰值集中重叠的区间中计数测序读数、以及基于结合亲和力（通过读数密度差异测量）的证据识别统计学上显著的差异结合位点。

DiffBind 的核心功能是差异结合亲和力分析，它可以识别在样本组之间统计学上显著差异结合的结合位点。核心分析程序默认使用 DESeq2 执行，并且可以选择使用 edgeR。每个工具都会为每个候选结合位点分配一个 p 值和 FDR，以指示它们是差异结合的置信度。

# 加载包
library(DiffBind)
library(tidyverse)

# 读取peaks集和相关元数据
samples 'meta/samplesheet.csv')
dbObj 
# 获取比对文件并计算共识集（在多个样本中都检测到的峰值区域）中的每个峰/区域的计数信息
dbObj 
# 构建比对信息
dbObj 
# 进行差异富集分析
dbObj 
# 查看结果摘要
dba.show(dbObj, bContrasts=T) 

# 提取结果
res_deseq out write.table(out, file="results/Nanog_vs_Pou5f1_deseq2.txt", sep="\t", quote=F, row.names=F)

# 为 DESeq2 识别出的每组显著区域创建 BED 文件，并根据富集的增加或减少进行分类
nanog_enrich % 
  filter(FDR  0) %>% 
  select(seqnames, start, end)
write.table(nanog_enrich, file="Nanog_enriched.bed", sep="\t", quote=F, row.names=F, col.names=F)

pou5f1_enrich % 
  filter(FDR % 
  select(seqnames, start, end)
write.table(pou5f1_enrich, file="Pou5f1_enriched.bed", sep="\t", quote=F, row.names=F, col.names=F)

我们先看一下差异富集结果

> res_deseq
GRanges object with 6957 ranges and 6 metadata columns:
           seqnames              ranges strand |      Conc Conc_Pou5f1 Conc_Nanog
                            |      
   519 NC_000001.11 214316564-214316964      * |   7.77216     3.25157    8.74039
   247 NC_000001.11   84479038-84479438      * |   6.97620     1.47896    7.96014
  2168 NC_000005.10   62043581-62043981      * |   7.43368     3.57583    8.38305
  3957 NC_000009.12 131531068-131531468      * |   6.59834     0.00000    7.59834
  3007 NC_000007.14   24427007-24427407      * |   7.60091     2.99133    8.57106
   ...          ...                 ...    ... .       ...         ...        ...
  6935  NC_012920.1          9675-10075      * |         0           0          0
  6936  NC_012920.1         10468-10868      * |         0           0          0
  6937  NC_012920.1         11068-11468      * |         0           0          0
  6938  NC_012920.1         12413-12813      * |         0           0          0
  6939  NC_012920.1         13758-14158      * |         0           0          0
            Fold     p-value         FDR
             
   519  -4.85541 5.36253e-10 3.00568e-06
   247  -5.58771 8.64943e-10 3.00568e-06
  2168  -4.36111 3.49062e-09 8.08661e-06
  3957  -6.64530 6.78297e-09 1.03517e-05
  3007  -4.86481 7.44727e-09 1.03517e-05
   ...       ...         ...         ...
  6935         0           1           1
  6936         0           1           1
  6937         0           1           1
  6938         0           1           1
  6939         0           1           1
  -------
  seqinfo: 34 sequences from an unspecified genome; no seqlengths

• Conc ：所有样本的平均读取浓度（默认计算使用 log2 标准化 ChIP 读取计数并减去控制读取计数）
• Conc_Nanog ：第一组（Nanog）的平均浓度
• Conc_Pou5f1 ：第二组（Pou5f1）的平均浓度
• Fold ：显示两组之间平均浓度的差异，正值表示 Nanog 组的结合亲和力增加，负值表示 Pou5f1 组的结合亲和力增加。

探索数据

样本 PCA

在这一步中我依旧可以对我们对数据做一些探索性分析。

比如根据识别到的所有共识位点绘制PCA图

dba.plotPCA(dbObj,  attributes=DBA_FACTOR, label=DBA_ID)

也可以根据差异分析识别出来的区域绘制 PCA

dba.plotPCA(dbObj, contrast=1, method=DBA_DESEQ2, attributes=DBA_FACTOR, label=DBA_ID)

相关性热图

plot(dbObj)

差异富集两工具交集

在进行差异分析的时候，我们使用了Deseq2和edgeR两种方法，我们可以查看一下有多少区域被两个工具都认定为差异区域。

这里我自己运行会报错，是因为在我的数据中edgeR没有发现差异区域，下图仅为示例

dba.plotVenn(dbObj,contrast=1,method=DBA_ALL_METHODS)

MA图

MA图用于展示在不同条件下结合位点的差异。MA图通过展示结合位点在不同条件下的平均信号强度（A值）和信号强度差异（M值），帮助识别显著差异结合的位点。

• 横轴（A值） ：表示结合位点在两个条件下的平均信号强度。
• 纵轴（M值） ：表示结合位点在两个条件下的信号强度差异（通常是log2 fold change）。
• 点的分布 ：每个点代表一个结合位点。大多数点应当集中在M值为0附近，表示没有显著差异结合。偏离0的点表示差异结合的位点。
• 显著差异结合位点 ：通常会用不同颜色标记显著差异结合的位点（例如，红色表示显著上调或下调的位点）。

dba.plotMA(dbObj, method=DBA_DESEQ2)

1. 横轴（log concentration） ：

• 表示结合位点在两个条件下的平均信号强度（A值）。
• 横轴上的数值是以对数刻度表示的结合位点的平均浓度。

• 表示结合位点在两个条件下的信号强度差异（M值）。
• 纵轴上的数值表示Pou5f1和Nanog之间的log2倍数变化（log2 fold change）。