peak注释

ChIPseeker里需要用到一个函数annotatePeak，它可以满足许多需求。

做注释需要注释信息，基因的起始终止，基因有那些内含子，外显子，以及它们的起始终止，非编码区的位置，功能元件的位置等各种信息。

ChIPseeker支持所有的物种，在进行注释时需要的是一个TxDb对象，这个TxDb就包含了我们需要的各种信息，ChIPseeker会把信息抽取出来，用于注释时使用。

require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)

看到这里有个参数tssRegion ，它指定了启动子区域（而启动子区域是没有明确定义的，需要自己指定，这里指定了上下游1kb）

输出

如果在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。

> x
Annotated peaks generated by ChIPseeker
1331/1331  peaks were annotated
Genomic Annotation Summary:
             Feature   Frequency
9           Promoter 48.15927874
4             5' UTR  0.75131480
3             3' UTR  4.20736289
1           1st Exon  0.75131480
7         Other Exon  3.90683696
2         1st Intron  6.53643877
8       Other Intron  4.88354621
6 Downstream (<=300)  0.07513148
5  Distal Intergenic 30.72877536

如果我想看具体的信息呢？你可以用as.GRanges方法，这里只打印前三行：

as.GRanges(x) |> head(3)

关于注释的类型

注释类型一：genomic annotation（annotation这一列）

指peak在基因组的位置：落在什么地方，例如外显子、内含子或是UTR。这一种策略是peak所在位置，可能就是调控的根本，比如你要做可变剪切的，最近基因的注释显然不是你关注的点。

注释类型二：nearest gene annotation（annotation后面的列）

指peak最近的基因：不管peak落在内含子、基因间区还是其他位置，按照peak相对于转录起始位点的距离，都能找到一个离它最近的基因一般做基因表达调控的，会关注promoter区域，离结合位点最近的基因更可能被调控

最近基因的注释信息虽然是以基因为单位给出，但我们针对的是转录起始位点来计算距离，针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，我们可以看到输出有一列是transcriptId.

注释类型三

两种注释有时候还不够，我想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因，annotatePeak也可以做到。
你只要传个参数说你要这个信息，还有什么区间内，就可以了。

f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=5000)

输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出。

基因注释

上面得到的结果都是以geneId（Entrez ID）给出，如果想要Symbol名称，可以再传参数annoDb

library(org.Hs.eg.db)
f = getSampleFiles()[[4]]
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb, 
                  addFlankGeneInfo=TRUE, flankDistance=5000,
                  annoDb = "org.Hs.eg.db")

同一物种的不同版本TxDb

例如TxDb.Hsapiens.UCSC.hg19.knownGene和TxDb.Hsapiens.UCSC.hg38.knownGene 的注释结果是不同的，不能混用。用哪个取决于上游分析比对使用的哪个版本的基因组

不同的版本中基因坐标是不一样的，如果硬要替换，可以使用liftOver将基因组版本坐标进行转换

ChIPseeker支持所有物种

Bioconductor上有30个左右TxDb，也只能覆盖一小部分物种（https://bioconductor.org/packages/3.11/data/annotation/），但UCSC和Ensemble的基因组都可以被ChIPseeker支持，因此所有物种都支持

因为我们可以使用GenomicFeatures包函数来制作TxDb对象：

• makeTxDbFromUCSC： 通过UCSC在线制作TxDb
• makeTxDbFromBiomart: 通过ensembl在线制作TxDb
• makeTxDbFromGRanges：通过GRanges对象制作TxDb
• makeTxDbFromGFF：通过解析GFF文件制作TxDb
  比如我想用人的参考基因信息来做注释，我们可以直接在线从UCSC生成TxDb

例如做裂殖酵母的ChIPseq注释

download.file("ftp://ftp.ebi.ac.uk/pub/databases/pombase/pombe/Chromosome_Dumps/gff3/schizosaccharomyces_pombe.chr.gff3", "schizosaccharomyces_pombe.chr.gff3")require(GenomicFeatures)
spombe <- makeTxDbFromGFF("schizosaccharomyces_pombe.chr.gff3")

有了TxDb怎么查看呢？
最详细的操作在官方文档：https://bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf

不管是从Bioconductor下载的还是自己制作的，都是一个GenomicFeatures对象

如果简单对名称操作，会返回这个注释文件的基本信息。要把TxDb当成一个数据库来对待，而不是一个简单的数据框或者矩阵。因此它的提取方法也会比较特别

如果想看其中包含的类目，可以用columns(txdb)

如果想指定提取转录本或外显子信息，可以：transcripts(txdb) 或者 exons(txdb)

如果想看全部的信息，可以：AnnotationDbi::select(glycine, columns=columns(glycine), keys=keys(glycine), keytype=c("GENEID"))

需要注意，如果使用这个select的时候，同时加载了tidyverse，那么同名的select就会发生冲突导致报错，这时可以用显式指定的形式来规范（如下图）

对自定义区域的注释

除了基因注释还能注释lincRNA
比如就可以利用：https://bioconductor.org/packages/3.11/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts.html

require("TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts")
linc_txdb=TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts
x=annotatePeak(peak, TxDb=linc_txdb)
as.GRanges(x)

如果你单纯只想看那些peak落在CpG上，或者说离CpG最近。annotatePeak的背景注释信息除了TxDb之外，其实它还可以是自定义的GRanges对象，这保证了用户各种各样的需求，因为TxDb也不是万能的，如果能自定义，比如说我就只想看蛋白的结合位点会不会在内含子和外显子的交界处，再比如说我做的并不是编码蛋白的基因表达调控，而是非编码RNA，那么我想要用lncRNA的位置信息来做注释。像这样的需求，ChIPseeker都是可以满足的。

按正负链分开注释

首先ChIPseq数据通常情况下是没有正负链信息的（有特殊的实验可以有），annotatePeak函数有参数是sameStrand，默认是FALSE，你可以给你的peak分别赋正负链，然后传入sameStrand=TRUE，分开做两次，你就可以分开拿到正链和负链的最近基因。

最近基因位置是相对于TSS的，能否相对于整个基因

首先如果peak和TSS有overlap，genomic annotation就是promoter，而最近基因也是同一个，所以在这种情况下，两种注释都指向同一个基因，可以说信息是冗余的，能不能不要冗余信息？这个是可以的，你可以传入参数ignoreOverlap=TRUE，那么最近基因就会去找不overlap的。
最近基因是相对于TSS，如果和TSS有overlap，距离是0, 必须是最近。回到标题的问题，如果我想说只要和基因有overlap就是最近基因，这种情况其实你的最近基因就是host gene，也就是annotation这个column给出来的是相对应的，我们就想找peak所在位置的基因信息，那么这当然也是可以的，默认参数overlap=”TSS”, 如果改为overlap=”all”，它看的就是整个基因而不是TSS，当然distanceToTSS也还是会计算，如果overlap的不是TSS，而是基因体里，并不会因为而设为0。

如果我只想注释上游或者是下游的基因呢

当然也可以，有ignoreUpstream和ignoreDownstream参数，默认都是FALSE。随便你想看上游还是下游，都可以。

注释信息可视化

f = getSampleFiles()[[4]]
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)
png(file="a.png")
plotAnnoPie(x)
dev.off()