使用ChIPSeeker进行ChIP-seq， ATAC-seq，cut&tag等富集峰的基因组注释

二代测序产生的数据类型

常规的下一代高通量测序（next generation sequencing, NGS）实验通常产生大量短片段(reads)，通常我们需要将这些reads比对到参考基因组/转录组上，即将它们置于生物学上有意义的基因背景下，才能获得有意义的结果。一般我们认为会产生两种类型的数据（当然两者并无严格意义上的区分）：

   1.表达类

一般为固定区域，关注于定量比较。例如转录组测序结果中，数据库中已有mRNA基因的表达，lncRNA基因的表达等，这类结果一般以矩阵形式存储，第一列是名字，其余列是表达值。CircRNA表达，miRNA表达等归于此类。

    2.基因组区域（富集峰)类

一般区域不固定，关注于定性。例如ChIP-seq、ATAC-seq、Cut&tag等比对后获得的富集峰。一般以bed格式存储，第一列是染色体，第二列是富集峰起始坐标，第三列是富集峰终止坐标（图1）。eccDNA，m6A，MeDIP等归于此类。

图1. 表达类（区域固定） vs 富集峰类（区域不固定）

什么是富集峰注释？

基于抗体富集的原理，众多reads片段比对到基因组上某区段，会形成一个类似山峰的富集区。由于我们是在基因组背景下进行生物医学研究的，因此需要将基因组区域（富集峰，peak）与基因联系起来，即确定峰落在哪个基因上，落在该基因的哪种基因组特征上，距离TSS的位置是多少bp等，然后才能进行后续的功能研究。这个过程叫做富集峰注释（peak annotation，图2）。如果仅关注某几个区域，就不需要用软件注释，建议直接用IGV或者UCSC Genome Browser查看。

富集峰注释的难点

虽然CHIP-Seq已经有近20年的历史，然而由于基因组上基因的结构非常复杂，不同注释软件在具体细节处理上往往不同，从而导致同一批数据，用不同的软件进行注释，获得的结果略有不同（大同小异）。

1. 注释到转录本还是注释到基因？

由于一个基因可能包含多个转录本，因此，我们在注释的时候，到底是注释到基因水平还是注释到转录本水平？

     2.基因位置重叠怎么处理？

由于同一个位置可能存在多个基因，如果两个基因的坐标有重叠，我们到底是注释到A基因还是注释到B基因？

    3.最邻近怎么判断？

如果富集峰的中点（或者顶点）正好落在两个基因的中间，那么这个富集峰是注释到A基因还是注释到B基因，如何定义最邻近？

    4.基因组特征如何分类？

不同文章、软件对基因组特征分类不同，例如有的分为promoter、intron、exon，5’UTR和3’UTR；有的分为：upstream、promoter、intron、exon、downstream等

    5.promoter如何定义？

不同文章、软件定义promoter区也不同，有的定义TSS上游3K到下游3K都是promoter，有的定义TSS上游200bp到下游800bp是启动子。

    6.不同基因组特征的注释优先级

当一个很长的富集峰横跨同一个基因的intron、exon、3’UTR时候，这个富集峰该分到什么特征中呢？

    7.注释数据库版本问题

同样是human hg38，如果注释库版本不同，那么注释结果也会有差异。原因是：虽然基因组序列不变，然而注释库却更新频繁，有基因会更新坐标，有基因会添加新的转录本，有基因会从非编码基因变成编码基因等。

图2. 富集峰注释及难点

凡此种种，给我们的注释工作带来了巨大困难。然而，作为用户（调包侠），我们基本不用深究注释背后的细节。我们需要做的就是：找一个引用比较多的注释工具，默认参数进行注释即可。

常见富集峰注释软件

软件	语言	默认支持物种
ChIPSeeker	R	有注释文件即支持
PAVIS	在线	人、大鼠、小鼠（细分编码和非编码等）
Homer	Perl	有注释文件即支持
GREAT	在线	human，mouse
ChIPpeakAnno	R	有注释文件即支持

表1. 常见富集峰注释工具

为什么要用新版注释？

由于注释数据库频繁更新，如果你使用的注释还是N年前的，那么reviewer在公共数据库（例如UCSC、Ensembl、NCBI）上使用网站默认版本查询时，就有可能查不到你的基因，或者你N年前的数据，与新的数据联合分析时，由于使用的注释数据库不同，取交集时，会漏掉一些基因。因此，我们强烈建议所有的测序数据，包括RNA-seq、ChIP-seq、m6A-seq等都使用同一套注释库进行注释分析，并在结果中明确说明所使用的注释库版本。这对于在不同公司，不同时间做的测序结果来说，是非常重要的。

