[13] 9 下游分析 & 9.1 基因组 contex

9 下游分析

---

• 本章介绍如何注释已鉴定的 Peak 的基因组序列，并注释到基因，然后对其进行功能特征分析。它还提出了解决所研究蛋白质的 DNA 序列特异性的初步步骤，并对如何将 ChIP-seq 与其他功能基因组学数据集成进行了展望。

9.1 基因组 context

---

• 转录因子 （ TF ）结合位点的基因组 context 可以告知其在细胞中的潜在功能。ChIP-seq Peak 的基因组分布可以在不同的间隔尺寸水平上进行评估，从全局对染色质类型（ 颜色 ）的分类到单个基因中的特定区域。

9.1.1 基因组位置

• 在大多数情况下，第一步是检查 ChIP-seq Peak 相对于注释基因的位置。然而，同样的方法也可以应用于其他基因组特征，例如重复区域，CpG 岛或增强子区域。

• 基因可分为编码蛋白基因、假基因和非编码 RNAs （ 称为基因生物型 ）。注释包括转录区域，但不包括前面的启动子，启动子通常被定义为转录起始位点（ TSS ）上游的 2kb 。基因本身分为内含子和外显子，如果是蛋白质编码基因，则进一步分为 5‘UTR、CDS 和 3’UTR 。基因组的其余部分被称为基因间区。

• 基因注释（ GTF 格式 ）可从 Ensembl、UCSC 或 NCBI 以及物种特定资源数据库（例如：flybase 、arabidopsis ）。根据来源的不同，注释文件在布局和信息内容方面可能会有所不同。例如，NCBI RefSeq 注释仅包括一组简明的手动整理的转录本，而 Ensembl 报告了潜在的异构体的图谱，包括自动注释的转录本，而没有实验支持。UCSC KnownGenes 是另一个广泛使用的注释源，具有相当数量的转录本。对 长非编码 RNA （ LncRNA ）基因最全面的分析可以从 PTANTOM 项目中获得（ FANTOM ）（ ref1 ）。

• 基因组特征也可以通过 R/Bioconductor 注释包检索（ 详情见 annotation ）。 我们在 R 中提供替代代码，用于在脚本中进行基因组位置分析，见附加在线文件。

• 基因注释可能很难处理，因为许多特征是重叠的。这发生在基因水平上，每个基因的多个转录异构体进一步扩增。重叠注释可以通过基于关于蛋白质功能的先验假设来定义符号的层次结构来解决（ 例如：exon > 5' UTR > 3' UTR > intron > promoter > intergenic ）或通过使用其他类别（ 例如：ambiguous ）。当使用层次结构时，需要注意确保相关分布不是由于强加的层次结构，而是反映了明确分配的 Peak 的分布。

