作了ChIP-seq试验，经过初步下机处理并比对参考基因组之后，获得了目标蛋白在基因组区域的结合峰位置。接下来，就需要对峰位置进行注释，获得这些峰结合在了哪些基因上，位于基因的上下游还是内部，以便了解目标蛋白的功能。
那么，怎样对ChIP-seq峰的结合位置进行注释呢？本篇简介一个非常好用的R包，ChIPseeker，能够快速高效地获得预期结果。
输入数据就是记录ChIP-seq峰位置的bed文件，这是ChIP-seq的常见数据格式，测序公司都会给提供的。如下示例，是在人类细胞中进行的ChIP-seq试验，经下机处理并和人参考基因组hg38比对后获得的bed文件。
bed文件中包含5列信息，以tab键分隔。第1列为峰值所处的染色体id，第2-3列为峰值在该染色体中的碱基位置，第4列为峰的id，第5列为信号强度。用Excel打开bed文件就是长这样。
ChIPseeker包可以使用Bioconductor进行安装。
为了对ChIP-seq的峰位置进行注释，查看蛋白主要结合到基因组的哪些区域，首先需要准备对应物种的基因组注释库。
例如，上述给定的示例bed文件来源物种是人，因此就需要准备人类参考基因组注释库。有些R包，例如org.Hs.eg.db中提供了已经构建好的人类参考基因组注释库，可以直接调用。此外，也可以通过GenomicFeatures包，读取基因组注释文件gff3本地构建。
推荐通过GenomicFeatures包本地构建注释库，因为绝大多数物种都是没有现成的库可用的，方法灵活。而且，像人类参考基因组也是有很多版本的，可能已知的库和自己使用的基因组版本不匹配，用错了就尴尬了，为了避免万一还是推荐从头构建一下。
好了，接下来就展示ChIPseeker注释ChIP-seq峰的全过程。
已知上述ChIP-seq的bed文件获取过程中，比对的参考基因组来自hg38，因此首先通过本地准备的人类参考基因组hg38的gff3注释文件，构建本地注释库，随后读取ChIP-seq的bed文件进行注释。
结果中，注释了ChIP-seq峰结合的染色体位置、区段、基因名称、基因位置等信息。
对于基因组区域，基因启动子区、外显子或内含子区均有。后续就可以根据这些信息作一些筛选，例如ChIP-seq试验使用的蛋白如果是转录因子，那么通常就是在启动子区发挥作用，只关注启动子区的结合位点即可。
同时，ChIPseeker包中自带了一些可视化函数，可以进行一些简单的统计帮助我们了解数据分布，例如ChIP-seq峰所在基因不同位置的数量占比等。
很多情况下，ChIP-seq试验不止作了一次，可能存在多个bed文件有待处理。如果这时，一个个分别注释就会略显繁琐。好在ChIPseeker也提供了批量注释的方法，如果有多个bed文件，可以放在一个list里面统一执行，更加高效。
参考以下示例，两个bed文件的注释，参考基因组仍然是人类hg38，因此继续使用上文构建好的库执行。
输出表格略，内容和单个运行的结果是一致的。对于输出结果图，则是将两组放在一起比较的图。
上一步骤用IDR对重复样本peaks的一致性进行了评估，同时得到了merge后的一致性的.narrowPeak文件，接下来就是对peaks的注释。
这篇主要用Y叔的R包ChIPseeker对peaks的位置（如peaks位置落在启动子、UTR、内含子等）以及peaks临近基因的注释。
.bed文件
由上游测序数据处理而来
TxDb文件
TxDb.Hsapiens.UCSC.hg19.knownGene
同名的包中含有相应的文件，直接引用即可，同样hg38也有同名的包，Bioconductor提供了30个TxDb包，可以供我们使用。要与测序数据比对一致
如果实在没有TxDb呢？
使用GenomicFeatures包函数来制作TxDb对象（这也是为什么说ChIPseeker支持所有物种）
用人的参考基因信息来做注释，从UCSC生成TxDb:
makeTxDbFromUCSC：通过UCSC在线制作TxDb
makeTxDbFromBiomart: 通过ensembl在线制作TxDb
makeTxDbFromGRanges：通过GRanges对象制作TxDb
makeTxDbFromGFF：通过解析GFF文件制作TxDb
这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。
有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。
1.peakAnnoList
2.查看具体信息，会显示出每个peak的位置，所在的区域
as.GRanges(peakAnnoList) %>% head(5)
输出文件解读:
genomic annotation注释：annotation列，注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子
nearest gene annotation（最近基因）注释：多出来的其它列，注释的是peak相对于转录起始位点（TSS）的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，都能够注释到一个离它最近的基因（即使它可能非常远）
针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，可以看到输出有一列是transcriptId
两种注释策略面对不同的研究对象，第一种策略，peak所在的位置可能就是调控本身（比如你要做可变剪切的，最近基因的注释显然不是你关注的点）；而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。
在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。
提取peakAnnolist中的基因，结合clusterProfiler包对peaks内的邻近基因进行富集注释。
Clustering with selected Principal ComponentsR语言：prcomp做主成分分析(PCA)
https://github.com/vallotlab/scChIPseq/blob/master/R_scChIP_seq_analysis.R
每R一点：层次聚类分析实例实战-dist、hclust、heatmap等
降维后进行层次聚类，通常会选择用dist 计算降维后样本距离，这篇文章用的 1-cor() 来代替距离计算，记录一下。
R函数：scale(data, center=T/F, scale=T/F)
center (中心化)：将数据减去均值
scale (标准化)：在中心化后的数据基础上再除以数据的标准差
不太确定哪一种比较好，欢迎交流讨论~
![[27] 《R数据科学》使用tibble实现简单数据框](/aiimages/%5B27%5D+%E3%80%8AR%E6%95%B0%E6%8D%AE%E7%A7%91%E5%AD%A6%E3%80%8B%E4%BD%BF%E7%94%A8tibble%E5%AE%9E%E7%8E%B0%E7%AE%80%E5%8D%95%E6%95%B0%E6%8D%AE%E6%A1%86.png)