由于上述所列在线工具都是N年前的，所以我们使用ChIPSeeker R包搭建了一个简易的在线peak注释工具，可以对人、大鼠、小鼠的ChIP-seq，ATAC-seq，cut&tag等富集峰进行一键注释。

1.打开绘图页面

首先，使用浏览器（推荐chrome或者edge）打开ChIP-Seq富集峰注释页面。左侧为常见作图导航，中间为数据输入框和可选参数，右侧为描述和结果示例。也可以在搜索框中搜索peak，找到注释页面。

http://www.bioinformatics.com.cn/basic_chipseq_atacseq_peak_annotation_by_chipseeker_t017

图3. 注释页面

2.示例数据

点击右侧“示例数据”链接下载excel格式的示例数据。

示例数据包括4列，分别为：chr，start，end，-LOG10(pvalue)。

注意：为了遵循各大数据库的使用，这里染色体必需使用chr+数字，即：chr1-22、chrX、chrY、chrM等。

图4. 输入数据示例

3.粘贴示例数据

直接拷贝示例数据中的ABCD四列数据，然后粘贴到输入框。注意必需带每一列的说明行（header），此行将用于最终的excel表头。

注意：不是拷贝excel文件，是拷贝excel文件里边的数据。另外粘贴到输入框后，格式乱了没关系，只要在excel中是整齐的就行。同时数据矩阵中不能有空的单元格，中文字符等。

图5. 必需输入

4.修改参数，并提交

我们设置了promoter区的范围选项，及注释所用的物种及注释版本选项（例如human所使用的是hg38基因组，注释版本为Ensembl v108），当前仅支持human，mouse和rat。后续将支持更多物种。

图6. 可选参数

5.提交出图

粘贴好输入数据，调整好参数（重点是物种及注释版本）后，点击提交按钮，约30秒钟后（取决于数据多少），会在页面右侧出现peak分类饼图及excel格式的数据下载链接，请下载后解压查看。

图7. 结果页面

图8. 注释结果

结果说明

以peak_class为界，结果包括两部分：左侧为输入的内容，其中start添加了1 bp（因为bed格式是0-based，这里变成了1-based），并添加了peak长度信息（end-start+1）；右侧为注释信息，包括：peak分类，基因位置，基因/转录本注释等信息。并基于peak_class的数据绘制了peak分布饼图。

注意：

     1.由于peak注释与注释库及优先级关系密切，因此最终放在paper里边的图，以IGV可视化结果为准。

     2.输入peak默认不考虑链，如需更精细地注释，请参考ChIPSeeker R包。

     3.默认一个peak仅注释到一个转录本

没有预览就是没有出图/结果，这时请参考示例数据，检查输入数据的格式。或者使用我们提供的小工具 pyinstaller打包python脚本为exe可执行文件实例：错误排查小脚本检查输入。

参考文献：

1. Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015 Jul 15;31(14):2382-3. doi: 10.1093/bioinformatics/btv145. Epub 2015 Mar 11. PMID: 25765347.

2. Huang W, Loganantharaj R, Schroeder B, Fargo D, Li L. PAVIS: a tool for Peak Annotation and Visualization. Bioinformatics. 2013 Dec 1;29(23):3097-9. doi: 10.1093/bioinformatics/btt520. Epub 2013 Sep 4. PMID: 24008416; PMCID: PMC3834791.

3. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010 May 28;38(4):576-89. doi: 10.1016/j.molcel.2010.05.004. PMID: 20513432; PMCID: PMC2898526.

4. McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, Wenger AM, Bejerano G. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010 May;28(5):495-501. doi: 10.1038/nbt.1630. Epub 2010 May 2. PMID: 20436461; PMCID: PMC4840234.

5. Zhu LJ, Gazin C, Lawson ND, Pagès H, Lin SM, Lapointe DS, Green MR. ChIPpeakAnno: a Bioconductor package to annotate ChIP-seq and ChIP-chip data. BMC Bioinformatics. 2010 May 11;11:237. doi: 10.1186/1471-2105-11-237. PMID: 20459804; PMCID: PMC3098059.

微生信助力高分文章，用户73000+，引用990+