• 重叠注释的问题由于 ChIP-seq Peak 可能非常宽而进一步恶化。解决这个问题的一种方法是只使用 TF Peak。对于较宽的区域，例如组蛋白修饰，可以考虑重叠的程度，使得例如需要 >50% 的 Peak 区域位于给定特征内。或者，可以通过与每个注释的重叠部分将区域指定给多个特征。以下代码显示了如何将 NRF1 和 H3K27ac Peak 分配给不同基因组特征的示例。它使用基于编码蛋白基因的Ensembl 注释的小鼠基因组预处理文件（ mm10）（ 作为附加在线文件提供 ）。

  # 使用以下层次结构的基因组位置：
  # exon(EXON) > 5’UTR (5UTR) > 3’UTR (3UTR) > intron (INTRON) > 2kb upstream promoter (P2000) > intergenic (INTER)
  
  head mm10_genomic_features.bed
  
  # 定义 TF Peak summit 区域
  
  awk -v OFS="\t" '{print $1,$4-1,$4}' peaks/NRF1_CHIP_WT_1_peaks_peakzilla.bed > peaks/NRF1_CHIP_WT_1_peaks_peakzilla_summits.bed
  
  # TF Peak summit 基因组中的位置
  intersectBed -wo -a peaks/NRF1_CHIP_WT_1_peaks_peakzilla_summits.bed -b mm10_genomic_features.bed |\
  awk -v OFS="\t" '{F[$7]++;t++}END{for(location in F){print location,F[location]*100/t}}' > peaks/NRF1_CHIP_WT_1_peaks_peakzilla_summits_location.txt
  
  cat peaks/NRF1_CHIP_WT_1_peaks_peakzilla_summits_location.txt
  
  

  • 组蛋白修饰

    # 组蛋白峰区域按基因组位置的重新划分
    intersectBed -wo -a peaks/H3K27AC_CHIP_WT_1_peaks_macs.bed -b mm10_genomic_features.bed |\
    awk -v OFS="\t" '{F[$10]+=$12;t+=$12}END{for(location in F){print location,F[location]*100/t}}' > peaks/H3K27AC_CHIP_WT_1_peaks_peakzilla_regions_location.txt
    
    cat
    peaks/H3K27AC_CHIP_WT_1_peaks_peakzilla_regions_location.txt
    

  • 基因组位置的整体重新划分

    awk -v OFS="\t" '{F[$4]+=$3-$2;t+=$3-$2}END{for(location in F){print location,F[location]*100/t}}' mm10_genomic_features.bed > peaks/genome_location.txt
    
    cat peaks/genome_location.txt
    

  • 绘图

    # 导入数据
    tf = read.table("peaks/
    NRF1_CHIP_WT_1_peaks_peakzilla_summits_location.txt")
    
    histone = read.table("peaks/H3K27AC_CHIP_WT_1_peaks_peakzilla_regions_location.txt")
    
    genome = read.table("peaks/genome_location.txt")
    
    # 合并三个文件
    d = merge(merge(tf,histone,by=1), genome, by=1)
    
    colnames(d) = c("Location","TF","Histone","Genome")
    
    # 基因组位置分布的饼图
    pdf("peaks/genomic_location_piechart.pdf",height=5,width=15)
    
    par(mfrow=c(1,3),bg="white")
    pie(d$TF,labels=d$Location,main="NRF1")
    pie(d$Histone,labels=d$Location,main="H3K27AC")
    pie(d$Genome,labels=d$Location,main="Genome")
    
    dev.off()
    
    # 基因组富集区的条形图
    pdf("peaks/genomic_location_barplot.pdf")
    
    par(mfrow=c(1,2),bg="white")
    barplot(log2(d$TF/d$Genome),
            names=d$Location,main="NRF1",
            las=2,ylim=c(-2,6))
    
    barplot(log2(d$Histone/d$Genome),
            names=d$Location,
            main="H3K27AC",
            las=2,ylim=c(-2,6))
    
    dev.off()
    
    q()
    

• 结果显示为饼图，以表示不同基因组区域中 Peak 的相对出现比例（ 图 9.1A ）。由于特征在基因组中以不同的频率出现，因此可以通过基因组中每种特征类型的总核苷酸数来标准化计数，显示为富集的条形图（ 图 9.1B ）。

9.1.2 距离基因的距离

• TF 结合位点可以发生在启动子区域内（ TSS 的近端 ）或基因间区位置（TSS的远端 ）。为了区分位于 TSS 近端或远端的 Peak ，检查每个 Peak 与最近的 TSS 的距离，而与特定的目标基因分配无关。由于许多 TF 既结合近端点又结合远端点，因此到 TSS 的距离通常呈双模态分布（ 图 9.1C ）。对于也可以具有位置偏好的组蛋白修饰，预期有不同的模型。例如，H3K27me3 修饰几乎只发生在启动子区，而 H3K4me3 修饰和 H3K4me1 修饰之间的平衡允许区分启动子和增强子区域（ ref2 ）。

• 下面的代码计算每个 Peak 到最接近的 TSS 的距离。它使用基于所有编码蛋白基因转录本的 Ensembl 注释的小鼠基因组预处理文件（ mm10 ）。

  # 所有基因的 TSS 位置
  ls mm10_tss.bed
  
  # TF Peak summits 距离 TSS 的距离
  awk -v OFS="\t" '{print $1,$4-1,$4}' peaks/NRF1_CHIP_WT_1_peaks_peakzilla.bed |\
  closestBed -d -t "first" -a stdin -b mm10_tss.bed |\
  awk '{print $NF}' > peaks/NRF1_CHIP_WT_1_peaks_peakzilla_summits_dist_tss.txt
  
  # 组蛋白 Peak 区域中心到 TSS 的距离
  awk -v OFS="\t" '{c=($2+$3)/2;print $1,c-1,c}' peaks/H3K27AC_CHIP_WT_1_peaks_macs.bed |\
  closestBed -d -t "first" -a stdin -b mm10_tss.bed |\
  awk '{print $NF}' > peaks/H3K27AC_CHIP_WT_1_peaks_macs_center_dist_tss.txt
  

  • R 中可视化

    # 导入距离 TSS 的矩阵
    tf = read.table("peaks/
    NRF1_CHIP_WT_1_peaks_peakzilla_summits_dist_tss.txt")
    
    histone = read.table("peaks/
    H3K27AC_CHIP_WT_1_peaks_macs_center_dist_tss.txt")
    
    # Peak 距最近TSS的距离直方图
    pdf("peaks/dist_tss_hist.pdf")
    
    par(mfrow=c(2,1),bg="white")
    hist(log10(tf[,1]),
         main="NRF1",
         xlab="Distance to closetgene TSS (log10)",
         breaks=seq(0,7,0.2))
    
    hist(log10(histone[,1]),
         main="H3K27AC",
         xlab="Distance to closet gene TSS (log10)",
         breaks=seq(0,7,0.2))
    
    dev.off()
    
    q()
    

• 结果显示为以 log10 比例显示的直方图，以覆盖广泛的距离（ 图 9.1C ）。这使得可以更好地可视化 NRF1 和 H3K27ac 结合特征的近端和远端 Peak 的双峰分布。近端和远端区域的自然分离出现在2kb 左右。NRF1 峰在代表富含 CpG 启动子区域结合的基因的近端略微富集，而远端结合位点具有较低的 CpG 含量（ ref3 ）。H3K27ac 峰在基因的远端更加富集，代表对远端活性增强子区域的预期偏好。

  图 9.1 Peak 区间的基因组序列

  （ A ）Peak 基因组中的分布比例饼图

  （ B ）柱状图

  （ C ）Peak 距离 TSS 位点距离的条形图

  image

9.2 功能分析

---

• 一种流行的下游分析是探索靶基因的功能。

9.2.1 注释到靶基因

• Peak 到基因的分配仍然是一项不平凡的任务，因为 TF 和增强子可以从非常长的差异激活它们的目标基因，小鼠中的基因被位于 1Mb 之外的增强子调控（ref4）。即使已经探索了几个概念来分配目标基因，最简单和最有效的方法是使用最近的 TSS（ ref5 ）。理想情况下，重新开发的技术，如 Capture Hi-C（ Chi-C ) ( ref6 )，可以用来推断可靠的关联，但数据的可用性和处理仍然是有限的。

  # 最接近 NRF1 Peak 的唯一基因列表（ 7167个 Peak 注释到 5595 个基因 ）
  closestBed -t "first" -a peaks/NRF1_CHIP_WT_1_peaks_peakzilla.bed -b mm10_tss.bed |\
  awk '{print $10}' | sort -u > peaks/NRF1_CHIP_WT_1_peaks_peakzilla_genes.txt
  

9.2.2 基因富集分析

• 以基因本体论（ GO ）的形式在许多物种上都可以获得对基因功能的全面描述（ ref7 ）。GO被组织成三个不重叠的本体，它们描述蛋白质的生理作用（ 生物学过程：Biological Process ），分子活性（ 分子功能：Molecular Function ）或在细胞内的位置（ 细胞成分：Cellular Component）。此外，分配给蛋白质的每个 GO 术语都与一个 GO 号相关联，指定所分配的功能是例如通过实验验证的，还是仅仅从正交学中推断出来的。

• 基于GO注释，可以检验一系列基因特定功能的富集。对于每个GO term，将列表中与该 term 相关联的基因的部分与其总体出现进行比较，以识别明显过度表达的 term。显著性通常使用超几何检验的 p 值来计算。值得注意的是，GO 富集可能受到 baseline 选择的强烈影响，即是否对基因组中的所有基因或一组特定的 control 基因（ 即背景文件 ）进行富集检验。通常应用的control集都是表达基因（ 例如根据 RNA-seq 数据 ）或具有共享的和差异的 ChIP-seq Peak 的基因。用于 GO 分析的流行在线工具包括 David 以及用于可视化结果的 REViGO。

  注意：与用于 Peak calling 的阈值选择类似（ 参见 第6.2.5章 ），应始终根据 p 值而不是变化倍数对富集的类别进行排序和选择。在报告或可视化围棋分析结果时，应避免任意选择GO terms。应提供完整的富集注释信息表作为补充信息。

9.2.3 其它类型的基因富集分析

• 富集的概念可以扩展到在研究上下文中感兴趣的任何预定义的基因列表。例如，可以对目标基因进行检验以富集发育调节基因或某一蛋白质的相互作用伙伴。可以从已发表或数据库中检索参考文献列表，也可以手动编辑参考文献列表。
• 另一个流行的功能注释来源是 KEGG 数据库，它收集手动整理的生物学途经。最初为酶和代谢过程设计的 KEGG 现在包含了数百张手工绘制的 map，包括人类疾病和药物设计（ ref8）。KEGG Mapper 工具允许将基因列表映射到通路上，通路图可以根据用户定义的信息进行着色。最后，像 g：profiler 这样的工具将广泛的不同功能注释集成到一个联合资源中，以便能够对基因列表进行全面的功能解释。

9.3 序列分析

---

• 分析 Peak 区域下的 DNA 序列提供了对所研究蛋白质的 DNA 结合偏好或在相邻位置重复结合的潜在协同因子的洞察。

9.3.1 Motif 分析

• De novo motif discovery motif 分析中的第一个策略是在没有先验假设的情况下搜索富含 Peak 区域的序列，也称为从头 motif 发现。搜索通常在围绕 TF Peak summits 或组蛋白修饰的整个区域的 50-200bp 的窗口中执行。大多数 Motif 发现工具都遵循基于 word-based或基于 profile-based 的方法（ ref9 ）。在例如在DREME（ ref10 ）中实现的基于 word-based 的方法中，所有可能的 k-mer（ 即长度为 k 的序列 ）都被穷举以生成在输入序列中以增加的频率出现的共识基序。相反，基于 Profile-based 的方法，如 MEME（ ref10 ），迭代地优化序列比对以获得最佳评分motif。最近，应用深度学习方法来发现 ChIP-seq 数据中的结合 motif（ ref11 ）。

• Motifs 在整个基因组中出现的频率很高。因此，任何富集的基序都应始终对照背景序列进行检验，要么由用户提供，要么由 randomisation 生成。这些背景序列的选择可能会强烈影响所发现的 motif。

• HOMER 是一种可以通过命令行运行的流行工具。它将目标区域和背景区域的基因组坐标作为输入，或者生成具有匹配目标区域的 GC 含量的可能性的随机背景区域。MEME-ChIP （ ref12 ）是一个所谓的集成工具，它结合了几种 Motif发现算法。它可以作为在线工具运行，将目标区域和背景区域的 FASTA序列作为输入，或使用随背景字母频率变化的随机控制。

  注意：Motif 表示为位置权重矩阵（ PWM ），这些矩阵由多序列比对构建而成。PWM 报告 motif 中每个位置的每个核苷酸出现的概率，这可以被可视化为 Sequence logo。

  # 定义 Peak summits 两侧 75 bp 区域
  awk -v OFS="\t" '{print $-1, $4-75, $4+75}' peaks/${sample}_peaks_peakzilla.bed >
  peaks/${sample}_peaks_peakzilla_151bp.bed
  
  # HOMER 进行 从头 motif 发现
  findMotifsGenome.pl
  peaks/${sample}_peaks_peakzilla_151bp.bed genomes/mm10/motifs -size given
  
  # 打开 html 结果文件
  see motifs/homerResults.html &
  

• HOMER 输出在目标序列中找到的 Motif 的排序列表（ 图 9.2A ）。对于每个 motif，它表示序列（以 logo 表示 ）与背景序列相比，靶标中该 motif 的富集相对应的 p 值，以及已知 motif 中该 motif 的最佳匹配。在 NRF1 中，如预期的那样，发现与已知的NRF1 motif 匹配的从头识别的 motif 在 Peak 区域中最富集，大约 64%。

• 已知 motif 搜索 motif 分析中的第二个策略是扫描已定义 motif 的Peak 区域，也称为已知 motif 搜索。许多 TF 的 motif 现在已经从体外（ 例如通过指数富集（ SELEX ）（ ref13 ）或蛋白质结合矩阵（ PBM ）或体内（ 例如使用 ChIP-seq ) 实验获得，并且可以在公共数据库中获得（例如：JASPAR （ ref14 ） 或者 HOCOMOCO （ ref15 ））。已知基序的 PWMs 可用于扫描感兴趣的基因组区域以识别 motif （ 例如：使用 MAST（ ref15 ） ）。为了选择有意义的 Motif 出现，需要应用 p 值阈值，我们建议根据 motif 的信息内容进行调整（ 例如： 根据 motif 的长度，相同的阈值将具有不同的严格性）。下面的代码显示了如何在我们的Peak 区域搜索已知的 NRF1 motif。

  # 从 JASPAR 下载 NRF1 motif
  # jaspar.genereg.net/matrix/MA0506.1
  wget http://jaspar.genereg.net/api/v1/matrix/MA0506.1.meme -O motifs/NRF1.meme
  
  #使用 p 值阈值 10ˆ-5 扫描基因组中出现的 motif，并将输出重新格式化到 BED 中
  mast -hit_list -mt 1e-04 motifs/NRF1.meme genomes/mm10.fa |
  awk ’($1!˜/#/){if($2=="+1"){s="+"}else{s="-"};print
  $1,$3-1,$4,"NRF1",$6,s}’ |\
  gzip > motifs/NRF1_mm10.bed.gz
  
  # Peak 区间数目：7167
  cat peaks/${sample}_peaks_peakzilla_151bp.bed |\
  wc -l
  
  
  # 有 NRF1 motif 的 Peak 区间数目：5245（ 73% ）
  intersectBed -u -sorted -a
  peaks/${sample}_peaks_peakzilla_151bp.bed -b
  motifs/NRF1_mm10.bed.gz | wc -l
  

• 使用已知的 NRF1 motif 在特定阈值下，我们发现 73% 的Peak 区域含有一个 motif。在 TF 的 ChIP-seq 数据中，带有 Motif 的 Peak 的比例通常在 60-80% 左右。一些非特异性峰可能是由实验偏差引起的，如 crosslinking artefacts 。可以将相同的代码调整为在 Control 区域上运行（ 使用命令 shuffledBed 生成 ）。可替换地，可以使用 Peak 的子选择，例如 TKO 特定的 Peak与共享的 Peak。最后，可以使用超几何检验来统计评估targets 区和control 区的富集程度的比较（例如：使用 R 中的函数 phyper ）。同样的分析可以运行更多的 motif，甚至所有可能的 k-mers。与从头开始的 motif 发现方法相比，使用已知 motif 扫描Peak区域的优点是，该信息可以用于进一步的分析，例如探索不同Motif 在特定区域中的组织和共生（ 例如，彼此之间的距离或方向 ）。此外，计算 metaplot 中的位置富集使我们能够可视化是否以及在何处在 Peak 周围富集了 motif。

9.3.2 序列保守性

• 当具有额外物种的多个比对可用时，可以探索 Peak 或 motif 的保守性水平。为此，可以从 UCSC 基因组浏览器以 bigwig 格式下载PhastCons 或 PhyloP 等保守性分数，并且可以使用 bwtool 或 bedtools 进行处理（ 见 章节 9.4.1）。

  图 9.2：Motif 和其他ChIP-seq 数据集集成

  （ A ）HOMER 从头搜索到的前三个 motif

  （ B ）JASPAR 数据库中已知的 NRF1 motif logo

  （ C ）在 WT 和 TKO 样本共有的 NRF1 Peak 中 NRF1 和 H3K27ac 的 reads 密度分布热图

  （ D ）来自图 C 的热图信号的 metaplot 图

  image

9.4 结合其他数据分析

• 基因组研究通常需要几种类型的实验来解决特定的生物学问题。此外，可以公开获得大量相关的基因组数据集。因此，ChIP-seq 数据与其他数据类型的结合分析是一种常见的分析。这种数据集成的一个示例可以在 NRF1 数据集的原始发布中找到。

9.4.1 额外的 ChIP-seq 数据集

• 第一步通常是与其他 ChIP-seq 数据集集成，这可能包括 TF 和组蛋白修饰的数据组合。

  注意：为了避免任何偏见和错误解释，强烈建议使用包括数据预处理（ 例如 reads 长度，修整 ）， reads 比对（ 例如索引，用于唯一 reads 的过滤阈值 ）和 Peak calling ( 例如算法、Peak 阈值 ）的类似流水线来处理每种类型的数据集（ 或重新处理公共数据）。

• 可视化和比较 TF 和组蛋白修饰的几个 ChIP-Sseq 数据集的流行方法是生成 Peak 区域中 reads 密度的热图。这种整合应该考虑到识别的 Peak 区域的不同性质：组蛋白修饰的信号通常较宽，并且 Peak 在 TF 信号周围。因此，建议对以特定位置为中心的区域进行比较分析，如 TF Peak summits 或 TSS，而不是合并所有富集区域。下面，我们提供代码为跨样本的 NRF1 共享和差异 Peak 区域生成这样的热图。有几种对用户友好的在线工具可用于根据测序数据生成热图和其他表达图（例如：deeptools2）。

  # 对于每个 Peak 选择，对于所有样本，从bigwig文件中提取 Peak 区域周围 5kb 内每个位置的 reads 密度
  
  for peaks in TKO_spec WT_spec shared
  do
    for sample in NRF1_CHIP_WT_1 NRF1_CHIP_TKO_1 H3K27AC_CHIP_WT_1 H3K27AC_CHIP_TKO_1
    do
        awk -v OFS="\t" '{center=int(($2+$3)/2);print $1,center-2500,center+2500,$2,$3}' changes/NRF1_all_regions_count_${peaks}_table.txt |\
        bwtools extract -tabs bed stdin
        tracks/${sample}.bw stdout >        changes/NRF1_all_regions_count_${peaks}_density_${sample}.txt   
    done
  done
  

  • R 可视化

    library(gplots)
    
    # 导入 RPKM 文件
    rpkm = read.table("changes/NRF1_all_regions_rpkm.txt")
    
    colnames(rpkm) = c("chr","start","end",
                   "NRF1_CHIP_WT_1","NRF1_CHIP_WT_2",              "NRF1_CHIP_TKO_1","NRF1_CHIP_TKO_2")
    
    # 选择每一个 Peak 
    for(peak in c("TKO_spec", "shared")){
        png(paste("changes/NRF1_all_regions_count_", peaks, "_density_heatmap.png", sep = ""))
        par(bg = "white", mfrow = c(1, 4))
        
        # For each sample
        for(sample in c("NRF1_CHIP_WT_1", "NRF1_CHIP_TKO_1", "H3K27AC_CHIP_WT_1","H3K27AC_CHIP_TKO_1")){
            # Load the density table
            x = read.table(paste("changes/NRF1_all_regions_count_",peaks,"_density_",sample,".txt",sep=""))
          
            # Add extra column with matching RPKM of region in NRF1_CHIP_WT_1 sample
            y=merge(x,             rpkm[,c("chr","start","end","NRF1_CHIP_WT_1")],
                    by.x=c(1,4,5),
                    by.y=c("chr","start","end"))
           # #绘制读取密度的图像，按 RPKM 对区域(行)进行排序，将色标从白色调整为黑色 0.1 (0到10之间)，再 one step 直到100 
            image(t(y[order(y$"NRF1_CHIP_WT_1"),7:5006]),
                  axes=F,
                 col=colorpanel(101,"white","black"),
    breaks=c(seq(0,10,0.1),100),main=sample)
        }
        dev.off()
    }
    
    # Plot colour scale
    pdf("changes/
    NRF1_all_regions_count_density_heatmap_scale.pdf",height=3)
    
    par(bg="white")
    
    plot(c(0,110),c(0,1),
         xlab="Readdensity",
         ylab="",
         pch="",axes=F)
    
    axis(1,at=seq(0,110,10),labels=c(0:10,100))
    
    cols=c(colorpanel(100,"white","black"),rep("black",10))
    
    rect(seq(0,109,1),0,seq(1,110,1),1,border=cols,col=cols)
    
    dev.off()
    
    q()
    

• 密度热图显示在共享 Peak 和 TKO 特定 Peak 中的 Peak 周围的 reads 密度的分布（图 9.2C ）。它也可以从 第 8.3 章 的整个 Peak 列表中生成。

  注意：重要的是要记住，尽管热图是很好的可视化工具，但它们不是表示数据的具体方式，因为颜色比例的细微变化可能会对人眼产生误导。在这里生成的密度热图的示例中，由于以相对较小的数字显示数千个区域，如果行不会按降序 RPKM 值排序，则具有低 reads 密度的一些区域在具有高 reads 密度的区域之间将不可见。此外，用户很容易以非线性步骤排列颜色标度以突出特定的特征，例如在我们的情况下，我们使用从白色到黑色的线性标度从0到10，并注释所有大于10直到100到黑色的附加值，因为密度值遵循下降的指数曲线。

9.4.2 表达数据

• 将 ChIP-seq 与RNA-seq 或芯片数据的基因表达信息结合，允许我们研究 TF 的结合或组蛋白修饰的存在是否与其目标基因的表达相关。为此，可以将 Peak 区域的信号与假定的目标基因的表达水平进行比较。如果有几个条件可用，在所谓的 delta-delta 散点图中比较结合和基因表达的变化可能更具信息性（ 见 章节 8.3.4 ）。请注意，由于可以将几个 Peak 分配给同一基因，因此某些基因表达值可能会多次出现。这可以通过取给定目标基因的所有相关 Peak 的最小、最大或平均信号来解决。

  注意：这种分析对假阳性目标基因引入的噪音很敏感，当被分配到最接近的 TSS 的基因时。由于远端 Peak 可能比近端 Peak 更经常被错误分配，因此分别对近端和远端 Peak 进行下游分析是有用的（ 例如，≤ 2kb vs > 2kb ）。

9.4.3 其他类型数据

• 最后，其他类型的基因组数据也可以整合到分析中，例如染色质可及性（ 例如 DNase-seq 或 ATAC-seq ）或DNA甲基化（ 例如 WGBS-seq 或 RRBS-seq)。这可以通过使用 bwtool 或 bedtools 在 Peak 区域上汇总信息来执行（ 见 章节 9.4.1 ）。例如，可以在 Peak 区域上计算信号平均值或在峰值区域的相关子集中进行比较。与基因表达的比较类似，delta-delta 散点图可用于将结合的变化与染色质可及性或DNA甲基化的变化进行比较。
• 除了用于靶基因分配，来自高分辨率基于 Hi-C 的方法的数据也可以通过比较 ChIP-seq 结合的变化与包含差异 Peak 的基因组区域的相互作用谱的变化来集成。

An atlas of human long non-coding RNAs with accurate 5' ends.

Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome

A long-range Shh enhancer regulates expression in the developing limb and fin and is associated with preaxial polydactyly

Assessing computational methods for transcription factor target gene identification based on ChIP-seq data.

Mapping long-range promoter contacts in human cells with high-resolution capture Hi-C.

Gene ontology: tool for the unification of biology

KEGG: new perspectives on genomes, pathways, diseases and drugs.

A review of ensemble methods for de novo motif discovery in ChIP-Seq data.

DREME: motif discovery in transcription factor ChIP-seq data.

Fitting a mixture model by expectation maximization to discover motifs in biopolymers.

Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning.

Motif-based analysis of large nucleotide data sets using MEME-ChIP.

DNA-binding specificities of human transcription factors

JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framewor

HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP-Seq analysis

Combining evidence using p-values: application to sequence homology searches.

Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes.

Detection of nonneutral substitution rates on mammalian phylogenies.